[cig-commits] commit: Now solCx has vertical gravity

[Mercurial]

changeset:   15:98bf302c0eb9
user:        Walter Landry <wlandry at caltech.edu>
date:        Sun Dec 04 00:32:52 2011 -0800
files:       initial.cxx main.cxx
description:
Now solCx has vertical gravity


diff -r b41ecaa82b4d -r 98bf302c0eb9 initial.cxx
--- a/initial.cxx	Sat Dec 03 20:47:20 2011 -0800
+++ b/initial.cxx	Sun Dec 04 00:32:52 2011 -0800
@@ -35,7 +35,12 @@ void initial(const Model &model, double 
                   const double csh(std::cosh(kx)), cs(std::cos(kx)),
                     snh(std::sinh(kx)), sn(std::sin(kx));
                 
-                  p[i][j]=-2*pi*((v0 + tau_0)*csh + (sn + snh)/4)*std::cos(pi*y);
+                  // double tau_xx=2*(csh*tau_0 - kx*snh*(v0+tau_0) - kx*snh/4);
+                  // double dux_x=csh*tau_0 - (v0+tau_0)*(csh + kx*snh)
+                  //   - (csh + kx*snh - cs)/4;
+                  // p[i][j]=pi*(2*dux_x - tau_xx);
+
+                  p[i][j]=-2*pi*((v0 + tau_0)*csh + (csh - cs)/4)*std::cos(pi*y);
                 }
                 break;
               default:
@@ -49,8 +54,8 @@ void initial(const Model &model, double 
 
             zx[i][j]=0;
 
-            // fx[i][j]=0;
-            fx[i][j]=pi*pi*cos(pi*y)*cos(pi*x);
+            fx[i][j]=0;
+            // fx[i][j]=pi*pi*cos(pi*y)*cos(pi*x);
             switch(model)
               {
               case Model::solKz:
@@ -70,8 +75,8 @@ void initial(const Model &model, double 
                   const double kx=pi*x;
                   const double csh(std::cosh(kx)), cs(std::cos(kx)),
                     snh(std::sinh(kx)), sn(std::sin(kx));
-                  zx[i][j]=(tau_0*snh - (v0+tau_0)*kx*csh
-                            - (cs - csh + kx*snh)/4)*std::cos(pi*y);
+                  zx[i][j]=(tau_0*snh - (v0+tau_0)*kx*csh - (kx*csh/2 - sn/2)/2)
+                    *std::cos(pi*y);
                 }
                 break;
               case Model::sinker:
@@ -118,8 +123,8 @@ void initial(const Model &model, double 
                 log_etay[i][j]=y*log_max_eta;
                 break;
               case Model::solCx:
-                fy[i][j]=0;
-                // fy[i][j]=sin(pi*y)*cos(pi*x);
+                // fy[i][j]=0;
+                fy[i][j]=pi*pi*sin(pi*y)*cos(pi*x);
                 if(x>middle)
                   {
                     log_etay[i][j]=log_max_eta;
@@ -133,7 +138,7 @@ void initial(const Model &model, double 
                   const double kx=pi*x;
                   const double csh(std::cosh(kx)), cs(std::cos(kx)),
                     snh(std::sinh(kx)), sn(std::sin(kx));
-                  zy[i][j]=(v0*csh + (v0+tau_0)*kx*snh + (kx*csh/2 - sn/2)/2)
+                  zy[i][j]=(v0*csh + (v0+tau_0)*kx*snh - (cs - csh - kx*snh)/4)
                     *std::sin(pi*y);
                 }
                 break;
diff -r b41ecaa82b4d -r 98bf302c0eb9 main.cxx
--- a/main.cxx	Sat Dec 03 20:47:20 2011 -0800
+++ b/main.cxx	Sun Dec 04 00:32:52 2011 -0800
@@ -170,17 +170,17 @@ int main()
                                   double dzx_xx_jump=-dvx_yy[j] + p_jump;
                                   dzx_xx_correction=eta_jump*(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
                                   // dzx_xx-=eta_jump*(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
-                                  if(n_sweeps==1)
-                                    std::cout << "dzx_xx "
-                                              << i << " "
-                                              << j << " "
-                                              << dzx_xx << " "
-                                              << dzx_xx_correction << " "
-                                              << eta_jump*zx_jump << " "
-                                              << zx[i][j] << " "
-                                              << zx[i+1][j] << " "
-                                              << zx[i-1][j] << " "
-                                              << "\n";
+                                  // if(n_sweeps==1)
+                                  //   std::cout << "dzx_xx "
+                                  //             << i << " "
+                                  //             << j << " "
+                                  //             << dzx_xx << " "
+                                  //             << dzx_xx_correction << " "
+                                  //             << eta_jump*zx_jump << " "
+                                  //             << zx[i][j] << " "
+                                  //             << zx[i+1][j] << " "
+                                  //             << zx[i-1][j] << " "
+                                  //             << "\n";
                                   if(dx>h/2)
                                     {
                                       dx-=h/2;
@@ -323,6 +323,9 @@ int main()
                           double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
 
                           Resid_y[i][j]=Ry;
+                          // Resid_y[i][j]=fy[i][j];
+                          if(i==N-1)
+                            Resid_y[i][j]=0;
                           // zy[i][j]-=theta_mom*Ry/dRy_dzy;
                         }
                     }