[cig-commits] commit: Fix interpolation of velocity to control points to be second order

[Mercurial]

changeset:   62:c996563dfb3d
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Wed Feb 29 14:54:14 2012 -0800
files:       compute_v_on_interface.cxx
description:
Fix interpolation of velocity to control points to be second order


diff -r 674bdfcc3cb8 -r c996563dfb3d compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx	Tue Feb 28 13:58:32 2012 -0800
+++ b/compute_v_on_interface.cxx	Wed Feb 29 14:54:14 2012 -0800
@@ -78,10 +78,7 @@ double compute_dv_yy(const double z[][ny
                      const double log_eta[][ny], const double &dx,
                      const int i, const int j)
 {
-  double dv_yy_p1=(vel(z,log_eta,i+1,j+1) - 2*vel(z,log_eta,i+1,j)
-                    + vel(z,log_eta,i+1,j-1))/(h*h);
-  double dv_yy_p2=(vel(z,log_eta,i+2,j+1) - 2*vel(z,log_eta,i+2,j)
-                    + vel(z,log_eta,i+2,j-1))/(h*h);
+  double dv_yy_p1, dv_yy_p2;
   if(j==0)
     {
       dv_yy_p1=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/(h*h);
@@ -92,12 +89,16 @@ double compute_dv_yy(const double z[][ny
       dv_yy_p1=(-vel(z,log_eta,i+1,j) + vel(z,log_eta,i+1,j-1))/(h*h);
       dv_yy_p2=(-vel(z,log_eta,i+2,j) + vel(z,log_eta,i+2,j-1))/(h*h);
     }
+  else
+    {
+      dv_yy_p1=(vel(z,log_eta,i+1,j+1) - 2*vel(z,log_eta,i+1,j)
+                + vel(z,log_eta,i+1,j-1))/(h*h);
+      dv_yy_p2=(vel(z,log_eta,i+2,j+1) - 2*vel(z,log_eta,i+2,j)
+                + vel(z,log_eta,i+2,j-1))/(h*h);
+    }
   const double dv_yy_p=(1-dx)*dv_yy_p1 + dx*dv_yy_p2;
 
-  double dv_yy_m0=(vel(z,log_eta,i,j+1) - 2*vel(z,log_eta,i,j)
-                    + vel(z,log_eta,i,j-1))/(h*h);
-  double dv_yy_m1=(vel(z,log_eta,i-1,j+1) - 2*vel(z,log_eta,i-1,j)
-                    + vel(z,log_eta,i-1,j-1))/(h*h);
+  double dv_yy_m0, dv_yy_m1;
   if(j==0)
     {
       dv_yy_m0=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/(h*h);
@@ -107,6 +108,13 @@ double compute_dv_yy(const double z[][ny
     {
       dv_yy_m0=(-vel(z,log_eta,i,j) + vel(z,log_eta,i,j-1))/(h*h);
       dv_yy_m1=(-vel(z,log_eta,i-1,j) + vel(z,log_eta,i-1,j-1))/(h*h);
+    }
+  else
+    {
+      dv_yy_m0=(vel(z,log_eta,i,j+1) - 2*vel(z,log_eta,i,j)
+                + vel(z,log_eta,i,j-1))/(h*h);
+      dv_yy_m1=(vel(z,log_eta,i-1,j+1) - 2*vel(z,log_eta,i-1,j)
+                + vel(z,log_eta,i-1,j-1))/(h*h);
     }
   double dv_yy_m=(1-dx)*dv_yy_m1 + dx*dv_yy_m0;
 
@@ -139,13 +147,21 @@ double compute_v(const double z[][ny],
                  const double &dx,
                  const int i, const int j)
 {
+  /* v=v0 + dv*x + ddv*h*h/2, where x is centered on v_p2.  The
+     distance between points is 1 here, which simplifies the
+     expressions. */
+
   double v_p1(vel(z,log_eta,i+1,j)), v_p2(vel(z,log_eta,i+2,j)),
     v_p3(vel(z,log_eta,i+3,j));
-  double v_p((1-dx)*v_p1 + dx*v_p2), v_pp((1-dx)*v_p2 + dx*v_p3);
+  double vp0(v_p2), dvp((v_p3-v_p1)/2), ddvp(v_p3 - 2*v_p2 + v_p1);
+  double v_p(vp0 - (1-dx)*dvp + (1-dx)*(1-dx)*ddvp/2),
+    v_pp(vp0 + dx*dvp + dx*dx*ddvp/2);
 
   double v_m0(vel(z,log_eta,i,j)), v_m1(vel(z,log_eta,i-1,j)),
     v_m2(vel(z,log_eta,i-2,j));
-  double v_m((1-dx)*v_m1 + dx*v_m0), v_mm((1-dx)*v_m2 + dx*v_m1);
+  double vm0(v_m1), dvm((v_m0-v_m2)/2), ddvm(v_m2 - 2*v_m1 + v_m0);
+  double v_m(vm0 + dx*dvm + dx*dx*ddvm/2),
+    v_mm(vm0 - (1-dx)*dvm + (1-dx)*(1-dx)*ddvm/2);
 
   return (4*(max_eta*v_p + min_eta*v_m) - (max_eta*v_pp + min_eta*v_mm)
           - 2*h*jump)