[cig-commits] commit: Set dvy_y correctly near the boundary

[Mercurial]

changeset:   30:2d4847e70d29
user:        Walter Landry <wlandry at caltech.edu>
date:        Tue Feb 07 11:30:31 2012 -0800
files:       compute_v_on_interface.cxx
description:
Set dvy_y correctly near the boundary


diff -r 7904eebd094d -r 2d4847e70d29 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx	Tue Feb 07 05:23:33 2012 -0800
+++ b/compute_v_on_interface.cxx	Tue Feb 07 11:30:31 2012 -0800
@@ -96,7 +96,7 @@ void compute_v_on_interface(double zx[N+
   /* Iterate to find v on the boundary */
 
   const double relax=0.25;
-  const double tolerance=1e-6;
+  const double tolerance=1e-8;
   double error=0;
   // for(int ii=0;ii<100;++ii)
   do
@@ -124,9 +124,7 @@ void compute_v_on_interface(double zx[N+
       gsl_spline *vy_spline=gsl_spline_alloc(gsl_interp_cspline, N+1);
       gsl_spline_init(vy_spline, y_y, vy, N+1);
 
-      dvy_y[0]=0;
-      dvy_y[N]=0;
-      for(int i=1;i<N+1;++i)
+      for(int i=0;i<N;++i)
         dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i],vy_accel);
 
       for(int i=0;i<N;++i)
@@ -145,19 +143,20 @@ void compute_v_on_interface(double zx[N+
             /(min_eta*a1*a2*(b1+b2) + max_eta*b1*b2*(a1+a2))
             - a1*a2*b1*b2*zx_jump/(min_eta*a1*a2*(b1+b2)
                                    + max_eta*b1*b2*(a1+a2));
-          if(i!=0 && i!=N-1)
-            {
-              error=
-                std::max(error,std::fabs((vx[i]-temp)
-                                         /(std::fabs(vx[i])+std::fabs(temp))));
-              if(i==100)
-              std::cout << "vx error "
-                        << i << " "
-                        << vx[i] << " "
-                        << temp << " "
-                        << error << " "
-                        << "\n";
-            }
+          if(vx[i]!=0 || temp!=0)
+            error=
+              std::max(error,std::fabs((vx[i]-temp)
+                                       /(std::fabs(vx[i])+std::fabs(temp))));
+          if(i<7)
+            std::cout << "vx error "
+                      << i << " "
+                      << vx[i] << " "
+                      << dvx_y[i] << " "
+                      << dvx_yy[i] << " "
+                      << temp << " "
+                      << error << " "
+                      << "\n";
+
           vx[i]=(1-relax)*vx[i] + relax*temp;
         }
 
@@ -170,21 +169,23 @@ void compute_v_on_interface(double zx[N+
             /(min_eta*a1*a2*(b1+b2) + max_eta*b1*b2*(a1+a2))
             - a1*a2*b1*b2*zy_jump/(min_eta*a1*a2*(b1+b2)
                                    + max_eta*b1*b2*(a1+a2));
-          if(i!=0 && i!=N)
-            {
-              error=
-                std::max(error,std::fabs((vy[i]-temp)
-                                         /(std::fabs(vy[i])+std::fabs(temp))));
-              if(i==100)
-              std::cout << "vy error "
-                        << i << " "
-                        << vy[i] << " "
-                        << temp << " "
-                        << std::fabs((vy[i]-temp)
-                                     /(std::fabs(vy[i])+std::fabs(temp))) << " "
-                        << error << " "
-                        << "\n";
-            }
+          if(vy[i]!=0 || temp!=0)
+            error=
+              std::max(error,std::fabs((vy[i]-temp)
+                                       /(std::fabs(vy[i])+std::fabs(temp))));
+
+          if(i<4)
+            std::cout << "vy error "
+                      << i << " "
+                      << vy[i] << " "
+                      << dvy_y[i] << " "
+                      << dvy_yy[i] << " "
+                      << temp << " "
+                      << std::fabs((vy[i]-temp)
+                                   /(std::fabs(vy[i])+std::fabs(temp))) << " "
+                      << error << " "
+                      << "\n";
+
           vy[i]=(1-relax)*vy[i] + relax*temp;
         }
       std::cout << "error "