[cig-commits] commit: Go back to computing correction in place
Mercurial
hg at geodynamics.org
Fri Mar 2 10:47:24 PST 2012
changeset: 68:0416056f296c
user: Walter Landry <wlandry at caltech.edu>
date: Fri Mar 02 09:45:24 2012 -0800
files: main.cxx
description:
Go back to computing correction in place
diff -r 3fed1c896c1f -r 0416056f296c main.cxx
--- a/main.cxx Thu Mar 01 03:24:29 2012 -0800
+++ b/main.cxx Fri Mar 02 09:45:24 2012 -0800
@@ -10,36 +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_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_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_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_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_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_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()
@@ -51,22 +51,22 @@ 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],
- 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_dzy[Nx][Ny+1], Cp[Nx][Ny];
+ p[Nx][Ny], dp[Nx][Ny], div[Nx][Ny], fx[Nx+1][Ny], fy[Nx][Ny+1];
+ double zx_new[Nx+1][Ny], zy_new[Nx][Ny+1], p_new[Nx][Ny];
+ // double 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);
+ const bool jacobi(true);
/* Initial conditions */
- const int max_iterations=100000;
+ const int max_iterations=1;
int n_sweeps=1000000;
- const double theta_mom=0.7;
- const double theta_continuity=1.0;
+ const double theta_mom=0.7/2;
+ const double theta_continuity=1.0/2;
const double tolerance=1.0e-6;
const double theta_correction=0.01;
- const int output_interval=1000000;
+ const int output_interval=10;
Model model(Model::solCx);
@@ -85,12 +85,12 @@ int main()
double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
- compute_correction(model,zx,zy,log_etax,log_etay,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,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 */
@@ -102,11 +102,8 @@ int main()
/* zx */
for(int rb=0;rb<2;++rb)
{
- // for(int j=0;j<Ny;++j)
- for(int jj=Ny/2;jj<Ny/2+Ny;++jj)
+ for(int j=0;j<Ny;++j)
{
- int j=jj%Ny;
-
int i_min=(j+rb)%2;
for(int i=i_min;i<Nx+1;i+=2)
{
@@ -191,18 +188,19 @@ 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[i][j];
+ // - dp_x - fx[i][j] + Cx[i][j];
+ - dp_x - fx[i][j] + Cx;
double dRx_dzx=
- // -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx[i][j];
- -6/(h*h) + 2*dlog_etaxx + dlog_etayy;
+ // -6/(h*h) + 2*dlog_etaxx + dlog_etayy;
+ -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
Resid_x[i][j]=Rx;
if(!jacobi)
@@ -211,15 +209,16 @@ int main()
zx_new[i][j]=zx[i][j]-theta_mom*Rx/dRx_dzx;
- // if(j==Ny/4 && Cx!=0)
- // std::cout << "Rx "
- // << i << " "
- // << j << " "
- // << middle << " "
- // << i*h << " "
- // << 6/(h*h) << " "
- // << dCx_dzx << " "
- // << "\n";
+ if(sweep%1000==0 && (j==Ny/4 || j==Ny/4+1) && Cx!=0)
+ std::cout << "Rx "
+ << i << " "
+ << j << " "
+ << middle << " "
+ << i*h << " "
+ << log_etax[i][j] << " "
+ << 6/(h*h) << " "
+ << dCx_dzx << " "
+ << "\n";
}
}
@@ -318,18 +317,19 @@ 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[i][j];
+ // - dp_y - fy[i][j] + Cy[i][j];
+ - dp_y - fy[i][j] + Cy;
double dRy_dzy=
- // -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy[i][j];
- -6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+ // -6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+ -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy;
// if(j==Ny/8 && i*h>middle-10*h && i*h<middle+10*h)
// std::cout << "Ry "
@@ -400,10 +400,11 @@ 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 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];
+ // + Cp[i][j];
+ + Cp;
double dRc_dzxp=1/h - dlog_etax/2;
double dRc_dzxm=-1/h - dlog_etax/2;
@@ -468,12 +469,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[i][j];
-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
+ // -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
if(!jacobi)
zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
@@ -500,12 +501,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[i][j];
-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+ // -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy;
if(!jacobi)
zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
@@ -529,7 +530,7 @@ int main()
p[i][j]=p_new[i][j];
}
- if(sweep%1000==0)
+ if(sweep%10==0)
{
std::cout << "sweep "
<< iteration << " "
More information about the CIG-COMMITS
mailing list