[cig-commits] commit: Compute dRx_dzx and dRy_dzy correctly near the boundary

[Mercurial]

changeset:   74:114d616b8587
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Mon Mar 05 17:40:26 2012 -0800
files:       compute_Cx.cxx main.cxx
description:
Compute dRx_dzx and dRy_dzy correctly near the boundary


diff -r fa26e2bf82f3 -r 114d616b8587 compute_Cx.cxx
--- a/compute_Cx.cxx	Fri Mar 02 16:27:50 2012 -0800
+++ b/compute_Cx.cxx	Mon Mar 05 17:40:26 2012 -0800
@@ -42,13 +42,16 @@ void compute_Cx(const Model &model, cons
 
           double eta=exp(log_etax[i][j]);
           double delta=(h-std::fabs(dx))/h;
-          double dvx_yy_dv=-2*eta*eta/((h*h)*(min_eta*min_eta + max_eta*max_eta));
+          double dvx_yy_dv=-2*eta*eta/(h*h*(min_eta*min_eta + max_eta*max_eta));
+          if(j==0 || j==Ny-1)
+            dvx_yy_dv/=2;
           double dzy_yx_dv=-eta_jump*dvx_yy_dv;
           double dvy_y_dv=-dzy_yx_dv*h/(min_eta+max_eta);
           double dzx_x_dv=eta_jump*dvy_y_dv;
           double dvp_dv=delta/2 + delta*delta/2;
           double dvpp_dv=-(1-delta)/2 + (1-delta)*(1-delta)/2;
-          double dvx_dv=(eta*(4*dvp_dv - dvpp_dv)  + 2*h*dzx_x_dv)/(3*(min_eta+max_eta));
+          double dvx_dv=(eta*(4*dvp_dv - dvpp_dv)  + 2*h*dzx_x_dv)
+            /(3*(min_eta+max_eta));
           double dp_x_dv=2*eta_jump*dvx_yy_dv;
           double dzx_xx_dv=-eta_jump*dvx_yy_dv + dp_x_dv;
           double dp_dv=-2*eta_jump*dvy_y_dv;
diff -r fa26e2bf82f3 -r 114d616b8587 main.cxx
--- a/main.cxx	Fri Mar 02 16:27:50 2012 -0800
+++ b/main.cxx	Mon Mar 05 17:40:26 2012 -0800
@@ -70,6 +70,8 @@ int main()
     for(int j=0;j<Ny;++j)
       p[i][j]=0;
 
+  zy[0][10]=1;
+
   double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
 
   for(int sweep=0;sweep<n_sweeps;++sweep)
@@ -100,17 +102,11 @@ int main()
                       double dzx_y, dzx_yy;
                       if(j==0)
                         {
-                          // dzx_yy=
-                          //   (-zx[i][Ny-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
-                          // dzx_y=(zx[i][j+1] + zx[i][Ny-1])/(2*h);
                           dzx_yy=(-zx[i][j] + zx[i][j+1])/(h*h);
                           dzx_y=(zx[i][j+1] - zx[i][j])/(2*h);
                         }
                       else if(j==Ny-1)
                         {
-                          // dzx_yy=
-                          //   (zx[i][j-1] - 2*zx[i][j] - zx[i][0])/(h*h);
-                          // dzx_y=(-zx[i][0] - zx[i][j-1])/(2*h);
                           dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
                           dzx_y=(zx[i][j] - zx[i][j-1])/(2*h);
                         }
@@ -166,14 +162,17 @@ int main()
                       double Cx, dCx_dzx;
                       compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,
                                  dCx_dzx);
+
                       double Rx=2*dzx_xx + dzx_yy + dzy_xy 
                         + 2*(dlog_etaxx*zx[i][j] + dlog_etax*dzx_x)
                         + dlog_etayy*zx[i][j] + dlog_etay*dzx_y
                         + dlog_etaxy*zy_avg + dlog_etax*dzy_y
                         - dp_x - fx[i][j] + Cx;
 
-                      double dRx_dzx=
-                        -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
+                      double dRx_dzx=-6/(h*h);
+                      if(j==0 || j==Ny-1)
+                        dRx_dzx=-5/(h*h);
+                      dRx_dzx+=2*dlog_etaxx + dlog_etayy + dCx_dzx;
 
                       Resid_x[i][j]=Rx;
                       if(!jacobi)
@@ -280,8 +279,10 @@ int main()
                         + dlog_etayx*zx_avg + dlog_etay*dzx_x
                         - dp_y - fy[i][j] + Cy;
 
-                      double dRy_dzy=
-                        -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy;
+                      double dRy_dzy=-6/(h*h);
+                      if(i==0 || i==Nx-1)
+                        dRy_dzy=-5/(h*h);
+                      dRy_dzy+=2*dlog_etayy + dlog_etaxx + dCy_dzy;
 
                       Resid_y[i][j]=Ry;
                       if(!jacobi)
@@ -392,8 +393,11 @@ int main()
                 double Cx, dCx_dzx;
                 compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,
                            dCx_dzx);
-                double dRx_dzx=
-                  -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
+
+                double dRx_dzx=-6/(h*h);
+                if(j==0 || j==Ny-1)
+                  dRx_dzx=-5/(h*h);
+                dRx_dzx+=2*dlog_etaxx + dlog_etayy + dCx_dzx;
 
                 if(!jacobi)
                   zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
@@ -423,8 +427,10 @@ int main()
                 double Cy, dCy_dzy;
                 compute_Cy(model,zx,zy,log_etax,log_etay,i,j,Cy,
                            dCy_dzy);
-                double dRy_dzy=
-                  -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy;
+                double dRy_dzy=-6/(h*h);
+                if(i==0 || i==Nx-1)
+                  dRy_dzy=-5/(h*h);
+                dRy_dzy+=2*dlog_etayy + dlog_etaxx + dCy_dzy;
 
                 if(!jacobi)
                   zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);