[cig-commits] commit: Move update of correction term to a separate function that operates outside of the individual sweep
Mercurial
hg at geodynamics.org
Tue Feb 21 16:47:29 PST 2012
changeset: 57:47c9472c5cd1
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Tue Feb 21 16:46:52 2012 -0800
files: compute_corrections.cxx constants.hxx main.cxx
description:
Move update of correction term to a separate function that operates outside of the individual sweep
diff -r fa2edb9def70 -r 47c9472c5cd1 compute_corrections.cxx
--- a/compute_corrections.cxx Tue Feb 21 16:46:04 2012 -0800
+++ b/compute_corrections.cxx Tue Feb 21 16:46:52 2012 -0800
@@ -52,7 +52,7 @@ void compute_corrections(const Model &mo
if(x>middle)
dzx_xx_correction*=-1;
- Cx_new=dzx_xx_correction;
+ Cx_new=2*dzx_xx_correction;
if(x+h/2>middle && x-h/2<middle)
{
diff -r fa2edb9def70 -r 47c9472c5cd1 constants.hxx
--- a/constants.hxx Tue Feb 21 16:46:04 2012 -0800
+++ b/constants.hxx Tue Feb 21 16:46:52 2012 -0800
@@ -1,16 +1,16 @@ const int N(64);
-const int N(64);
+const int N(512);
const int Nx(N);
const int Ny(N);
const double min_eta=1;
-const double max_eta=1;
+const double max_eta=2;
const double eta_jump=max_eta-min_eta;
#include <cmath>
const double log_max_eta=std::log(max_eta);
const double log_min_eta=std::log(min_eta);
// const double middle=(25 + 0.00000001)/64;
// const double middle=(25 + 32/Nx + 0.00000001)/64;
-// const double middle=(25 - 32/2048.0)/64;
-const double middle=0.4;
+const double middle=(25 - 1/32.0)/64;
+// const double middle=0.4;
const double h(1.0/N);
const double pi(atan(1.0)*4);
diff -r fa2edb9def70 -r 47c9472c5cd1 main.cxx
--- a/main.cxx Tue Feb 21 16:46:04 2012 -0800
+++ b/main.cxx Tue Feb 21 16:46:52 2012 -0800
@@ -35,10 +35,10 @@ int main()
/* Initial conditions */
const int max_iterations=1;
- int n_sweeps=10000;
+ int n_sweeps=1;
const double theta_mom=0.7;
const double theta_continuity=1.0;
- const double theta_correction=0.01;
+ const double theta_correction=1.0;
const double tolerance=1.0e-6;
Model model(Model::solCx);
@@ -142,6 +142,27 @@ int main()
/* Compute the residual and update zx */
+
+ // if(j==Ny/8 && i*h>middle-3*h && i*h<middle+3*h)
+ // std::cout << "Rx "
+ // << i << " "
+ // << j << " "
+ // << i*h << " "
+ // << middle << " "
+ // << dzx_xx << " "
+ // << Cx[i][j] << " "
+ // << (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]) << " "
+ // << (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]) << " "
+ // << "\n";
+
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
@@ -151,17 +172,20 @@ int main()
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;
}
}
}
}
- if(sweep%1000==0 || sweep>236000)
+ if(sweep%1000==0)
{
std::stringstream ss;
ss << "zx_resid" << sweep;
write_vtk(ss.str(),Nx+1,N,Resid_x);
+ ss.str("");
+ ss << "Cx" << sweep;
+ write_vtk(ss.str(),Nx+1,N,Cx);
}
/* zy */
@@ -242,28 +266,42 @@ int main()
dlog_etayx=0;
}
+ // if(j==Ny/8 && i*h>middle-3*h && i*h<middle+3*h)
+ // std::cout << "Ry "
+ // << i << " "
+ // << j << " "
+ // << i*h+h/2 << " "
+ // << middle << " "
+ // << dzy_xx << " "
+ // << Cy[i][j] << " "
+ // << dzy_xx+Cy[i][j] << " "
+ // << "\n";
+
/* Compute the residual and update zy */
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];
+ - dp_y - fy[i][j] + Cy[i][j];
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;
}
}
}
}
- if(sweep%1000==0 || sweep>236000)
+ if(sweep%1000==0)
{
std::stringstream ss;
ss << "zy_resid" << sweep;
write_vtk(ss.str(),Nx,N,Resid_y);
+ ss.str("");
+ ss << "Cy" << sweep;
+ write_vtk(ss.str(),Nx,N,Cy);
}
/* Pressure update */
@@ -315,14 +353,17 @@ 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)
+ if(sweep%1000==0)
{
std::stringstream ss;
ss << "p_resid" << sweep;
write_vtk(ss.str(),Nx,N,Resid_p);
+ ss.str("");
+ ss << "Cp" << sweep;
+ write_vtk(ss.str(),Nx,N,Cp);
}
/* Velocity Fix */
@@ -352,7 +393,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);
}
@@ -376,14 +417,14 @@ 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);
}
max_p_resid=std::max(fabs(Resid_p[i][j]),max_p_resid);
}
- if(sweep%1000==0 || sweep>236000)
+ if(sweep%1000==0)
{
std::cout << iteration << " "
<< sweep << " "
@@ -413,9 +454,12 @@ int main()
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
- && fabs(max_p_resid-max_p_resid_previous)/max_p_resid < tolerance)
+ // 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)
+ if(max_x_resid < tolerance
+ && max_y_resid < tolerance
+ && max_p_resid < tolerance)
{
std::cout << "Solved "
<< sweep << " "
More information about the CIG-COMMITS
mailing list