[cig-commits] commit: Have two ghost points for the spline interpolation
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:33 PST 2012
changeset: 33:39d15ef50cb2
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 08 01:18:14 2012 -0800
files: compute_v_on_interface.cxx main.cxx
description:
Have two ghost points for the spline interpolation
diff -r ea7929ba3ac3 -r 39d15ef50cb2 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Tue Feb 07 12:59:32 2012 -0800
+++ b/compute_v_on_interface.cxx Wed Feb 08 01:18:14 2012 -0800
@@ -21,12 +21,12 @@ double interpolate_v(T z, T log_eta, con
void compute_v_on_interface(double zx[N+1][N], double zy[N][N+1],
double log_etax[N+1][N], double log_etay[N][N+1],
- double vx[N+2], double dvx_y[N+1], double dvx_yy[N],
+ double vx[N+4], double dvx_y[N+1], double dvx_yy[N],
double vy[N+1], double dvy_y[N], double dvy_yy[N+1])
{
- double y_x[N+2], y_y[N+1];
- for(int i=0;i<N+2;++i)
- y_x[i]=h*(i-0.5);
+ double y_x[N+4], y_y[N+1];
+ for(int i=0;i<N+4;++i)
+ y_x[i]=h*(i-1.5);
for(int i=0;i<N+1;++i)
y_y[i]=h*i;
@@ -50,7 +50,7 @@ void compute_v_on_interface(double zx[N+
double low=interpolate_v(zx,log_etax,i-1,j,h,(1+dx)*h);
double high=interpolate_v(zx,log_etax,i+2,j,h,(-2+dx)*h);
- vx[j+1]=low*(1-dx) + high*dx;
+ vx[j+2]=low*(1-dx) + high*dx;
ux1_p[j]=zx[i+1][j]*exp(-log_etax[i+1][j])*(1-dx)
+ zx[i+2][j]*exp(-log_etax[i+2][j])*dx;
@@ -62,8 +62,15 @@ void compute_v_on_interface(double zx[N+
ux2_m[j]=zx[i-2][j]*exp(-log_etax[i-2][j])*(1-dx)
+ zx[i-1][j]*exp(-log_etax[i-1][j])*dx;
}
- vx[0]=vx[1];
- vx[N+1]=vx[N];
+ // vx[0]=vx[3];
+ // vx[N+3]=vx[N];
+ // vx[1]=vx[2];
+ // vx[N+2]=vx[N+1];
+
+ vx[0]=-vx[N];
+ vx[N+3]=-vx[3];
+ vx[1]=-vx[N+1];
+ vx[N+2]=-vx[2];
}
{
int i(middle/h-0.5);
@@ -107,8 +114,8 @@ void compute_v_on_interface(double zx[N+
/* Compute the derivatives from the current guess */
/* dvx_dy, dvx_dyy */
gsl_interp_accel *vx_accel=gsl_interp_accel_alloc ();
- gsl_spline *vx_spline=gsl_spline_alloc(gsl_interp_cspline, N+2);
- gsl_spline_init(vx_spline, y_x, vx, N+2);
+ gsl_spline *vx_spline=gsl_spline_alloc(gsl_interp_cspline, N+4);
+ gsl_spline_init(vx_spline, y_x, vx, N+4);
dvx_y[0]=0;
dvx_y[N]=0;
@@ -116,7 +123,7 @@ void compute_v_on_interface(double zx[N+
dvx_y[i]=gsl_spline_eval_deriv(vx_spline,y_y[i],vx_accel);
for(int i=0;i<N;++i)
- dvx_yy[i]=gsl_spline_eval_deriv2(vx_spline,y_x[i+1],vx_accel);
+ dvx_yy[i]=gsl_spline_eval_deriv2(vx_spline,y_x[i+2],vx_accel);
gsl_spline_free(vx_spline);
gsl_interp_accel_free(vx_accel);
@@ -127,7 +134,7 @@ void compute_v_on_interface(double zx[N+
gsl_spline_init(vy_spline, y_y, vy, N+1);
for(int i=0;i<N;++i)
- dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i+1],vy_accel);
+ dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i+2],vy_accel);
for(int i=0;i<N;++i)
dvy_yy[i]=gsl_spline_eval_deriv2(vy_spline,y_y[i],vy_accel);
@@ -145,24 +152,31 @@ void compute_v_on_interface(double zx[N+
/(min_eta*a1*a2*(b1+b2) + max_eta*b1*b2*(a1+a2))
- a1*a2*b1*b2*zx_jump/(min_eta*a1*a2*(b1+b2)
+ max_eta*b1*b2*(a1+a2));
- if(vx[i+1]!=0 || temp!=0)
+ if(vx[i+2]!=0 || temp!=0)
error=
- std::max(error,std::fabs((vx[i+1]-temp)
- /(std::fabs(vx[i+1])+std::fabs(temp))));
+ std::max(error,std::fabs((vx[i+2]-temp)
+ /(std::fabs(vx[i+2])+std::fabs(temp))));
if(i<7)
std::cout << "vx error "
<< i << " "
- << vx[i+1] << " "
+ << vx[i+2] << " "
<< dvx_y[i] << " "
<< dvx_yy[i] << " "
<< temp << " "
<< error << " "
<< "\n";
- vx[i+1]=(1-relax)*vx[i+1] + relax*temp;
+ vx[i+2]=(1-relax)*vx[i+2] + relax*temp;
}
- vx[0]=vx[1];
- vx[N+1]=vx[N];
+ // vx[0]=vx[3];
+ // vx[N+3]=vx[N];
+ // vx[1]=vx[2];
+ // vx[N+2]=vx[N+1];
+
+ vx[0]=-vx[N];
+ vx[N+3]=-vx[3];
+ vx[1]=-vx[N+1];
+ vx[N+2]=-vx[2];
for(int i=0;i<N+1;++i)
{
diff -r ea7929ba3ac3 -r 39d15ef50cb2 main.cxx
--- a/main.cxx Tue Feb 07 12:59:32 2012 -0800
+++ b/main.cxx Wed Feb 08 01:18:14 2012 -0800
@@ -13,7 +13,7 @@ extern void compute_v_on_interface(doubl
extern void compute_v_on_interface(double zx[N+1][N], double zy[N][N+1],
double log_etax[N+1][N],
double log_etay[N][N+1],
- double vx[N+2], double dvx_y[N+1],
+ double vx[N+4], double dvx_y[N+1],
double dvx_yy[N],
double vy[N+1], double dvy_y[N],
double dvy_yy[N+1]);
@@ -50,7 +50,7 @@ int main()
{
/* Compute the boundary jumps */
- double vx[N+2], dvx_y[N+1], dvx_yy[N],
+ double vx[N+4], dvx_y[N+1], dvx_yy[N],
vy[N+1], dvy_y[N], dvy_yy[N+1];
if(model==Model::solCx)
@@ -160,7 +160,7 @@ int main()
if(x-h<middle && x+h>middle)
{
double dx=(x>middle) ? (middle-(x-h)) : (middle-(x+h));
- double zx_jump=eta_jump*vx[j+1];
+ double zx_jump=eta_jump*vx[j+2];
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;
@@ -450,7 +450,7 @@ int main()
if(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+1];
+ double zx_jump=vx[j+2];
double dzx_x_jump=dvy_y[j];
dzx_x_correction=eta_jump*(zx_jump + dx*dzx_x_jump)/h;
More information about the CIG-COMMITS
mailing list