[cig-commits] commit: Seems to have the right initial conditions for solCx with constant viscosity and sideways gravity
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:13 PST 2012
changeset: 14:b41ecaa82b4d
user: Walter Landry <wlandry at caltech.edu>
date: Sat Dec 03 20:47:20 2011 -0800
files: constants.hxx initial.cxx main.cxx
description:
Seems to have the right initial conditions for solCx with constant viscosity and sideways gravity
diff -r a48ee28dbd74 -r b41ecaa82b4d constants.hxx
--- a/constants.hxx Mon Nov 28 14:32:56 2011 -0800
+++ b/constants.hxx Sat Dec 03 20:47:20 2011 -0800
@@ -1,6 +1,6 @@ const int N(64);
const int N(64);
// const double max_eta=1e10;
-const double max_eta=1.1;
+const double max_eta=1;
const double min_eta=1;
const double eta_jump=max_eta-min_eta;
#include <cmath>
diff -r a48ee28dbd74 -r b41ecaa82b4d initial.cxx
--- a/initial.cxx Mon Nov 28 14:32:56 2011 -0800
+++ b/initial.cxx Sat Dec 03 20:47:20 2011 -0800
@@ -2,6 +2,8 @@
#include "Model.hxx"
#include "constants.hxx"
#include "write_vtk.hxx"
+
+#include <cmath>
void initial(const Model &model, double zx[N+1][N], double zy[N][N+1],
double log_etax[N+1][N], double log_etay[N][N+1],
@@ -9,40 +11,52 @@ void initial(const Model &model, double
{
const double A(min_eta*(max_eta-min_eta)/(max_eta+min_eta));
std::complex<double> phi, psi, dphi, v;
+ const double v0(0.0), tau_0(0.0);
for(int i=0;i<N+1;++i)
for(int j=0;j<N+1;++j)
{
if(i<N && j<N)
{
- p[i][j]=0;
- if(model==Model::inclusion)
+ double x((i+0.5)*h), y((j+0.5)*h), r2(x*x+y*y);
+ switch(model)
{
- double x((i+0.5)*h), y((j+0.5)*h), r2(x*x+y*y);
+ case Model::inclusion:
+ p[i][j]=0;
if(r2>=r2_inclusion)
{
p[i][j]=-2*A*(r2_inclusion/r2)*2*epsilon*(2*x*x/r2 -1);
}
+ break;
+ case Model::solCx:
+ {
+ const double pi(atan(1.0)*4);
+ const double kx=pi*x;
+ const double csh(std::cosh(kx)), cs(std::cos(kx)),
+ snh(std::sinh(kx)), sn(std::sin(kx));
+
+ p[i][j]=-2*pi*((v0 + tau_0)*csh + (sn + snh)/4)*std::cos(pi*y);
+ }
+ break;
+ default:
+ abort();
+ break;
}
}
if(j<N)
{
double x(i*h), y((j+0.5)*h), r2(x*x+y*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;
- fx[i][j]=0;
- if(model==Model::solKz)
+ // fx[i][j]=0;
+ fx[i][j]=pi*pi*cos(pi*y)*cos(pi*x);
+ switch(model)
{
+ case Model::solKz:
log_etax[i][j]=y*log_max_eta;
- }
- else if(model==Model::solCx)
- {
+ break;
+ case Model::solCx:
if(x>middle)
{
log_etax[i][j]=log_max_eta;
@@ -51,9 +65,16 @@ void initial(const Model &model, double
{
log_etax[i][j]=log_min_eta;
}
- }
- else if(model==Model::sinker)
- {
+ {
+ const double pi(atan(1.0)*4);
+ const double kx=pi*x;
+ const double csh(std::cosh(kx)), cs(std::cos(kx)),
+ snh(std::sinh(kx)), sn(std::sin(kx));
+ zx[i][j]=(tau_0*snh - (v0+tau_0)*kx*csh
+ - (cs - csh + kx*snh)/4)*std::cos(pi*y);
+ }
+ break;
+ case Model::sinker:
if(x>.4 && x < .6 && y>.2 && y<.4)
{
log_etax[i][j]=log_max_eta;
@@ -62,10 +83,8 @@ void initial(const Model &model, double
{
log_etax[i][j]=log_min_eta;
}
- }
- else if(model==Model::inclusion)
- {
- const double r2=x*x+y*y;
+ break;
+ case Model::inclusion:
std::complex<double> z(x,y);
if(r2<r2_inclusion)
{
@@ -83,41 +102,42 @@ void initial(const Model &model, double
}
v=(phi - z*conj(dphi) - conj(psi))/2.0;
zx[i][j]=v.real();
+ break;
}
}
if(i<N)
{
double x((i+0.5)*h), y(j*h), r2(x*x+y*y);
-
- // if(j==0 || j==N)
- // zy[i][j]=0;
- // else
- // zy[i][j]=5.6 + 7.8*x + 9.0*y + 2.1*x*y;
zy[i][j]=0;
- if(model==Model::solKz || model==Model::solCx)
+ switch(model)
{
- // fy[i][j]=0;
+ case Model::solKz:
fy[i][j]=sin(pi*y)*cos(pi*x);
- if(model==Model::solKz)
+ log_etay[i][j]=y*log_max_eta;
+ break;
+ case Model::solCx:
+ fy[i][j]=0;
+ // fy[i][j]=sin(pi*y)*cos(pi*x);
+ if(x>middle)
{
- log_etay[i][j]=y*log_max_eta;
+ log_etay[i][j]=log_max_eta;
}
else
{
- if(x>middle)
- {
- log_etay[i][j]=log_max_eta;
- }
- else
- {
- log_etay[i][j]=log_min_eta;
- }
+ log_etay[i][j]=log_min_eta;
}
- }
- else if(model==Model::sinker)
- {
+ {
+ const double pi(atan(1.0)*4);
+ const double kx=pi*x;
+ const double csh(std::cosh(kx)), cs(std::cos(kx)),
+ snh(std::sinh(kx)), sn(std::sin(kx));
+ zy[i][j]=(v0*csh + (v0+tau_0)*kx*snh + (kx*csh/2 - sn/2)/2)
+ *std::sin(pi*y);
+ }
+ break;
+ case Model::sinker:
if(x>.4 && x < .6 && y>.2 && y<.4)
{
log_etay[i][j]=log_max_eta;
@@ -128,10 +148,8 @@ void initial(const Model &model, double
log_etay[i][j]=log_min_eta;
fy[i][j]=0;
}
- }
- else if(model==Model::inclusion)
- {
- const double r2=x*x+y*y;
+ break;
+ case Model::inclusion:
std::complex<double> z(x,y);
if(r2<r2_inclusion)
{
@@ -149,10 +167,13 @@ void initial(const Model &model, double
}
v=(phi - z*conj(dphi) - conj(psi))/2.0;
zy[i][j]=v.imag();
+ break;
}
}
}
write_vtk("p_initial",N,N,p);
write_vtk("zx_initial",N+1,N,zx);
write_vtk("zy_initial",N,N,zy);
+ write_vtk("log_etax",N+1,N,log_etax);
+ write_vtk("log_etay",N,N,log_etay);
}
diff -r a48ee28dbd74 -r b41ecaa82b4d main.cxx
--- a/main.cxx Mon Nov 28 14:32:56 2011 -0800
+++ b/main.cxx Sat Dec 03 20:47:20 2011 -0800
@@ -1,5 +1,6 @@
#include <iostream>
#include <fstream>
+#include <sstream>
#include "Model.hxx"
#include "constants.hxx"
@@ -18,78 +19,69 @@ int main()
with z, it is easy to compute v=z/eta. */
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];
+ p[N][N], dp[N][N], div[N][N], fx[N+1][N], fy[N][N+1];
/* Initial conditions */
- const int max_iterations=1;
- const int n_sweeps=1;
- // const int n_sweeps=10000;
+ const int max_iterations=2;
+ int n_sweeps=1;
const double theta_mom=1.0;
const double theta_c=1.0;
const double tolerance=1.0e-6;
- Model model(Model::inclusion);
+ Model model(Model::solCx);
initial(model,zx,zy,log_etax,log_etay,p,fx,fy);
double Resid_p[N][N], Resid_x[N+1][N], Resid_y[N][N+1];
+ double dzx_x_correction, dzx_xx_correction, dzx_yx_correction,
+ dzy_xx_correction, dzy_xy_correction, dp_x_correction;
+
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];
+ double vx[N], dvx_y[N+1], dvx_yy[N],
+ vy[N+1], dvy_y[N], dvy_yy[N+1];
if(model==Model::solCx)
{
int i(middle/h);
double dx=middle/h-i;
+ /* vx, dvx_y, dvx_yy */
+ for(int j=0;j<N;++j)
+ vx[j]=zx[i][j]*(1-dx)*exp(-log_etax[i][j])
+ + zx[i+1][j]*dx*exp(-log_etax[i+1][j]);
+
+ dvx_y[0]=dvx_y[N]=0;
+ for(int j=1;j<N+1;++j)
+ dvx_y[j]=(vx[j]-vx[j-1])/h;
+
+ for(int j=0;j<N;++j)
+ dvx_yy[j]=(dvx_y[j+1]-dvx_y[j])/h;
+
+ /* vy, dvy_y, dvy_yy */
for(int j=0;j<N+1;++j)
- vy_jump[j]=zy[i][j]*(1-dx)*exp(-log_etay[i][j])
+ vy[j]=zy[i][j]*(1-dx)*exp(-log_etay[i][j])
+ zy[i+1][j]*dx*exp(-log_etay[i+1][j]);
for(int j=0;j<N;++j)
- {
- 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[j]=(vy[j+1]-vy[j])/h;
- 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)
- {
- 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;
- }
+ dvy_yy[0]=dvy_yy[N]=0;
+ for(int j=1;j<N;++j)
+ dvy_yy[j]=(dvy_y[j+1]-dvy_y[j])/h;
}
/* Smoothing sweeps */
+ // if(iteration==max_iterations-1)
+ // n_sweeps=1;
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)
{
@@ -135,6 +127,10 @@ int main()
double dp_x=(p[i][j]-p[i-1][j])/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]
@@ -160,40 +156,41 @@ int main()
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)
+ if(x-h<middle && x+h>middle)
{
- double dx=(x>middle) ? (x+h/2-middle) : (x-h/2-middle);
+ double dx=(x>middle) ? (x+h-middle) : (middle-(x-h));
+ double zx_jump=vx[j];
+ double dzx_x_jump=dvy_y[j];
+ double p_jump=-2*dvy_y[j];
+ double dzx_xx_jump=-dvx_yy[j] + p_jump;
+ dzx_xx_correction=eta_jump*(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
+ // dzx_xx-=eta_jump*(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
+ if(n_sweeps==1)
+ std::cout << "dzx_xx "
+ << i << " "
+ << j << " "
+ << dzx_xx << " "
+ << dzx_xx_correction << " "
+ << eta_jump*zx_jump << " "
+ << zx[i][j] << " "
+ << zx[i+1][j] << " "
+ << zx[i-1][j] << " "
+ << "\n";
+ if(dx>h/2)
+ {
+ dx-=h/2;
+ double dp_x_jump=2*dvx_yy[j];
+ dp_x_correction=eta_jump*(p_jump + dx*dp_x_jump)/h;
+ // dp_x-=eta_jump*(p_jump + dx*dp_x_jump)/h;
- /* 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(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);
+ dzy_xy_correction=eta_jump*h*(dvy_y[j] - dx*dvx_yy[j])/(h*h);
+ // dzy_xy-=eta_jump*h*(dvy_y[j] - dx*dvx_yy[j])/(h*h);
+ }
}
}
@@ -207,6 +204,8 @@ int main()
double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
+ // Resid_x[i][j]=dzx_xx + dzx_yy - dp_x -fx[i][j];
+ // Resid_x[i][j]=fx[i][j];
Resid_x[i][j]=Rx;
// zx[i][j]-=theta_mom*Rx/dRx_dzx;
}
@@ -214,8 +213,12 @@ int main()
}
}
- write_vtk("zx_resid",N+1,N,Resid_x);
-
+ if(sweep%100==0)
+ {
+ std::stringstream ss;
+ ss << "zx_resid" << sweep;
+ write_vtk(ss.str(),N+1,N,Resid_x);
+ }
/* zy */
for(int rb=0;rb<2;++rb)
@@ -285,37 +288,29 @@ int main()
zx[i+1][j] + zx[i+1][j-1])/4;
/* Now do the jump conditions */
-
- double x((i+0.5)*h), y(j*h);
-
if(model==Model::solCx)
{
+ double x((i+0.5)*h), y(j*h);
dlog_etayy=dlog_etay=dlog_etaxx=dlog_etax=
dlog_etayx=0;
- if(x+h/2>middle && x-h/2>middle)
+ if(x-h<middle && x+h>middle)
{
- double dx=(x>middle) ? (x+h/2-middle) : (x-h/2-middle);
+ double dx=(x>middle) ? (x+h-middle) : (middle-(x-h));
+ double zy_jump=vy[j];
+ double dzy_x_jump=-dvx_y[j];
+ double dzy_xx_jump=-dvy_yy[j];
+ dzy_xx_correction=eta_jump*(zy_jump + dx*dzy_x_jump + dx*dx*dzy_xx_jump/2)/(h*h);
+ // dzy_xx-=eta_jump*(zy_jump + dx*dzy_x_jump + dx*dx*dzy_xx_jump/2)/(h*h);
- 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);
+ if(dx>h/2)
+ {
+ dx-=h/2;
+ dzx_yx_correction=eta_jump*(dvx_y[j] + dx*dvy_yy[j])/h;
+ // dzx_yx-=eta_jump*(dvx_y[j] + dx*dvy_yy[j])/h;
+ }
}
}
-
/* Compute the residual and update zy */
@@ -334,7 +329,12 @@ int main()
}
}
- write_vtk("zy_resid",N,N,Resid_y);
+ if(sweep%100==0)
+ {
+ std::stringstream ss;
+ ss << "zy_resid" << sweep;
+ write_vtk(ss.str(),N,N,Resid_y);
+ }
/* Pressure update */
@@ -366,11 +366,11 @@ int main()
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];
+ double zx_jump=vx[j];
+ double dzx_x_jump=dvy_y[j];
- dzx_x-=eta_jump*(zx_jump + dx*dzx_jump + dx*dx*ddzx_jump/2)/h;
+ dzx_x_correction=eta_jump*(zx_jump + dx*dzx_x_jump)/h;
+ // dzx_x-=eta_jump*(zx_jump + dx*dzx_x_jump)/h;
}
}
@@ -397,8 +397,12 @@ int main()
// p[i][j]+=dp[i][j];
}
- write_vtk("p_resid",N,N,Resid_p);
-
+ if(sweep%100==0)
+ {
+ std::stringstream ss;
+ ss << "p_resid" << sweep;
+ write_vtk(ss.str(),N,N,Resid_p);
+ }
/* Velocity Fix */
double max_x(0), max_y(0), max_p(0);
@@ -423,7 +427,7 @@ int main()
zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
- max_x=std::max(std::fabs(zx[i][j]),max_x);
+ max_x=std::max(std::fabs(Resid_x[i][j]),max_x);
}
/* Fix vy */
if(j>0)
@@ -442,17 +446,39 @@ int main()
zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
- max_y=std::max(fabs(zy[i][j]),max_y);
+ max_y=std::max(fabs(Resid_y[i][j]),max_y);
}
- max_p=std::max(fabs(p[i][j]),max_p);
+ max_p=std::max(fabs(Resid_p[i][j]),max_p);
}
if(sweep%100==0)
+ {
std::cout << sweep << " "
<< max_x << " "
<< max_y << " "
<< max_p << " "
<< "\n";
+
+ std::stringstream ss;
+ ss << "zx" << sweep;
+ write_vtk(ss.str(),N+1,N,zx);
+ ss.str("");
+ ss << "zy" << sweep;
+ write_vtk(ss.str(),N,N,zy);
+ ss.str("");
+ ss << "p" << sweep;
+ write_vtk(ss.str(),N,N,p);
+
+ for(int i=0;i<N;++i)
+ for(int j=0;j<N;++j)
+ {
+ 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(),N,N,div);
+ }
if(fabs(max_x-max_x_previous)/max_x < tolerance
&& fabs(max_y-max_y_previous)/max_y < tolerance
@@ -472,4 +498,14 @@ int main()
max_p_previous=max_p;
}
}
+ for(int i=0;i<N+1;++i)
+ for(int j=0;j<N+1;++j)
+ {
+ if(i<N)
+ zx[i][j]/=exp(log_etax[i][j]);
+ if(j<N)
+ zy[i][j]/=exp(log_etay[i][j]);
+ }
+ write_vtk("vx",N+1,N,zx);
+ write_vtk("vy",N,N,zy);
}
More information about the CIG-COMMITS
mailing list