[cig-commits] commit: Use cubic interpolation for computing v on the interface
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:24 PST 2012
changeset: 25:b28d5b89ebe7
user: Walter Landry <wlandry at caltech.edu>
date: Mon Feb 06 15:33:20 2012 -0800
files: main.cxx
description:
Use cubic interpolation for computing v on the interface
diff -r 81927ce18594 -r b28d5b89ebe7 main.cxx
--- a/main.cxx Fri Feb 03 13:40:17 2012 -0800
+++ b/main.cxx Mon Feb 06 15:33:20 2012 -0800
@@ -9,6 +9,21 @@ extern void initial(const Model &model,
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],
double p[N][N], double fx[N+1][N], double fy[N][N+1]);
+
+
+template<class T>
+double interpolate_v(T z, T log_eta,
+ const int &i, const int &j, const double &h, const double &x)
+{
+ double u0=z[i][j]*exp(-log_eta[i][j]);
+ double up=z[i+1][j]*exp(-log_eta[i+1][j]);
+ double um=z[i-1][j]*exp(-log_eta[i-1][j]);
+
+ double du=(up-um)/(2*h);
+ double ddu=(up - 2*u0 + um)/(h*h);
+
+ return u0 + x*du + x*x*ddu/2;
+}
int main()
{
@@ -54,11 +69,15 @@ int main()
/* vx, dvx_y, dvx_yy */
for(int j=0;j<N;++j)
{
- double low=(zx[i][j]*exp(-log_etax[i][j])*(dx+1)
+ double low1=(zx[i][j]*exp(-log_etax[i][j])*(dx+1)
- zx[i-1][j]*exp(-log_etax[i-1][j])*dx);
- double high=(zx[i+1][j]*exp(-log_etax[i+1][j])*(2-dx)
+ double high1=(zx[i+1][j]*exp(-log_etax[i+1][j])*(2-dx)
- zx[i+2][j]*exp(-log_etax[i+2][j])*(1-dx));
+
+
+ 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;
}
@@ -77,29 +96,37 @@ int main()
/* vy, dvy_y, dvy_yy */
for(int j=0;j<N+1;++j)
{
- double low=(zy[i][j]*exp(-log_etay[i][j])*(dx+1)
+ double low1=(zy[i][j]*exp(-log_etay[i][j])*(dx+1)
- zy[i-1][j]*exp(-log_etay[i-1][j])*dx);
- double high=(zy[i+1][j]*exp(-log_etay[i+1][j])*(2-dx)
+ double high1=(zy[i+1][j]*exp(-log_etay[i+1][j])*(2-dx)
- zy[i+2][j]*exp(-log_etay[i+2][j])*(1-dx));
- vy[j]=low*(1-dx) + high*dx;
+ double low=interpolate_v(zy,log_etay,i-1,j,h,h*(dx+1));
+ double high=interpolate_v(zy,log_etay,i+2,j,h,h*(-2+dx));
- // if(j>=N/4-2 && j<=N/4+2)
- // // if(j==N/4)
- // std::cout << "vy "
- // << j << " "
- // << vy[j] << " "
- // << low << " "
- // << high << " "
- // // << dx << " "
- // // << 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]) << " "
- // // << zy[i+2][j]*exp(-log_etay[i+2][j]) << " "
- // // << zy[i+3][j]*exp(-log_etay[i+3][j]) << " "
- // << "\n";
+ // vy[j]=low*(1-dx) + high*dx;
+
+ vy[j]=high;
+
+ if(j>=N/4-2 && j<=N/4+2)
+ // if(j==N/4)
+ std::cout << "vy "
+ << j << " "
+ << vy[j] << " "
+ << low << " "
+ << high << " "
+ << log_etay[i][j] << " "
+ << log_etay[i+1][j] << " "
+ << log_etay[i-1][j] << " "
+ // << dx << " "
+ // << 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]) << " "
+ // << zy[i+2][j]*exp(-log_etay[i+2][j]) << " "
+ // << zy[i+3][j]*exp(-log_etay[i+3][j]) << " "
+ << "\n";
}
for(int j=0;j<N;++j)
@@ -296,91 +323,31 @@ int main()
// << "\n";
if(x+h/2>middle && x-h/2<middle)
- // if(dx>h/2 || dx<-h/2)
{
dx=(x>middle) ? ((x-h/2)-middle) : ((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;
+ dp_x-=(p_jump + dx*dp_x_jump)/h;
if(j==N/4)
std::cout << "p "
<< i << " "
<< j << " "
<< x << " "
+ << dp_x << " "
+ << p_jump << " "
+ << dp_x_jump << " "
+ // << dp_fixed << " "
+ << p[i+1][j] << " "
+ << p[i][j] << " "
+ << p[i-1][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+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 << " "
- // << (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) << " "
+ << (p[i-2][j]-p[i-3][j])/h << " "
<< "\n";
More information about the CIG-COMMITS
mailing list