[cig-commits] commit: eta derivatives in pressure update seem ok
Mercurial
hg at geodynamics.org
Fri Feb 10 15:59:59 PST 2012
changeset: 3:a5b02f86d9f8
user: Walter Landry <wlandry at caltech.edu>
date: Thu Aug 04 01:02:43 2011 -0700
files: main.cxx
description:
eta derivatives in pressure update seem ok
diff -r c44834c77baa -r a5b02f86d9f8 main.cxx
--- a/main.cxx Thu Aug 04 00:22:36 2011 -0700
+++ b/main.cxx Thu Aug 04 01:02:43 2011 -0700
@@ -8,7 +8,7 @@ int main()
const int N(64);
double zx[N+1][N], zy[N][N+1], log_etax[N+1][N], log_etay[N][N+1],
- p[N][N], fx[N+1][N], fy[N][N+1];
+ p[N][N], dp[N][N], fx[N+1][N], fy[N][N+1];
/* Initial conditions */
@@ -16,11 +16,12 @@ int main()
const double max_eta=2;
const double log_max_eta=std::log(max_eta);
const double middle=0.4;
- const int n_sweeps=10000;
+ const int n_sweeps=1;
+ // const int n_sweeps=10000;
const double h(1.0/N);
const double pi(atan(1.0)*4);
const double theta_mom=1.0;
-
+ const double theta_c=1.0;
Model model(Model::solCx);
@@ -116,8 +117,10 @@ int main()
p[i][j]=0;
}
+
for(int sweep=0;sweep<n_sweeps;++sweep)
{
+ // /* zx */
// for(int rb=0;rb<2;++rb)
// {
// for(int j=0;j<N;++j)
@@ -213,102 +216,154 @@ int main()
// }
- for(int rb=0;rb<2;++rb)
- {
- for(int j=0;j<N+1;++j)
- {
- int i_min=(j+rb)%2;
- for(int i=i_min;i<N;i+=2)
- {
- if(j==0 || j==N)
- {
- zy[i][j]=0;
- }
- else
- {
- /* Do the finite difference parts */
+ // /* zy */
- /* Derivatives of z */
- double dzy_yy=(zy[i][j-1] - 2*zy[i][j] + zy[i][j+1])/(h*h);
+ // for(int rb=0;rb<2;++rb)
+ // {
+ // for(int j=0;j<N+1;++j)
+ // {
+ // int i_min=(j+rb)%2;
+ // for(int i=i_min;i<N;i+=2)
+ // {
+ // if(j==0 || j==N)
+ // {
+ // zy[i][j]=0;
+ // }
+ // else
+ // {
+ // /* Do the finite difference parts */
- double dzx_yx=((zx[i+1][j] - zx[i+1][j-1])
- - (zx[i][j] - zx[i][j-1]))/(h*h);
+ // /* Derivatives of z */
+ // double dzy_yy=(zy[i][j-1] - 2*zy[i][j] + zy[i][j+1])/(h*h);
- double dzy_y=(zy[i][j+1] - zy[i][j-1])/(2*h);
- double dzx_x=(zx[i+1][j] - zx[i][j]
- + zx[i+1][j-1] - zx[i][j-1])/(2*h);
+ // double dzx_yx=((zx[i+1][j] - zx[i+1][j-1])
+ // - (zx[i][j] - zx[i][j-1]))/(h*h);
- double dzy_x, dzy_xx;
- if(i==0)
- {
- dzy_xx=(-zy[i][j] + zy[i+1][j])/(h*h);
- dzy_x=(zy[i+1][j] - zy[i][j])/(2*h);
- }
- else if(i==N-1)
- {
- dzy_xx=(zy[i-1][j] - zy[i][j])/(h*h);
- dzy_x=(zy[i][j] - zy[i-1][j])/(2*h);
- }
- else
- {
- dzy_xx=(zy[i-1][j] - 2*zy[i][j] + zy[i+1][j])/(h*h);
- dzy_x=(zy[i+1][j] - zy[i-1][j])/(2*h);
- }
+ // double dzy_y=(zy[i][j+1] - zy[i][j-1])/(2*h);
+ // double dzx_x=(zx[i+1][j] - zx[i][j]
+ // + zx[i+1][j-1] - zx[i][j-1])/(2*h);
- /* Derivatives of p and eta. Eta is collocated
- with z. Maxbe eta should be cell & vertey
- centered? It would simplifx most of these
- differences, but it would make dlog_etayx
- more messx. Also, with it collocated with z,
- it is easx to compute v=z/eta. */
+ // double dzy_x, dzy_xx;
+ // if(i==0)
+ // {
+ // dzy_xx=(-zy[i][j] + zy[i+1][j])/(h*h);
+ // dzy_x=(zy[i+1][j] - zy[i][j])/(2*h);
+ // }
+ // else if(i==N-1)
+ // {
+ // dzy_xx=(zy[i-1][j] - zy[i][j])/(h*h);
+ // dzy_x=(zy[i][j] - zy[i-1][j])/(2*h);
+ // }
+ // else
+ // {
+ // dzy_xx=(zy[i-1][j] - 2*zy[i][j] + zy[i+1][j])/(h*h);
+ // dzy_x=(zy[i+1][j] - zy[i-1][j])/(2*h);
+ // }
- double dp_y=(p[i][j]-p[i][j-1])/h;
+ // /* Derivatives of p and eta. Eta is collocated
+ // with z. Maxbe eta should be cell & vertey
+ // centered? It would simplifx most of these
+ // differences, but it would make dlog_etayx
+ // more messx. Also, with it collocated with z,
+ // it is easx to compute v=z/eta. */
- double dlog_etayy=
- ((log_etax[i+1][j] + log_etax[i][j])/2
- - 2*log_etay[i][j]
- + (log_etax[i+1][j-1] + log_etax[i][j-1])/2)/(h*h/4);
- double dlog_etay=
- ((log_etax[i+1][j] + log_etax[i][j])/2
- - (log_etax[i+1][j-1] + log_etax[i][j-1])/2)/h;
+ // double dp_y=(p[i][j]-p[i][j-1])/h;
- double dlog_etaxx=
- ((log_etax[i+1][j] + log_etax[i+1][j-1])/2
- - 2*log_etay[i][j]
- + (log_etax[i][j] + log_etax[i][j-1])/2)/(h*h/4);
+ // double dlog_etayy=
+ // ((log_etax[i+1][j] + log_etax[i][j])/2
+ // - 2*log_etay[i][j]
+ // + (log_etax[i+1][j-1] + log_etax[i][j-1])/2)/(h*h/4);
+ // double dlog_etay=
+ // ((log_etax[i+1][j] + log_etax[i][j])/2
+ // - (log_etax[i+1][j-1] + log_etax[i][j-1])/2)/h;
- double dlog_etax=
- ((log_etax[i+1][j] + log_etax[i+1][j-1])
- - (log_etax[i][j] + log_etax[i][j-1]))/(2*h);
+ // double dlog_etaxx=
+ // ((log_etax[i+1][j] + log_etax[i+1][j-1])/2
+ // - 2*log_etay[i][j]
+ // + (log_etax[i][j] + log_etax[i][j-1])/2)/(h*h/4);
- double dlog_etayx=
- ((log_etax[i+1][j] - log_etax[i+1][j-1])
- - (log_etax[i][j] - log_etax[i][j-1]))/(h*h);
+ // double dlog_etax=
+ // ((log_etax[i+1][j] + log_etax[i+1][j-1])
+ // - (log_etax[i][j] + log_etax[i][j-1]))/(2*h);
- const double zx_avg=(zx[i][j] + zx[i][j-1] +
- zx[i+1][j] + zx[i+1][j-1])/4;
+ // double dlog_etayx=
+ // ((log_etax[i+1][j] - log_etax[i+1][j-1])
+ // - (log_etax[i][j] - log_etax[i][j-1]))/(h*h);
- /* Now do the jump conditions */
+ // const double zx_avg=(zx[i][j] + zx[i][j-1] +
+ // zx[i+1][j] + zx[i+1][j-1])/4;
+ // /* Now do the jump conditions */
- /* 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];
+ // /* Compute the residual and update zy */
- double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+ // 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];
- zy[i][j]-=theta_mom*Ry/dRy_dzy;
+ // double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
- }
- }
- }
- }
+ // zy[i][j]-=theta_mom*Ry/dRy_dzy;
+ // }
+ // }
+ // }
+ // }
+ /* Pressure update */
+
+ for(int j=0;j<N;++j)
+ for(int i=0;i<N;++i)
+ {
+ double dzx_x=(zx[i+1][j] - zx[i][j])/h;
+ double dzy_y=(zy[i][j+1] - zy[i][j])/h;
+
+ double dlog_etax=(log_etax[i+1][j] - log_etax[i][j])/h;
+ double dlog_etay=(log_etay[i][j+1] - log_etay[i][j])/h;
+
+ double dlog_etaxx=
+ (log_etax[i+1][j] - (log_etay[i][j+1] + log_etay[i][j])
+ + log_etax[i][j])/(h*h/4);
+ double dlog_etayy=
+ (log_etay[i][j+1] - (log_etax[i+1][j] + log_etax[i][j])
+ + log_etay[i][j])/(h*h/4);
+
+ double zx_avg=(zx[i+1][j] + zx[i][j])/2;
+ double zy_avg=(zy[i][j+1] + zy[i][j])/2;
+
+ double Rc=dzx_x + dzy_y + dlog_etax*zx_avg + dlog_etay*zy_avg;
+
+ double dRc_dzxp=1/h - dlog_etax/2;
+ double dRc_dzxm=-1/h - dlog_etax/2;
+ double dRc_dzyp=1/h - dlog_etay/2;
+ double dRc_dzym=-1/h - dlog_etay/2;
+
+ double dRmp_dp=-1/h;
+ double dRmm_dp=1/h;
+
+ double dRm_dz=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+
+ double dRc_dp=dRc_dzxp*dRmp_dp/dRm_dz
+ + dRc_dzxm*dRmm_dp/dRm_dz
+ + dRc_dzyp*dRmp_dp/dRm_dz
+ + dRc_dzym*dRmm_dp/dRm_dz;
+
+ dp[i][j]=-theta_c*Rc/dRc_dp;
+ p[i][j]+=dp[i][j];
+
+
+ std::cout << i << " " << j << " "
+ << dlog_etax << " "
+ << dlog_etay << " "
+ << dlog_etaxx << " "
+ << dlog_etayy << " "
+ << Rc << " "
+ << dp[i][j] << " "
+ << "\n";
+ }
std::cout << sweep << " "
<< zy[1][1] << " "
More information about the CIG-COMMITS
mailing list