[cig-commits] commit: Add waf and inclusion
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:08 PST 2012
changeset: 10:d2666ac3ef8e
user: Walter Landry <wlandry at caltech.edu>
date: Mon Nov 28 12:07:54 2011 -0800
files: Model.hxx constants.hxx initial.cxx main.cxx waf-1.6.3 wscript
description:
Add waf and inclusion
diff -r 35ce498dd25d -r d2666ac3ef8e Model.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Model.hxx Mon Nov 28 12:07:54 2011 -0800
@@ -0,0 +1,2 @@
+enum class Model {solKz, solCx, sinker, inclusion};
+
diff -r 35ce498dd25d -r d2666ac3ef8e constants.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/constants.hxx Mon Nov 28 12:07:54 2011 -0800
@@ -0,0 +1,14 @@
+const int N(64);
+// const double max_eta=1e10;
+const double max_eta=1.1;
+const double min_eta=1;
+const double eta_jump=max_eta-min_eta;
+#include <cmath>
+const double log_max_eta=std::log(max_eta);
+const double log_min_eta=std::log(min_eta);
+const double middle=0.4;
+const double h(1.0/N);
+const double pi(atan(1.0)*4);
+
+const double r2_inclusion=0.1;
+const double epsilon=1;
diff -r 35ce498dd25d -r d2666ac3ef8e initial.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/initial.cxx Mon Nov 28 12:07:54 2011 -0800
@@ -0,0 +1,154 @@
+#include <complex>
+#include "Model.hxx"
+#include "constants.hxx"
+
+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],
+ double p[N][N], double fx[N+1][N], double fy[N][N+1])
+{
+ const double A(min_eta*(max_eta-min_eta)/(max_eta+min_eta));
+ std::complex<double> phi, psi, dphi, v;
+
+ 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);
+ if(r2>=r2_inclusion)
+ {
+ p[i][j]=-2*A*(r2_inclusion/r2)*2*epsilon*(2*x*x/r2 -1);
+ }
+ }
+ }
+ 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)
+ {
+ log_etax[i][j]=y*log_max_eta;
+ }
+ else if(model==Model::solCx)
+ {
+ if(x>middle)
+ {
+ log_etax[i][j]=log_max_eta;
+ }
+ else
+ {
+ log_etax[i][j]=log_min_eta;
+ }
+ }
+ else if(model==Model::sinker)
+ {
+ if(x>.4 && x < .6 && y>.2 && y<.4)
+ {
+ log_etax[i][j]=log_max_eta;
+ }
+ else
+ {
+ log_etax[i][j]=log_min_eta;
+ }
+ }
+ else if(model==Model::inclusion)
+ {
+ const double r2=x*x+y*y;
+ std::complex<double> z(x,y);
+ if(r2<r2_inclusion)
+ {
+ log_etax[i][j]=log_max_eta;
+ phi=0;
+ dphi=0;
+ psi=-4*epsilon*(max_eta*min_eta/(min_eta+max_eta))*z;
+ }
+ else
+ {
+ log_etax[i][j]=log_min_eta;
+ phi=-2*epsilon*A*r2_inclusion/z;
+ dphi=-phi/z;
+ psi=-2*epsilon*(min_eta*z+A*r2_inclusion*r2_inclusion/(z*z*z));
+ }
+ v=(phi - z*conj(dphi) - conj(psi))/2.0;
+ zx[i][j]=v.real();
+ }
+ }
+
+ 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)
+ {
+ // fy[i][j]=0;
+ fy[i][j]=sin(pi*y)*cos(pi*x);
+ if(model==Model::solKz)
+ {
+ log_etay[i][j]=y*log_max_eta;
+ }
+ else
+ {
+ if(x>middle)
+ {
+ log_etay[i][j]=log_max_eta;
+ }
+ else
+ {
+ log_etay[i][j]=log_min_eta;
+ }
+ }
+ }
+ else if(model==Model::sinker)
+ {
+ if(x>.4 && x < .6 && y>.2 && y<.4)
+ {
+ log_etay[i][j]=log_max_eta;
+ fy[i][j]=1;
+ }
+ else
+ {
+ log_etay[i][j]=log_min_eta;
+ fy[i][j]=0;
+ }
+ }
+ else if(model==Model::inclusion)
+ {
+ const double r2=x*x+y*y;
+ std::complex<double> z(x,y);
+ if(r2<r2_inclusion)
+ {
+ log_etay[i][j]=log_max_eta;
+ phi=0;
+ dphi=0;
+ psi=-4*epsilon*(max_eta*min_eta/(min_eta+max_eta))*z;
+ }
+ else
+ {
+ log_etay[i][j]=log_min_eta;
+ phi=-2*epsilon*A*r2_inclusion/z;
+ dphi=-phi/z;
+ psi=-2*epsilon*(min_eta*z+A*r2_inclusion*r2_inclusion/(z*z*z));
+ }
+ v=(phi - z*conj(dphi) - conj(psi))/2.0;
+ zy[i][j]=v.imag();
+ }
+ }
+ }
+}
diff -r 35ce498dd25d -r d2666ac3ef8e main.cxx
--- a/main.cxx Thu Aug 04 18:55:11 2011 -0700
+++ b/main.cxx Mon Nov 28 12:07:54 2011 -0800
@@ -1,13 +1,15 @@
-#include <cmath>
#include <iostream>
#include <fstream>
-enum class Model {solKz, solCx, sinker};
+#include "Model.hxx"
+#include "constants.hxx"
+
+extern 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],
+ double p[N][N], double fx[N+1][N], double fy[N][N+1]);
int main()
{
- const int N(64);
-
/* 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 for the momentum residual and dlog_etax for
@@ -19,115 +21,16 @@ int main()
/* Initial conditions */
- // 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);
- const double pi(atan(1.0)*4);
+ const int max_iterations=1;
+ const int n_sweeps=1;
+ // const int n_sweeps=10000;
const double theta_mom=1.0;
const double theta_c=1.0;
const double tolerance=1.0e-6;
- Model model(Model::solCx);
+ Model model(Model::inclusion);
- for(int i=0;i<N+1;++i)
- for(int j=0;j<N+1;++j)
- {
- if(j<N)
- {
- 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;
-
- zx[i][j]=0;
-
- fx[i][j]=0;
- if(model==Model::solKz)
- {
- log_etax[i][j]=y*log_max_eta;
- }
- else if(model==Model::solCx)
- {
- if(x>middle)
- {
- log_etax[i][j]=log_max_eta;
- }
- else
- {
- log_etax[i][j]=log_min_eta;
- }
- }
- else if(model==Model::sinker)
- {
- if(x>.4 && x < .6 && y>.2 && y<.4)
- {
- log_etax[i][j]=log_max_eta;
- }
- else
- {
- log_etax[i][j]=log_min_eta;
- }
- }
- }
-
- if(i<N)
- {
- double x((i+0.5)*h), y(j*h);
-
- // 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)
- {
- // fy[i][j]=0;
- fy[i][j]=sin(pi*y)*cos(pi*x);
- if(model==Model::solKz)
- {
- log_etay[i][j]=y*log_max_eta;
- }
- else
- {
- if(x>middle)
- {
- log_etay[i][j]=log_max_eta;
- }
- else
- {
- log_etay[i][j]=log_min_eta;
- }
- }
- }
- else if(model==Model::sinker)
- {
- if(x>.4 && x < .6 && y>.2 && y<.4)
- {
- log_etay[i][j]=log_max_eta;
- fy[i][j]=1;
- }
- else
- {
- log_etay[i][j]=log_min_eta;
- fy[i][j]=0;
- }
- }
- }
-
- if(i<N && j<N)
- p[i][j]=0;
- }
+ initial(model,zx,zy,log_etax,log_etay,p,fx,fy);
for(int iteration=0;iteration<max_iterations;++iteration)
{
@@ -185,6 +88,7 @@ int main()
{
/* zx */
+ std::ofstream outfilex("zx");
for(int rb=0;rb<2;++rb)
{
for(int j=0;j<N;++j)
@@ -192,15 +96,11 @@ int main()
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
+ if(i!=0 && i!=N)
{
/* Do the finite difference parts */
- const double x(i*h), y((j+0/5)*h);
+ 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);
@@ -281,14 +181,6 @@ int main()
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)
{
@@ -301,14 +193,6 @@ int main()
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 */
@@ -321,7 +205,12 @@ int main()
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;
+
+ outfilex << x << " "
+ << y << " "
+ << Rx << " "
+ << "\n";
}
}
}
@@ -330,6 +219,7 @@ int main()
/* zy */
+ std::ofstream outfiley("zy");
for(int rb=0;rb<2;++rb)
{
for(int j=0;j<N+1;++j)
@@ -337,11 +227,7 @@ int main()
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
+ if(j!=0 && j!=N)
{
/* Do the finite difference parts */
@@ -402,12 +288,12 @@ int main()
/* Now do the jump conditions */
+ double x((i+0.5)*h), y(j*h);
+
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)
{
@@ -443,7 +329,12 @@ int main()
double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
- zy[i][j]-=theta_mom*Ry/dRy_dzy;
+ // zy[i][j]-=theta_mom*Ry/dRy_dzy;
+
+ outfiley << x << " "
+ << y << " "
+ << Ry << " "
+ << "\n";
}
}
}
@@ -451,6 +342,7 @@ int main()
/* Pressure update */
+ std::ofstream outfilep("p");
for(int j=0;j<N;++j)
for(int i=0;i<N;++i)
{
@@ -470,9 +362,10 @@ int main()
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 */
+ double x((i+0.5)*h), y((j+0.5)*h);
if(model==Model::solCx)
{
- double x((i+0.5)*h);
dlog_etax=dlog_etay=dlog_etaxx=dlog_etayy=0;
if(x+h/2>middle && x-h/2<middle)
@@ -504,7 +397,12 @@ int main()
+ dRc_dzym*dRmm_dp/dRm_dz;
dp[i][j]=-theta_c*Rc/dRc_dp;
- p[i][j]+=dp[i][j];
+ // p[i][j]+=dp[i][j];
+
+ outfilep << x << " "
+ << y << " "
+ << Rc << " "
+ << "\n";
}
/* Velocity Fix */
@@ -580,28 +478,28 @@ int main()
max_p_previous=max_p;
}
}
- {
- std::ofstream outfile("zx");
- for(int j=0;j<N;++j)
- for(int i=0;i<N+1;++i)
- outfile << i << " " << j << " "
- << zx[i][j] << " "
- << "\n";
- }
- {
- std::ofstream outfile("zy");
- for(int j=0;j<N+1;++j)
- for(int i=0;i<N;++i)
- outfile << i << " " << j << " "
- << zy[i][j] << " "
- << "\n";
- }
- {
- std::ofstream outfile("p");
- for(int j=0;j<N;++j)
- for(int i=0;i<N;++i)
- outfile << i << " " << j << " "
- << p[i][j] << " "
- << "\n";
- }
+ // {
+ // std::ofstream outfile("zx");
+ // for(int j=0;j<N;++j)
+ // for(int i=0;i<N+1;++i)
+ // outfile << i << " " << j << " "
+ // << zx[i][j] << " "
+ // << "\n";
+ // }
+ // {
+ // std::ofstream outfile("zy");
+ // for(int j=0;j<N+1;++j)
+ // for(int i=0;i<N;++i)
+ // outfile << i << " " << j << " "
+ // << zy[i][j] << " "
+ // << "\n";
+ // }
+ // {
+ // std::ofstream outfile("p");
+ // for(int j=0;j<N;++j)
+ // for(int i=0;i<N;++i)
+ // outfile << i << " " << j << " "
+ // << p[i][j] << " "
+ // << "\n";
+ // }
}
diff -r 35ce498dd25d -r d2666ac3ef8e waf-1.6.3
Binary file waf-1.6.3 has changed
diff -r 35ce498dd25d -r d2666ac3ef8e wscript
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/wscript Mon Nov 28 12:07:54 2011 -0800
@@ -0,0 +1,18 @@
+# top = '.'
+# out = 'build'
+
+def options(opt):
+ opt.load('compiler_cxx')
+
+def configure(conf):
+ conf.load('compiler_cxx')
+def build(bld):
+ bld.program(
+ features = ['cxx','cprogram'],
+ source = ['main.cxx',
+ 'initial.cxx'],
+
+ target = 'Gamr_disc',
+ # cxxflags = ['-std=c++0x','-g','-D_GLIBCXX_DEBUG'],
+ cxxflags = ['-std=c++0x','-O3'],
+ )
More information about the CIG-COMMITS
mailing list