[cig-commits] commit: Go back to Gauss-Seidel. Scale the relaxation parameter by the jump in viscosity. Seems to work for high viscosity variations, though it is too slow to actually get an answer
Mercurial
hg at geodynamics.org
Fri Feb 17 09:47:22 PST 2012
changeset: 52:c91ba26ea86d
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Feb 17 09:47:09 2012 -0800
files: main.cxx
description:
Go back to Gauss-Seidel. Scale the relaxation parameter by the jump in viscosity. Seems to work for high viscosity variations, though it is too slow to actually get an answer
diff -r d8119d27f1e7 -r c91ba26ea86d main.cxx
--- a/main.cxx Fri Feb 17 09:45:39 2012 -0800
+++ b/main.cxx Fri Feb 17 09:47:09 2012 -0800
@@ -29,14 +29,12 @@ 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];
- double zx_new[Nx+1][Ny], zy_new[Nx][Ny+1], p_new[Nx][Ny];
-
/* Initial conditions */
- const int max_iterations=10000;
+ const int max_iterations=1;
int n_sweeps=10000000;
- const double theta_mom=0.7/100;
- const double theta_c=1.0/100;
+ const double theta_mom=0.7/eta_jump;
+ const double theta_c=1.0/eta_jump;
const double tolerance=1.0e-6;
Model model(Model::solCx);
@@ -130,6 +128,7 @@ int main()
const double zy_avg=(zy[i][j] + zy[i-1][j] +
zy[i][j+1] + zy[i-1][j+1])/4;
+ double theta(theta_mom);
/* Now do the jump conditions */
if(model==Model::solCx)
{
@@ -138,6 +137,7 @@ int main()
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
@@ -210,13 +210,13 @@ int main()
double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
Resid_x[i][j]=Rx;
- zx_new[i][j]=zx[i][j]-theta_mom*Rx/dRx_dzx;
+ zx[i][j]-=theta*Rx/dRx_dzx;
}
}
}
}
- if(sweep%100==0)
+ if(sweep%1000==0 || sweep>236000)
{
std::stringstream ss;
ss << "zx_resid" << sweep;
@@ -293,6 +293,7 @@ int main()
const double zx_avg=(zx[i][j] + zx[i][j-1] +
zx[i+1][j] + zx[i+1][j-1])/4;
+ double theta(theta_mom);
/* Now do the jump conditions */
if(model==Model::solCx)
{
@@ -302,6 +303,7 @@ int main()
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
@@ -391,13 +393,13 @@ int main()
double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
Resid_y[i][j]=Ry;
- zy_new[i][j]=zy[i][j]-theta_mom*Ry/dRy_dzy;
+ zy[i][j]-=theta*Ry/dRy_dzy;
}
}
}
}
- if(sweep%100==0)
+ if(sweep%1000==0 || sweep>236000)
{
std::stringstream ss;
ss << "zy_resid" << sweep;
@@ -425,6 +427,8 @@ 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);
+
/* Apply the jump condition */
double x((i+0.5)*h);
if(model==Model::solCx)
@@ -433,6 +437,7 @@ int main()
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
@@ -471,11 +476,11 @@ int main()
Resid_p[i][j]=Rc;
- dp[i][j]=-theta_c*Rc/dRc_dp;
- p_new[i][j]=p[i][j]+dp[i][j];
+ dp[i][j]=-theta*Rc/dRc_dp;
+ p[i][j]+=dp[i][j];
}
- if(sweep%100==0)
+ if(sweep%1000==0 || sweep>236000)
{
std::stringstream ss;
ss << "p_resid" << sweep;
@@ -509,7 +514,7 @@ int main()
double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
- zx_new[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);
}
@@ -533,14 +538,14 @@ int main()
double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
- zy_new[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);
}
max_p_resid=std::max(fabs(Resid_p[i][j]),max_p_resid);
}
- if(sweep%100==0)
+ if(sweep%1000==0 || sweep>236000)
{
std::cout << iteration << " "
<< sweep << " "
@@ -570,34 +575,22 @@ int main()
write_vtk(ss.str(),Nx,NN,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
- // && fabs(max_p_resid-max_p_resid_previous)/max_p_resid < tolerance)
- // {
- // std::cout << "Solved "
- // << sweep << " "
- // << tolerance << " "
- // << fabs(max_x_resid-max_x_resid_previous)/max_x_resid << " "
- // << fabs(max_y_resid-max_y_resid_previous)/max_y_resid << " "
- // << fabs(max_p_resid-max_p_resid_previous)/max_p_resid << " "
- // << "\n";
- // break;
- // }
+ 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
+ && fabs(max_p_resid-max_p_resid_previous)/max_p_resid < tolerance)
+ {
+ std::cout << "Solved "
+ << sweep << " "
+ << tolerance << " "
+ << fabs(max_x_resid-max_x_resid_previous)/max_x_resid << " "
+ << fabs(max_y_resid-max_y_resid_previous)/max_y_resid << " "
+ << fabs(max_p_resid-max_p_resid_previous)/max_p_resid << " "
+ << "\n";
+ break;
+ }
max_x_resid_previous=max_x_resid;
max_y_resid_previous=max_y_resid;
max_p_resid_previous=max_p_resid;
-
- 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];
- }
-
}
}
for(int i=0;i<Nx+1;++i)
More information about the CIG-COMMITS
mailing list