[cig-commits] commit: Uses the average of extrapolation on each side of the interface. Seems to work for fixing the pressure derivatives
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:23 PST 2012
changeset: 24:81927ce18594
user: Walter Landry <wlandry at caltech.edu>
date: Fri Feb 03 13:40:17 2012 -0800
files: constants.hxx initial.cxx main.cxx
description:
Uses the average of extrapolation on each side of the interface. Seems to work for fixing the pressure derivatives
diff -r 567f1fbb28be -r 81927ce18594 constants.hxx
--- a/constants.hxx Fri Feb 03 11:32:24 2012 -0800
+++ b/constants.hxx Fri Feb 03 13:40:17 2012 -0800
@@ -6,7 +6,7 @@ const double log_max_eta=std::log(max_et
const double log_max_eta=std::log(max_eta);
const double log_min_eta=std::log(min_eta);
// const double middle=25.0/64 + 0.0000000001;
-const double middle=400.00000001/1024;
+const double middle=0.4;
const double h(1.0/N);
const double pi(atan(1.0)*4);
diff -r 567f1fbb28be -r 81927ce18594 initial.cxx
--- a/initial.cxx Fri Feb 03 11:32:24 2012 -0800
+++ b/initial.cxx Fri Feb 03 13:40:17 2012 -0800
@@ -175,24 +175,24 @@ void initial(const Model &model, double
- 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";
+ // if(j==N/4 && i>middle*N-5 && i<middle*N+5)
+ // std::cout << "zx "
+ // << 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:
@@ -272,6 +272,16 @@ void initial(const Model &model, double
zy[i][j]=(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-2 && i<middle*N+2)
+ // std::cout << "zy "
+ // << i << " "
+ // << j << " "
+ // << x << " "
+ // << y << " "
+ // << zy[i][j]/eta << " "
+ // << "\n";
+
}
break;
case Model::sinker:
diff -r 567f1fbb28be -r 81927ce18594 main.cxx
--- a/main.cxx Fri Feb 03 11:32:24 2012 -0800
+++ b/main.cxx Fri Feb 03 13:40:17 2012 -0800
@@ -5,7 +5,6 @@
#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],
@@ -48,27 +47,20 @@ int main()
if(model==Model::solCx)
{
- 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;
+ double dx=middle/h-i;
/* 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]);
+ double low=(zx[i][j]*exp(-log_etax[i][j])*(dx+1)
+ - zx[i-1][j]*exp(-log_etax[i-1][j])*dx);
- gsl_spline_init(spline,xa,ya,4);
- vx[j]=gsl_spline_eval(spline,h+dx,accel);
+ double high=(zx[i+1][j]*exp(-log_etax[i+1][j])*(2-dx)
+ - zx[i+2][j]*exp(-log_etax[i+2][j])*(1-dx));
- // vx[j]=zx[i][j]*(1-dx)*exp(-log_etax[i][j])
- // + zx[i+1][j]*dx*exp(-log_etax[i+1][j]);
+ vx[j]=low*(1-dx) + high*dx;
}
dvx_y[0]=dvx_y[N]=0;
@@ -80,40 +72,35 @@ int main()
}
{
int i(middle/h-0.5);
- double dx=h*(middle/h-0.5-i);
-
- // const gsl_interp_type *t=gsl_interp_polynomial;
- // gsl_interp* interp=gsl_interp_alloc(t,4);
+ double dx=middle/h-0.5-i;
/* 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]);
+ double low=(zy[i][j]*exp(-log_etay[i][j])*(dx+1)
+ - zy[i-1][j]*exp(-log_etay[i-1][j])*dx);
- gsl_spline_init(spline,xa,ya,4);
- vy[j]=gsl_spline_eval(spline,h+dx,accel);
+ double high=(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]=zy[i][j]*(1-dx)*exp(-log_etay[i][j])
- // + zy[i+1][j]*dx*exp(-log_etay[i+1][j]);
+ vy[j]=low*(1-dx) + high*dx;
- 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";
+ // 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";
}
- // gsl_interp_free(interp);
for(int j=0;j<N;++j)
dvy_y[j]=(vy[j+1]-vy[j])/h;
@@ -122,8 +109,6 @@ int main()
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 */
@@ -237,83 +222,83 @@ int main()
// dzx_xx-=eta_jump*(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
- 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] << " "
+ // 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] << " "
+ // // << 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 << " "
+ // // << (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+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 << " "
+ // // << (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] << " "
+ // // << 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";
+ // // << 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=(x>middle) ? (middle-(x-h/2)) : ((x+h/2)-middle);
+ 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;
@@ -328,20 +313,21 @@ int main()
std::cout << "p "
<< i << " "
<< j << " "
+ << x << " "
- << p_jump << " "
- << dp_x_jump << " "
+ // << 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-2][j])/h << " "
- << (p[i-2][j]-p[i-3][j])/h << " "
+ // << (p[i-2][j]-p[i-3][j])/h << " "
// << dp_numerical << " "
- << std::fabs((dp_fixed-dp_numerical)/dp_numerical) << " "
+ // << std::fabs((dp_fixed-dp_numerical)/dp_numerical) << " "
// // << x << " "
More information about the CIG-COMMITS
mailing list