[cig-commits] commit: Add a ghost region for vx on the interface to enforce the zero derivative
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:31 PST 2012
changeset: 31:0fd4b7aff763
user: Walter Landry <wlandry at caltech.edu>
date: Tue Feb 07 12:36:14 2012 -0800
files: compute_v_on_interface.cxx main.cxx
description:
Add a ghost region for vx on the interface to enforce the zero derivative
diff -r 2d4847e70d29 -r 0fd4b7aff763 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Tue Feb 07 11:30:31 2012 -0800
+++ b/compute_v_on_interface.cxx Tue Feb 07 12:36: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], double dvx_y[N+1], double dvx_yy[N],
+ double vx[N+2], 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], y_y[N+1];
- for(int i=0;i<N;++i)
- y_x[i]=h*(i+0.5);
+ double y_x[N+2], y_y[N+1];
+ for(int i=0;i<N+2;++i)
+ y_x[i]=h*(i-0.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]=low*(1-dx) + high*dx;
+ vx[j+1]=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,6 +62,8 @@ 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];
}
{
int i(middle/h-0.5);
@@ -105,8 +107,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);
- gsl_spline_init(vx_spline, y_x, vx, N);
+ gsl_spline *vx_spline=gsl_spline_alloc(gsl_interp_cspline, N+2);
+ gsl_spline_init(vx_spline, y_x, vx, N+2);
dvx_y[0]=0;
dvx_y[N]=0;
@@ -114,7 +116,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],vx_accel);
+ dvx_yy[i]=gsl_spline_eval_deriv2(vx_spline,y_x[i+1],vx_accel);
gsl_spline_free(vx_spline);
gsl_interp_accel_free(vx_accel);
@@ -125,7 +127,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],vy_accel);
+ dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i+1],vy_accel);
for(int i=0;i<N;++i)
dvy_yy[i]=gsl_spline_eval_deriv2(vy_spline,y_y[i],vy_accel);
@@ -143,21 +145,21 @@ 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]!=0 || temp!=0)
+ if(vx[i+1]!=0 || temp!=0)
error=
- std::max(error,std::fabs((vx[i]-temp)
- /(std::fabs(vx[i])+std::fabs(temp))));
+ std::max(error,std::fabs((vx[i+1]-temp)
+ /(std::fabs(vx[i+1])+std::fabs(temp))));
if(i<7)
std::cout << "vx error "
<< i << " "
- << vx[i] << " "
+ << vx[i+1] << " "
<< dvx_y[i] << " "
<< dvx_yy[i] << " "
<< temp << " "
<< error << " "
<< "\n";
- vx[i]=(1-relax)*vx[i] + relax*temp;
+ vx[i+1]=(1-relax)*vx[i+1] + relax*temp;
}
for(int i=0;i<N+1;++i)
diff -r 2d4847e70d29 -r 0fd4b7aff763 main.cxx
--- a/main.cxx Tue Feb 07 11:30:31 2012 -0800
+++ b/main.cxx Tue Feb 07 12:36: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], double dvx_y[N+1],
+ double vx[N+2], 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], dvx_y[N+1], dvx_yy[N],
+ double vx[N+2], 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];
+ double zx_jump=eta_jump*vx[j+1];
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;
@@ -250,7 +250,7 @@ int main()
dp_x_correction=(p_jump + dx*dp_x_jump)/h;
dp_x-=(p_jump + dx*dp_x_jump)/h;
- if(j==N/4)
+ if(j<3)
std::cout << "p "
<< i << " "
<< j << " "
@@ -290,6 +290,8 @@ int main()
Resid_x[i][j]=Rx;
// zx[i][j]-=theta_mom*Rx/dRx_dzx;
+
+ Resid_x[i][j]=dp_x;
}
}
}
@@ -448,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];
+ double zx_jump=vx[j+1];
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