[cig-commits] commit: Added the velocity fix, which seems to work.

[Mercurial]

changeset:   5:135b6a4f6c31
user:        Walter Landry <wlandry at caltech.edu>
date:        Thu Aug 04 02:28:30 2011 -0700
files:       main.cxx
description:
Added the velocity fix, which seems to work.


diff -r 3d2e8ad57043 -r 135b6a4f6c31 main.cxx
--- a/main.cxx	Thu Aug 04 02:06:39 2011 -0700
+++ b/main.cxx	Thu Aug 04 02:28:30 2011 -0700
@@ -118,6 +118,8 @@ int main()
       }
 
 
+  /* Smoothing sweeps */
+
   for(int sweep=0;sweep<n_sweeps;++sweep)
     {
       /* zx */
@@ -306,26 +308,7 @@ int main()
 
                       double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
 
-                      // if(i>30 && i<33 && j==63)
-                      //   std::cout << "Ry "
-                      //             << i << " " 
-                      //             << j << " "
-                      //             << zy[i][j] << " "
-                      //             << Ry/dRy_dzy << " "
-                      //             << Ry << " "
-                      //             << dRy_dzy << " "
-                      //             << 2*dzy_yy << " "
-                      //             << dzy_xx << " "
-                      //             << dzx_yx << " "
-                      //             << 2*(dlog_etayy*zy[i][j] + dlog_etay*dzy_y) << " "
-                      //             << dlog_etaxx*zy[i][j] + dlog_etax*dzy_x << " "
-                      //             << dlog_etayx*zx_avg + dlog_etay*dzx_x << " "
-                      //             << dp_y << " "
-                      //             << fy[i][j] << " "
-                      //             << "\n";
-
                       zy[i][j]-=theta_mom*Ry/dRy_dzy;
-
                     }
                 }
             }
@@ -372,59 +355,62 @@ int main()
             dp[i][j]=-theta_c*Rc/dRc_dp;
             p[i][j]+=dp[i][j];
 
-
-            // std::cout << i << " " << j << " "
-            //           << dzx_x << " "
-            //           << dzy_y << " "
-            //           << zx_avg << " "
-            //           << zy_avg << " "
-            //           << dlog_etax << " "
-            //           << dlog_etay << " "
-            //           << dlog_etaxx << " "
-            //           << dlog_etayy << " "
-            //           << Rc << " "
-            //           << dp[i][j] << " "
-            //           << "\n";
-
-            // std::cout << "sweep "
-            //           << sweep << " "
-            //           << i << " "
-            //           << j << " "
-            //           << zx[i][j] << " "
-            //           << zy[i][j] << " "
-            //           << p[i][j] << " "
-            //           << "\n";
           }
 
-      std::cout << sweep << " "
-                << zx[32][1] << " "
-                << zx[32][32] << " "
-                << zx[32][62] << " "
-                << zx[32][63] << " "
-                << zy[32][1] << " "
-                << zy[32][32] << " "
-                << zy[32][62] << " "
-                << zy[32][63] << " "
-                << p[32][1] << " "
-                << p[32][32] << " "
-                << p[32][62] << " "
-                << p[32][63] << " "
-                << "\n";
+      /* Velocity Fix */
 
+      double max_x(0), max_y(0), max_p(0);
 
-      // std::cout << sweep << " "
-      //           << zx[1][1] << " "
-      //           << zx[10][10] << " "
-      //           << zx[32][32] << " "
-      //           << zx[63][63] << " "
-      //           << zy[1][1] << " "
-      //           << zy[10][10] << " "
-      //           << zy[32][32] << " "
-      //           << zy[63][63] << " "
-      //           << p[1][1] << " "
-      //           << p[10][10] << " "
-      //           << p[32][32] << " "
-      //           << p[63][63] << " "
-      //           << "\n";
+      for(int j=0;j<N;++j)
+        for(int i=0;i<N;++i)
+          {
+            /* Fix vx */
+            if(i>0)
+              {
+                double dlog_etaxx=
+                  ((log_etay[i][j+1] + log_etay[i][j])/2
+                   - 2*log_etax[i][j]
+                   + (log_etay[i-1][j+1] + log_etay[i-1][j])/2)/(h*h/4);
+
+                double dlog_etayy=
+                  ((log_etay[i][j+1] + log_etay[i-1][j+1])/2
+                   - 2*log_etax[i][j]
+                   + (log_etay[i][j] + log_etay[i-1][j])/2)/(h*h/4);
+
+                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(fabs(zx[i][j]),max_x);
+              }
+            /* Fix vy */
+            if(j>0)
+              {
+                double dlog_etayy=
+                  ((log_etax[i+1][j] + log_etax[i][j])/2
+                   - 2*log_etay[i][j]
+                   + (log_etax[i+1][j-1] + log_etax[i][j-1])/2)/(h*h/4);
+
+                double dlog_etaxx=
+                  ((log_etax[i+1][j] + log_etax[i+1][j-1])/2
+                   - 2*log_etay[i][j]
+                   + (log_etax[i][j] + log_etax[i][j-1])/2)/(h*h/4);
+
+                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(zy[i][j]),max_y);
+              }
+            max_p=std::max(fabs(p[i][j]),max_p);
+          }
+
+      if(sweep%100==0)
+        std::cout << sweep << " "
+                  << max_x << " "
+                  << max_y << " "
+                  << max_p << " "
+                  << "\n";
+
     }
 }