[cig-commits] commit: Do Jacobi or Gauss-Seidel with a compile time flag
Mercurial
hg at geodynamics.org
Thu Mar 1 03:24:37 PST 2012
changeset: 66:e7f5c2b63f02
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 29 20:39:46 2012 -0800
files: main.cxx
description:
Do Jacobi or Gauss-Seidel with a compile time flag
diff -r 4e071e769780 -r e7f5c2b63f02 main.cxx
--- a/main.cxx Wed Feb 29 19:20:03 2012 -0800
+++ b/main.cxx Wed Feb 29 20:39:46 2012 -0800
@@ -41,17 +41,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];
+ 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];
+
+ const bool jacobi(false);
/* Initial conditions */
const int max_iterations=1;
- int n_sweeps=10000000;
+ int n_sweeps=100000;
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 int output_interval=100000000;
+ const int output_interval=1000000;
Model model(Model::solCx);
@@ -59,8 +64,13 @@ 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,i,j,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);
+
/* Smoothing sweeps */
// if(iteration==max_iterations-1)
@@ -154,7 +164,6 @@ int main()
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
@@ -165,7 +174,22 @@ int main()
-6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
Resid_x[i][j]=Rx;
- zx[i][j]-=theta_mom*Rx/dRx_dzx;
+ if(!jacobi)
+ zx[i][j]-=theta_mom*Rx/dRx_dzx;
+ else
+ 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";
+
}
}
}
@@ -175,7 +199,7 @@ int main()
{
std::stringstream ss;
ss << "zx_resid" << sweep;
- // write_vtk(ss.str(),Nx+1,N,Resid_x);
+ write_vtk(ss.str(),Nx+1,N,Resid_x);
// if(sweep==0)
// {
// ss.str("");
@@ -295,7 +319,10 @@ int main()
Resid_y[i][j]=Ry;
- zy[i][j]-=theta_mom*Ry/dRy_dzy;
+ if(!jacobi)
+ zy[i][j]-=theta_mom*Ry/dRy_dzy;
+ else
+ zy_new[i][j]=zy[i][j]-theta_mom*Ry/dRy_dzy;
}
}
}
@@ -305,7 +332,7 @@ int main()
{
std::stringstream ss;
ss << "zy_resid" << sweep;
- // write_vtk(ss.str(),Nx,N,Resid_y);
+ write_vtk(ss.str(),Nx,N,Resid_y);
// if(sweep==0)
// {
// ss.str("");
@@ -364,14 +391,18 @@ int main()
Resid_p[i][j]=Rc;
dp[i][j]=-theta_continuity*Rc/dRc_dp;
- p[i][j]+=dp[i][j];
+
+ if(!jacobi)
+ p[i][j]+=dp[i][j];
+ else
+ p_new[i][j]=p[i][j]+dp[i][j];
}
if(sweep%output_interval==0)
{
std::stringstream ss;
ss << "p_resid" << sweep;
- // write_vtk(ss.str(),Nx,N,Resid_p);
+ write_vtk(ss.str(),Nx,N,Resid_p);
// if(sweep==0)
// {
// ss.str("");
@@ -405,9 +436,16 @@ int main()
dlog_etaxx=dlog_etayy=0;
}
- double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
+ 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;
- zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
+ if(!jacobi)
+ zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
+ else
+ zx_new[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);
}
@@ -429,14 +467,33 @@ int main()
dlog_etaxx=dlog_etayy=0;
}
- double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+ 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;
- zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
+ if(!jacobi)
+ zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
+ else
+ zy_new[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);
}
max_p_resid=std::max(fabs(Resid_p[i][j]),max_p_resid);
}
+
+ if(jacobi)
+ for(int i=0;i<Nx+1;++i)
+ for(int j=0;j<Ny+1;++j)
+ {
+ if(j<Ny)
+ zx[i][j]=zx_new[i][j];
+ if(i<Nx)
+ zy[i][j]=zy_new[i][j];
+ if(j<Ny && i<Nx)
+ p[i][j]=p_new[i][j];
+ }
if(sweep%100==0)
{
More information about the CIG-COMMITS
mailing list