[cig-commits] commit: Move compute_Cx to compute_Cxyz
Mercurial
hg at geodynamics.org
Sun Mar 11 04:58:07 PDT 2012
changeset: 83:5fe6aa6142b8
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Sun Mar 11 04:57:54 2012 -0700
files: compute_Cx.cxx compute_Cxyz.cxx main.cxx wscript
description:
Move compute_Cx to compute_Cxyz
diff -r dec30a218834 -r 5fe6aa6142b8 compute_Cx.cxx
--- a/compute_Cx.cxx Sun Mar 11 04:46:02 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,100 +0,0 @@
-#include "constants.hxx"
-#include "Model.hxx"
-#include <iostream>
-#include "compute_v_on_interface.hxx"
-#include "FTensor.hpp"
-
-void compute_Cx(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 int &i, const int &j,
- double &Cx, double &dCx_dzx)
-{
- Cx=0;
- dCx_dzx=0;
- if(model==Model::solCx)
- {
- const double x(i*h);
- if(j<Ny && i!=0 && i!=Nx && x-h<middle && x+h>middle)
- {
- double dCx_dvx(0);
- /* xyz= the component we are correcting, Cx, Cy, or Cz.
- dir=direction of derivative
- dir2, dir3=components of mixed derivative we are correcting.
- So for zy,xy with boundary in x, dir2=y, dir3=y. */
- FTensor::Tensor1<double,2> v, dv, ddv, v_dv, dv_dv, ddv_dv,
- xyz(1,0), dir(1,0), dir2(0,1), dir3(0,1), dir4(1,0);
- compute_v_on_interface(zx,zy,log_etax,log_etay,middle,i,j,0,
- v,dv,ddv,v_dv,dv_dv,ddv_dv);
-
- FTensor::Tensor2_symmetric<double,2> ddel(0,1,0);
- const FTensor::Index<'a',2> a;
- const FTensor::Index<'b',2> b;
- const FTensor::Index<'c',2> c;
-
- double dx=(x>middle) ? (middle-(x-h))
- : (middle-(x+h));
- double z_jump=eta_jump*v(a)*xyz(a);
-
- FTensor::Tensor2<double,2,2> dz_jump;
- FTensor::Number<0> n;
- FTensor::Number<1> t;
-
- dz_jump(a,n)=-eta_jump*dv(b)*ddel(a,b);
- dz_jump(a,t)=eta_jump*dv(a);
- double dz_d_jump=dz_jump(a,b)*xyz(a)*dir(b);
-
- double p_jump=-2*eta_jump*dv(1);
-
- FTensor::Tensor1<double,2> dp_jump(2*eta_jump*ddv(n),
- -2*eta_jump*ddv(t));
-
- FTensor::Tensor3_christof<double,2,2> ddz_jump;
- ddz_jump(a,t,t)=eta_jump*ddv(a);
- ddz_jump(a,n,n)=-ddz_jump(a,t,t) + dp_jump(a);
- ddz_jump(a,n,t)=-eta_jump*ddv(b)*ddel(a,b);
-
- double ddz_dd_jump=ddz_jump(a,b,c)*xyz(a)*dir(b)*dir(c);
-
- double dz_dd_correction=
- -(z_jump - dx*dz_d_jump
- + dx*dx*ddz_dd_jump/2)/(h*h);
- const int dz_dd_sign(x>middle ? -1 : 1);
- Cx=dz_dd_sign*2*dz_dd_correction;
-
- FTensor::Tensor1<double,2> dp_dv(2*eta_jump*ddv_dv(n),
- -2*eta_jump*ddv_dv(t));
- double p_dv=-2*eta_jump*dv_dv(t);
-
- FTensor::Tensor2<double,2,2> dz_dv;
- dz_dv(a,n)=-eta_jump*dv_dv(b)*ddel(a,b);
- dz_dv(a,t)=eta_jump*dv_dv(a);
-
- FTensor::Tensor3_christof<double,2,2> ddz_dv;
- ddz_dv(a,t,t)=eta_jump*ddv_dv(a);
- ddz_dv(a,n,n)=-ddz_dv(a,t,t) + dp_dv(a);
- ddz_dv(a,n,t)=-eta_jump*ddv_dv(b)*ddel(a,b);
-
- dCx_dvx=-dz_dd_sign*2*(eta_jump*v_dv(a)*xyz(a)
- - dx*dz_dv(a,b)*xyz(a)*dir(b)
- + dx*dx*ddz_dv(a,b,c)*xyz(a)*dir(b)*dir(c)/2)
- /(h*h);
-
- if(x+h/2>middle && x-h/2<middle)
- {
- dx=(x>middle) ? ((x-h/2)-middle)
- : ((x+h/2)-middle);
-
- Cx+=(p_jump + dx*dp_jump(a)*dir(a))/h
- - (dz_jump(a,b)*dir2(a)*dir3(b)
- + dx*ddz_jump(a,b,c)*dir2(a)*dir3(b)*dir4(c))/h;
-
- dCx_dvx+=(p_dv + dx*dp_dv(a)*dir(a))/h
- - (dz_dv(a,b)*dir2(a)*dir3(b)
- + dx*ddz_dv(a,b,c)*dir2(a)*dir3(b)*dir4(c))/h;
- }
- dCx_dzx=dCx_dvx*exp(-log_etax[i][j]);
- }
- }
-}
diff -r dec30a218834 -r 5fe6aa6142b8 compute_Cxyz.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_Cxyz.cxx Sun Mar 11 04:57:54 2012 -0700
@@ -0,0 +1,100 @@
+#include "constants.hxx"
+#include "Model.hxx"
+#include <iostream>
+#include "compute_v_on_interface.hxx"
+#include "FTensor.hpp"
+
+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 int &i, const int &j,
+ double &C, double &dC_dz)
+{
+ C=0;
+ dC_dz=0;
+ if(model==Model::solCx)
+ {
+ const double x(i*h);
+ if(j<Ny && i!=0 && i!=Nx && x-h<middle && x+h>middle)
+ {
+ double dC_dv(0);
+ /* xyz= the component we are correcting, Cx, Cy, or Cz.
+ dir=direction of derivative
+ dir2, dir3=components of mixed derivative we are correcting.
+ So for zy,xy with boundary in x, dir2=y, dir3=y. */
+ FTensor::Tensor1<double,2> v, dv, ddv, v_dv, dv_dv, ddv_dv,
+ xyz(1,0), dir(1,0), dir2(0,1), dir3(0,1), dir4(1,0);
+ compute_v_on_interface(zx,zy,log_etax,log_etay,middle,i,j,0,
+ v,dv,ddv,v_dv,dv_dv,ddv_dv);
+
+ FTensor::Tensor2_symmetric<double,2> ddel(0,1,0);
+ const FTensor::Index<'a',2> a;
+ const FTensor::Index<'b',2> b;
+ const FTensor::Index<'c',2> c;
+
+ double dx=(x>middle) ? (middle-(x-h))
+ : (middle-(x+h));
+ double z_jump=eta_jump*v(a)*xyz(a);
+
+ FTensor::Tensor2<double,2,2> dz_jump;
+ FTensor::Number<0> n;
+ FTensor::Number<1> t;
+
+ dz_jump(a,n)=-eta_jump*dv(b)*ddel(a,b);
+ dz_jump(a,t)=eta_jump*dv(a);
+ double dz_d_jump=dz_jump(a,b)*xyz(a)*dir(b);
+
+ double p_jump=-2*eta_jump*dv(1);
+
+ FTensor::Tensor1<double,2> dp_jump(2*eta_jump*ddv(n),
+ -2*eta_jump*ddv(t));
+
+ FTensor::Tensor3_christof<double,2,2> ddz_jump;
+ ddz_jump(a,t,t)=eta_jump*ddv(a);
+ ddz_jump(a,n,n)=-ddz_jump(a,t,t) + dp_jump(a);
+ ddz_jump(a,n,t)=-eta_jump*ddv(b)*ddel(a,b);
+
+ double ddz_dd_jump=ddz_jump(a,b,c)*xyz(a)*dir(b)*dir(c);
+
+ double dz_dd_correction=
+ -(z_jump - dx*dz_d_jump
+ + dx*dx*ddz_dd_jump/2)/(h*h);
+ const int dz_dd_sign(x>middle ? -1 : 1);
+ C=dz_dd_sign*2*dz_dd_correction;
+
+ FTensor::Tensor1<double,2> dp_dv(2*eta_jump*ddv_dv(n),
+ -2*eta_jump*ddv_dv(t));
+ double p_dv=-2*eta_jump*dv_dv(t);
+
+ FTensor::Tensor2<double,2,2> dz_dv;
+ dz_dv(a,n)=-eta_jump*dv_dv(b)*ddel(a,b);
+ dz_dv(a,t)=eta_jump*dv_dv(a);
+
+ FTensor::Tensor3_christof<double,2,2> ddz_dv;
+ ddz_dv(a,t,t)=eta_jump*ddv_dv(a);
+ ddz_dv(a,n,n)=-ddz_dv(a,t,t) + dp_dv(a);
+ ddz_dv(a,n,t)=-eta_jump*ddv_dv(b)*ddel(a,b);
+
+ dC_dv=-dz_dd_sign*2*(eta_jump*v_dv(a)*xyz(a)
+ - dx*dz_dv(a,b)*xyz(a)*dir(b)
+ + dx*dx*ddz_dv(a,b,c)*xyz(a)*dir(b)*dir(c)/2)
+ /(h*h);
+
+ if(x+h/2>middle && x-h/2<middle)
+ {
+ dx=(x>middle) ? ((x-h/2)-middle)
+ : ((x+h/2)-middle);
+
+ C+=(p_jump + dx*dp_jump(a)*dir(a))/h
+ - (dz_jump(a,b)*dir2(a)*dir3(b)
+ + dx*ddz_jump(a,b,c)*dir2(a)*dir3(b)*dir4(c))/h;
+
+ dC_dv+=(p_dv + dx*dp_dv(a)*dir(a))/h
+ - (dz_dv(a,b)*dir2(a)*dir3(b)
+ + dx*ddz_dv(a,b,c)*dir2(a)*dir3(b)*dir4(c))/h;
+ }
+ dC_dz=dC_dv*exp(-log_etax[i][j]);
+ }
+ }
+}
diff -r dec30a218834 -r 5fe6aa6142b8 main.cxx
--- a/main.cxx Sun Mar 11 04:46:02 2012 -0700
+++ b/main.cxx Sun Mar 11 04:57:54 2012 -0700
@@ -10,12 +10,12 @@ extern void initial(const Model &model,
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]);
-extern void compute_Cx(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 int &i, const int &j,
- double &Cx, double &dCx_dzx);
+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 int &i, const int &j,
+ double &C, double &dCx_dz);
extern void compute_Cy(const Model &model, const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -160,8 +160,8 @@ int main()
/* Compute the residual and update zx */
double Cx, dCx_dzx;
- compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,
- dCx_dzx);
+ compute_Cxyz(model,zx,zy,log_etax,log_etay,i,j,Cx,
+ dCx_dzx);
double Rx=2*dzx_xx + dzx_yy + dzy_xy
+ 2*(dlog_etaxx*zx[i][j] + dlog_etax*dzx_x)
@@ -391,8 +391,8 @@ int main()
}
double Cx, dCx_dzx;
- compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,
- dCx_dzx);
+ compute_Cxyz(model,zx,zy,log_etax,log_etay,i,j,Cx,
+ dCx_dzx);
double dRx_dzx=-6/(h*h);
if(j==0 || j==Ny-1)
diff -r dec30a218834 -r 5fe6aa6142b8 wscript
--- a/wscript Sun Mar 11 04:46:02 2012 -0700
+++ b/wscript Sun Mar 11 04:57:54 2012 -0700
@@ -9,7 +9,7 @@ def build(bld):
source = ['main.cxx',
'initial.cxx',
'compute_v_on_interface.cxx',
- 'compute_Cx.cxx',
+ 'compute_Cxyz.cxx',
'compute_Cy.cxx',
'compute_Cp.cxx'],
More information about the CIG-COMMITS
mailing list