[cig-commits] commit: zx smoothing seems to work

[Mercurial]

changeset:   1:49a5cfe71a1e
user:        Walter Landry <wlandry at caltech.edu>
date:        Wed Aug 03 23:57:31 2011 -0700
files:       main.cxx
description:
zx smoothing seems to work


diff -r 52dc23ba6718 -r 49a5cfe71a1e main.cxx
--- a/main.cxx	Wed Aug 03 23:29:26 2011 -0700
+++ b/main.cxx	Wed Aug 03 23:57:31 2011 -0700
@@ -16,9 +16,11 @@ int main()
   const double max_eta=2;
   const double log_max_eta=std::log(max_eta);
   const double middle=0.4;
-  const int n_sweeps=1;
+  const int n_sweeps=10000;
   const double h(1.0/N);
   const double pi(atan(1.0)*4);
+  const double theta_mom=1.0;
+
 
   Model model(Model::solCx);
 
@@ -28,7 +30,12 @@ int main()
         if(j<N)
           {
             double x(i*h), y((j+0.5)*h);
-            zx[i][j]=1.2 + 3.4*x + 4.5*y;
+
+            if(i==0 || i==N)
+              zx[i][j]=0;
+            else
+              zx[i][j]=1.2 + 3.4*x + 4.5*y;
+
             // zx[i][j]=0;
             fx[i][j]=0;
             if(model==Model::solCx)
@@ -62,8 +69,8 @@ int main()
         if(i<N)
           {
             double x((i+0.5)*h), y(j*h);
-            zy[i][j]=5.6 + 7.8*x + 9.0*y + 2.1*x*y;
-            // zy[i][j]=0;
+            // zy[i][j]=5.6 + 7.8*x + 9.0*y + 2.1*x*y;
+            zy[i][j]=0;
 
             if(model==Model::solCx || model==Model::solKz)
               {
@@ -176,7 +183,24 @@ int main()
                       double dlog_etaxy=
                         ((log_etay[i][j+1] - log_etay[i-1][j+1])
                          - (log_etay[i][j] - log_etay[i-1][j]))/(h*h);
+
+                      const double zy_avg=(zy[i][j] + zy[i-1][j] +
+                                           zy[i][j+1] + zy[i-1][j+1])/4;
+
                       /* Now do the jump conditions */
+
+
+                      /* Compute the residual and update zx */
+
+                      double Rx=2*dzx_xx + dzx_yy + dzy_xy 
+                        + 2*(dlog_etaxx*zx[i][j] + dlog_etax*dzx_x)
+                        + dlog_etayy*zx[i][j] + dlog_etay*dzx_y
+                        + dlog_etaxy*zy_avg + dlog_etax*dzy_y
+                        + dp_x - fx[i][j];
+
+                      double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
+
+                      zx[i][j]-=theta_mom*Rx/dRx_dzx;
 
                       // const double x(i*h), y((j+0.5)*h);
                       // std::cout << i << " "
@@ -205,5 +229,12 @@ int main()
                 }
             }
         }
+
+      std::cout << sweep << " "
+                << zx[1][1] << " "
+                << zx[10][10] << " "
+                << zx[32][32] << " "
+                << zx[63][63] << " "
+                << "\n";
     }
 }