[cig-commits] commit: Compute dC_dv ahead of time instead of getting it from compute_Cxyz.
Mercurial
hg at geodynamics.org
Thu Mar 22 11:44:58 PDT 2012
changeset: 108:115c2a1443ff
user: Walter Landry <wlandry at caltech.edu>
date: Thu Mar 22 11:36:26 2012 -0700
files: compute_Cxyz.cxx compute_Cxyz.hxx compute_dC_dv.cxx main.cxx wscript
description:
Compute dC_dv ahead of time instead of getting it from compute_Cxyz.
diff -r 199e688f48c2 -r 115c2a1443ff compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Thu Mar 22 11:35:43 2012 -0700
+++ b/compute_Cxyz.cxx Thu Mar 22 11:36:26 2012 -0700
@@ -11,10 +11,10 @@ void compute_Cxyz(const double zx[Nx+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)
+ double &C)
{
C=0;
- dC_dz=0;
+ double dC_dz=0;
const Interface interface(d,i,j,distx,disty);
if(!interface.intersects())
diff -r 199e688f48c2 -r 115c2a1443ff compute_Cxyz.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_Cxyz.hxx Thu Mar 22 11:36:26 2012 -0700
@@ -0,0 +1,15 @@
+#ifndef GAMR_COMPUTE_CXYZ_HXX
+#define GAMR_COMPUTE_CXYZ_HXX
+
+#include "constants.hxx"
+void compute_Cxyz(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);
+
+
+#endif
diff -r 199e688f48c2 -r 115c2a1443ff compute_dC_dv.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_dC_dv.cxx Thu Mar 22 11:36:26 2012 -0700
@@ -0,0 +1,41 @@
+#include "constants.hxx"
+#include "compute_Cxyz.hxx"
+
+void compute_dC_dv(double zx[Nx+1][Ny],
+ 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],
+ double dCx[Nx+1][Ny],
+ double dCy[Nx][Ny+1])
+{
+ for(int i=0;i<Nx+1;++i)
+ for(int j=0;j<Ny;++j)
+ zx[i][j]=0;
+
+ for(int i=0;i<Nx;++i)
+ for(int j=0;j<Ny+1;++j)
+ zy[i][j]=0;
+
+ for(int j=0;j<Ny;++j)
+ for(int i=1;i<Nx;++i)
+ {
+ zx[i][j]=1;
+ compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
+ 0,i,j,dCx[i][j]);
+ zx[i][j]=0;
+ }
+
+ for(int j=1;j<Ny;++j)
+ for(int i=0;i<Nx;++i)
+ {
+ zy[i][j]=1;
+ compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
+ 1,i,j,dCy[i][j]);
+ zy[i][j]=0;
+ }
+
+}
+
+
diff -r 199e688f48c2 -r 115c2a1443ff main.cxx
--- a/main.cxx Thu Mar 22 11:35:43 2012 -0700
+++ b/main.cxx Thu Mar 22 11:36:26 2012 -0700
@@ -5,21 +5,12 @@
#include "Model.hxx"
#include "constants.hxx"
#include "write_vtk.hxx"
+#include "compute_Cxyz.hxx"
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 distx[Nx+1][Ny], double disty[Nx][Ny+1]);
-
-extern void compute_Cxyz(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);
extern void compute_Cp(const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -29,6 +20,15 @@ extern void compute_Cp(const double zx[N
const double disty[Nx][Ny+1],
const int &i, const int &j,
double &Cp);
+
+extern void compute_dC_dv(double zx[Nx+1][Ny],
+ 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],
+ double dCx[Nx+1][Ny],
+ double dCy[Nx][Ny+1]);
int main()
@@ -41,13 +41,13 @@ int main()
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],
- distx[Nx+1][Ny], disty[Nx][Ny+1];
+ distx[Nx+1][Ny], disty[Nx][Ny+1], dCx[Nx+1][Ny], dCy[Nx][Ny+1];
double zx_new[Nx+1][Ny], zy_new[Nx][Ny+1], p_new[Nx][Ny];
const bool jacobi(false);
/* Initial conditions */
- int n_sweeps=1000000;
+ int n_sweeps=1001;
const double theta_mom=0.7;
const double theta_continuity=1.0;
const double tolerance=1.0e-6;
@@ -57,6 +57,10 @@ int main()
Model model(Model::solCx);
initial(model,zx,zy,log_etax,log_etay,p,fx,fy,distx,disty);
+
+ compute_dC_dv(zx,zy,log_etax,log_etay,distx,disty,dCx,dCy);
+
+ double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
for(int i=0;i<Nx+1;++i)
for(int j=0;j<Ny;++j)
@@ -69,8 +73,6 @@ int main()
for(int i=0;i<Nx;++i)
for(int j=0;j<Ny;++j)
p[i][j]=0;
-
- double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
for(int sweep=0;sweep<n_sweeps;++sweep)
{
@@ -157,9 +159,9 @@ int main()
/* Compute the residual and update zx */
- double Cx, dCx_dzx;
+ double Cx;
compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
- 0,i,j,Cx,dCx_dzx);
+ 0,i,j,Cx);
double Rx=2*dzx_xx + dzx_yy + dzy_xy
+ 2*(dlog_etaxx*zx[i][j] + dlog_etax*dzx_x)
@@ -170,7 +172,7 @@ int main()
double dRx_dzx=-6/(h*h);
if(j==0 || j==Ny-1)
dRx_dzx=-5/(h*h);
- dRx_dzx+=2*dlog_etaxx + dlog_etayy + dCx_dzx;
+ dRx_dzx+=2*dlog_etaxx + dlog_etayy + dCx[i][j];
Resid_x[i][j]=Rx;
if(!jacobi)
@@ -268,9 +270,10 @@ int main()
/* Compute the residual and update zy */
- double Cy, dCy_dzy;
+ double Cy;
compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
- 1,i,j,Cy,dCy_dzy);
+ 1,i,j,Cy);
+
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
@@ -280,7 +283,7 @@ int main()
double dRy_dzy=-6/(h*h);
if(i==0 || i==Nx-1)
dRy_dzy=-5/(h*h);
- dRy_dzy+=2*dlog_etayy + dlog_etaxx + dCy_dzy;
+ dRy_dzy+=2*dlog_etayy + dlog_etaxx + dCy[i][j];
Resid_y[i][j]=Ry;
if(!jacobi)
@@ -388,14 +391,14 @@ int main()
dlog_etaxx=dlog_etayy=0;
}
- double Cx, dCx_dzx;
+ double Cx;
compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
- 0,i,j,Cx,dCx_dzx);
+ 0,i,j,Cx);
double dRx_dzx=-6/(h*h);
if(j==0 || j==Ny-1)
dRx_dzx=-5/(h*h);
- dRx_dzx+=2*dlog_etaxx + dlog_etayy + dCx_dzx;
+ dRx_dzx+=2*dlog_etaxx + dlog_etayy + dCx[i][j];
if(!jacobi)
zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
@@ -422,14 +425,14 @@ int main()
dlog_etaxx=dlog_etayy=0;
}
- double Cy, dCy_dzy;
+ double Cy;
compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
- 1,i,j,Cy,dCy_dzy);
+ 1,i,j,Cy);
double dRy_dzy=-6/(h*h);
if(i==0 || i==Nx-1)
dRy_dzy=-5/(h*h);
- dRy_dzy+=2*dlog_etayy + dlog_etaxx + dCy_dzy;
+ dRy_dzy+=2*dlog_etayy + dlog_etaxx + dCy[i][j];
if(!jacobi)
zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
diff -r 199e688f48c2 -r 115c2a1443ff wscript
--- a/wscript Thu Mar 22 11:35:43 2012 -0700
+++ b/wscript Thu Mar 22 11:36:26 2012 -0700
@@ -11,7 +11,8 @@ def build(bld):
'compute_v_on_interface/compute_v_on_interface.cxx',
'compute_v_on_interface/compute_dv_dtt.cxx',
'compute_Cxyz.cxx',
- 'compute_Cp.cxx'],
+ 'compute_Cp.cxx',
+ 'compute_dC_dv.cxx'],
target = 'Gamr_disc',
# cxxflags = ['-std=c++0x','-g','-Wall','-Drestrict='],
More information about the CIG-COMMITS
mailing list