[cig-commits] commit: Use periodic spline interpolation
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:38 PST 2012
changeset: 37:c408f265d389
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 08 03:46:42 2012 -0800
files: compute_v_on_interface.cxx constants.hxx main.cxx
description:
Use periodic spline interpolation
diff -r 127bfe48757b -r c408f265d389 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Wed Feb 08 03:27:21 2012 -0800
+++ b/compute_v_on_interface.cxx Wed Feb 08 03:46:42 2012 -0800
@@ -21,14 +21,14 @@ double interpolate_v(T z, T log_eta, con
void compute_v_on_interface(double zx[Nx+1][Ny], double zy[Nx][Ny+1],
double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
- double vx[Ny+4], double dvx_y[Ny+1],
+ double vx[Ny+1], double dvx_y[Ny+1],
double dvx_yy[Ny],
double vy[Ny+1], double dvy_y[Ny],
double dvy_yy[Ny+1])
{
- double y_x[Ny+4], y_y[Ny+1];
- for(int i=0;i<Ny+4;++i)
- y_x[i]=h*(i-1.5);
+ double y_x[Ny+1], y_y[Ny+1];
+ for(int i=0;i<Ny+1;++i)
+ y_x[i]=h*(i-0.5);
for(int i=0;i<Ny+1;++i)
y_y[i]=h*i;
@@ -52,7 +52,7 @@ void compute_v_on_interface(double zx[Nx
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+2]=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;
@@ -64,15 +64,7 @@ void compute_v_on_interface(double zx[Nx
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[3];
- // vx[Ny+3]=vx[Ny];
- // vx[1]=vx[2];
- // vx[Ny+2]=vx[Ny+1];
-
- vx[0]=-vx[Ny];
- vx[Ny+3]=-vx[3];
- vx[1]=-vx[Ny+1];
- vx[Ny+2]=-vx[2];
+ vx[0]=vx[Ny];
}
{
int i(middle/h-0.5);
@@ -116,8 +108,8 @@ void compute_v_on_interface(double zx[Nx
/* 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, Ny+4);
- gsl_spline_init(vx_spline, y_x, vx, Ny+4);
+ gsl_spline *vx_spline=gsl_spline_alloc(gsl_interp_cspline_periodic, Ny+1);
+ gsl_spline_init(vx_spline, y_x, vx, Ny+1);
dvx_y[0]=0;
dvx_y[Ny]=0;
@@ -125,28 +117,24 @@ void compute_v_on_interface(double zx[Nx
dvx_y[i]=gsl_spline_eval_deriv(vx_spline,y_y[i],vx_accel);
for(int i=0;i<Ny;++i)
- dvx_yy[i]=gsl_spline_eval_deriv2(vx_spline,y_x[i+2],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);
/* dvy_dy, dvy_dyy */
gsl_interp_accel *vy_accel=gsl_interp_accel_alloc ();
- gsl_spline *vy_spline=gsl_spline_alloc(gsl_interp_cspline, Ny+1);
+ gsl_spline *vy_spline=gsl_spline_alloc(gsl_interp_cspline_periodic, Ny+1);
gsl_spline_init(vy_spline, y_y, vy, Ny+1);
for(int i=0;i<Ny;++i)
- dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i+2],vy_accel);
+ dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i+1],vy_accel);
for(int i=0;i<Ny;++i)
dvy_yy[i]=gsl_spline_eval_deriv2(vy_spline,y_y[i],vy_accel);
gsl_spline_free(vy_spline);
gsl_interp_accel_free(vy_accel);
-
- // std::cout << "vx "
- // << vx[Ny/4+2] << " "
- // << "\n";
/* Update the current guess based on the derivatives */
for(int i=0;i<Ny;++i)
@@ -158,21 +146,13 @@ void compute_v_on_interface(double zx[Nx
/(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+2]!=0 || temp!=0)
+ if(vx[i+1]!=0 || temp!=0)
error=
- std::max(error,std::fabs((vx[i+2]-temp)
- /(std::fabs(vx[i+2])+std::fabs(temp))));
- vx[i+2]=(1-relax)*vx[i+2] + relax*temp;
+ std::max(error,std::fabs((vx[i+1]-temp)
+ /(std::fabs(vx[i+1])+std::fabs(temp))));
+ vx[i+1]=(1-relax)*vx[i+1] + relax*temp;
}
- // vx[0]=vx[3];
- // vx[Ny+3]=vx[Ny];
- // vx[1]=vx[2];
- // vx[Ny+2]=vx[Ny+1];
-
- vx[0]=-vx[Ny];
- vx[Ny+3]=-vx[3];
- vx[1]=-vx[Ny+1];
- vx[Ny+2]=-vx[2];
+ vx[0]=vx[Ny];
for(int i=0;i<Ny+1;++i)
{
diff -r 127bfe48757b -r c408f265d389 constants.hxx
--- a/constants.hxx Wed Feb 08 03:27:21 2012 -0800
+++ b/constants.hxx Wed Feb 08 03:46:42 2012 -0800
@@ -1,6 +1,6 @@ const int NN(1024);
-const int NN(1024);
+const int NN(256);
const int Nx(NN);
-const int Ny(NN);
+const int Ny(2*NN);
const double min_eta=1;
const double max_eta=1e6;
const double eta_jump=max_eta-min_eta;
diff -r 127bfe48757b -r c408f265d389 main.cxx
--- a/main.cxx Wed Feb 08 03:27:21 2012 -0800
+++ b/main.cxx Wed Feb 08 03:46:42 2012 -0800
@@ -13,7 +13,7 @@ extern void compute_v_on_interface(doubl
extern void compute_v_on_interface(double zx[Nx+1][Ny], double zy[Nx][Ny+1],
double log_etax[Nx+1][Ny],
double log_etay[Nx][Ny+1],
- double vx[Ny+4], double dvx_y[Ny+1],
+ double vx[Ny+1], double dvx_y[Ny+1],
double dvx_yy[Ny],
double vy[Ny+1], double dvy_y[Ny],
double dvy_yy[Ny+1]);
@@ -50,7 +50,7 @@ int main()
{
/* Compute the boundary jumps */
- double vx[Ny+4], dvx_y[Ny+1], dvx_yy[Ny],
+ double vx[Ny+1], dvx_y[Ny+1], dvx_yy[Ny],
vy[Ny+1], dvy_y[Ny], dvy_yy[Ny+1];
if(model==Model::solCx)
@@ -150,7 +150,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+2];
+ 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 dp_x_jump=2*dvx_yy[j]*eta_jump;
@@ -163,7 +163,7 @@ int main()
dzx_xx+=dzx_xx_correction;
// dzx_xx-=eta_jump*(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
- if(j==Ny/4)
+ if(j==Ny/8)
std::cout << "dzx_xx "
<< i << " "
<< j << " "
@@ -374,7 +374,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+2];
+ 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