[cig-commits] commit: Implement Jacobi iterations. Does not improve things

[Mercurial]

changeset:   49:5354d53bfb2b
user:        Walter Landry <wlandry at caltech.edu>
date:        Wed Feb 15 11:24:28 2012 -0800
files:       main.cxx
description:
Implement Jacobi iterations.  Does not improve things


diff -r 85bab1d05473 -r 5354d53bfb2b main.cxx
--- a/main.cxx	Wed Feb 15 06:34:25 2012 -0800
+++ b/main.cxx	Wed Feb 15 11:24:28 2012 -0800
@@ -29,12 +29,14 @@ int main()
   double zx[Nx+1][Ny], zy[Nx][Ny+1], log_etax[Nx+1][Ny], log_etay[Nx][Ny+1],
     p[Nx][Ny], dp[Nx][Ny], div[Nx][Ny], fx[Nx+1][Ny], fy[Nx][Ny+1];
 
+  double zx_new[Nx+1][Ny], zy_new[Nx][Ny+1], p_new[Nx][Ny];
+
   /* Initial conditions */
 
-  const int max_iterations=1;
-  int n_sweeps=1;
-  const double theta_mom=0.7;
-  const double theta_c=1.0;
+  const int max_iterations=10000;
+  int n_sweeps=10000000;
+  const double theta_mom=0.7/100;
+  const double theta_c=1.0/100;
   const double tolerance=1.0e-6;
 
   Model model(Model::solCx);
@@ -169,7 +171,10 @@ int main()
                                   //             << middle << " "
                                   //             // << zx_jump << " "
                                   //             // << dzx_x_jump << " "
-                                  //             // << dzx_xx_jump << " "
+                                  //             << dzx_xx_jump << " "
+                                  //             << pi*pi*zx[i+1][j]*exp(-log_etax[i+1][j])*eta_jump << " "
+                                  //             << pi*pi*zx[i][j]*exp(-log_etax[i][j])*eta_jump << " "
+                                  //             << pi*pi*zx[i-1][j]*exp(-log_etax[i-1][j])*eta_jump << " "
                                   //             // << dzx_xx << " "
                                   //             << (zx[i+3][j]-2*zx[i+2][j]+zx[i+1][j])/(h*h) << " "
                                   //             << (zx[i+2][j]-2*zx[i+1][j]+zx[i][j])/(h*h) << " "
@@ -205,9 +210,7 @@ int main()
                           double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
 
                           Resid_x[i][j]=Rx;
-                          // zx[i][j]-=theta_mom*Rx/dRx_dzx;
-
-                          // Resid_x[i][j]=dzx_xx;
+                          zx_new[i][j]=zx[i][j]-theta_mom*Rx/dRx_dzx;
                         }
                     }
                 }
@@ -388,7 +391,7 @@ int main()
                           double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
 
                           Resid_y[i][j]=Ry;
-                          // zy[i][j]-=theta_mom*Ry/dRy_dzy;
+                          zy_new[i][j]=zy[i][j]-theta_mom*Ry/dRy_dzy;
                         }
                     }
                 }
@@ -469,7 +472,7 @@ int main()
                 Resid_p[i][j]=Rc;
 
                 dp[i][j]=-theta_c*Rc/dRc_dp;
-                // p[i][j]+=dp[i][j];
+                p_new[i][j]=p[i][j]+dp[i][j];
               }
 
           if(sweep%100==0)
@@ -506,7 +509,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_new[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);
                   }
@@ -530,7 +533,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_new[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);
                   }
@@ -539,11 +542,12 @@ int main()
 
           if(sweep%100==0)
           {
-            // std::cout << sweep << " "
-            //           << max_x_resid << " "
-            //           << max_y_resid << " "
-            //           << max_p_resid << " "
-            //           << "\n";
+            std::cout << iteration << " "
+                      << sweep << " "
+                      << max_x_resid << " "
+                      << max_y_resid << " "
+                      << max_p_resid << " "
+                      << "\n";
 
             std::stringstream ss;
             ss << "zx" << sweep;
@@ -566,22 +570,34 @@ int main()
             write_vtk(ss.str(),Nx,NN,div);
           }
 
-          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_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;
-            }
+          // 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_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_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)
+            for(int j=0;j<Ny+1;++j)
+              {
+                if(j<Ny)
+                  zx[i][j]=zx_new[i][j];
+                if(i<Nx)
+                  zy[i][j]=zy_new[i][j];
+                if(j<Ny && i<Nx)
+                  p[i][j]=p_new[i][j];
+              }
+
         }
     }
   for(int i=0;i<Nx+1;++i)