[cig-commits] commit: An in-progress version that uses cubic interpolation onto the interface. Does not work because the derivatives are not continuous at the interface.
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:22 PST 2012
changeset: 23:567f1fbb28be
user: Walter Landry <wlandry at caltech.edu>
date: Fri Feb 03 11:32:24 2012 -0800
files: constants.hxx initial.cxx main.cxx wscript
description:
An in-progress version that uses cubic interpolation onto the interface. Does not work because the derivatives are not continuous at the interface.
diff -r abd75122918b -r 567f1fbb28be constants.hxx
--- a/constants.hxx Sun Dec 04 16:11:05 2011 -0800
+++ b/constants.hxx Fri Feb 03 11:32:24 2012 -0800
@@ -1,11 +1,12 @@ const int N(64);
-const int N(64);
-const double min_eta=.0001;
-const double max_eta=1000;
+const int N(1024);
+const double min_eta=1;
+const double max_eta=2;
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=25.0/64 + 0.0000000001;
+const double middle=400.00000001/1024;
const double h(1.0/N);
const double pi(atan(1.0)*4);
diff -r abd75122918b -r 567f1fbb28be initial.cxx
--- a/initial.cxx Sun Dec 04 16:11:05 2011 -0800
+++ b/initial.cxx Fri Feb 03 11:32:24 2012 -0800
@@ -25,7 +25,8 @@ 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*eta0, (cosh0-kx0*sinh0)*eta0, kx1*sinh1*eta1, -(cosh1-kx1*sinh1)*eta1,
+ -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};
@@ -35,30 +36,36 @@ void initial(const Model &model, double
(kx0*sinh0 + kx1*sinh1)*rho/2,
((sinh0 + kx0*cosh0) + (sinh1 + kx1*cosh1))*rho/2};
+ // double m_temp[]={m[0],m[1],m[2],m[3],m[4],m[5],m[6],m[7],m[8],m[9],m[10],m[11],m[12],m[13],m[14],m[15],rhs[3]};
+
+ // std::cout << "rhs "
+ // << rhs[0] << " "
+ // << rhs[1] << " "
+ // << rhs[2] << " "
+ // << rhs[3] << " "
+ // << (sinh0 + kx0*cosh0)*rho/2 << " "
+ // << (sinh1 + kx1*cosh1)*rho/2 << " "
+ // << "\n";
+
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);
-
- 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);
gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,x);
- gsl_permutation_free(perm);
-
- std::cout << "soln\n";
- gsl_vector_fprintf(stdout,x,"%g ");
-
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);
+ gsl_permutation_free(perm);
+
+ // std::cout << "v tau "
+ // << v0 << " "
+ // << tau_eta_0 << " "
+ // << v1 << " "
+ // << tau_eta_1 << " "
+ // << "\n";
for(int i=0;i<N+1;++i)
for(int j=0;j<N+1;++j)
@@ -97,16 +104,13 @@ 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;
+ double tau_xx=2*pi*(-v*eta*kx*snh + tau_eta*eta*(csh-kx*snh)
+ - sign*kx*snh*rho/2);
- // p[i][j]=v*(-snh-kx*csh)*eta - tau_eta*kx*csh*eta
- // -sign*(snh+kx*csh)*rho/2;
+ double dvx_x=(-v*(csh + kx*snh) + tau_eta*(-kx*snh)
+ - sign*(kx*snh + csh - cs)*rho/(2*eta));
- p[i][j]=-2*pi*((v + tau_eta)*eta*csh + sign*(csh - cs)*rho/2)
- *std::cos(pi*y);
+ p[i][j]=(2*pi*eta*dvx_x - tau_xx)*std::cos(pi*y);
}
break;
default:
@@ -156,6 +160,39 @@ void initial(const Model &model, double
zx[i][j]=(tau_eta*snh - (v+tau_eta)*kx*csh
- sign*rho*(kx*csh/2 - sn/2)/eta)
*std::cos(pi*y)*eta;
+
+
+ double tau_xz=std::sin(pi*y)*(-v*eta*(snh + kx*csh) - tau_eta*eta*kx*csh
+ - sign*(snh + kx*csh)*rho/2);
+
+ double dvx_z=pi*std::sin(pi*y)*(eta*(v*kx*csh + tau_eta*(-snh + kx*csh))
+ + sign*(kx*csh - sn)*rho/2);
+
+ double dvz_x=pi*std::sin(pi*y)*(eta*(v*(2*snh + kx*csh) + tau_eta*(snh + kx*csh))
+ + sign*(sn + 2*snh + kx*csh)*rho/2);
+
+ double zy=(v*csh + (v+tau_eta)*kx*snh
+ - sign*rho*(cs - csh - kx*snh)/(2*eta))
+ *std::sin(pi*y)*eta;
+
+ if(j==N/4 && i>middle*N-5 && i<middle*N+5)
+ std::cout << "tau_xx "
+ << i << " "
+ << j << " "
+ << x << " "
+ << y << " "
+ // << tau_xx << " "
+ // << dvx_z << " "
+ // << dvz_x << " "
+ // << dvx_z + dvz_x << " "
+ // << tau_xz << " "
+ // << zx[i][j] << " "
+ << zx[i][j]*2*pi*pi*eta_jump/eta << " "
+ // << 2*zy*pi*pi*eta_jump/eta << " "
+ << eta_jump*2*zy*pi/(tan(pi*y)*eta) << " "
+ // << p[i][j] << " "
+ << zy/eta << " "
+ << "\n";
}
break;
case Model::sinker:
diff -r abd75122918b -r 567f1fbb28be main.cxx
--- a/main.cxx Sun Dec 04 16:11:05 2011 -0800
+++ b/main.cxx Fri Feb 03 11:32:24 2012 -0800
@@ -5,6 +5,7 @@
#include "Model.hxx"
#include "constants.hxx"
#include "write_vtk.hxx"
+#include <gsl/gsl_spline.h>
extern 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],
@@ -47,32 +48,82 @@ int main()
if(model==Model::solCx)
{
- int i(middle/h);
- double dx=middle/h-i;
+ double xa[]={0,h,2*h,3*h};
+ double ya[]={12,15,18,22};
+ gsl_interp_accel *accel=gsl_interp_accel_alloc();
+ gsl_spline *spline=gsl_spline_alloc(gsl_interp_cspline,4);
+ {
+ int i(middle/h);
+ double dx=middle-h*i;
- /* vx, dvx_y, dvx_yy */
- for(int j=0;j<N;++j)
- vx[j]=zx[i][j]*(1-dx)*exp(-log_etax[i][j])
- + zx[i+1][j]*dx*exp(-log_etax[i+1][j]);
+ /* vx, dvx_y, dvx_yy */
+ for(int j=0;j<N;++j)
+ {
+ ya[0]=zx[i-1][j]*exp(-log_etax[i-1][j]);
+ ya[1]=zx[i][j]*exp(-log_etax[i][j]);
+ ya[2]=zx[i+1][j]*exp(-log_etax[i+1][j]);
+ ya[3]=zx[i+2][j]*exp(-log_etax[i+2][j]);
- dvx_y[0]=dvx_y[N]=0;
- for(int j=1;j<N+1;++j)
- dvx_y[j]=(vx[j]-vx[j-1])/h;
+ gsl_spline_init(spline,xa,ya,4);
+ vx[j]=gsl_spline_eval(spline,h+dx,accel);
- for(int j=0;j<N;++j)
- dvx_yy[j]=(dvx_y[j+1]-dvx_y[j])/h;
+ // vx[j]=zx[i][j]*(1-dx)*exp(-log_etax[i][j])
+ // + zx[i+1][j]*dx*exp(-log_etax[i+1][j]);
+ }
- /* vy, dvy_y, dvy_yy */
- for(int j=0;j<N+1;++j)
- vy[j]=zy[i][j]*(1-dx)*exp(-log_etay[i][j])
- + zy[i+1][j]*dx*exp(-log_etay[i+1][j]);
+ dvx_y[0]=dvx_y[N]=0;
+ for(int j=1;j<N;++j)
+ dvx_y[j]=(vx[j]-vx[j-1])/h;
- for(int j=0;j<N;++j)
- dvy_y[j]=(vy[j+1]-vy[j])/h;
+ for(int j=0;j<N;++j)
+ dvx_yy[j]=(dvx_y[j+1]-dvx_y[j])/h;
+ }
+ {
+ int i(middle/h-0.5);
+ double dx=h*(middle/h-0.5-i);
- dvy_yy[0]=dvy_yy[N]=0;
- for(int j=1;j<N;++j)
- dvy_yy[j]=(dvy_y[j+1]-dvy_y[j])/h;
+ // const gsl_interp_type *t=gsl_interp_polynomial;
+ // gsl_interp* interp=gsl_interp_alloc(t,4);
+
+ /* vy, dvy_y, dvy_yy */
+ for(int j=0;j<N+1;++j)
+ {
+ ya[0]=zy[i-1][j]*exp(-log_etay[i-1][j]);
+ ya[1]=zy[i][j]*exp(-log_etay[i][j]);
+ ya[2]=zy[i+1][j]*exp(-log_etay[i+1][j]);
+ ya[3]=zy[i+2][j]*exp(-log_etay[i+2][j]);
+
+ gsl_spline_init(spline,xa,ya,4);
+ vy[j]=gsl_spline_eval(spline,h+dx,accel);
+
+ // vy[j]=zy[i][j]*(1-dx)*exp(-log_etay[i][j])
+ // + zy[i+1][j]*dx*exp(-log_etay[i+1][j]);
+
+ if(j==N/4)
+ std::cout << "vy "
+ << j << " "
+ << xa[0] << " "
+ << xa[1] << " "
+ << xa[2] << " "
+ << xa[3] << " "
+ << ya[0] << " "
+ << ya[1] << " "
+ << ya[2] << " "
+ << ya[3] << " "
+ << vy[j] << " "
+ << "\n";
+ }
+ // gsl_interp_free(interp);
+
+ for(int j=0;j<N;++j)
+ dvy_y[j]=(vy[j+1]-vy[j])/h;
+
+ dvy_yy[0]=dvy_yy[N]=0;
+ for(int j=1;j<N;++j)
+ dvy_yy[j]=(dvy_y[j+1]-dvy_y[j])/h;
+ }
+ gsl_spline_free(spline);
+ gsl_interp_accel_free(accel);
}
/* Smoothing sweeps */
@@ -161,32 +212,191 @@ int main()
dlog_etaxx=dlog_etax=dlog_etayy=dlog_etay=
dlog_etaxy=0;
+ // if(j==N/4)
+ // std::cout << "dzx_xx "
+ // << i << " "
+ // << j << " "
+ // << x << " "
+ // << y << " "
+ // // << zx[i][j] << " "
+ // << zy[i][j] << " "
+ // // << log_etax[i][j] << " "
+ // // << log_etay[i][j] << " "
+ // << "\n";
+
if(x-h<middle && x+h>middle)
{
- double dx=(x>middle) ? (x+h-middle) : (middle-(x-h));
- double zx_jump=vx[j];
- double dzx_x_jump=dvy_y[j];
- double p_jump=-2*dvy_y[j];
- double dzx_xx_jump=-dvx_yy[j] + p_jump;
- dzx_xx_correction=eta_jump*(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
+ double dx=(x>middle) ? (middle-(x-h)) : (middle-(x+h));
+ double zx_jump=eta_jump*vx[j];
+ double dzx_x_jump=eta_jump*dvy_y[j];
+ double p_jump=-2*dvy_y[j]*eta_jump;
+ double dzx_xx_jump=-dvx_yy[j]*eta_jump + p_jump;
+ dzx_xx_correction=-(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
+ if(x>middle)
+ dzx_xx_correction*=-1;
+
// dzx_xx-=eta_jump*(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
- // if(n_sweeps==1)
- // std::cout << "dzx_xx "
- // << i << " "
- // << j << " "
- // << dzx_xx << " "
- // << dzx_xx_correction << " "
- // << eta_jump*zx_jump << " "
- // << zx[i][j] << " "
- // << zx[i+1][j] << " "
- // << zx[i-1][j] << " "
- // << "\n";
- if(dx>h/2)
+
+ if(j==N/4)
+ std::cout << "dzx_xx "
+ << i << " "
+ << j << " "
+ << x << " "
+ << y << " "
+ // << zx[i+2][j]*exp(-log_etax[i+2][j]) << " "
+ // << zx[i+1][j]*exp(-log_etax[i+1][j]) << " "
+ // << zx[i][j]*exp(-log_etax[i][j]) << " "
+ // << zx[i-1][j]*exp(-log_etax[i-1][j]) << " "
+ // << zx[i-2][j]*exp(-log_etax[i-2][j]) << " "
+ // << zx[i-3][j]*exp(-log_etax[i-3][j]) << " "
+ // << zy[i+2][j]*exp(-log_etay[i+2][j]) << " "
+ // << zy[i+1][j]*exp(-log_etay[i+1][j]) << " "
+ // << zy[i][j]*exp(-log_etay[i][j]) << " "
+ // << zy[i-1][j]*exp(-log_etay[i-1][j]) << " "
+ // << p[i+1][j] << " "
+ // << p[i][j] << " "
+ // << p[i-1][j] << " "
+ // << p[i-2][j] << " "
+
+ // << 2*(zx[i+2][j]-zx[i+1][j])/h - p[i+1][j] << " "
+ // << 2*(zx[i+1][j]-zx[i][j])/h - p[i][j] << " "
+ // << 2*(zx[i][j]-zx[i-1][j])/h - p[i-1][j] << " "
+ // << 2*(zx[i-1][j]-zx[i-2][j])/h - p[i-2][j] << " "
+
+ // << (zx[i+2][j] - zx[i+1][j])/h + (zy[i+1][j+1] - zy[i+1][j])/h << " "
+ // << (zx[i+1][j] - zx[i][j])/h + (zy[i][j+1] - zy[i][j])/h << " "
+ // << (zx[i][j] - zx[i-1][j])/h + (zy[i-1][j+1] - zy[i-1][j])/h << " "
+ // << (zx[i-1][j] - zx[i-2][j])/h + (zy[i-2][j+1] - zy[i-2][j])/h << " "
+
+ // << (zy[i+2][j+1]-zy[i+1][j+1])/h + (zx[i+2][j+1]-zx[i+2][j])/h << " "
+ // << (zy[i+1][j+1]-zy[i ][j+1])/h + (zx[i+1][j+1]-zx[i+1][j])/h << " "
+ // << (zy[i ][j+1]-zy[i-1][j+1])/h + (zx[i ][j+1]-zx[i ][j])/h << " "
+ // << (zy[i-1][j+1]-zy[i-2][j+1])/h + (zx[i-1][j+1]-zx[i-1][j])/h << " "
+ // << (zy[i-2][j+1]-zy[i-3][j+1])/h + (zx[i-2][j+1]-zx[i-2][j])/h << " "
+
+
+ // << (zy[i+2][j]-zy[i+1][j])/h << " "
+ // << (zx[i+2][j+1]-zx[i+2][j])/h << " "
+ // << (zy[i+1][j]-zy[i][j])/h << " "
+ // << (zx[i+1][j+1]-zx[i+1][j])/h << " "
+ // << (zy[i][j]-zy[i-1][j])/h << " "
+ // << (zx[i][j+1]-zx[i][j])/h << " "
+ // << (zy[i-1][j]-zy[i-2][j])/h << " "
+ // << (zx[i-1][j+1]-zx[i-1][j])/h << " "
+ // << (zy[i-2][j]-zy[i-3][j])/h << " "
+ // << (zx[i-2][j+1]-zx[i-2][j])/h << " "
+
+ // << zx[i+2][j] << " "
+ // << zx[i+1][j] << " "
+ // << zx[i][j] << " "
+ // << zx[i-1][j] << " "
+ // << zx[i-2][j] << " "
+ // << zy[i+2][j] << " "
+ // << zy[i+1][j] << " "
+ // << zy[i][j] << " "
+ // << zy[i-1][j] << " "
+ // << zy[i-2][j] << " "
+
+ // << dx << " "
+ // << dzx_xx << " "
+ // << dzx_xx_correction << " "
+ // << dzx_xx + dzx_xx_correction << " "
+ // << eta_jump*zx_jump/(h*h) << " "
+ // << dx*eta_jump*dzx_x_jump/(h*h) << " "
+ // << dx*dx*dzx_xx_jump/(2*h*h) << " "
+ // << (zx[i+2][j]-2*zx[i+1][j]+zx[i][j])/(h*h) << " "
+ // << (zx[i+3][j]-2*zx[i+2][j]+zx[i+1][j])/(h*h) << " "
+ // << (zx[i+4][j]-2*zx[i+3][j]+zx[i+2][j])/(h*h) << " "
+ // << (zx[i-3][j]-2*zx[i-2][j]+zx[i-1][j])/(h*h) << " "
+ << "\n";
+
+ if(x+h/2>middle && x-h/2<middle)
+ // if(dx>h/2 || dx<-h/2)
{
- dx-=h/2;
- double dp_x_jump=2*dvx_yy[j];
- dp_x_correction=eta_jump*(p_jump + dx*dp_x_jump)/h;
+ dx=(x>middle) ? (middle-(x-h/2)) : ((x+h/2)-middle);
+ double dp_x_jump=2*dvx_yy[j]*eta_jump;
+ dp_x_correction=(p_jump + dx*dp_x_jump)/h;
// dp_x-=eta_jump*(p_jump + dx*dp_x_jump)/h;
+
+ double dp_fixed=
+ ((dp_x - p_jump/h) - dp_x_jump*dx/h);
+
+ double dp_numerical=
+ (p[i-1][j]-p[i-2][j])/h;
+
+ if(j==N/4)
+ std::cout << "p "
+ << i << " "
+ << j << " "
+
+ << p_jump << " "
+ << dp_x_jump << " "
+
+ << dp_x << " "
+ << (dp_x - p_jump/h) << " "
+
+ << dp_fixed << " "
+ << (p[i+3][j]-p[i+2][j])/h << " "
+ << (p[i+2][j]-p[i+1][j])/h << " "
+ << (p[i-1][j]-p[i-2][j])/h << " "
+ << (p[i-2][j]-p[i-3][j])/h << " "
+ // << dp_numerical << " "
+ << std::fabs((dp_fixed-dp_numerical)/dp_numerical) << " "
+
+
+ // // << x << " "
+ // // << dx << " "
+ // // << h << " "
+ // // << dx/h << " "
+ // << (dp_x - eta_jump*p_jump/h)*2
+ // + dp_x_jump*2*dx/(h) << " "
+ // // << dp_x << " "
+ // // << dp_x - eta_jump*p_jump/h << " "
+ // // << dp_x_correction << " "
+ // // << dp_x - dp_x_correction << " "
+ // // << p_jump << " "
+ // // << p[i][j] - p[i-1][j] << " "
+
+ // // << p[i+2][j] << " "
+ // // << p[i+1][j] << " "
+ // // << p[i][j] << " "
+ // // << p[i-1][j] << " "
+ // // << p[i-2][j] << " "
+
+ // // << zy[i+1][j] << " "
+ // // << zy[i][j] << " "
+ // // << zy[i-1][j] << " "
+ // // << 2*(zy[i+2][j+1]-zy[i+2][j])/h << " "
+ // // << 2*(zy[i+1][j+1]-zy[i+1][j])/h << " "
+ // // << 2*(zy[i][j+1]-zy[i][j])/h << " "
+ // // << 2*(zy[i-1][j+1]-zy[i-1][j])/h << " "
+
+ // // << 2*(zy[i+2][j+1]*exp(-log_etay[i+2][j+1])-zy[i+2][j]*exp(-log_etay[i+2][j]))/h << " "
+ // // << 2*(zy[i+1][j+1]*exp(-log_etay[i+1][j+1])-zy[i+1][j]*exp(-log_etay[i+1][j]))/h << " "
+ // // << 2*(zy[i][j+1]*exp(-log_etay[i][j+1])-zy[i][j]*exp(-log_etay[i][j]))/h << " "
+ // // << 2*(zy[i-1][j+1]*exp(-log_etay[i-1][j+1])-zy[i-1][j]*exp(-log_etay[i-1][j]))/h << " "
+
+ // // << dp_x_jump << " "
+ // // << ((p[i+1][j]-p[i][j])/h
+ // // - (p[i-1][j]-p[i-2][j])/h) << " "
+ // // << dx/h << " "
+ // // << (p[i+3][j]-p[i+2][j])/h << " "
+ // // << (p[i+2][j]-p[i+1][j])/h << " "
+ // // << (p[i+1][j]-p[i][j])/h << " "
+ // // << (p[i][j]-p[i-1][j])/h << " "
+ // << (p[i-1][j]-p[i-2][j])/h << " "
+ // << ((dp_x - eta_jump*p_jump/h)*2
+ // + dp_x_jump*2*dx/h)/((p[i-1][j]-p[i-2][j])/h) << " "
+ // // << (p[i-2][j]-p[i-3][j])/h << " "
+ // // << (p[i-3][j]-p[i-4][j])/h << " "
+ // // << (dp_x)
+ // // /((p[i-1][j]-p[i-2][j])/h) << " "
+ // // << (dp_x - eta_jump*(p_jump)/h)
+ // // /((p[i-1][j]-p[i-2][j])/h) << " "
+ // // << (dp_x - dp_x_correction)
+ // // /((p[i-1][j]-p[i-2][j])/h) << " "
+ << "\n";
+
dzy_xy_correction=eta_jump*h*(dvy_y[j] - dx*dvx_yy[j])/(h*h);
// dzy_xy-=eta_jump*h*(dvy_y[j] - dx*dvx_yy[j])/(h*h);
@@ -363,12 +573,31 @@ int main()
if(x+h/2>middle && x-h/2<middle)
{
- double dx=(x>middle) ? (x+h/2-middle) : (x-h/2-middle);
+ double dx=(x>middle) ? (middle-(x-h/2)) : (middle-(x+h/2));
double zx_jump=vx[j];
double dzx_x_jump=dvy_y[j];
dzx_x_correction=eta_jump*(zx_jump + dx*dzx_x_jump)/h;
// dzx_x-=eta_jump*(zx_jump + dx*dzx_x_jump)/h;
+
+ // if(j==N/4)
+ // std::cout << "dzx_x "
+ // << i << " "
+ // << j << " "
+ // << x << " "
+ // << dzx_x << " "
+ // << dzx_x_correction << " "
+ // << dzx_x-dzx_x_correction << " "
+ // << zx_jump/h << " "
+ // << dx*dzx_x_jump/h << " "
+ // << (zx[i+3][j]-zx[i+2][j])/h << " "
+ // << (zx[i+2][j]-zx[i+1][j])/h << " "
+ // << (zx[i+1][j]-zx[i][j])/h << " "
+ // << (zx[i][j]-zx[i-1][j])/h << " "
+ // << (zx[i-1][j]-zx[i-2][j])/h << " "
+ // << (zx[i-2][j]-zx[i-3][j])/h << " "
+ // << "\n";
+
}
}
@@ -451,11 +680,11 @@ int main()
if(sweep%100==0)
{
- std::cout << sweep << " "
- << max_x << " "
- << max_y << " "
- << max_p << " "
- << "\n";
+ // std::cout << sweep << " "
+ // << max_x << " "
+ // << max_y << " "
+ // << max_p << " "
+ // << "\n";
std::stringstream ss;
ss << "zx" << sweep;
diff -r abd75122918b -r 567f1fbb28be wscript
--- a/wscript Sun Dec 04 16:11:05 2011 -0800
+++ b/wscript Fri Feb 03 11:32:24 2012 -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', 'gsl', 'cblas'],
+ lib = ['boost_filesystem', 'boost_system', 'gsl', 'gslcblas'],
)
More information about the CIG-COMMITS
mailing list