[cig-commits] commit: Had to fix a sign error in the pressure term. Now the complete thing seems to work.
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:00 PST 2012
changeset: 4:3d2e8ad57043
user: Walter Landry <wlandry at caltech.edu>
date: Thu Aug 04 02:06:39 2011 -0700
files: main.cxx
description:
Had to fix a sign error in the pressure term. Now the complete thing seems to work.
diff -r a5b02f86d9f8 -r 3d2e8ad57043 main.cxx
--- a/main.cxx Thu Aug 04 01:02:43 2011 -0700
+++ b/main.cxx Thu Aug 04 02:06:39 2011 -0700
@@ -16,8 +16,8 @@ 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=1;
- // 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;
@@ -32,12 +32,12 @@ int main()
{
double x(i*h), y((j+0.5)*h);
- // if(i==0 || i==N)
- // zx[i][j]=0;
- // else
- // zx[i][j]=1.2 + 3.4*x + 4.5*y;
+ if(i==0 || i==N)
+ zx[i][j]=0;
+ else
+ zx[i][j]=1.2 + 3.4*x + 4.5*y;
- zx[i][j]=0;
+ // zx[i][j]=0;
fx[i][j]=0;
if(model==Model::solCx)
@@ -120,198 +120,216 @@ int main()
for(int sweep=0;sweep<n_sweeps;++sweep)
{
- // /* zx */
- // for(int rb=0;rb<2;++rb)
- // {
- // for(int j=0;j<N;++j)
- // {
- // int i_min=(j+rb)%2;
- // for(int i=i_min;i<N+1;i+=2)
- // {
- // if(i==0 || i==N)
- // {
- // zx[i][j]=0;
- // }
- // else
- // {
- // /* Do the finite difference parts */
+ /* zx */
+ for(int rb=0;rb<2;++rb)
+ {
+ for(int j=0;j<N;++j)
+ {
+ int i_min=(j+rb)%2;
+ for(int i=i_min;i<N+1;i+=2)
+ {
+ if(i==0 || i==N)
+ {
+ zx[i][j]=0;
+ }
+ else
+ {
+ /* Do the finite difference parts */
- // /* Derivatives of z */
- // double dzx_xx=(zx[i-1][j] - 2*zx[i][j] + zx[i+1][j])/(h*h);
+ /* Derivatives of z */
+ double dzx_xx=(zx[i-1][j] - 2*zx[i][j] + zx[i+1][j])/(h*h);
- // double dzy_xy=((zy[i][j+1] - zy[i-1][j+1])
- // - (zy[i][j] - zy[i-1][j]))/(h*h);
+ double dzy_xy=((zy[i][j+1] - zy[i-1][j+1])
+ - (zy[i][j] - zy[i-1][j]))/(h*h);
- // double dzx_x=(zx[i+1][j] - zx[i-1][j])/(2*h);
- // double dzy_y=(zy[i][j+1] - zy[i][j]
- // + zy[i-1][j+1] - zy[i-1][j])/(2*h);
+ double dzx_x=(zx[i+1][j] - zx[i-1][j])/(2*h);
+ double dzy_y=(zy[i][j+1] - zy[i][j]
+ + zy[i-1][j+1] - zy[i-1][j])/(2*h);
- // double dzx_y, dzx_yy;
- // if(j==0)
- // {
- // dzx_yy=(-zx[i][j] + zx[i][j+1])/(h*h);
- // dzx_y=(zx[i][j+1] - zx[i][j])/(2*h);
- // }
- // else if(j==N-1)
- // {
- // dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
- // dzx_y=(zx[i][j] - zx[i][j-1])/(2*h);
- // }
- // else
- // {
- // dzx_yy=(zx[i][j-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
- // dzx_y=(zx[i][j+1] - zx[i][j-1])/(2*h);
- // }
+ double dzx_y, dzx_yy;
+ if(j==0)
+ {
+ dzx_yy=(-zx[i][j] + zx[i][j+1])/(h*h);
+ dzx_y=(zx[i][j+1] - zx[i][j])/(2*h);
+ }
+ else if(j==N-1)
+ {
+ dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
+ dzx_y=(zx[i][j] - zx[i][j-1])/(2*h);
+ }
+ else
+ {
+ dzx_yy=(zx[i][j-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
+ dzx_y=(zx[i][j+1] - zx[i][j-1])/(2*h);
+ }
- // /* Derivatives of p and eta. Eta is collocated
- // with z. Maybe eta should be cell & vertex
- // centered? It would simplify most of these
- // differences, but it would make dlog_etaxy
- // more messy. Also, with it collocated with z,
- // it is easy to compute v=z/eta. */
+ /* Derivatives of p and eta. Eta is collocated
+ with z. Maybe eta should be cell & vertex
+ centered? It would simplify most of these
+ differences, but it would make dlog_etaxy
+ more messy. Also, with it collocated with z,
+ it is easy to compute v=z/eta. */
- // double dp_x=(p[i][j]-p[i-1][j])/h;
+ double dp_x=(p[i][j]-p[i-1][j])/h;
- // double dlog_etaxx=
- // ((log_etay[i][j+1] + log_etay[i][j])/2
- // - 2*log_etax[i][j]
- // + (log_etay[i-1][j+1] + log_etay[i-1][j])/2)/(h*h/4);
- // double dlog_etax=
- // ((log_etay[i][j+1] + log_etay[i][j])/2
- // - (log_etay[i-1][j+1] + log_etay[i-1][j])/2)/h;
+ double dlog_etaxx=
+ ((log_etay[i][j+1] + log_etay[i][j])/2
+ - 2*log_etax[i][j]
+ + (log_etay[i-1][j+1] + log_etay[i-1][j])/2)/(h*h/4);
+ double dlog_etax=
+ ((log_etay[i][j+1] + log_etay[i][j])/2
+ - (log_etay[i-1][j+1] + log_etay[i-1][j])/2)/h;
- // double dlog_etayy=
- // ((log_etay[i][j+1] + log_etay[i-1][j+1])/2
- // - 2*log_etax[i][j]
- // + (log_etay[i][j] + log_etay[i-1][j])/2)/(h*h/4);
+ double dlog_etayy=
+ ((log_etay[i][j+1] + log_etay[i-1][j+1])/2
+ - 2*log_etax[i][j]
+ + (log_etay[i][j] + log_etay[i-1][j])/2)/(h*h/4);
- // double dlog_etay=
- // ((log_etay[i][j+1] + log_etay[i-1][j+1])
- // - (log_etay[i][j] + log_etay[i-1][j]))/(2*h);
+ double dlog_etay=
+ ((log_etay[i][j+1] + log_etay[i-1][j+1])
+ - (log_etay[i][j] + log_etay[i-1][j]))/(2*h);
- // double dlog_etaxy=
- // ((log_etay[i][j+1] - log_etay[i-1][j+1])
- // - (log_etay[i][j] - log_etay[i-1][j]))/(h*h);
+ double dlog_etaxy=
+ ((log_etay[i][j+1] - log_etay[i-1][j+1])
+ - (log_etay[i][j] - log_etay[i-1][j]))/(h*h);
- // const double zy_avg=(zy[i][j] + zy[i-1][j] +
- // zy[i][j+1] + zy[i-1][j+1])/4;
+ const double zy_avg=(zy[i][j] + zy[i-1][j] +
+ zy[i][j+1] + zy[i-1][j+1])/4;
- // /* Now do the jump conditions */
+ /* Now do the jump conditions */
- // /* Compute the residual and update zx */
+ /* Compute the residual and update zx */
- // 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
- // + dlog_etaxy*zy_avg + dlog_etax*dzy_y
- // + dp_x - fx[i][j];
+ 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
+ + dlog_etaxy*zy_avg + dlog_etax*dzy_y
+ - dp_x - fx[i][j];
- // double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
+ double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
- // zx[i][j]-=theta_mom*Rx/dRx_dzx;
- // }
- // }
- // }
- // }
+ zx[i][j]-=theta_mom*Rx/dRx_dzx;
+ }
+ }
+ }
+ }
- // /* zy */
+ /* zy */
- // 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 */
+ 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 */
- // /* Derivatives of z */
- // double dzy_yy=(zy[i][j-1] - 2*zy[i][j] + zy[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 dzx_yx=((zx[i+1][j] - zx[i+1][j-1])
- // - (zx[i][j] - zx[i][j-1]))/(h*h);
+ double dzx_yx=((zx[i+1][j] - zx[i+1][j-1])
+ - (zx[i][j] - zx[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 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 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_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);
+ }
- // /* 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. */
+ /* 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 dp_y=(p[i][j]-p[i][j-1])/h;
+ double dp_y=(p[i][j]-p[i][j-1])/h;
- // 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_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_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_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_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_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_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_etayx=
+ ((log_etax[i+1][j] - log_etax[i+1][j-1])
+ - (log_etax[i][j] - log_etax[i][j-1]))/(h*h);
- // const double zx_avg=(zx[i][j] + zx[i][j-1] +
- // zx[i+1][j] + zx[i+1][j-1])/4;
+ 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 */
+ /* Now do the jump conditions */
- // /* Compute the residual and update zy */
+ /* 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];
+ 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];
- // double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+ double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
- // zy[i][j]-=theta_mom*Ry/dRy_dzy;
+ // if(i>30 && i<33 && j==63)
+ // std::cout << "Ry "
+ // << i << " "
+ // << j << " "
+ // << zy[i][j] << " "
+ // << Ry/dRy_dzy << " "
+ // << Ry << " "
+ // << dRy_dzy << " "
+ // << 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] << " "
+ // << "\n";
- // }
- // }
- // }
- // }
+ zy[i][j]-=theta_mom*Ry/dRy_dzy;
+
+ }
+ }
+ }
+ }
/* Pressure update */
@@ -355,21 +373,58 @@ int main()
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 << i << " " << j << " "
+ // << dzx_x << " "
+ // << dzy_y << " "
+ // << zx_avg << " "
+ // << zy_avg << " "
+ // << dlog_etax << " "
+ // << dlog_etay << " "
+ // << dlog_etaxx << " "
+ // << dlog_etayy << " "
+ // << Rc << " "
+ // << dp[i][j] << " "
+ // << "\n";
+
+ // std::cout << "sweep "
+ // << sweep << " "
+ // << i << " "
+ // << j << " "
+ // << zx[i][j] << " "
+ // << zy[i][j] << " "
+ // << p[i][j] << " "
+ // << "\n";
}
std::cout << sweep << " "
- << zy[1][1] << " "
- << zy[10][10] << " "
+ << zx[32][1] << " "
+ << zx[32][32] << " "
+ << zx[32][62] << " "
+ << zx[32][63] << " "
+ << zy[32][1] << " "
<< zy[32][32] << " "
- << zy[63][63] << " "
+ << zy[32][62] << " "
+ << zy[32][63] << " "
+ << p[32][1] << " "
+ << p[32][32] << " "
+ << p[32][62] << " "
+ << p[32][63] << " "
<< "\n";
+
+
+ // std::cout << sweep << " "
+ // << zx[1][1] << " "
+ // << zx[10][10] << " "
+ // << zx[32][32] << " "
+ // << zx[63][63] << " "
+ // << zy[1][1] << " "
+ // << zy[10][10] << " "
+ // << zy[32][32] << " "
+ // << zy[63][63] << " "
+ // << p[1][1] << " "
+ // << p[10][10] << " "
+ // << p[32][32] << " "
+ // << p[63][63] << " "
+ // << "\n";
}
}
More information about the CIG-COMMITS
mailing list