[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