[cig-commits] commit: When computing dv_dyy at the interface, weight the high viscosity side more than the low viscosity side. Seems to work up to 10^10 viscosity contrasts
Mercurial
hg at geodynamics.org
Fri Mar 2 10:47:25 PST 2012
changeset: 69:a3ec1ee1fbd1
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Mar 02 10:47:18 2012 -0800
files: compute_Cx.cxx compute_Cy.cxx compute_v_on_interface.cxx constants.hxx main.cxx
description:
When computing dv_dyy at the interface, weight the high viscosity side more than the low viscosity side. Seems to work up to 10^10 viscosity contrasts
diff -r 0416056f296c -r a3ec1ee1fbd1 compute_Cx.cxx
--- a/compute_Cx.cxx Fri Mar 02 09:45:24 2012 -0800
+++ b/compute_Cx.cxx Fri Mar 02 10:47:18 2012 -0800
@@ -42,7 +42,7 @@ void compute_Cx(const Model &model, cons
double eta=exp(log_etax[i][j]);
double delta=(h-std::fabs(dx))/h;
- double dvx_yy_dv=-delta/(h*h);
+ double dvx_yy_dv=-2*eta*eta/((h*h)*(min_eta*min_eta + max_eta*max_eta));
double dzy_yx_dv=-eta_jump*dvx_yy_dv;
double dvy_y_dv=-dzy_yx_dv*h/(min_eta+max_eta);
double dzx_x_dv=eta_jump*dvy_y_dv;
diff -r 0416056f296c -r a3ec1ee1fbd1 compute_Cy.cxx
--- a/compute_Cy.cxx Fri Mar 02 09:45:24 2012 -0800
+++ b/compute_Cy.cxx Fri Mar 02 10:47:18 2012 -0800
@@ -39,7 +39,7 @@ void compute_Cy(const Model &model, cons
double eta=exp(log_etay[i][j]);
double delta=(h-std::fabs(dx))/h;
- double dvy_yy_dv=-delta/(h*h);
+ double dvy_yy_dv=-2*eta*eta/((h*h)*(min_eta*min_eta + max_eta*max_eta));
double dzx_yx_dv=-eta_jump*dvy_yy_dv;
double dvx_y_dv=(-h/(min_eta+max_eta))*dzx_yx_dv;
double dzy_x_dv=-eta_jump*dvx_y_dv;
diff -r 0416056f296c -r a3ec1ee1fbd1 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Fri Mar 02 09:45:24 2012 -0800
+++ b/compute_v_on_interface.cxx Fri Mar 02 10:47:18 2012 -0800
@@ -78,63 +78,47 @@ double compute_dv_yy(const double z[][ny
const double log_eta[][ny], const double &dx,
const int i, const int j)
{
- double dv_yy_p1, dv_yy_p2;
+ double dv_yy_p;
if(j==0)
{
- 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);
+ dv_yy_p=(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_p=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/(h*h);
}
else if(j==ny-1)
{
- 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);
+ dv_yy_p=(-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_p=(-vel(z,log_eta,i+1,j) + vel(z,log_eta,i+1,j-1))/(h*h);
}
else
{
- 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,j-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,j-1))/(h*h);
+ dv_yy_p=(vel(z,log_eta,i+1,j+1) - 2*vel(z,log_eta,i+1,j)
+ + vel(z,log_eta,i+1,j-1))/(h*h);
}
- const double dv_yy_p=(1-dx)*dv_yy_p1 + dx*dv_yy_p2;
- double dv_yy_m0, dv_yy_m1;
+ double dv_yy_m;
if(j==0)
{
- 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);
+ dv_yy_m=(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_m=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/(h*h);
}
else if(j==ny-1)
{
- 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);
+ dv_yy_m=(-vel(z,log_eta,i,0) - 2*vel(z,log_eta,i,j)
+ + vel(z,log_eta,i,j-1))/(h*h);
+ // dv_yy_m=(-vel(z,log_eta,i,j) + vel(z,log_eta,i,j-1))/(h*h);
}
else
{
- dv_yy_m0=(vel(z,log_eta,i,j+1) - 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,j+1) - 2*vel(z,log_eta,i-1,j)
- + vel(z,log_eta,i-1,j-1))/(h*h);
+ dv_yy_m=(vel(z,log_eta,i,j+1) - 2*vel(z,log_eta,i,j)
+ + vel(z,log_eta,i,j-1))/(h*h);
}
- double dv_yy_m=(1-dx)*dv_yy_m1 + dx*dv_yy_m0;
- return (dv_yy_p + dv_yy_m)/2;
+ double eta_p(exp(log_eta[i+1][j])), eta_m(exp(log_eta[i][j]));
+
+ return (dv_yy_p*eta_p*eta_p + dv_yy_m*eta_m*eta_m)/(eta_m*eta_m + eta_p*eta_p);
}
template<int ny>
diff -r 0416056f296c -r a3ec1ee1fbd1 constants.hxx
--- a/constants.hxx Fri Mar 02 09:45:24 2012 -0800
+++ b/constants.hxx Fri Mar 02 10:47:18 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=100;
+const double max_eta=1e10;
const double eta_jump=max_eta-min_eta;
#include <cmath>
const double log_max_eta=std::log(max_eta);
diff -r 0416056f296c -r a3ec1ee1fbd1 main.cxx
--- a/main.cxx Fri Mar 02 09:45:24 2012 -0800
+++ b/main.cxx Fri Mar 02 10:47:18 2012 -0800
@@ -60,28 +60,31 @@ int main()
/* Initial conditions */
const int max_iterations=1;
- int n_sweeps=1000000;
+ int n_sweeps=100000;
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=10;
+ const int output_interval=1000000000;
Model model(Model::solCx);
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)
+ // for(int i=0;i<Nx+1;++i)
+ // for(int j=0;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)
+ // zx[i][j]=0;
+ // if(i<Nx)
+ // zy[i][j]=0;
// if(j<Ny && i<Nx)
- // p[i][j]=p_new[i][j];
+ // p[i][j]=0;
// }
+ // // zx[static_cast<int>(middle/h)+1][Ny/8]=1;
+ // zy[static_cast<int>(middle/h)+1][Ny/8]=1;
+
double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
@@ -195,11 +198,9 @@ int main()
+ 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;
double dRx_dzx=
- // -6/(h*h) + 2*dlog_etaxx + dlog_etayy;
-6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
Resid_x[i][j]=Rx;
@@ -208,18 +209,16 @@ int main()
else
zx_new[i][j]=zx[i][j]-theta_mom*Rx/dRx_dzx;
-
- 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";
-
+ // if(zx[i][j]!=0)
+ // std::cout << "Rx "
+ // << i << " "
+ // << j << " "
+ // << middle << " "
+ // << i*h << " "
+ // << zx[i][j] << " "
+ // << Cx << " "
+ // << dCx_dzx << " "
+ // << "\n";
}
}
}
@@ -324,11 +323,9 @@ int main()
+ 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;
double dRy_dzy=
- // -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)
@@ -355,6 +352,18 @@ int main()
zy[i][j]-=theta_mom*Ry/dRy_dzy;
else
zy_new[i][j]=zy[i][j]-theta_mom*Ry/dRy_dzy;
+
+ // if(zy[i][j]!=0)
+ // std::cout << "Rx "
+ // << i << " "
+ // << j << " "
+ // << middle << " "
+ // << i*h << " "
+ // << zy[i][j] << " "
+ // << Cy << " "
+ // << dCy_dzy << " "
+ // << "\n";
+
}
}
}
@@ -403,7 +412,6 @@ int main()
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;
double dRc_dzxp=1/h - dlog_etax/2;
@@ -473,8 +481,7 @@ int main()
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;
- // -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
+ -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);
@@ -505,8 +512,7 @@ int main()
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;
- // -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy;
+ -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);
@@ -530,7 +536,7 @@ int main()
p[i][j]=p_new[i][j];
}
- if(sweep%10==0)
+ if(sweep%1000==0)
{
std::cout << "sweep "
<< iteration << " "
More information about the CIG-COMMITS
mailing list