[cig-commits] commit: Clean up
Mercurial
hg at geodynamics.org
Fri Mar 2 10:54:23 PST 2012
changeset: 70:cc41ec7c81ee
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Mar 02 10:54:16 2012 -0800
files: compute_correction.cxx main.cxx wscript
description:
Clean up
diff -r a3ec1ee1fbd1 -r cc41ec7c81ee compute_correction.cxx
--- a/compute_correction.cxx Fri Mar 02 10:47:18 2012 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,80 +0,0 @@
-#include "constants.hxx"
-#include "Model.hxx"
-
-extern void compute_Cx(const Model &model, const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const int &i, const int &j,
- double &Cx, double &dCx_dzx);
-
-extern void compute_Cy(const Model &model, const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const int &i, const int &j,
- double &Cy, double &dCy_dzy);
-
-extern void compute_Cp(const Model &model, const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const int &i, const int &j,
- double &Cp);
-
-void compute_correction(const Model &model, const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- double Cx[Nx+1][Ny], double dCx[Nx+1][Ny],
- double Cy[Nx][Ny+1], double dCy[Nx][Ny+1],
- double Cp[Nx][Ny],
- const double &theta_correction)
-{
- if(model==Model::solCx)
- {
- double C_new, dC_new;
- for(int j=0;j<Ny;++j)
- for(int i=0;i<Nx+1;++i)
- {
- compute_Cx(model,zx,zy,log_etax,log_etay,i,j,C_new,dC_new);
- if(theta_correction==1.0)
- {
- Cx[i][j]=C_new;
- dCx[i][j]=dC_new;
- }
- else
- {
- Cx[i][j]+=theta_correction*(C_new-Cx[i][j]);
- dCx[i][j]+=theta_correction*(dC_new-dCx[i][j]);
- }
- }
-
- for(int j=0;j<Ny+1;++j)
- for(int i=0;i<Nx;++i)
- {
- compute_Cy(model,zx,zy,log_etax,log_etay,i,j,C_new,dC_new);
- if(theta_correction==1.0)
- {
- Cy[i][j]=C_new;
- dCy[i][j]=dC_new;
- }
- else
- {
- Cy[i][j]+=theta_correction*(C_new-Cy[i][j]);
- dCy[i][j]+=theta_correction*(dC_new-dCy[i][j]);
- }
- }
-
- for(int j=0;j<Ny;++j)
- for(int i=0;i<Nx;++i)
- {
- compute_Cp(model,zx,zy,log_etax,log_etay,i,j,C_new);
- if(theta_correction==1.0)
- Cp[i][j]=C_new;
- else
- Cp[i][j]+=theta_correction*(C_new-Cp[i][j]);
- }
- }
-}
-
diff -r a3ec1ee1fbd1 -r cc41ec7c81ee main.cxx
--- a/main.cxx Fri Mar 02 10:47:18 2012 -0800
+++ b/main.cxx Fri Mar 02 10:54:16 2012 -0800
@@ -9,16 +9,6 @@ extern void initial(const Model &model,
extern void initial(const Model &model, double zx[Nx+1][Ny], double zy[Nx][Ny+1],
double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1]);
-
-// extern void compute_correction(const Model &model, const double zx[Nx+1][Ny],
-// const double zy[Nx][Ny+1],
-// const double log_etax[Nx+1][Ny],
-// const double log_etay[Nx][Ny+1],
-// double Cx[Nx+1][Ny], double dCx[Nx+1][Ny],
-// double Cy[Nx][Ny+1], double dCy[Nx][Ny+1],
-// double Cp[Nx][Ny],
-// const double &theta_correction);
-
extern void compute_Cx(const Model &model, const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -51,20 +41,16 @@ int main()
with z, it is easy to compute v=z/eta. */
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];
+ p[Nx][Ny], dp[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];
- // double Cx[Nx+1][Ny], dCx_dzx[Nx+1][Ny], Cy[Nx][Ny+1], dCy_dzy[Nx][Ny+1], Cp[Nx][Ny];
-
const bool jacobi(true);
/* Initial conditions */
- const int max_iterations=1;
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=1000000000;
@@ -72,525 +58,419 @@ int main()
initial(model,zx,zy,log_etax,log_etay,p,fx,fy);
- // for(int i=0;i<Nx+1;++i)
- // for(int j=0;j<Ny+1;++j)
- // {
- // if(j<Ny)
- // zx[i][j]=0;
- // if(i<Nx)
- // zy[i][j]=0;
- // if(j<Ny && i<Nx)
- // 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];
- // compute_correction(model,zx,zy,log_etax,log_etay,Cx,dCx_dzx,
- // Cy,dCy_dzy,Cp,1.0);
- for(int iteration=0;iteration<max_iterations;++iteration)
+ for(int sweep=0;sweep<n_sweeps;++sweep)
{
- // compute_correction(model,zx,zy,log_etax,log_etay,Cx,dCx_dzx,
- // Cy,dCy_dzy,Cp,theta_correction);
+ /* zx */
+ for(int rb=0;rb<2;++rb)
+ {
+ for(int j=0;j<Ny;++j)
+ {
+ int i_min=(j+rb)%2;
+ for(int i=i_min;i<Nx+1;i+=2)
+ {
+ if(i!=0 && i!=Nx)
+ {
+ /* Do the finite difference parts */
- /* Smoothing sweeps */
+ /* Derivatives of z */
+ double dzx_xx=
+ (zx[i-1][j] - 2*zx[i][j] + zx[i+1][j])/(h*h);
- // if(iteration==max_iterations-1)
- // n_sweeps=1;
- int sweep;
- for(sweep=0;sweep<n_sweeps;++sweep)
- {
- /* zx */
- for(int rb=0;rb<2;++rb)
- {
- for(int j=0;j<Ny;++j)
- {
- int i_min=(j+rb)%2;
- for(int i=i_min;i<Nx+1;i+=2)
- {
- if(i!=0 && i!=Nx)
+ 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_y, dzx_yy;
+ if(j==0)
{
- /* Do the finite difference parts */
+ dzx_yy=
+ (-zx[i][Ny-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
+ dzx_y=(zx[i][j+1] + zx[i][Ny-1])/(2*h);
+ // 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==Ny-1)
+ {
+ dzx_yy=
+ (zx[i][j-1] - 2*zx[i][j] - zx[i][0])/(h*h);
+ dzx_y=(-zx[i][0] - zx[i][j-1])/(2*h);
+ // 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 z */
- double dzx_xx=
- (zx[i-1][j] - 2*zx[i][j] + zx[i+1][j])/(h*h);
+ /* Derivatives of p and eta. */
- double dzy_xy=((zy[i][j+1] - zy[i-1][j+1])
- - (zy[i][j] - zy[i-1][j]))/(h*h);
+ double dp_x=(p[i][j]-p[i-1][j])/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);
+ /* This plays a lot of games to reduce the
+ size of the stencil. It may not be worth
+ it, and it makes it more difficult to
+ correct for the jumps. */
+ 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 dzx_y, dzx_yy;
- if(j==0)
- {
- dzx_yy=
- (-zx[i][Ny-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
- dzx_y=(zx[i][j+1] + zx[i][Ny-1])/(2*h);
- // 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==Ny-1)
- {
- dzx_yy=
- (zx[i][j-1] - 2*zx[i][j] - zx[i][0])/(h*h);
- dzx_y=(-zx[i][0] - zx[i][j-1])/(2*h);
- // 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 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);
- /* Derivatives of p and eta. */
+ 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 dp_x=(p[i][j]-p[i-1][j])/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);
- /* This plays a lot of games to reduce the
- size of the stencil. It may not be worth
- it, and it makes it more difficult to
- correct for the jumps. */
- 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;
+ const double zy_avg=(zy[i][j] + zy[i-1][j] +
+ zy[i][j+1] + zy[i-1][j+1])/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);
+ /* Now do the jump conditions */
+ if(model==Model::solCx)
+ {
+ dlog_etaxx=dlog_etax=dlog_etayy=dlog_etay=
+ dlog_etaxy=0;
+ }
- 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);
+ /* Compute the residual and update zx */
- 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 Cx, dCx_dzx;
+ compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,
+ dCx_dzx);
+ 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] + Cx;
- const double zy_avg=(zy[i][j] + zy[i-1][j] +
- zy[i][j+1] + zy[i-1][j+1])/4;
+ double dRx_dzx=
+ -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
- /* Now do the jump conditions */
- if(model==Model::solCx)
- {
- dlog_etaxx=dlog_etax=dlog_etayy=dlog_etay=
- dlog_etaxy=0;
- }
-
- /* Compute the residual and update zx */
-
- double Cx, dCx_dzx;
- compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,
- dCx_dzx);
- 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] + Cx;
-
- double dRx_dzx=
- -6/(h*h) + 2*dlog_etaxx + dlog_etayy + dCx_dzx;
-
- Resid_x[i][j]=Rx;
- if(!jacobi)
- zx[i][j]-=theta_mom*Rx/dRx_dzx;
- else
- zx_new[i][j]=zx[i][j]-theta_mom*Rx/dRx_dzx;
-
- // if(zx[i][j]!=0)
- // std::cout << "Rx "
- // << i << " "
- // << j << " "
- // << middle << " "
- // << i*h << " "
- // << zx[i][j] << " "
- // << Cx << " "
- // << dCx_dzx << " "
- // << "\n";
- }
+ Resid_x[i][j]=Rx;
+ if(!jacobi)
+ zx[i][j]-=theta_mom*Rx/dRx_dzx;
+ else
+ zx_new[i][j]=zx[i][j]-theta_mom*Rx/dRx_dzx;
}
}
}
+ }
- if(sweep%output_interval==0)
- {
- std::stringstream ss;
- ss << "zx_resid" << sweep;
- write_vtk(ss.str(),Nx+1,N,Resid_x);
- // if(sweep==0)
- // {
- // ss.str("");
- // ss << "Cx" << iteration;
- // write_vtk(ss.str(),Nx+1,N,Cx);
- // }
- }
- /* zy */
+ if(sweep%output_interval==0)
+ {
+ std::stringstream ss;
+ ss << "zx_resid" << sweep;
+ write_vtk(ss.str(),Nx+1,N,Resid_x);
+ }
+ /* zy */
- for(int rb=0;rb<2;++rb)
+ for(int rb=0;rb<2;++rb)
+ {
+ for(int j=0;j<Ny+1;++j)
{
- for(int j=0;j<Ny+1;++j)
+ int i_min=(j+rb)%2;
+ for(int i=i_min;i<Nx;i+=2)
{
- int i_min=(j+rb)%2;
- for(int i=i_min;i<Nx;i+=2)
+ if(j!=0 && j!=Ny)
{
- if(j!=0 && j!=Ny)
+ /* 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);
+
+ 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_x, dzy_xx;
+ if(i==0)
{
- /* Do the finite difference parts */
+ 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==Nx-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 z */
- double dzy_yy=
- (zy[i][j-1] - 2*zy[i][j] + zy[i][j+1])/(h*h);
+ /* Derivatives of p and eta. */
- double dzx_yx=((zx[i+1][j] - zx[i+1][j-1])
- - (zx[i][j] - zx[i][j-1]))/(h*h);
+ double dp_y=(p[i][j]-p[i][j-1])/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 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 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==Nx-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 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);
- /* Derivatives of p and eta. */
+ 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 dp_y=(p[i][j]-p[i][j-1])/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_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;
+ const double zx_avg=(zx[i][j] + zx[i][j-1] +
+ zx[i+1][j] + zx[i+1][j-1])/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);
+ /* Now do the jump conditions */
+ if(model==Model::solCx)
+ {
+ dlog_etayy=dlog_etay=dlog_etaxx=dlog_etax=
+ dlog_etayx=0;
+ }
- 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);
+ /* Compute the residual and update zy */
- 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 Cy, dCy_dzy;
+ compute_Cy(model,zx,zy,log_etax,log_etay,i,j,Cy,
+ dCy_dzy);
+ 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] + Cy;
- const double zx_avg=(zx[i][j] + zx[i][j-1] +
- zx[i+1][j] + zx[i+1][j-1])/4;
+ double dRy_dzy=
+ -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy;
- /* Now do the jump conditions */
- if(model==Model::solCx)
- {
- dlog_etayy=dlog_etay=dlog_etaxx=dlog_etax=
- dlog_etayx=0;
- }
-
- /* Compute the residual and update zy */
-
- double Cy, dCy_dzy;
- compute_Cy(model,zx,zy,log_etax,log_etay,i,j,Cy,
- dCy_dzy);
- 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] + Cy;
-
- double dRy_dzy=
- -6/(h*h) + 2*dlog_etayy + dlog_etaxx + dCy_dzy;
-
- // if(j==Ny/8 && i*h>middle-10*h && i*h<middle+10*h)
- // std::cout << "Ry "
- // << i << " "
- // << j << " "
- // << i*h+h/2 << " "
- // << middle << " "
-
- // // << 2*dzy_yy << " "
- // // << dzy_xx << " "
- // // << dzx_yx << " "
- // // << Ry << " "
- // // << Cy << " "
-
- // << dzy_xx + Cy << " "
-
- // << "\n";
-
-
-
- Resid_y[i][j]=Ry;
- if(!jacobi)
- 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";
-
- }
+ Resid_y[i][j]=Ry;
+ if(!jacobi)
+ zy[i][j]-=theta_mom*Ry/dRy_dzy;
+ else
+ zy_new[i][j]=zy[i][j]-theta_mom*Ry/dRy_dzy;
}
}
}
+ }
- if(sweep%output_interval==0)
+ if(sweep%output_interval==0)
+ {
+ std::stringstream ss;
+ ss << "zy_resid" << sweep;
+ write_vtk(ss.str(),Nx,N,Resid_y);
+ }
+
+ /* Pressure update */
+
+ for(int j=0;j<Ny;++j)
+ for(int i=0;i<Nx;++i)
{
- std::stringstream ss;
- ss << "zy_resid" << sweep;
- write_vtk(ss.str(),Nx,N,Resid_y);
- // if(sweep==0)
- // {
- // ss.str("");
- // ss << "Cy" << iteration;
- // write_vtk(ss.str(),Nx,N,Cy);
- // }
+ 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;
+
+ /* Apply the jump condition */
+ if(model==Model::solCx)
+ {
+ dlog_etax=dlog_etay=dlog_etaxx=dlog_etayy=0;
+ }
+
+ 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;
+
+ 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;
+
+ Resid_p[i][j]=Rc;
+
+ dp[i][j]=-theta_continuity*Rc/dRc_dp;
+
+ if(!jacobi)
+ p[i][j]+=dp[i][j];
+ else
+ p_new[i][j]=p[i][j]+dp[i][j];
}
- /* Pressure update */
+ if(sweep%output_interval==0)
+ {
+ std::stringstream ss;
+ ss << "p_resid" << sweep;
+ write_vtk(ss.str(),Nx,N,Resid_p);
+ }
- for(int j=0;j<Ny;++j)
- for(int i=0;i<Nx;++i)
+ /* Velocity Fix */
+
+ double max_x_resid(0), max_y_resid(0), max_p_resid(0);
+
+ for(int j=0;j<Ny;++j)
+ for(int i=0;i<Nx;++i)
+ {
+ /* Fix vx */
+ if(i>0)
{
- double dzx_x=(zx[i+1][j] - zx[i][j])/h;
- double dzy_y=(zy[i][j+1] - zy[i][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_etax[i+1][j] - log_etax[i][j])/h;
- double dlog_etay=(log_etay[i][j+1] - log_etay[i][j])/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);
+
+ if(model==Model::solCx)
+ {
+ dlog_etaxx=dlog_etayy=0;
+ }
+
+ double Cx, dCx_dzx;
+ 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 + dCx_dzx;
+
+ if(!jacobi)
+ zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
+ else
+ zx_new[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);
+ }
+ /* Fix vy */
+ if(j>0)
+ {
+ 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_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);
+ ((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 zx_avg=(zx[i+1][j] + zx[i][j])/2;
- double zy_avg=(zy[i][j+1] + zy[i][j])/2;
-
- /* Apply the jump condition */
if(model==Model::solCx)
{
- dlog_etax=dlog_etay=dlog_etaxx=dlog_etayy=0;
+ dlog_etaxx=dlog_etayy=0;
}
- 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;
-
- 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;
-
- Resid_p[i][j]=Rc;
-
- dp[i][j]=-theta_continuity*Rc/dRc_dp;
+ double Cy, dCy_dzy;
+ 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 + dCy_dzy;
if(!jacobi)
- p[i][j]+=dp[i][j];
+ zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
else
- p_new[i][j]=p[i][j]+dp[i][j];
+ zy_new[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);
}
-
- if(sweep%output_interval==0)
- {
- std::stringstream ss;
- ss << "p_resid" << sweep;
- write_vtk(ss.str(),Nx,N,Resid_p);
- // if(sweep==0)
- // {
- // ss.str("");
- // ss << "Cp" << iteration;
- // write_vtk(ss.str(),Nx,N,Cp);
- // }
+ max_p_resid=std::max(fabs(Resid_p[i][j]),max_p_resid);
}
- /* Velocity Fix */
+ if(jacobi)
+ 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];
+ }
- double max_x_resid(0), max_y_resid(0), max_p_resid(0);
+ if(sweep%1000==0)
+ {
+ std::cout << "sweep "
+ << sweep << " "
+ << max_x_resid << " "
+ << max_y_resid << " "
+ << max_p_resid << " "
+ << "\n";
+ }
+ if(sweep%output_interval==0)
+ {
+ std::stringstream ss;
+ ss << "zx" << sweep;
+ write_vtk(ss.str(),Nx+1,N,zx);
+ ss.str("");
+ ss << "zy" << sweep;
+ write_vtk(ss.str(),Nx,N,zy);
+ ss.str("");
+ ss << "p" << sweep;
+ write_vtk(ss.str(),Nx,N,p);
+ }
- for(int j=0;j<Ny;++j)
- for(int i=0;i<Nx;++i)
- {
- /* Fix vx */
- if(i>0)
- {
- 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_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);
-
- if(model==Model::solCx)
- {
- dlog_etaxx=dlog_etayy=0;
- }
-
- double Cx, dCx_dzx;
- 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 + dCx_dzx;
-
- if(!jacobi)
- zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
- else
- zx_new[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);
- }
- /* Fix vy */
- if(j>0)
- {
- 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_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);
-
- if(model==Model::solCx)
- {
- dlog_etaxx=dlog_etayy=0;
- }
-
- double Cy, dCy_dzy;
- 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 + dCy_dzy;
-
- if(!jacobi)
- zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
- else
- zy_new[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(jacobi)
- 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];
- }
-
- if(sweep%1000==0)
- {
- std::cout << "sweep "
- << iteration << " "
- << sweep << " "
- << max_x_resid << " "
- << max_y_resid << " "
- << max_p_resid << " "
- << "\n";
- }
- if(sweep%output_interval==0)
- {
- std::stringstream ss;
- ss << "zx" << sweep;
- write_vtk(ss.str(),Nx+1,N,zx);
- ss.str("");
- ss << "zy" << sweep;
- write_vtk(ss.str(),Nx,N,zy);
- ss.str("");
- ss << "p" << sweep;
- write_vtk(ss.str(),Nx,N,p);
-
- // for(int j=0;j<Ny;++j)
- // for(int i=0;i<Nx;++i)
- // {
- // div[i][j]=(zx[i+1][j]/exp(log_etax[i+1][j])-zx[i][j]/exp(log_etax[i][j]))
- // + (zy[i][j+1]/exp(log_etay[i][j+1])-zy[i][j]/exp(log_etay[i][j]));
- // }
- // ss.str("");
- // ss << "div" << sweep;
- // 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(max_x_resid < tolerance
- && max_y_resid < tolerance
- && max_p_resid < tolerance)
- {
- std::cout << "Solved "
- << sweep << " "
- << tolerance << " "
- << max_x_resid << " "
- << max_y_resid << " "
- << max_p_resid << " "
- // << 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(max_x_resid < tolerance
+ && max_y_resid < tolerance
+ && max_p_resid < tolerance)
+ {
+ std::cout << "Solved "
+ << sweep << " "
+ << tolerance << " "
+ << max_x_resid << " "
+ << max_y_resid << " "
+ << max_p_resid << " "
+ << "\n";
+ break;
}
- // if(sweep==0)
- // break;
}
write_vtk("zx",Nx+1,N,zx);
diff -r a3ec1ee1fbd1 -r cc41ec7c81ee wscript
--- a/wscript Fri Mar 02 10:47:18 2012 -0800
+++ b/wscript Fri Mar 02 10:54:16 2012 -0800
@@ -12,7 +12,6 @@ def build(bld):
source = ['main.cxx',
'initial.cxx',
'compute_v_on_interface.cxx',
- 'compute_correction.cxx',
'compute_Cx.cxx',
'compute_Cy.cxx',
'compute_Cp.cxx'],
More information about the CIG-COMMITS
mailing list