[cig-commits] commit: Solve using the corrections to the residual. Only works for small variations in viscosity.

[Mercurial]

changeset:   65:4e071e769780
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Wed Feb 29 19:20:03 2012 -0800
files:       constants.hxx main.cxx
description:
Solve using the corrections to the residual.  Only works for small variations in viscosity.


diff -r 57ee741e9644 -r 4e071e769780 constants.hxx
--- a/constants.hxx	Wed Feb 29 15:31:48 2012 -0800
+++ b/constants.hxx	Wed Feb 29 19:20:03 2012 -0800
@@ -1,16 +1,16 @@ const int N(1024);
-const int N(1024);
+const int N(64);
 const int Nx(N);
 const int Ny(N);
 const double min_eta=1;
-const double max_eta=1e2;
+const double max_eta=10;
 const double eta_jump=max_eta-min_eta;
 #include <cmath>
 const double log_max_eta=std::log(max_eta);
 const double log_min_eta=std::log(min_eta);
 // const double middle=(25 + 0.00000001)/64;
-const double middle=(25 + 0.00000001)/64 + 1/(2*Nx);
+// const double middle=(25 + 0.00000001)/64 + 1/(2*Nx);
 // const double middle=(25 - 1/32.0)/64;
-// const double middle=0.4;
+const double middle=0.4;
 const double h(1.0/N);
 const double pi(atan(1.0)*4);
 
diff -r 57ee741e9644 -r 4e071e769780 main.cxx
--- a/main.cxx	Wed Feb 29 15:31:48 2012 -0800
+++ b/main.cxx	Wed Feb 29 19:20:03 2012 -0800
@@ -46,7 +46,7 @@ int main()
   /* Initial conditions */
 
   const int max_iterations=1;
-  int n_sweeps=1;
+  int n_sweeps=10000000;
   const double theta_mom=0.7;
   const double theta_continuity=1.0;
   const double tolerance=1.0e-6;
@@ -152,7 +152,8 @@ int main()
                           /* Compute the residual and update zx */
 
                           double Cx, dCx_dzx;
-                          compute_Cx(model,zx,zy,log_etax,log_etay,i,j,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)
@@ -160,10 +161,11 @@ int main()
                             + 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;
+                          double dRx_dzx=
+                            -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
 
                           Resid_x[i][j]=Rx;
-                          // zx[i][j]-=theta_mom*Rx/dRx_dzx;
+                          zx[i][j]-=theta_mom*Rx/dRx_dzx;
                         }
                     }
                 }
@@ -270,7 +272,8 @@ 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;
+                          double dRy_dzy=
+                            -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy;
 
                           // if(j==Ny/8 && i*h>middle-10*h && i*h<middle+10*h)
                           //   std::cout << "Ry "
@@ -292,7 +295,7 @@ int main()
 
 
                           Resid_y[i][j]=Ry;
-                          // zy[i][j]-=theta_mom*Ry/dRy_dzy;
+                          zy[i][j]-=theta_mom*Ry/dRy_dzy;
                         }
                     }
                 }
@@ -361,7 +364,7 @@ int main()
                 Resid_p[i][j]=Rc;
 
                 dp[i][j]=-theta_continuity*Rc/dRc_dp;
-                // p[i][j]+=dp[i][j];
+                p[i][j]+=dp[i][j];
               }
 
           if(sweep%output_interval==0)
@@ -404,7 +407,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);
                   }
@@ -428,14 +431,14 @@ 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);
                   }
                 max_p_resid=std::max(fabs(Resid_p[i][j]),max_p_resid);
               }
 
-          if(sweep%1000==0)
+          if(sweep%100==0)
           {
             std::cout << "sweep "
                       << iteration << " "