[cig-commits] commit: solCX is set correctly when eta!=1

[Mercurial]

changeset:   20:83d2c60ea56b
user:        Walter Landry <wlandry at caltech.edu>
date:        Sun Dec 04 14:12:07 2011 -0800
files:       constants.hxx initial.cxx
description:
solCX is set correctly when eta!=1


diff -r e54175da78d8 -r 83d2c60ea56b constants.hxx
--- a/constants.hxx	Sun Dec 04 12:09:42 2011 -0800
+++ b/constants.hxx	Sun Dec 04 14:12:07 2011 -0800
@@ -1,7 +1,7 @@ const int N(64);
 const int N(64);
 // const double max_eta=1e10;
-const double max_eta=1;
-const double min_eta=1;
+const double max_eta=2;
+const double min_eta=2;
 const double eta_jump=max_eta-min_eta;
 #include <cmath>
 const double log_max_eta=std::log(max_eta);
diff -r e54175da78d8 -r 83d2c60ea56b initial.cxx
--- a/initial.cxx	Sun Dec 04 12:09:42 2011 -0800
+++ b/initial.cxx	Sun Dec 04 14:12:07 2011 -0800
@@ -29,10 +29,10 @@ void initial(const Model &model, double 
     -sinh0-kx0*cosh0, -kx0*cosh0, (sinh1+kx1*cosh1), kx1*cosh1};
 
   double rhs[]={
-    ((kx0*cosh0 - sin0)*eta0 + (kx1*cosh1 - sin1)*eta1)*rho/2,
-    ((cos0 - cosh0 - kx0*sinh0)*eta0 + (cos1 - cosh1 - kx1*sinh1)*eta1)*rho/2,
-    (kx0*sinh0*eta0 + kx1*sinh1*eta1)*rho/2,
-    ((sinh0 + kx0*cosh0)*eta0 + (sinh1 + kx1*cosh1)*eta1)*rho/2};
+    ((kx0*cosh0 - sin0)/eta0 + (kx1*cosh1 - sin1)/eta1)*rho/2,
+    ((cos0 - cosh0 - kx0*sinh0)/eta0 + (cos1 - cosh1 - kx1*sinh1)/eta1)*rho/2,
+    (kx0*sinh0/eta0 + kx1*sinh1/eta1)*rho/2,
+    ((sinh0 + kx0*cosh0)/eta0 + (sinh1 + kx1*cosh1)/eta1)*rho/2};
 
   gsl_matrix_view mv=gsl_matrix_view_array(m,4,4);
 
@@ -56,8 +56,8 @@ void initial(const Model &model, double 
   std::cout << "soln\n";
   gsl_vector_fprintf(stdout,x,"%g ");
 
-  double v0(gsl_vector_get(x,0)), tau_0(gsl_vector_get(x,1)),
-    v1(gsl_vector_get(x,2)), tau_1(gsl_vector_get(x,3));
+  double v0(gsl_vector_get(x,0)), tau_eta_0(gsl_vector_get(x,1)),
+    v1(gsl_vector_get(x,2)), tau_eta_1(gsl_vector_get(x,3));
   gsl_vector_free(x);
 
   for(int i=0;i<N+1;++i)
@@ -77,12 +77,12 @@ void initial(const Model &model, double 
                 break;
               case Model::solCx:
                 {
-                  double kx, v, tau, eta, sign;
+                  double kx, v, tau_eta, eta, sign;
                   if(x<middle)
                     {
                       kx=pi*x;
                       v=v0;
-                      tau=tau_0;
+                      tau_eta=tau_eta_0;
                       eta=min_eta;
                       sign=1;
                     }
@@ -90,13 +90,13 @@ void initial(const Model &model, double 
                     {
                       kx=pi*(x-1);
                       v=v1;
-                      tau=tau_1;
+                      tau_eta=tau_eta_1;
                       eta=max_eta;
                       sign=-1;
                     }                      
                   const double csh(std::cosh(kx)), cs(std::cos(kx)),
                     snh(std::sinh(kx)), sn(std::sin(kx));
-                  p[i][j]=-2*pi*((v/eta + tau)*csh + sign*(csh - cs)*rho/2)
+                  p[i][j]=-2*pi*((v + tau_eta)*eta*csh + sign*(csh - cs)*rho/2)
                     *std::cos(pi*y);
                 }
                 break;
@@ -125,12 +125,12 @@ void initial(const Model &model, double 
                     log_etax[i][j]=log_min_eta;
                   }
                 {
-                  double kx, v, tau, eta, sign;
+                  double kx, v, tau_eta, eta, sign;
                   if(x<middle)
                     {
                       kx=pi*x;
                       v=v0;
-                      tau=tau_0;
+                      tau_eta=tau_eta_0;
                       eta=min_eta;
                       sign=1;
                     }
@@ -138,15 +138,15 @@ void initial(const Model &model, double 
                     {
                       kx=pi*(x-1);
                       v=v1;
-                      tau=tau_1;
+                      tau_eta=tau_eta_1;
                       eta=max_eta;
                       sign=-1;
                     }                      
                   const double csh(std::cosh(kx)), cs(std::cos(kx)),
                     snh(std::sinh(kx)), sn(std::sin(kx));
-                  zx[i][j]=(eta*tau*snh - (v+tau*eta)*kx*csh
-                            - sign*rho*eta*(kx*csh/2 - sn/2))
-                    *std::cos(pi*y);
+                  zx[i][j]=(tau_eta*snh - (v+tau_eta)*kx*csh
+                            - sign*rho*(kx*csh/2 - sn/2)/eta)
+                    *std::cos(pi*y)*eta;
                 }
                 break;
               case Model::sinker:
@@ -204,12 +204,12 @@ void initial(const Model &model, double 
                   }
                 {
                   const double pi(atan(1.0)*4);
-                  double kx, v, tau, eta, sign;
+                  double kx, v, tau_eta, eta, sign;
                   if(x<middle)
                     {
                       kx=pi*x;
                       v=v0;
-                      tau=tau_0;
+                      tau_eta=tau_eta_0;
                       eta=min_eta;
                       sign=1;
                     }
@@ -217,15 +217,15 @@ void initial(const Model &model, double 
                     {
                       kx=pi*(x-1);
                       v=v1;
-                      tau=tau_1;
+                      tau_eta=tau_eta_1;
                       eta=max_eta;
                       sign=-1;
                     }                      
                   const double csh(std::cosh(kx)), cs(std::cos(kx)),
                     snh(std::sinh(kx)), sn(std::sin(kx));
-                  zy[i][j]=(v*csh + (v+tau*eta)*kx*snh
-                            - sign*eta*rho*(cs - csh - kx*snh)/2)
-                    *std::sin(pi*y);
+                  zy[i][j]=(v*csh + (v+tau_eta)*kx*snh
+                            - sign*rho*(cs - csh - kx*snh)/(2*eta))
+                    *std::sin(pi*y)*eta;
                 }
                 break;
               case Model::sinker: