[cig-commits] commit: Add in an explicit correction term
Mercurial
hg at geodynamics.org
Tue Feb 21 13:00:16 PST 2012
changeset: 55:bc2368bc255d
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Tue Feb 21 13:00:03 2012 -0800
files: compute_corrections.cxx compute_v_on_interface.cxx constants.hxx initial.cxx main.cxx wscript
description:
Add in an explicit correction term
diff -r f5cd216d64c9 -r bc2368bc255d compute_corrections.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_corrections.cxx Tue Feb 21 13:00:03 2012 -0800
@@ -0,0 +1,125 @@
+#include "constants.hxx"
+#include "Model.hxx"
+
+extern void compute_v_on_interface(const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double &x, const int &j,
+ const int &xyz,
+ double &vx, double &vy,
+ double &dvx_y, double &dvy_y,
+ double &dvx_yy, double &dvy_yy);
+
+void compute_corrections(const Model &model, const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ double Cx[Nx+1][Ny], double Cy[Nx][Ny+1],
+ double Cp[Nx][Ny],
+ const double &theta_correction)
+{
+ if(model==Model::solCx)
+ {
+ for(int j=0;j<Ny+1;++j)
+ for(int i=0;i<Nx+1;++i)
+ {
+ /* Cx */
+ double x(i*h);
+ if(j<Ny && i!=0 && i!=Nx && x-h<middle && x+h>middle)
+ {
+ double Cx_new(0);
+
+ double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
+
+ compute_v_on_interface(zx,zy,log_etax,
+ log_etay,
+ middle,j,0,
+ vx,vy,dvx_y,dvy_y,
+ dvx_yy,dvy_yy);
+
+ double dx=(x>middle) ? (middle-(x-h))
+ : (middle-(x+h));
+ double zx_jump=eta_jump*vx;
+ double dzx_x_jump=eta_jump*dvy_y;
+ double p_jump=-2*dvy_y*eta_jump;
+ double dp_x_jump=2*dvx_yy*eta_jump;
+ double dzx_xx_jump=
+ -dvx_yy*eta_jump + dp_x_jump;
+ double dzx_xx_correction=
+ -(zx_jump + dx*dzx_x_jump
+ + dx*dx*dzx_xx_jump/2)/(h*h);
+ if(x>middle)
+ dzx_xx_correction*=-1;
+
+ Cx_new=dzx_xx_correction;
+
+ if(x+h/2>middle && x-h/2<middle)
+ {
+ dx=(x>middle) ? ((x-h/2)-middle)
+ : ((x+h/2)-middle);
+
+ Cx_new+=(p_jump + dx*dp_x_jump)/h;
+ Cx_new-=eta_jump*(dvy_y - dx*dvx_yy)/h;
+ }
+ Cx[i][j]+=theta_correction*(Cx_new-Cx[i][j]);
+ }
+
+ /* Cy */
+ x=(i+0.5)*h;
+ if(i<Nx && j!=0 && j!=Ny && x-h<middle && x+h>middle)
+ {
+ double Cy_new(0);
+ double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
+
+ compute_v_on_interface(zx,zy,log_etax,
+ log_etay,
+ middle,j,1,
+ vx,vy,dvx_y,dvy_y,
+ dvx_yy,dvy_yy);
+
+ double dx=(x>middle) ? (middle-(x-h))
+ : (middle-(x+h));
+ double zy_jump=eta_jump*vy;
+ double dzy_x_jump=-eta_jump*dvx_y;
+ double dzy_xx_jump=-3*eta_jump*dvy_yy;
+ double dzy_xx_correction=
+ -(zy_jump - dx*dzy_x_jump + dx*dx*dzy_xx_jump/2)/(h*h);
+ if(x>middle)
+ dzy_xx_correction*=-1;
+ Cy_new=dzy_xx_correction;
+
+ if(x+h/2>middle && x-h/2<middle)
+ {
+ dx=(x>middle) ? ((x-h/2)-middle)
+ : ((x+h/2)-middle);
+ double dzx_yx_correction=
+ eta_jump*(dvx_y - dx*dvy_yy)/h;
+ Cy_new-=dzx_yx_correction;
+ }
+ Cy[i][j]+=theta_correction*(Cy_new-Cy[i][j]);
+ }
+
+ /* Cp */
+ if(i<Nx && j<Ny && x+h/2>middle && x-h/2<middle)
+ {
+ double Cp_new(0);
+ double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
+
+ compute_v_on_interface(zx,zy,log_etax,log_etay,
+ middle,j,0,
+ vx,vy,dvx_y,dvy_y,
+ dvx_yy,dvy_yy);
+
+ double dx=(x>middle) ? (middle-(x-h/2)) : (middle-(x+h/2));
+ double zx_jump=eta_jump*vx;
+ double dzx_x_jump=eta_jump*dvy_y;
+
+ double dzx_x_correction=(zx_jump + dx*dzx_x_jump)/h;
+ Cp_new=-dzx_x_correction;
+
+ Cp[i][j]+=theta_correction*(Cp_new-Cp[i][j]);
+ }
+ }
+ }
+}
diff -r f5cd216d64c9 -r bc2368bc255d compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Tue Feb 21 11:21:25 2012 -0800
+++ b/compute_v_on_interface.cxx Tue Feb 21 13:00:03 2012 -0800
@@ -48,33 +48,34 @@ double analytic(const double x, const do
template<int ny>
-double vel(double z[][ny],double log_eta[][ny], const int i, const int j)
+double vel(const double z[][ny], const double log_eta[][ny],
+ const int i, const int j)
{
- if(j>ny/8-4 && j<ny/8+4)
- {
- double x,y;
- if(ny%2==0)
- {
- x=i*h-middle;
- y=(j+0.5)*h-1.0/8;
- // y=j*h-1.0/8;
- return analytic(x,y,true);
- }
- else
- {
- x=(i+0.5)*h-middle;
- y=j*h-1.0/8;
- // y=(j-0.5)*h-1.0/8;
- return analytic(x,y,false);
- }
- }
+ // if(j>ny/8-4 && j<ny/8+4)
+ // {
+ // double x,y;
+ // if(ny%2==0)
+ // {
+ // x=i*h-middle;
+ // y=(j+0.5)*h-1.0/8;
+ // // y=j*h-1.0/8;
+ // return analytic(x,y,true);
+ // }
+ // else
+ // {
+ // x=(i+0.5)*h-middle;
+ // y=j*h-1.0/8;
+ // // y=(j-0.5)*h-1.0/8;
+ // return analytic(x,y,false);
+ // }
+ // }
return z[i][j]*std::exp(-log_eta[i][j]);
}
template<int ny>
-double compute_dv_yy(double z[][ny],
- double log_eta[][ny], const double &dx,
+double compute_dv_yy(const double z[][ny],
+ const double log_eta[][ny], const double &dx,
const int i, const int j)
{
double dv_yy_p1=(vel(z,log_eta,i+1,j+1) - 2*vel(z,log_eta,i+1,j)
@@ -109,24 +110,12 @@ double compute_dv_yy(double z[][ny],
}
double dv_yy_m=(1-dx)*dv_yy_m1 + dx*dv_yy_m0;
- // if(j==ny/8)
- // std::cout << "dvy_yy "
- // << ny << " "
- // << dv_yy_m << " "
- // << dv_yy_m0 << " "
- // << dv_yy_m1 << " "
- // << dv_yy_p << " "
- // << dv_yy_p1 << " "
- // << dv_yy_p2 << " "
- // << (dv_yy_p + dv_yy_m)/2 << " "
- // << "\n";
-
return (dv_yy_p + dv_yy_m)/2;
}
template<int ny>
-double compute_dv_y(double z[][ny],
- double log_eta[][ny],
+double compute_dv_y(const double z[][ny],
+ const double log_eta[][ny],
const double &jump,
const double &dx,
const int i, const int j)
@@ -139,33 +128,13 @@ double compute_dv_y(double z[][ny],
double dv_y_m1=(vel(z,log_eta,i-1,j+1) - vel(z,log_eta,i-1,j))/h;
double dv_y_m=(1-dx)*dv_y_m1 + dx*dv_y_m0;
- if(j==ny/8 && ny%2==1)
- std::cout << "dvy_y "
- << ny << " "
- << i*h-middle << " "
- << (j+0.5)*h-1.0/8 << " "
- << (j+1.5)*h-1.0/8 << " "
- << vel(z,log_eta,i+1,j+1) << " "
- << vel(z,log_eta,i+1,j) << " "
- << vel(z,log_eta,i+2,j+1) << " "
- << vel(z,log_eta,i+2,j) << " "
- << dv_y_m << " "
- << dv_y_m0 << " "
- << dv_y_m1 << " "
- << dv_y_p << " "
- << dv_y_p1 << " "
- << dv_y_p2 << " "
- << (max_eta*dv_y_p + min_eta*dv_y_m - h*jump)
- /(max_eta + min_eta) << " "
- << "\n";
-
return (max_eta*dv_y_p + min_eta*dv_y_m - h*jump)
/(max_eta + min_eta);
}
template<int ny>
-double compute_v(double z[][ny],
- double log_eta[][ny],
+double compute_v(const double z[][ny],
+ const double log_eta[][ny],
const double &jump,
const double &dx,
const int i, const int j)
@@ -204,8 +173,10 @@ double dzy_x_jump(const double &dvx_y)
return -eta_jump*dvx_y;
}
-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],
+void compute_v_on_interface(const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
const double &x, const int &j,
const int &xyz,
double &vx, double &vy,
@@ -220,14 +191,14 @@ void compute_v_on_interface(double zx[Nx
dvy_y=compute_dv_y(zy,log_etay,dzy_yx_jump(dvx_yy),dy,iy,j);
vx=compute_v(zx,log_etax,dzx_x_jump(dvy_y),dx,ix,j);
- if(j==Ny/8)
- std::cout << "compute x "
- << x << " "
- << j << " "
- << vx << " "
- << dvy_y << " "
- << dvx_yy << " "
- << "\n";
+ // if(j==Ny/8)
+ // std::cout << "compute x "
+ // << x << " "
+ // << j << " "
+ // << vx << " "
+ // << dvy_y << " "
+ // << dvx_yy << " "
+ // << "\n";
vy=std::numeric_limits<double>::infinity();
dvx_y=std::numeric_limits<double>::infinity();
@@ -239,14 +210,14 @@ void compute_v_on_interface(double zx[Nx
dvx_y=compute_dv_y(zx,log_etax,dzx_yx_jump(dvy_yy),dx,ix,j-1);
vy=compute_v(zy,log_etay,dzy_x_jump(dvx_y),dy,iy,j);
- if(j==Ny/8)
- std::cout << "compute y "
- << x << " "
- << j << " "
- << vy << " "
- << dvx_y << " "
- << dvy_yy << " "
- << "\n";
+ // if(j==Ny/8)
+ // std::cout << "compute y "
+ // << x << " "
+ // << j << " "
+ // << vy << " "
+ // << dvx_y << " "
+ // << dvy_yy << " "
+ // << "\n";
vx=std::numeric_limits<double>::infinity();
dvy_y=std::numeric_limits<double>::infinity();
diff -r f5cd216d64c9 -r bc2368bc255d constants.hxx
--- a/constants.hxx Tue Feb 21 11:21:25 2012 -0800
+++ b/constants.hxx Tue Feb 21 13:00:03 2012 -0800
@@ -1,8 +1,8 @@ const int N(256);
-const int N(256);
+const int N(64);
const int Nx(N);
const int Ny(N);
const double min_eta=1;
-const double max_eta=1e4;
+const double max_eta=1;
const double eta_jump=max_eta-min_eta;
#include <cmath>
const double log_max_eta=std::log(max_eta);
diff -r f5cd216d64c9 -r bc2368bc255d initial.cxx
--- a/initial.cxx Tue Feb 21 11:21:25 2012 -0800
+++ b/initial.cxx Tue Feb 21 13:00:03 2012 -0800
@@ -8,7 +8,8 @@
void initial(const Model &model, double zx[Nx+1][Ny], double zy[Nx][Ny+1],
double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
- double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1])
+ double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1],
+ double Cx[Nx+1][Ny], double Cy[Nx][Ny+1], double Cp[Nx][Ny])
{
const double A(min_eta*(max_eta-min_eta)/(max_eta+min_eta));
std::complex<double> phi, psi, dphi, v;
@@ -54,6 +55,7 @@ void initial(const Model &model, double
{
if(i<Nx && j<Ny)
{
+ Cp[i][j]=0;
double x((i+0.5)*h), y((j+0.5)*h), r2(x*x+y*y);
switch(model)
{
@@ -102,6 +104,7 @@ void initial(const Model &model, double
}
if(j<Ny)
{
+ Cx[i][j]=0;
double x(i*h), y((j+0.5)*h), r2(x*x+y*y);
zx[i][j]=0;
fx[i][j]=0;
@@ -178,6 +181,7 @@ void initial(const Model &model, double
if(i<Nx)
{
+ Cy[i][j]=0;
double x((i+0.5)*h), y(j*h), r2(x*x+y*y);
zy[i][j]=0;
diff -r f5cd216d64c9 -r bc2368bc255d main.cxx
--- a/main.cxx Tue Feb 21 11:21:25 2012 -0800
+++ b/main.cxx Tue Feb 21 13:00:03 2012 -0800
@@ -8,15 +8,17 @@
extern void initial(const Model &model, double zx[Nx+1][Ny], double zy[Nx][Ny+1],
double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
- double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1]);
+ double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1],
+ double Cx[Nx+1][Ny], double Cy[Nx][Ny+1],
+ double Cp[Nx][Ny]);
-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],
- const double &x, const int &j,
- const int &xyz,
- double &vx, double &vy,
- double &dvx_y, double &dvy_y,
- double &dvx_yy, double &dvy_yy);
+extern void compute_corrections(const Model &model, const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ double Cx[Nx+1][Ny], double Cy[Nx][Ny+1],
+ double Cp[Nx][Ny],
+ const double &theta_correction);
int main()
{
@@ -27,19 +29,21 @@ int main()
with z, it is easy to compute v=z/eta. */
double zx[Nx+1][Ny], zy[Nx][Ny+1], log_etax[Nx+1][Ny], log_etay[Nx][Ny+1],
- p[Nx][Ny], dp[Nx][Ny], div[Nx][Ny], fx[Nx+1][Ny], fy[Nx][Ny+1];
+ p[Nx][Ny], dp[Nx][Ny], div[Nx][Ny], fx[Nx+1][Ny], fy[Nx][Ny+1],
+ Cx[Nx+1][Ny], Cy[Nx][Ny+1], Cp[Nx][Ny];
/* Initial conditions */
const int max_iterations=1;
- int n_sweeps=1;
- const double theta_mom=0.7/eta_jump;
- const double theta_c=1.0/eta_jump;
+ int n_sweeps=10000;
+ const double theta_mom=0.7;
+ const double theta_continuity=1.0;
+ const double theta_correction=0.01;
const double tolerance=1.0e-6;
Model model(Model::solCx);
- initial(model,zx,zy,log_etax,log_etay,p,fx,fy);
+ initial(model,zx,zy,log_etax,log_etay,p,fx,fy,Cx,Cy,Cp);
double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
@@ -52,6 +56,8 @@ int main()
double max_x_resid_previous(0), max_y_resid_previous(0), max_p_resid_previous(0);
for(int sweep=0;sweep<n_sweeps;++sweep)
{
+ compute_corrections(model,zx,zy,log_etax,log_etay,Cx,Cy,Cp,theta_correction);
+
/* zx */
for(int rb=0;rb<2;++rb)
{
@@ -63,8 +69,6 @@ int main()
if(i!=0 && i!=Nx)
{
/* Do the finite difference parts */
-
- const double x(i*h);
/* Derivatives of z */
double dzx_xx=
@@ -134,103 +138,6 @@ int main()
{
dlog_etaxx=dlog_etax=dlog_etayy=dlog_etay=
dlog_etaxy=0;
-
- if(x-h<middle && x+h>middle)
- {
- // theta/=eta_jump;
- double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
-
- /* We waste some effort by computing
- vx etc. once for each side,
- though it does not change. */
- compute_v_on_interface(zx,zy,log_etax,
- log_etay,
- middle,j,0,
- vx,vy,dvx_y,dvy_y,
- dvx_yy,dvy_yy);
-
- double dx=(x>middle) ? (middle-(x-h))
- : (middle-(x+h));
- double zx_jump=eta_jump*vx;
- double dzx_x_jump=eta_jump*dvy_y;
- double p_jump=-2*dvy_y*eta_jump;
- double dp_x_jump=2*dvx_yy*eta_jump;
- double dzx_xx_jump=
- -dvx_yy*eta_jump + dp_x_jump;
- double dzx_xx_correction=
- -(zx_jump + dx*dzx_x_jump
- + dx*dx*dzx_xx_jump/2)/(h*h);
- if(x>middle)
- dzx_xx_correction*=-1;
-
- // if(j==Ny/8)
- // std::cout << "vx "
- // << i << " "
- // << j << " "
- // << x << " "
- // << middle << " "
- // // << zx_jump << " "
- // // << dzx_x_jump << " "
- // // << dzx_xx_jump << " "
- // // << pi*pi*zx[i+1][j]*exp(-log_etax[i+1][j])*eta_jump << " "
- // // << pi*pi*zx[i][j]*exp(-log_etax[i][j])*eta_jump << " "
- // // << pi*pi*zx[i-1][j]*exp(-log_etax[i-1][j])*eta_jump << " "
- // // << dzx_xx << " "
- // // << (zx[i+3][j]-2*zx[i+2][j]+zx[i+1][j])/(h*h) << " "
- // // << (zx[i+2][j]-2*zx[i+1][j]+zx[i][j])/(h*h) << " "
- // // << (zx[i+1][j]-2*zx[i][j]+zx[i-1][j])/(h*h) << " "
- // // << (zx[i][j]-2*zx[i-1][j]+zx[i-2][j])/(h*h) << " "
- // // << (zx[i-1][j]-2*zx[i-2][j]+zx[i-3][j])/(h*h) << " "
- // // << dzx_xx+dzx_xx_correction << " "
-
- // // << vx << " "
- // // << exp(-log_etax[i+2][j])*zx[i+2][j] << " "
- // // << exp(-log_etax[i+1][j])*zx[i+1][j] << " "
- // // << exp(-log_etax[i][j])*zx[i][j] << " "
- // // << exp(-log_etax[i-1][j])*zx[i-1][j] << " "
- // // << exp(-log_etax[i-2][j])*zx[i-2][j] << " "
-
- // // << dvy_y << " "
- // // << exp(-log_etay[i+2][j])*(zy[i+2][j+1]-zy[i+2][j])/(h) << " "
- // // << exp(-log_etay[i+1][j])*(zy[i+1][j+1]-zy[i+1][j])/(h) << " "
- // // << exp(-log_etay[i][j])*(zy[i][j+1]-zy[i][j])/(h) << " "
- // // << exp(-log_etay[i-1][j])*(zy[i-1][j+1]-zy[i-1][j])/(h) << " "
-
-
- // // << dvx_yy << " "
- // // << exp(-log_etax[i+2][j])*(zx[i+2][j+1]-2*zx[i+2][j]+zx[i+2][j-1])/(h*h) << " "
- // // << exp(-log_etax[i+1][j])*(zx[i+1][j+1]-2*zx[i+1][j]+zx[i+1][j-1])/(h*h) << " "
- // // << exp(-log_etax[i][j])*(zx[i][j+1]-2*zx[i][j]+zx[i][j-1])/(h*h) << " "
- // // << exp(-log_etax[i-1][j])*(zx[i-1][j+1]-2*zx[i-1][j]+zx[i-1][j-1])/(h*h) << " "
- // << "\n";
-
-
- dzx_xx+=dzx_xx_correction;
-
- if(x+h/2>middle && x-h/2<middle)
- {
- dx=(x>middle) ? ((x-h/2)-middle)
- : ((x+h/2)-middle);
-
- // if(j==Ny/8)
- // std::cout << "dp "
- // << i << " "
- // << j << " "
- // << x << " "
- // << middle << " "
- // << dp_x-(p_jump + dx*dp_x_jump)/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 << " "
- // << "\n";
-
- dp_x-=(p_jump + dx*dp_x_jump)/h;
-
- dzy_xy-=eta_jump*(dvy_y - dx*dvx_yy)/h;
- }
- }
}
/* Compute the residual and update zx */
@@ -239,12 +146,12 @@ int main()
+ 2*(dlog_etaxx*zx[i][j] + dlog_etax*dzx_x)
+ dlog_etayy*zx[i][j] + dlog_etay*dzx_y
+ dlog_etaxy*zy_avg + dlog_etax*dzy_y
- - dp_x - fx[i][j];
+ - dp_x - fx[i][j] + Cx[i][j];
double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
Resid_x[i][j]=Rx;
- // zx[i][j]-=theta*Rx/dRx_dzx;
+ zx[i][j]-=theta*Rx/dRx_dzx;
}
}
}
@@ -331,111 +238,8 @@ int main()
/* Now do the jump conditions */
if(model==Model::solCx)
{
- double x((i+0.5)*h);
dlog_etayy=dlog_etay=dlog_etaxx=dlog_etax=
dlog_etayx=0;
-
- if(x-h<middle && x+h>middle)
- {
- // theta/=eta_jump;
- double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
-
- /* We waste some effort by computing
- vx etc. once for each side,
- though it does not change. */
- compute_v_on_interface(zx,zy,log_etax,
- log_etay,
- middle,j,1,
- vx,vy,dvx_y,dvy_y,
- dvx_yy,dvy_yy);
-
- double dx=(x>middle) ? (middle-(x-h))
- : (middle-(x+h));
- double zy_jump=eta_jump*vy;
- double dzy_x_jump=-eta_jump*dvx_y;
- double dzy_xx_jump=-3*eta_jump*dvy_yy;
- double dzy_xx_correction=
- -(zy_jump - dx*dzy_x_jump + dx*dx*dzy_xx_jump/2)/(h*h);
- if(x>middle)
- dzy_xx_correction*=-1;
- dzy_xx+=dzy_xx_correction;
-
- if(x+h/2>middle && x-h/2<middle)
- {
- dx=(x>middle) ? ((x-h/2)-middle)
- : ((x+h/2)-middle);
- double dzx_yx_correction=
- eta_jump*(dvx_y - dx*dvy_yy)/h;
- if(j==Ny/8)
- std::cout << "dzx_yx "
- << i << " "
- << j << " "
- << dx << " "
- << middle << " "
- << x << " "
- // << ((x-h/2)-middle)/dx << " "
- // << ((x+h/2)-middle)/dx << " "
- // << dzx_yx << " "
- // << dzx_yx_correction << " "
- // << eta_jump*dvx_y[j]/h << " "
- // << eta_jump*dx*dvy_yy[j]/h << " "
- // << dzx_yx << " "
- // << eta_jump*(dvx_y)/h << " "
- // << eta_jump*(-dx*dvy_yy)/h << " "
- << dzx_yx - dzx_yx_correction << " "
-
- << ((zx[i+3][j] - zx[i+3][j-1])
- - (zx[i+2][j] - zx[i+2][j-1]))/(h*h) << " "
- << ((zx[i+2][j] - zx[i+2][j-1])
- - (zx[i+1][j] - zx[i+1][j-1]))/(h*h) << " "
- << ((zx[i+1][j] - zx[i+1][j-1])
- - (zx[i][j] - zx[i][j-1]))/(h*h) << " "
- << ((zx[i][j] - zx[i][j-1])
- - (zx[i-1][j] - zx[i-1][j-1]))/(h*h) << " "
- << ((zx[i-1][j] - zx[i-1][j-1])
- - (zx[i-2][j] - zx[i-2][j-1]))/(h*h) << " "
- << ((zx[i-2][j] - zx[i-2][j-1])
- - (zx[i-3][j] - zx[i-3][j-1]))/(h*h) << " "
-
- << dvx_y << " "
- << exp(-log_etax[i+2][j])*(zx[i+2][j]-zx[i+2][j-1])/h << " "
- << exp(-log_etax[i+1][j])*(zx[i+1][j]-zx[i+1][j-1])/h << " "
- << exp(-log_etax[i][j])*(zx[i][j]-zx[i][j-1])/h << " "
- << exp(-log_etax[i-1][j])*(zx[i-1][j]-zx[i-1][j-1])/h << " "
- << exp(-log_etax[i-2][j])*(zx[i-2][j]-zx[i-2][j-1])/h << " "
-
-
- // << dvy_yy << " "
- // << exp(-log_etay[i+2][j])*(zy[i+2][j+1]-2*zy[i+2][j]+zy[i+2][j-1])/(h*h) << " "
- // << exp(-log_etay[i+1][j])*(zy[i+1][j+1]-2*zy[i+1][j]+zy[i+1][j-1])/(h*h) << " "
- // << exp(-log_etay[i][j])*(zy[i][j+1]-2*zy[i][j]+zy[i][j-1])/(h*h) << " "
- // << exp(-log_etay[i-1][j])*(zy[i-1][j+1]-2*zy[i-1][j]+zy[i-1][j-1])/(h*h) << " "
- // << exp(-log_etay[i-2][j])*(zy[i-2][j+1]-2*zy[i-2][j]+zy[i-2][j-1])/(h*h) << " "
-
-
- << "\n";
- dzx_yx-=dzx_yx_correction;
- }
- }
-
- // if(j==Ny/8 && (x-2*h<middle && x+2*h>middle))
- // std::cout << "Ry "
- // << i << " "
- // << j << " "
- // << x << " "
- // << middle << " "
- // << 2*dzy_yy + dzy_xx + dzx_yx
- // + 2*(dlog_etayy*zy[i][j] + dlog_etay*dzy_y)
- // + dlog_etaxx*zy[i][j] + dlog_etax*dzy_x
- // + dlog_etayx*zx_avg + dlog_etay*dzx_x
- // - dp_y - fy[i][j] << " "
- // << 2*dzy_yy << " "
- // << dzy_xx << " "
- // << dzx_yx << " "
- // << dp_y << " "
- // << fy[i][j] << " "
- // << "\n";
-
}
/* Compute the residual and update zy */
@@ -449,7 +253,7 @@ int main()
double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
Resid_y[i][j]=Ry;
- // zy[i][j]-=theta*Ry/dRy_dzy;
+ zy[i][j]-=theta*Ry/dRy_dzy;
}
}
}
@@ -483,37 +287,15 @@ int main()
double zx_avg=(zx[i+1][j] + zx[i][j])/2;
double zy_avg=(zy[i][j+1] + zy[i][j])/2;
- double theta(theta_c);
+ double theta(theta_continuity);
/* Apply the jump condition */
- double x((i+0.5)*h);
if(model==Model::solCx)
{
dlog_etax=dlog_etay=dlog_etaxx=dlog_etayy=0;
-
- if(x+h/2>middle && x-h/2<middle)
- {
- // theta/=eta_jump;
- double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
-
- /* We waste some effort by computing vx
- etc. once for each side, though it does not
- change. */
- compute_v_on_interface(zx,zy,log_etax,log_etay,
- middle,j,0,
- vx,vy,dvx_y,dvy_y,
- dvx_yy,dvy_yy);
-
- double dx=(x>middle) ? (middle-(x-h/2)) : (middle-(x+h/2));
- double zx_jump=eta_jump*vx;
- double dzx_x_jump=eta_jump*dvy_y;
-
- double dzx_x_correction=(zx_jump + dx*dzx_x_jump)/h;
- dzx_x-=dzx_x_correction;
- }
}
- double Rc=dzx_x + dzy_y + dlog_etax*zx_avg + dlog_etay*zy_avg;
+ double Rc=dzx_x + dzy_y + dlog_etax*zx_avg + dlog_etay*zy_avg + Cp[i][j];
double dRc_dzxp=1/h - dlog_etax/2;
double dRc_dzxm=-1/h - dlog_etax/2;
@@ -533,7 +315,7 @@ int main()
Resid_p[i][j]=Rc;
dp[i][j]=-theta*Rc/dRc_dp;
- // p[i][j]+=dp[i][j];
+ p[i][j]+=dp[i][j];
}
if(sweep%1000==0 || sweep>236000)
@@ -570,7 +352,7 @@ int main()
double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
- // zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
+ zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
max_x_resid=std::max(std::fabs(Resid_x[i][j]),max_x_resid);
}
@@ -594,7 +376,7 @@ int main()
double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
- // zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
+ zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
max_y_resid=std::max(fabs(Resid_y[i][j]),max_y_resid);
}
diff -r f5cd216d64c9 -r bc2368bc255d wscript
--- a/wscript Tue Feb 21 11:21:25 2012 -0800
+++ b/wscript Tue Feb 21 13:00:03 2012 -0800
@@ -11,7 +11,8 @@ def build(bld):
features = ['cxx','cprogram'],
source = ['main.cxx',
'initial.cxx',
- 'compute_v_on_interface.cxx'],
+ 'compute_v_on_interface.cxx',
+ 'compute_corrections.cxx'],
target = 'Gamr_disc',
# cxxflags = ['-std=c++0x','-g','-D_GLIBCXX_DEBUG','-Wall'],
More information about the CIG-COMMITS
mailing list