[cig-commits] commit: Enable smoothing and move correction updates out of the regular sweeps and into the iterations

[Mercurial]

changeset:   58:32ba145fcf28
user:        Walter Landry <wlandry at caltech.edu>
date:        Wed Feb 22 05:22:31 2012 -0800
files:       main.cxx
description:
Enable smoothing and move correction updates out of the regular sweeps and into the iterations


diff -r 47c9472c5cd1 -r 32ba145fcf28 main.cxx
--- a/main.cxx	Tue Feb 21 16:46:52 2012 -0800
+++ b/main.cxx	Wed Feb 22 05:22:31 2012 -0800
@@ -34,11 +34,11 @@ int main()
 
   /* Initial conditions */
 
-  const int max_iterations=1;
-  int n_sweeps=1;
+  const int max_iterations=100000;
+  int n_sweeps=10000000;
   const double theta_mom=0.7;
   const double theta_continuity=1.0;
-  const double theta_correction=1.0;
+  const double theta_correction=0.00001;
   const double tolerance=1.0e-6;
 
   Model model(Model::solCx);
@@ -47,8 +47,10 @@ int main()
 
   double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
 
+  compute_corrections(model,zx,zy,log_etax,log_etay,Cx,Cy,Cp,1.0);
   for(int iteration=0;iteration<max_iterations;++iteration)
     {
+      compute_corrections(model,zx,zy,log_etax,log_etay,Cx,Cy,Cp,theta_correction);
       /* Smoothing sweeps */
 
       // if(iteration==max_iterations-1)
@@ -56,7 +58,7 @@ int main()
       double max_x_resid_previous(0), max_y_resid_previous(0), max_p_resid_previous(0);
       for(int sweep=0;sweep<n_sweeps;++sweep)
         {
-          compute_corrections(model,zx,zy,log_etax,log_etay,Cx,Cy,Cp,theta_correction);
+          // compute_corrections(model,zx,zy,log_etax,log_etay,Cx,Cy,Cp,theta_correction);
 
           /* zx */
           for(int rb=0;rb<2;++rb)
@@ -132,7 +134,6 @@ int main()
                           const double zy_avg=(zy[i][j] + zy[i-1][j] +
                                                zy[i][j+1] + zy[i-1][j+1])/4;
 
-                          double theta(theta_mom);
                           /* Now do the jump conditions */
                           if(model==Model::solCx)
                             {
@@ -172,7 +173,7 @@ int main()
                           double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
 
                           Resid_x[i][j]=Rx;
-                          // zx[i][j]-=theta*Rx/dRx_dzx;
+                          zx[i][j]-=theta_mom*Rx/dRx_dzx;
                         }
                     }
                 }
@@ -258,7 +259,6 @@ int main()
                           const double zx_avg=(zx[i][j] + zx[i][j-1] +
                                                zx[i+1][j] + zx[i+1][j-1])/4;
 
-                          double theta(theta_mom);
                           /* Now do the jump conditions */
                           if(model==Model::solCx)
                             {
@@ -288,7 +288,7 @@ int main()
                           double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
 
                           Resid_y[i][j]=Ry;
-                          // zy[i][j]-=theta*Ry/dRy_dzy;
+                          zy[i][j]-=theta_mom*Ry/dRy_dzy;
                         }
                     }
                 }
@@ -325,8 +325,6 @@ int main()
                 double zx_avg=(zx[i+1][j] + zx[i][j])/2;
                 double zy_avg=(zy[i][j+1] + zy[i][j])/2;
 
-                double theta(theta_continuity);
-
                 /* Apply the jump condition */
                 if(model==Model::solCx)
                   {
@@ -352,8 +350,8 @@ int main()
 
                 Resid_p[i][j]=Rc;
 
-                dp[i][j]=-theta*Rc/dRc_dp;
-                // p[i][j]+=dp[i][j];
+                dp[i][j]=-theta_continuity*Rc/dRc_dp;
+                p[i][j]+=dp[i][j];
               }
 
           if(sweep%1000==0)
@@ -393,7 +391,7 @@ int main()
 
                     double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
 
-                    // zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
+                    zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
 
                     max_x_resid=std::max(std::fabs(Resid_x[i][j]),max_x_resid);
                   }
@@ -417,7 +415,7 @@ int main()
 
                     double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
 
-                    // zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
+                    zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
 
                     max_y_resid=std::max(fabs(Resid_y[i][j]),max_y_resid);
                   }