[cig-commits] commit: solCx seems to work with viscosity jump=0
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:06 PST 2012
changeset: 9:35ce498dd25d
user: Walter Landry <wlandry at caltech.edu>
date: Thu Aug 04 18:55:11 2011 -0700
files: main.cxx
description:
solCx seems to work with viscosity jump=0
diff -r e97eff318ea1 -r 35ce498dd25d main.cxx
--- a/main.cxx Thu Aug 04 07:08:18 2011 -0700
+++ b/main.cxx Thu Aug 04 18:55:11 2011 -0700
@@ -17,12 +17,16 @@ int main()
double zx[N+1][N], zy[N][N+1], log_etax[N+1][N], log_etay[N][N+1],
p[N][N], dp[N][N], fx[N+1][N], fy[N][N+1];
-
/* Initial conditions */
- const double max_eta=1e10;
+ // const double max_eta=1e10;
+ const double max_eta=1;
+ const double min_eta=1;
+ const double eta_jump=max_eta-min_eta;
const double log_max_eta=std::log(max_eta);
+ const double log_min_eta=std::log(min_eta);
const double middle=0.4;
+ const int max_iterations=10;
// const int n_sweeps=1;
const int n_sweeps=10000;
const double h(1.0/N);
@@ -31,7 +35,7 @@ int main()
const double theta_c=1.0;
const double tolerance=1.0e-6;
- Model model(Model::solKz);
+ Model model(Model::solCx);
for(int i=0;i<N+1;++i)
for(int j=0;j<N+1;++j)
@@ -60,7 +64,7 @@ int main()
}
else
{
- log_etax[i][j]=0;
+ log_etax[i][j]=log_min_eta;
}
}
else if(model==Model::sinker)
@@ -71,7 +75,7 @@ int main()
}
else
{
- log_etax[i][j]=0;
+ log_etax[i][j]=log_min_eta;
}
}
}
@@ -102,7 +106,7 @@ int main()
}
else
{
- log_etay[i][j]=0;
+ log_etay[i][j]=log_min_eta;
}
}
}
@@ -115,7 +119,7 @@ int main()
}
else
{
- log_etay[i][j]=0;
+ log_etay[i][j]=log_min_eta;
fy[i][j]=0;
}
}
@@ -125,311 +129,457 @@ int main()
p[i][j]=0;
}
+ for(int iteration=0;iteration<max_iterations;++iteration)
+ {
+ /* Compute the boundary jumps */
+
+ double vx_jump[N], dvx_y_jump[N+1], dvx_yy_jump[N],
+ vy_jump[N+1], dvy_y_jump[N], dvy_yy_jump[N+1], dvy_yyy_jump[N];
- /* Smoothing sweeps */
+ if(model==Model::solCx)
+ {
+ int i(middle/h);
+ double dx=middle/h-i;
- double max_x_previous(0), max_y_previous(0), max_p_previous(0);
- for(int sweep=0;sweep<n_sweeps;++sweep)
- {
+ for(int j=0;j<N+1;++j)
+ vy_jump[j]=zy[i][j]*(1-dx)*exp(-log_etay[i][j])
+ + zy[i+1][j]*dx*exp(-log_etay[i+1][j]);
- /* 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)
+ vx_jump[j]=zx[i][j]*(1-dx)*exp(-log_etax[i][j])
+ + zx[i+1][j]*dx*exp(-log_etax[i+1][j]);
+
+ dvy_y_jump[j]=(vy_jump[j+1]-vy_jump[j])/h;
+ }
+
+ for(int j=0;j<N+1;++j)
+ {
+ if(j==0 || j==N)
{
- if(i==0 || i==N)
+ dvx_y_jump[j]=0;
+ /* Extrapolate the derivative to the boundary */
+ if(j==0)
+ dvy_yy_jump[j]=(vy_jump[j+2]-2*vy_jump[j+1])/(h*h);
+ else
+ dvy_yy_jump[j]=(vy_jump[j-2]-2*vy_jump[j-1])/(h*h);
+ }
+ else
+ {
+ dvx_y_jump[j]=(vx_jump[j]-vx_jump[j-1])/h;
+ dvy_yy_jump[j]=(dvy_y_jump[j]-dvy_y_jump[j-1])/h;
+ }
+ }
+
+ for(int j=0;j<N;++j)
+ {
+ dvx_yy_jump[j]=(dvx_y_jump[j+1]-dvx_y_jump[j])/h;
+ dvy_yyy_jump[j]=(dvy_yy_jump[j+1]-dvy_yy_jump[j])/h;
+ }
+ }
+
+ /* Smoothing sweeps */
+
+ double max_x_previous(0), max_y_previous(0), max_p_previous(0);
+ 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)
{
- 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);
-
- 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)
+ if(i==0 || i==N)
{
- 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);
+ zx[i][j]=0;
}
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);
+ /* Do the finite difference parts */
+
+ const double x(i*h), y((j+0/5)*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 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);
+ }
+
+ /* Derivatives of p and eta. */
+
+ 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_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_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;
+
+ /* Now do the jump conditions */
+
+ if(model==Model::solCx)
+ {
+ dlog_etaxx=dlog_etax=dlog_etayy=dlog_etay=
+ dlog_etaxy=0;
+
+ if(x+h/2>middle && x-h/2>middle)
+ {
+ double dx=(x>middle) ? (x+h/2-middle) : (x-h/2-middle);
+
+ /* dzy_xy */
+ double zy_jump=vy_jump[j+1]-vy_jump[j];
+ double dzy_jump=dvx_y_jump[j+1]-dvx_y_jump[j];
+ double ddzy_jump=-3*(dvy_yy_jump[j+1]-dvy_yy_jump[j]);
+ dzy_xy-=
+ eta_jump*(zy_jump + dx*dzy_jump + dx*dx*ddzy_jump/2)/(h*h);
+
+ /* p */
+ double p_jump=-2*dvy_y_jump[j];
+ double dp_jump=2*dvx_yy_jump[j];
+ double ddp_jump=2*dvy_yyy_jump[j];
+
+ dp_x-=eta_jump*(p_jump + dx*dp_jump + dx*dx*ddp_jump/2)/(h*h);
+
+
+ if(dvx_yy_jump[j]>100000000)
+ std::cout << " ";
+ if(dp_jump>100000000)
+ std::cout << " ";
+ if(ddp_jump>100000000)
+ std::cout << " ";
+ }
+ if(x+h>middle && x-h<middle)
+ {
+ double dx=(x>middle) ? (x+h-middle) : (x-h-middle);
+
+ double zx_jump=vx_jump[j];
+ double dzx_jump=-dvy_y_jump[j];
+ double ddzx_jump=dvx_yy_jump[j];
+
+ dzx_xx-=
+ eta_jump*(zx_jump + dx*dzx_jump + dx*dx*ddzx_jump/2)/(h*h);
+ }
+
+ if(dzx_xx>100000000)
+ std::cout << " ";
+ if(dzy_xy>100000000)
+ std::cout << " ";
+ if(dp_x>100000000)
+ std::cout << " ";
+
+ }
+
+ /* 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 dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
+
+ zx[i][j]-=theta_mom*Rx/dRx_dzx;
}
-
- /* Derivatives of p and eta. */
-
- 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_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_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;
-
- /* Now do the jump conditions */
-
-
- /* 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 dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
-
- 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)
+ for(int rb=0;rb<2;++rb)
{
- int i_min=(j+rb)%2;
- for(int i=i_min;i<N;i+=2)
+ for(int j=0;j<N+1;++j)
{
- if(j==0 || j==N)
+ int i_min=(j+rb)%2;
+ for(int i=i_min;i<N;i+=2)
{
- 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);
-
- 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)
+ if(j==0 || j==N)
{
- 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);
+ zy[i][j]=0;
}
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);
+ /* 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)
+ {
+ 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. */
+
+ 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_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_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;
+
+ /* Now do the jump conditions */
+
+ if(model==Model::solCx)
+ {
+ dlog_etayy=dlog_etay=dlog_etaxx=dlog_etax=
+ dlog_etayx=0;
+
+ double x((i+0.5)*h);
+
+ if(x+h/2>middle && x-h/2>middle)
+ {
+ double dx=(x>middle) ? (x+h/2-middle) : (x-h/2-middle);
+
+ double zx_jump=vx_jump[j+1]-vx_jump[j];
+ double dzx_jump=-(dvy_y_jump[j+1]-dvy_y_jump[j]);
+ double ddzx_jump=dvx_yy_jump[j+1]-dvx_yy_jump[j];
+ dzx_yx-=
+ eta_jump*(zx_jump + dx*dzx_jump + dx*dx*ddzx_jump/2)/(h*h);
+ }
+ if(x+h>middle && x-h<middle)
+ {
+ double dx=(x>middle) ? (x+h-middle) : (x-h-middle);
+
+ double zy_jump=vy_jump[j];
+ double dzy_jump=-dvx_y_jump[j];
+ double ddzy_jump=-3*dvy_yy_jump[j];
+
+ dzy_xx-=
+ eta_jump*(zy_jump + dx*dzy_jump + dx*dx*ddzy_jump/2)/(h*h);
+ }
+ }
+
+
+ /* 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 dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+
+ zy[i][j]-=theta_mom*Ry/dRy_dzy;
}
-
- /* Derivatives of p and eta. */
-
- 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_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_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;
-
- /* 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];
-
- double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
-
- zy[i][j]-=theta_mom*Ry/dRy_dzy;
}
}
}
- }
- /* Pressure update */
+ /* 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;
+ 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];
-
- }
-
- /* Velocity Fix */
-
- double max_x(0), max_y(0), max_p(0);
-
- for(int j=0;j<N;++j)
- for(int i=0;i<N;++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);
-
- double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
-
- zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
-
- max_x=std::max(fabs(zx[i][j]),max_x);
- }
- /* 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_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_etax[i+1][j-1])/2
- - 2*log_etay[i][j]
- + (log_etax[i][j] + log_etax[i][j-1])/2)/(h*h/4);
+ (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 dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+ double zx_avg=(zx[i+1][j] + zx[i][j])/2;
+ double zy_avg=(zy[i][j+1] + zy[i][j])/2;
- zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
+ if(model==Model::solCx)
+ {
+ double x((i+0.5)*h);
+ dlog_etax=dlog_etay=dlog_etaxx=dlog_etayy=0;
- max_y=std::max(fabs(zy[i][j]),max_y);
+ if(x+h/2>middle && x-h/2<middle)
+ {
+ double dx=(x>middle) ? (x+h/2-middle) : (x-h/2-middle);
+ double zx_jump=vx_jump[j];
+ double dzx_jump=-dvy_y_jump[j];
+ double ddzx_jump=dvx_yy_jump[j];
+
+ dzx_x-=eta_jump*(zx_jump + dx*dzx_jump + dx*dx*ddzx_jump/2)/h;
+ }
+ }
+
+ 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];
}
- max_p=std::max(fabs(p[i][j]),max_p);
- }
- if(sweep%100==0)
- std::cout << sweep << " "
- << max_x << " "
- << max_y << " "
- << max_p << " "
- << "\n";
+ /* Velocity Fix */
- if(fabs(max_x-max_x_previous)/max_x < tolerance
- && fabs(max_y-max_y_previous)/max_y < tolerance
- && fabs(max_p-max_p_previous)/max_p < tolerance)
- {
- std::cout << "Solved "
- << sweep << " "
- << tolerance << " "
- << fabs(max_x-max_x_previous)/max_x << " "
- << fabs(max_y-max_y_previous)/max_y << " "
- << fabs(max_p-max_p_previous)/max_p << " "
- << "\n";
- break;
+ double max_x(0), max_y(0), max_p(0);
+
+ for(int j=0;j<N;++j)
+ for(int i=0;i<N;++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);
+
+ double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
+
+ zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
+
+ max_x=std::max(std::fabs(zx[i][j]),max_x);
+ }
+ /* 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);
+
+ double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+
+ zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
+
+ max_y=std::max(fabs(zy[i][j]),max_y);
+ }
+ max_p=std::max(fabs(p[i][j]),max_p);
+ }
+
+ if(sweep%100==0)
+ std::cout << sweep << " "
+ << max_x << " "
+ << max_y << " "
+ << max_p << " "
+ << "\n";
+
+ if(fabs(max_x-max_x_previous)/max_x < tolerance
+ && fabs(max_y-max_y_previous)/max_y < tolerance
+ && fabs(max_p-max_p_previous)/max_p < tolerance)
+ {
+ std::cout << "Solved "
+ << sweep << " "
+ << tolerance << " "
+ << fabs(max_x-max_x_previous)/max_x << " "
+ << fabs(max_y-max_y_previous)/max_y << " "
+ << fabs(max_p-max_p_previous)/max_p << " "
+ << "\n";
+ break;
+ }
+ max_x_previous=max_x;
+ max_y_previous=max_y;
+ max_p_previous=max_p;
}
- max_x_previous=max_x;
- max_y_previous=max_y;
- max_p_previous=max_p;
}
-
{
std::ofstream outfile("zx");
for(int j=0;j<N;++j)
More information about the CIG-COMMITS
mailing list