[cig-commits] commit: Make the stresses match for solCx

Mercurial hg at geodynamics.org
Fri Feb 10 16:00:20 PST 2012


changeset:   22:abd75122918b
user:        Walter Landry <wlandry at caltech.edu>
date:        Sun Dec 04 16:11:05 2011 -0800
files:       constants.hxx initial.cxx
description:
Make the stresses match for solCx


diff -r ff2dd6947300 -r abd75122918b constants.hxx
--- a/constants.hxx	Sun Dec 04 15:17:42 2011 -0800
+++ b/constants.hxx	Sun Dec 04 16:11:05 2011 -0800
@@ -1,12 +1,11 @@ const int N(64);
 const int N(64);
-// const double max_eta=1e10;
-const double max_eta=2;
-const double min_eta=2;
+const double min_eta=.0001;
+const double max_eta=1000;
 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.4;
+const double middle=0.5;
 const double h(1.0/N);
 const double pi(atan(1.0)*4);
 
diff -r ff2dd6947300 -r abd75122918b initial.cxx
--- a/initial.cxx	Sun Dec 04 15:17:42 2011 -0800
+++ b/initial.cxx	Sun Dec 04 16:11:05 2011 -0800
@@ -25,14 +25,15 @@ void initial(const Model &model, double 
   double m[]={
     -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};
+    -kx0*sinh0*eta0, (cosh0-kx0*sinh0)*eta0, kx1*sinh1*eta1, -(cosh1-kx1*sinh1)*eta1,
+    (-sinh0-kx0*cosh0)*eta0, -kx0*cosh0*eta0,
+    (sinh1+kx1*cosh1)*eta1, kx1*cosh1*eta1};
 
   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*sinh0 + kx1*sinh1)*rho/2,
+    ((sinh0 + kx0*cosh0) + (sinh1 + kx1*cosh1))*rho/2};
 
   gsl_matrix_view mv=gsl_matrix_view_array(m,4,4);
 
@@ -40,6 +41,10 @@ void initial(const Model &model, double 
   gsl_matrix_fprintf(stdout,&mv.matrix,"%g ");
 
   gsl_vector_view rhsv=gsl_vector_view_array(rhs,4);
+
+  std::cout << "rhs\n";
+  gsl_vector_fprintf(stdout,&rhsv.vector,"%g ");
+
   gsl_vector *x=gsl_vector_alloc(4);
   int s;
   gsl_permutation *perm=gsl_permutation_alloc(4);
@@ -48,11 +53,6 @@ void initial(const Model &model, double 
 
   gsl_permutation_free(perm);
 
-
-  // std::cout << "inverse\n";
-  // gsl_matrix_fprintf(stdout,&mv.matrix,"%g ");
-  std::cout << "rhs\n";
-  gsl_vector_fprintf(stdout,&rhsv.vector,"%g ");
   std::cout << "soln\n";
   gsl_vector_fprintf(stdout,x,"%g ");
 
@@ -96,6 +96,15 @@ void initial(const Model &model, double 
                     }                      
                   const double csh(std::cosh(kx)), cs(std::cos(kx)),
                     snh(std::sinh(kx)), sn(std::sin(kx));
+
+                  /* set pressure to stresses to make sure that they
+                     are continuous */
+                  // p[i][j]=v*(-kx*snh*eta) + tau_eta*(csh-kx*snh)*eta
+                  //   -sign*(kx*snh)*rho/2;
+
+                  // p[i][j]=v*(-snh-kx*csh)*eta - tau_eta*kx*csh*eta
+                  //   -sign*(snh+kx*csh)*rho/2;
+
                   p[i][j]=-2*pi*((v + tau_eta)*eta*csh + sign*(csh - cs)*rho/2)
                     *std::cos(pi*y);
                 }



More information about the CIG-COMMITS mailing list