[cig-commits] commit: Seems to compute v0 etc correctly, but not zx
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:16 PST 2012
changeset: 17:d933ef1bd792
user: Walter Landry <wlandry at caltech.edu>
date: Sun Dec 04 12:00:01 2011 -0800
files: constants.hxx initial.cxx wscript
description:
Seems to compute v0 etc correctly, but not zx
diff -r 5a513d0f7ec3 -r d933ef1bd792 constants.hxx
--- a/constants.hxx Sun Dec 04 00:39:47 2011 -0800
+++ b/constants.hxx Sun Dec 04 12:00:01 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.4;
+const double middle=0.5;
const double h(1.0/N);
const double pi(atan(1.0)*4);
diff -r 5a513d0f7ec3 -r d933ef1bd792 initial.cxx
--- a/initial.cxx Sun Dec 04 00:39:47 2011 -0800
+++ b/initial.cxx Sun Dec 04 12:00:01 2011 -0800
@@ -4,6 +4,7 @@
#include "write_vtk.hxx"
#include <cmath>
+#include <gsl/gsl_linalg.h>
void initial(const Model &model, double zx[N+1][N], double zy[N][N+1],
double log_etax[N+1][N], double log_etay[N][N+1],
@@ -11,7 +12,53 @@ void initial(const Model &model, double
{
const double A(min_eta*(max_eta-min_eta)/(max_eta+min_eta));
std::complex<double> phi, psi, dphi, v;
- const double v0(1.0), tau_0(1.0);
+
+ const double pi(std::atan(1.0)*4);
+ const double kx0(pi*middle), kx1(pi*(1-middle));
+ 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)),
+ cos1(std::cos(kx1)), sin1(std::sin(kx1));
+
+ const double eta0(min_eta), eta1(max_eta), rho(0.5);
+
+ 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};
+
+ 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};
+
+ gsl_matrix_view mv=gsl_matrix_view_array(m,4,4);
+
+ std::cout << "m\n";
+ gsl_matrix_fprintf(stdout,&mv.matrix,"%g ");
+
+ gsl_vector_view rhsv=gsl_vector_view_array(rhs,4);
+ gsl_vector *x=gsl_vector_alloc(4);
+ int s;
+ gsl_permutation *perm=gsl_permutation_alloc(4);
+ gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
+ gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,x);
+
+ 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 ");
+
+ 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));
+ gsl_vector_free(x);
for(int i=0;i<N+1;++i)
for(int j=0;j<N+1;++j)
@@ -30,11 +77,25 @@ void initial(const Model &model, double
break;
case Model::solCx:
{
- const double pi(atan(1.0)*4);
- const double kx=pi*x;
+ double kx, v, tau, eta;
+ if(x<middle)
+ {
+ kx=pi*x;
+ v=v0;
+ tau=tau_0;
+ eta=min_eta;
+ }
+ else
+ {
+ kx=pi*(x-1);
+ v=v1;
+ tau=tau_1;
+ eta=max_eta;
+ }
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 + (csh - cs)/4)*std::cos(pi*y);
+ p[i][j]=-2*pi*((v/eta + tau)*csh + (csh - cs)*rho/2)
+ *std::cos(pi*y);
}
break;
default:
@@ -62,11 +123,25 @@ void initial(const Model &model, double
log_etax[i][j]=log_min_eta;
}
{
- const double pi(atan(1.0)*4);
- const double kx=pi*x;
+ double kx, v, tau, eta;
+ if(x<middle)
+ {
+ kx=pi*x;
+ v=v0;
+ tau=tau_0;
+ eta=min_eta;
+ }
+ else
+ {
+ kx=pi*(1-x);
+ v=v1;
+ tau=tau_1;
+ eta=max_eta;
+ }
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 - (kx*csh/2 - sn/2)/2)
+ zx[i][j]=
+ (eta*tau*snh - (v+tau*eta)*kx*csh - rho*eta*(kx*csh/2 - sn/2))
*std::cos(pi*y);
}
break;
@@ -125,10 +200,25 @@ void initial(const Model &model, double
}
{
const double pi(atan(1.0)*4);
- const double kx=pi*x;
+ double kx, v, tau, eta;
+ if(x<middle)
+ {
+ kx=pi*x;
+ v=v0;
+ tau=tau_0;
+ eta=min_eta;
+ }
+ else
+ {
+ kx=pi*(x-1);
+ v=v1;
+ tau=tau_1;
+ eta=max_eta;
+ }
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 - (cs - csh - kx*snh)/4)
+ zy[i][j]=
+ (v*csh + (v+tau*eta)*kx*snh - eta*rho*(cs - csh - kx*snh)/2)
*std::sin(pi*y);
}
break;
diff -r 5a513d0f7ec3 -r d933ef1bd792 wscript
--- a/wscript Sun Dec 04 00:39:47 2011 -0800
+++ b/wscript Sun Dec 04 12:00:01 2011 -0800
@@ -15,5 +15,5 @@ def build(bld):
target = 'Gamr_disc',
# cxxflags = ['-std=c++0x','-g','-D_GLIBCXX_DEBUG'],
cxxflags = ['-std=c++0x','-O3'],
- lib = ['boost_filesystem', 'boost_system'],
+ lib = ['boost_filesystem', 'boost_system', 'gsl', 'cblas'],
)
More information about the CIG-COMMITS
mailing list