[cig-commits] commit: Use periodic BC on top and bottom. Also go back to computing correction once per iteration, rather than once per sweep
Mercurial
hg at geodynamics.org
Thu Mar 1 03:24:38 PST 2012
changeset: 67:3fed1c896c1f
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Thu Mar 01 03:24:29 2012 -0800
files: compute_correction.cxx compute_v_on_interface.cxx constants.hxx main.cxx wscript
description:
Use periodic BC on top and bottom. Also go back to computing correction once per iteration, rather than once per sweep
diff -r e7f5c2b63f02 -r 3fed1c896c1f compute_correction.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_correction.cxx Thu Mar 01 03:24:29 2012 -0800
@@ -0,0 +1,80 @@
+#include "constants.hxx"
+#include "Model.hxx"
+
+extern void compute_Cx(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],
+ const int &i, const int &j,
+ double &Cx, double &dCx_dzx);
+
+extern void compute_Cy(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],
+ const int &i, const int &j,
+ double &Cy, double &dCy_dzy);
+
+extern void compute_Cp(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],
+ const int &i, const int &j,
+ double &Cp);
+
+void compute_correction(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 dCx[Nx+1][Ny],
+ double Cy[Nx][Ny+1], double dCy[Nx][Ny+1],
+ double Cp[Nx][Ny],
+ const double &theta_correction)
+{
+ if(model==Model::solCx)
+ {
+ double C_new, dC_new;
+ for(int j=0;j<Ny;++j)
+ for(int i=0;i<Nx+1;++i)
+ {
+ compute_Cx(model,zx,zy,log_etax,log_etay,i,j,C_new,dC_new);
+ if(theta_correction==1.0)
+ {
+ Cx[i][j]=C_new;
+ dCx[i][j]=dC_new;
+ }
+ else
+ {
+ Cx[i][j]+=theta_correction*(C_new-Cx[i][j]);
+ dCx[i][j]+=theta_correction*(dC_new-dCx[i][j]);
+ }
+ }
+
+ for(int j=0;j<Ny+1;++j)
+ for(int i=0;i<Nx;++i)
+ {
+ compute_Cy(model,zx,zy,log_etax,log_etay,i,j,C_new,dC_new);
+ if(theta_correction==1.0)
+ {
+ Cy[i][j]=C_new;
+ dCy[i][j]=dC_new;
+ }
+ else
+ {
+ Cy[i][j]+=theta_correction*(C_new-Cy[i][j]);
+ dCy[i][j]+=theta_correction*(dC_new-dCy[i][j]);
+ }
+ }
+
+ for(int j=0;j<Ny;++j)
+ for(int i=0;i<Nx;++i)
+ {
+ compute_Cp(model,zx,zy,log_etax,log_etay,i,j,C_new);
+ if(theta_correction==1.0)
+ Cp[i][j]=C_new;
+ else
+ Cp[i][j]+=theta_correction*(C_new-Cp[i][j]);
+ }
+ }
+}
+
diff -r e7f5c2b63f02 -r 3fed1c896c1f compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Wed Feb 29 20:39:46 2012 -0800
+++ b/compute_v_on_interface.cxx Thu Mar 01 03:24:29 2012 -0800
@@ -81,13 +81,21 @@ double compute_dv_yy(const double z[][ny
double dv_yy_p1, dv_yy_p2;
if(j==0)
{
- dv_yy_p1=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/(h*h);
- dv_yy_p2=(vel(z,log_eta,i+2,j+1) - vel(z,log_eta,i+2,j))/(h*h);
+ dv_yy_p1=(vel(z,log_eta,i+1,j+1) - 2*vel(z,log_eta,i+1,j)
+ - vel(z,log_eta,i+1,ny-1))/(h*h);
+ dv_yy_p2=(vel(z,log_eta,i+2,j+1) - 2*vel(z,log_eta,i+2,j)
+ - vel(z,log_eta,i+2,ny-1))/(h*h);
+ // dv_yy_p1=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/(h*h);
+ // dv_yy_p2=(vel(z,log_eta,i+2,j+1) - vel(z,log_eta,i+2,j))/(h*h);
}
else if(j==ny-1)
{
- dv_yy_p1=(-vel(z,log_eta,i+1,j) + vel(z,log_eta,i+1,j-1))/(h*h);
- dv_yy_p2=(-vel(z,log_eta,i+2,j) + vel(z,log_eta,i+2,j-1))/(h*h);
+ dv_yy_p1=(-vel(z,log_eta,i+1,0) - 2*vel(z,log_eta,i+1,j)
+ + vel(z,log_eta,i+1,j-1))/(h*h);
+ dv_yy_p2=(-vel(z,log_eta,i+2,0) - 2*vel(z,log_eta,i+2,j)
+ + vel(z,log_eta,i+2,j-1))/(h*h);
+ // dv_yy_p1=(-vel(z,log_eta,i+1,j) + vel(z,log_eta,i+1,j-1))/(h*h);
+ // dv_yy_p2=(-vel(z,log_eta,i+2,j) + vel(z,log_eta,i+2,j-1))/(h*h);
}
else
{
@@ -101,13 +109,21 @@ double compute_dv_yy(const double z[][ny
double dv_yy_m0, dv_yy_m1;
if(j==0)
{
- dv_yy_m0=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/(h*h);
- dv_yy_m1=(vel(z,log_eta,i-1,j+1) - vel(z,log_eta,i-1,j))/(h*h);
+ dv_yy_m0=(vel(z,log_eta,i,j+1) - 2*vel(z,log_eta,i,j)
+ - vel(z,log_eta,i,ny-1))/(h*h);
+ dv_yy_m1=(vel(z,log_eta,i-1,j+1) - 2*vel(z,log_eta,i-1,j)
+ - vel(z,log_eta,i-1,ny-1))/(h*h);
+ // dv_yy_m0=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/(h*h);
+ // dv_yy_m1=(vel(z,log_eta,i-1,j+1) - vel(z,log_eta,i-1,j))/(h*h);
}
else if(j==ny-1)
{
- dv_yy_m0=(-vel(z,log_eta,i,j) + vel(z,log_eta,i,j-1))/(h*h);
- dv_yy_m1=(-vel(z,log_eta,i-1,j) + vel(z,log_eta,i-1,j-1))/(h*h);
+ dv_yy_m0=(-vel(z,log_eta,i,0) - 2*vel(z,log_eta,i,j)
+ + vel(z,log_eta,i,j-1))/(h*h);
+ dv_yy_m1=(-vel(z,log_eta,i-1,0) - 2*vel(z,log_eta,i-1,j)
+ + vel(z,log_eta,i-1,j-1))/(h*h);
+ // dv_yy_m0=(-vel(z,log_eta,i,j) + vel(z,log_eta,i,j-1))/(h*h);
+ // dv_yy_m1=(-vel(z,log_eta,i-1,j) + vel(z,log_eta,i-1,j-1))/(h*h);
}
else
{
diff -r e7f5c2b63f02 -r 3fed1c896c1f constants.hxx
--- a/constants.hxx Wed Feb 29 20:39:46 2012 -0800
+++ b/constants.hxx Thu Mar 01 03:24:29 2012 -0800
@@ -2,7 +2,7 @@ const int Nx(N);
const int Nx(N);
const int Ny(N);
const double min_eta=1;
-const double max_eta=10;
+const double max_eta=100;
const double eta_jump=max_eta-min_eta;
#include <cmath>
const double log_max_eta=std::log(max_eta);
diff -r e7f5c2b63f02 -r 3fed1c896c1f main.cxx
--- a/main.cxx Wed Feb 29 20:39:46 2012 -0800
+++ b/main.cxx Thu Mar 01 03:24:29 2012 -0800
@@ -10,26 +10,36 @@ extern void initial(const Model &model,
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]);
-extern void compute_Cx(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],
- const int &i, const int &j,
- double &Cx, double &dCx_dzx);
+extern void compute_correction(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 dCx[Nx+1][Ny],
+ double Cy[Nx][Ny+1], double dCy[Nx][Ny+1],
+ double Cp[Nx][Ny],
+ const double &theta_correction);
-extern void compute_Cy(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],
- const int &i, const int &j,
- double &Cy, double &dCy_dzy);
-extern void compute_Cp(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],
- const int &i, const int &j,
- double &Cp);
+// extern void compute_Cx(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],
+// const int &i, const int &j,
+// double &Cx, double &dCx_dzx);
+
+// extern void compute_Cy(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],
+// const int &i, const int &j,
+// double &Cy, double &dCy_dzy);
+
+// extern void compute_Cp(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],
+// const int &i, const int &j,
+// double &Cp);
int main()
@@ -43,18 +53,18 @@ int main()
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],
zx_new[Nx+1][Ny], zy_new[Nx][Ny+1], p_new[Nx][Ny],
- Cx[Nx+1][Ny], dCx_dzx[Nx+1][Ny], Cy[Nx][Ny+1], dCy_zy[Nx][Ny+1], Cp[Nx][Ny];
+ Cx[Nx+1][Ny], dCx_dzx[Nx+1][Ny], Cy[Nx][Ny+1], dCy_dzy[Nx][Ny+1], Cp[Nx][Ny];
const bool jacobi(false);
/* Initial conditions */
- const int max_iterations=1;
- int n_sweeps=100000;
+ const int max_iterations=100000;
+ int n_sweeps=1000000;
const double theta_mom=0.7;
const double theta_continuity=1.0;
const double tolerance=1.0e-6;
- // const double theta_correction=0.01;
+ const double theta_correction=0.01;
const int output_interval=1000000;
@@ -62,14 +72,25 @@ int main()
initial(model,zx,zy,log_etax,log_etay,p,fx,fy);
+ // for(int i=1;i<Nx+1;++i)
+ // for(int j=1;j<Ny+1;++j)
+ // {
+ // if(j<Ny-1)
+ // zx[i][j]=zx_new[i][j];
+ // if(i<Nx-1)
+ // zy[i][j]=zy_new[i][j];
+ // if(j<Ny && i<Nx)
+ // p[i][j]=p_new[i][j];
+ // }
+
double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
- // compute_correction(model,zx,zy,log_etax,log_etay,i,j,Cx,dCx_dzx,
- // Cy,dCy_dzy,Cp,1.0);
+ compute_correction(model,zx,zy,log_etax,log_etay,Cx,dCx_dzx,
+ Cy,dCy_dzy,Cp,1.0);
for(int iteration=0;iteration<max_iterations;++iteration)
{
- // compute_correction(model,zx,zy,log_etax,log_etay,i,j,Cx,dCx_dzx,
- // Cy,dCy_dzy,Cp,theta_correction);
+ compute_correction(model,zx,zy,log_etax,log_etay,Cx,dCx_dzx,
+ Cy,dCy_dzy,Cp,theta_correction);
/* Smoothing sweeps */
@@ -81,8 +102,11 @@ int main()
/* zx */
for(int rb=0;rb<2;++rb)
{
- for(int j=0;j<Ny;++j)
+ // for(int j=0;j<Ny;++j)
+ for(int jj=Ny/2;jj<Ny/2+Ny;++jj)
{
+ int j=jj%Ny;
+
int i_min=(j+rb)%2;
for(int i=i_min;i<Nx+1;i+=2)
{
@@ -104,13 +128,19 @@ int main()
double dzx_y, dzx_yy;
if(j==0)
{
- dzx_yy=(-zx[i][j] + zx[i][j+1])/(h*h);
- dzx_y=(zx[i][j+1] - zx[i][j])/(2*h);
+ dzx_yy=
+ (-zx[i][Ny-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
+ dzx_y=(zx[i][j+1] + zx[i][Ny-1])/(2*h);
+ // dzx_yy=(-zx[i][j] + zx[i][j+1])/(h*h);
+ // dzx_y=(zx[i][j+1] - zx[i][j])/(2*h);
}
else if(j==Ny-1)
{
- dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
- dzx_y=(zx[i][j] - zx[i][j-1])/(2*h);
+ dzx_yy=
+ (zx[i][j-1] - 2*zx[i][j] - zx[i][0])/(h*h);
+ dzx_y=(-zx[i][0] - zx[i][j-1])/(2*h);
+ // dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
+ // dzx_y=(zx[i][j] - zx[i][j-1])/(2*h);
}
else
{
@@ -161,17 +191,18 @@ int main()
/* Compute the residual and update zx */
- double Cx, dCx_dzx;
- compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,
- dCx_dzx);
+ // double Cx, dCx_dzx;
+ // compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,
+ // dCx_dzx);
double Rx=2*dzx_xx + dzx_yy + dzy_xy
+ 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] + Cx;
+ - dp_x - fx[i][j] + Cx[i][j];
double dRx_dzx=
- -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
+ // -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx[i][j];
+ -6/(h*h) + 2*dlog_etaxx + dlog_etayy;
Resid_x[i][j]=Rx;
if(!jacobi)
@@ -287,17 +318,18 @@ int main()
/* Compute the residual and update zy */
- double Cy, dCy_dzy;
- compute_Cy(model,zx,zy,log_etax,log_etay,i,j,Cy,
- dCy_dzy);
+ // double Cy, dCy_dzy;
+ // compute_Cy(model,zx,zy,log_etax,log_etay,i,j,Cy,
+ // dCy_dzy);
double Ry=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] + Cy;
+ - dp_y - fy[i][j] + Cy[i][j];
double dRy_dzy=
- -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy;
+ // -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy[i][j];
+ -6/(h*h) + 2*dlog_etayy + dlog_etaxx;
// if(j==Ny/8 && i*h>middle-10*h && i*h<middle+10*h)
// std::cout << "Ry "
@@ -368,10 +400,10 @@ int main()
dlog_etax=dlog_etay=dlog_etaxx=dlog_etayy=0;
}
- double Cp;
- compute_Cp(model,zx,zy,log_etax,log_etay,i,j,Cp);
- double Rc=
- dzx_x + dzy_y + dlog_etax*zx_avg + dlog_etay*zy_avg + Cp;
+ // double Cp;
+ // compute_Cp(model,zx,zy,log_etax,log_etay,i,j,Cp);
+ 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;
@@ -436,11 +468,12 @@ int main()
dlog_etaxx=dlog_etayy=0;
}
- double Cx, dCx_dzx;
- compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,
- dCx_dzx);
+ // double Cx, dCx_dzx;
+ // compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,
+ // dCx_dzx);
double dRx_dzx=
- -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
+ // -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx[i][j];
+ -6/(h*h) + 2*dlog_etaxx + dlog_etayy;
if(!jacobi)
zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
@@ -467,11 +500,12 @@ int main()
dlog_etaxx=dlog_etayy=0;
}
- double Cy, dCy_dzy;
- compute_Cy(model,zx,zy,log_etax,log_etay,i,j,Cy,
- dCy_dzy);
+ // double Cy, dCy_dzy;
+ // compute_Cy(model,zx,zy,log_etax,log_etay,i,j,Cy,
+ // dCy_dzy);
double dRy_dzy=
- -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy;
+ // -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy[i][j];
+ -6/(h*h) + 2*dlog_etayy + dlog_etaxx;
if(!jacobi)
zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
@@ -495,7 +529,7 @@ int main()
p[i][j]=p_new[i][j];
}
- if(sweep%100==0)
+ if(sweep%1000==0)
{
std::cout << "sweep "
<< iteration << " "
@@ -505,17 +539,17 @@ int main()
<< max_p_resid << " "
<< "\n";
}
- // if(sweep%output_interval==0)
- // {
- // std::stringstream ss;
- // ss << "zx" << sweep;
- // write_vtk(ss.str(),Nx+1,N,zx);
- // ss.str("");
- // ss << "zy" << sweep;
- // write_vtk(ss.str(),Nx,N,zy);
- // ss.str("");
- // ss << "p" << sweep;
- // write_vtk(ss.str(),Nx,N,p);
+ if(sweep%output_interval==0)
+ {
+ std::stringstream ss;
+ ss << "zx" << sweep;
+ write_vtk(ss.str(),Nx+1,N,zx);
+ ss.str("");
+ ss << "zy" << sweep;
+ write_vtk(ss.str(),Nx,N,zy);
+ ss.str("");
+ ss << "p" << sweep;
+ write_vtk(ss.str(),Nx,N,p);
// for(int j=0;j<Ny;++j)
// for(int i=0;i<Nx;++i)
@@ -526,7 +560,7 @@ int main()
// ss.str("");
// ss << "div" << sweep;
// write_vtk(ss.str(),Nx,N,div);
- // }
+ }
// if(fabs(max_x_resid-max_x_resid_previous)/max_x_resid < tolerance
// && fabs(max_y_resid-max_y_resid_previous)/max_y_resid < tolerance
@@ -548,8 +582,8 @@ int main()
break;
}
}
- if(sweep==0)
- break;
+ // if(sweep==0)
+ // break;
}
write_vtk("zx",Nx+1,N,zx);
diff -r e7f5c2b63f02 -r 3fed1c896c1f wscript
--- a/wscript Wed Feb 29 20:39:46 2012 -0800
+++ b/wscript Thu Mar 01 03:24:29 2012 -0800
@@ -12,6 +12,7 @@ def build(bld):
source = ['main.cxx',
'initial.cxx',
'compute_v_on_interface.cxx',
+ 'compute_correction.cxx',
'compute_Cx.cxx',
'compute_Cy.cxx',
'compute_Cp.cxx'],
More information about the CIG-COMMITS
mailing list