[cig-commits] commit: Change max_x to max_x_resid etc., theta_momentum to 0.7

[Mercurial]

changeset:   44:2ceeeafb95a9
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Fri Feb 10 14:11:18 2012 -0800
files:       main.cxx
description:
Change max_x to max_x_resid etc., theta_momentum to 0.7


diff -r f82ffcbde572 -r 2ceeeafb95a9 main.cxx
--- a/main.cxx	Wed Feb 08 13:37:12 2012 -0800
+++ b/main.cxx	Fri Feb 10 14:11:18 2012 -0800
@@ -33,7 +33,7 @@ int main()
 
   const int max_iterations=1;
   int n_sweeps=1;
-  const double theta_mom=1.0;
+  const double theta_mom=0.7;
   const double theta_c=1.0;
   const double tolerance=1.0e-6;
 
@@ -58,7 +58,7 @@ int main()
 
       // if(iteration==max_iterations-1)
       //   n_sweeps=1;
-      double max_x_previous(0), max_y_previous(0), max_p_previous(0);
+      double max_x_resid_previous(0), max_y_resid_previous(0), max_p_resid_previous(0);
       for(int sweep=0;sweep<n_sweeps;++sweep)
         {
           /* zx */
@@ -158,6 +158,7 @@ int main()
                                       + dx*dx*dzx_xx_jump/2)/(h*h);
                                   if(x>middle)
                                     dzx_xx_correction*=-1;
+
                                   dzx_xx+=dzx_xx_correction;
 
                                   if(x+h/2>middle && x-h/2<middle)
@@ -184,6 +185,8 @@ int main()
 
                           Resid_x[i][j]=Rx;
                           // zx[i][j]-=theta_mom*Rx/dRx_dzx;
+
+                          // Resid_x[i][j]=dzx_xx;
                         }
                     }
                 }
@@ -401,25 +404,6 @@ int main()
 
                         double dzx_x_correction=(zx_jump + dx*dzx_x_jump)/h;
                         dzx_x-=dzx_x_correction;
-
-                        // if(j==Ny/8)
-                        // std::cout << "dzx_x "
-                        //           << i << " "
-                        //           << j << " "
-                        //           << x << " "
-                        //           << dzx_x << " "
-                        //           // << dzx_x_correction << " "
-                        //           // << dzx_x-dzx_x_correction << " "
-                        //           // << zx_jump/h << " "
-                        //           // << dx*dzx_x_jump/h << " "
-                        //           << (zx[i+3][j]-zx[i+2][j])/h << " "
-                        //           << (zx[i+2][j]-zx[i+1][j])/h << " "
-                        //           << (zx[i+1][j]-zx[i][j])/h << " "
-                        //           << (zx[i][j]-zx[i-1][j])/h << " "
-                        //           << (zx[i-1][j]-zx[i-2][j])/h << " "
-                        //           << (zx[i-2][j]-zx[i-3][j])/h << " "
-                        //           << "\n";
-
                       }
                   }
 
@@ -452,9 +436,10 @@ int main()
             ss << "p_resid" << sweep;
             write_vtk(ss.str(),Nx,NN,Resid_p);
           }
+
           /* Velocity Fix */
 
-          double max_x(0), max_y(0), max_p(0);
+          double max_x_resid(0), max_y_resid(0), max_p_resid(0);
 
           for(int j=0;j<Ny;++j)
             for(int i=0;i<Nx;++i)
@@ -472,11 +457,16 @@ int main()
                        - 2*log_etax[i][j]
                        + (log_etay[i][j] + log_etay[i-1][j])/2)/(h*h/4);
 
+                    if(model==Model::solCx)
+                      {
+                        dlog_etaxx=dlog_etayy=0;
+                      }
+
                     double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
 
                     // zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
 
-                    max_x=std::max(std::fabs(Resid_x[i][j]),max_x);
+                    max_x_resid=std::max(std::fabs(Resid_x[i][j]),max_x_resid);
                   }
                 /* Fix vy */
                 if(j>0)
@@ -491,21 +481,26 @@ int main()
                        - 2*log_etay[i][j]
                        + (log_etax[i][j] + log_etax[i][j-1])/2)/(h*h/4);
 
+                    if(model==Model::solCx)
+                      {
+                        dlog_etaxx=dlog_etayy=0;
+                      }
+
                     double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
 
                     // zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
 
-                    max_y=std::max(fabs(Resid_y[i][j]),max_y);
+                    max_y_resid=std::max(fabs(Resid_y[i][j]),max_y_resid);
                   }
-                max_p=std::max(fabs(Resid_p[i][j]),max_p);
+                max_p_resid=std::max(fabs(Resid_p[i][j]),max_p_resid);
               }
 
           if(sweep%100==0)
           {
             // std::cout << sweep << " "
-            //           << max_x << " "
-            //           << max_y << " "
-            //           << max_p << " "
+            //           << max_x_resid << " "
+            //           << max_y_resid << " "
+            //           << max_p_resid << " "
             //           << "\n";
 
             std::stringstream ss;
@@ -529,22 +524,22 @@ int main()
             write_vtk(ss.str(),Nx,NN,div);
           }
 
-          if(fabs(max_x-max_x_previous)/max_x < tolerance
-             && fabs(max_y-max_y_previous)/max_y < tolerance
-             && fabs(max_p-max_p_previous)/max_p < tolerance)
+          if(fabs(max_x_resid-max_x_resid_previous)/max_x_resid < tolerance
+             && fabs(max_y_resid-max_y_resid_previous)/max_y_resid < tolerance
+             && fabs(max_p_resid-max_p_resid_previous)/max_p_resid < tolerance)
             {
               std::cout << "Solved "
                         << sweep << " "
                         << tolerance << " "
-                        << fabs(max_x-max_x_previous)/max_x << " "
-                        << fabs(max_y-max_y_previous)/max_y << " "
-                        << fabs(max_p-max_p_previous)/max_p << " "
+                        << fabs(max_x_resid-max_x_resid_previous)/max_x_resid << " "
+                        << fabs(max_y_resid-max_y_resid_previous)/max_y_resid << " "
+                        << fabs(max_p_resid-max_p_resid_previous)/max_p_resid << " "
                         << "\n";
               break;
             }
-          max_x_previous=max_x;
-          max_y_previous=max_y;
-          max_p_previous=max_p;
+          max_x_resid_previous=max_x_resid;
+          max_y_resid_previous=max_y_resid;
+          max_p_resid_previous=max_p_resid;
         }
     }
   for(int i=0;i<Nx+1;++i)