[cig-commits] commit: Add a level set function, but do not use it yet. Also make the initial guess really zero
Mercurial
hg at geodynamics.org
Tue Mar 13 06:44:19 PDT 2012
changeset: 90:3b23d7928e38
user: Walter Landry <wlandry at caltech.edu>
date: Tue Mar 13 06:11:34 2012 -0700
files: compute_Cp.cxx compute_Cxyz.cxx initial.cxx main.cxx
description:
Add a level set function, but do not use it yet. Also make the initial guess really zero
diff -r b3ac581506e0 -r 3b23d7928e38 compute_Cp.cxx
--- a/compute_Cp.cxx Tue Mar 13 05:40:49 2012 -0700
+++ b/compute_Cp.cxx Tue Mar 13 06:11:34 2012 -0700
@@ -8,6 +8,7 @@ void compute_Cp(const Model &model, cons
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
const int &i, const int &j,
double &Cp)
{
diff -r b3ac581506e0 -r 3b23d7928e38 compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Tue Mar 13 05:40:49 2012 -0700
+++ b/compute_Cxyz.cxx Tue Mar 13 06:11:34 2012 -0700
@@ -9,6 +9,7 @@ void compute_Cxyz(const Model &model, co
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
const int &d,
const int &i, const int &j,
double &C, double &dC_dz)
diff -r b3ac581506e0 -r 3b23d7928e38 initial.cxx
--- a/initial.cxx Tue Mar 13 05:40:49 2012 -0700
+++ b/initial.cxx Tue Mar 13 06:11:34 2012 -0700
@@ -8,7 +8,8 @@
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])
+ double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1],
+ double distx[Nx+1][Ny], double disty[Nx][Ny+1])
{
const double A(min_eta*(max_eta-min_eta)/(max_eta+min_eta));
std::complex<double> phi, psi, dphi, v;
@@ -111,6 +112,7 @@ void initial(const Model &model, double
log_etax[i][j]=y*log_max_eta;
break;
case Model::solCx:
+ distx[i][j]=x-middle;
if(x>middle)
{
log_etax[i][j]=log_max_eta;
@@ -189,6 +191,7 @@ void initial(const Model &model, double
break;
case Model::solCx:
fy[i][j]=pi*pi*sin(pi*y)*cos(pi*x);
+ disty[i][j]=x-middle;
if(x>middle)
{
log_etay[i][j]=log_max_eta;
@@ -260,6 +263,8 @@ void initial(const Model &model, double
write_vtk("p_initial",Nx,Ny,p);
write_vtk("zx_initial",Nx+1,Ny,zx);
write_vtk("zy_initial",Nx,Ny,zy);
+ write_vtk("distx",Nx+1,Ny,zx);
+ write_vtk("disty",Nx,Ny,zy);
write_vtk("log_etax",Nx+1,Ny,log_etax);
write_vtk("log_etay",Nx,Ny,log_etay);
}
diff -r b3ac581506e0 -r 3b23d7928e38 main.cxx
--- a/main.cxx Tue Mar 13 05:40:49 2012 -0700
+++ b/main.cxx Tue Mar 13 06:11:34 2012 -0700
@@ -8,12 +8,15 @@
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]);
+ double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1],
+ double distx[Nx+1][Ny], double disty[Nx][Ny+1]);
extern void compute_Cxyz(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 double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
const int &d,
const int &i, const int &j,
double &C, double &dCx_dz);
@@ -22,6 +25,8 @@ extern void compute_Cp(const Model &mode
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
const int &i, const int &j,
double &Cp);
@@ -35,7 +40,8 @@ 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], fx[Nx+1][Ny], fy[Nx][Ny+1];
+ p[Nx][Ny], dp[Nx][Ny], fx[Nx+1][Ny], fy[Nx][Ny+1],
+ distx[Nx+1][Ny], disty[Nx][Ny+1];
double zx_new[Nx+1][Ny], zy_new[Nx][Ny+1], p_new[Nx][Ny];
const bool jacobi(false);
@@ -50,7 +56,7 @@ int main()
Model model(Model::solCx);
- initial(model,zx,zy,log_etax,log_etay,p,fx,fy);
+ initial(model,zx,zy,log_etax,log_etay,p,fx,fy,distx,disty);
for(int i=0;i<Nx+1;++i)
for(int j=0;j<Ny;++j)
@@ -63,8 +69,6 @@ int main()
for(int i=0;i<Nx;++i)
for(int j=0;j<Ny;++j)
p[i][j]=0;
-
- zy[0][10]=1;
double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
@@ -154,8 +158,8 @@ int main()
/* Compute the residual and update zx */
double Cx, dCx_dzx;
- compute_Cxyz(model,zx,zy,log_etax,log_etay,0,i,j,Cx,
- dCx_dzx);
+ compute_Cxyz(model,zx,zy,log_etax,log_etay,distx,disty,
+ 0,i,j,Cx,dCx_dzx);
double Rx=2*dzx_xx + dzx_yy + dzy_xy
+ 2*(dlog_etaxx*zx[i][j] + dlog_etax*dzx_x)
@@ -265,8 +269,8 @@ int main()
/* Compute the residual and update zy */
double Cy, dCy_dzy;
- compute_Cxyz(model,zx,zy,log_etax,log_etay,1,i,j,
- Cy,dCy_dzy);
+ compute_Cxyz(model,zx,zy,log_etax,log_etay,distx,disty,
+ 1,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
@@ -323,7 +327,7 @@ int main()
}
double Cp;
- compute_Cp(model,zx,zy,log_etax,log_etay,i,j,Cp);
+ compute_Cp(model,zx,zy,log_etax,log_etay,distx,disty,i,j,Cp);
double Rc=dzx_x + dzy_y + dlog_etax*zx_avg + dlog_etay*zy_avg
+ Cp;
@@ -385,7 +389,8 @@ int main()
}
double Cx, dCx_dzx;
- compute_Cxyz(model,zx,zy,log_etax,log_etay,0,i,j,Cx,dCx_dzx);
+ compute_Cxyz(model,zx,zy,log_etax,log_etay,distx,disty,
+ 0,i,j,Cx,dCx_dzx);
double dRx_dzx=-6/(h*h);
if(j==0 || j==Ny-1)
@@ -418,7 +423,8 @@ int main()
}
double Cy, dCy_dzy;
- compute_Cxyz(model,zx,zy,log_etax,log_etay,1,i,j,Cy,dCy_dzy);
+ compute_Cxyz(model,zx,zy,log_etax,log_etay,distx,disty,
+ 1,i,j,Cy,dCy_dzy);
double dRy_dzy=-6/(h*h);
if(i==0 || i==Nx-1)
More information about the CIG-COMMITS
mailing list