[cig-commits] commit: Computes solCx correctly for constant viscosity

[Mercurial]

changeset:   18:868bc23b3588
user:        Walter Landry <wlandry at caltech.edu>
date:        Sun Dec 04 12:09:25 2011 -0800
files:       constants.hxx initial.cxx
description:
Computes solCx correctly for constant viscosity


diff -r d933ef1bd792 -r 868bc23b3588 constants.hxx
--- a/constants.hxx	Sun Dec 04 12:00:01 2011 -0800
+++ b/constants.hxx	Sun Dec 04 12:09:25 2011 -0800
@@ -6,7 +6,7 @@ const double eta_jump=max_eta-min_eta;
 #include <cmath>
 const double log_max_eta=std::log(max_eta);
 const double log_min_eta=std::log(min_eta);
-const double middle=0.5;
+const double middle=0.4;
 const double h(1.0/N);
 const double pi(atan(1.0)*4);
 
diff -r d933ef1bd792 -r 868bc23b3588 initial.cxx
--- a/initial.cxx	Sun Dec 04 12:00:01 2011 -0800
+++ b/initial.cxx	Sun Dec 04 12:09:25 2011 -0800
@@ -14,7 +14,7 @@ void initial(const Model &model, double 
   std::complex<double> phi, psi, dphi, v;
 
   const double pi(std::atan(1.0)*4);
-  const double kx0(pi*middle), kx1(pi*(1-middle));
+  const double kx0(pi*middle), kx1(pi*(middle-1));
   const double cosh0(std::cosh(kx0)), sinh0(std::sinh(kx0)),
     cos0(std::cos(kx0)), sin0(std::sin(kx0)),
     cosh1(std::cosh(kx1)), sinh1(std::sinh(kx1)),
@@ -23,16 +23,16 @@ void initial(const Model &model, double 
   const double eta0(min_eta), eta1(max_eta), rho(0.5);
 
   double m[]={
-    -kx0*cosh0, sinh0-kx0*cosh0, -kx1*cosh1, (sinh1-kx1*cosh1),
+    -kx0*cosh0, sinh0-kx0*cosh0, kx1*cosh1, -(sinh1-kx1*cosh1),
     cosh0+kx0*sinh0, kx0*sinh0, -(cosh1+kx1*sinh1), -kx1*sinh1,
     -kx0*sinh0, cosh0-kx0*sinh0, kx1*sinh1, -(cosh1-kx1*sinh1),
-    -sinh0-kx0*cosh0, -kx0*cosh0, -(sinh1+kx1*cosh1), -kx1*cosh1};
+    -sinh0-kx0*cosh0, -kx0*cosh0, (sinh1+kx1*cosh1), kx1*cosh1};
 
   double rhs[]={
-    ((kx0*cosh0 - sin0)*eta0 - (kx1*cosh1 - sin1)*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};
+    ((sinh0 + kx0*cosh0)*eta0 + (sinh1 + kx1*cosh1)*eta1)*rho/2};
 
   gsl_matrix_view mv=gsl_matrix_view_array(m,4,4);
 
@@ -77,13 +77,14 @@ void initial(const Model &model, double 
                 break;
               case Model::solCx:
                 {
-                  double kx, v, tau, eta;
+                  double kx, v, tau, eta, sign;
                   if(x<middle)
                     {
                       kx=pi*x;
                       v=v0;
                       tau=tau_0;
                       eta=min_eta;
+                      sign=1;
                     }
                   else
                     {
@@ -91,10 +92,11 @@ void initial(const Model &model, double 
                       v=v1;
                       tau=tau_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 + (csh - cs)*rho/2)
+                  p[i][j]=-2*pi*((v/eta + tau)*csh + sign*(csh - cs)*rho/2)
                     *std::cos(pi*y);
                 }
                 break;
@@ -123,25 +125,27 @@ void initial(const Model &model, double 
                     log_etax[i][j]=log_min_eta;
                   }
                 {
-                  double kx, v, tau, eta;
+                  double kx, v, tau, eta, sign;
                   if(x<middle)
                     {
                       kx=pi*x;
                       v=v0;
                       tau=tau_0;
                       eta=min_eta;
+                      sign=1;
                     }
                   else
                     {
-                      kx=pi*(1-x);
+                      kx=pi*(x-1);
                       v=v1;
                       tau=tau_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 - rho*eta*(kx*csh/2 - sn/2))
+                  zx[i][j]=(eta*tau*snh - (v+tau*eta)*kx*csh
+                            - sign*rho*eta*(kx*csh/2 - sn/2))
                     *std::cos(pi*y);
                 }
                 break;
@@ -200,13 +204,14 @@ void initial(const Model &model, double 
                   }
                 {
                   const double pi(atan(1.0)*4);
-                  double kx, v, tau, eta;
+                  double kx, v, tau, eta, sign;
                   if(x<middle)
                     {
                       kx=pi*x;
                       v=v0;
                       tau=tau_0;
                       eta=min_eta;
+                      sign=1;
                     }
                   else
                     {
@@ -214,11 +219,12 @@ void initial(const Model &model, double 
                       v=v1;
                       tau=tau_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 - eta*rho*(cs - csh - kx*snh)/2)
+                  zy[i][j]=(v*csh + (v+tau*eta)*kx*snh
+                            - sign*eta*rho*(cs - csh - kx*snh)/2)
                     *std::sin(pi*y);
                 }
                 break;