[cig-commits] commit: Fix Cp. This involved getting rid of the eta_jump global, fixing
Mercurial
hg at geodynamics.org
Wed Mar 28 02:44:27 PDT 2012
changeset: 124:a659d92e837f
user: Walter Landry <wlandry at caltech.edu>
date: Tue Mar 27 16:14:49 2012 -0700
files: compute_Cp.cxx compute_Cxyz.cxx compute_jumps.cxx compute_v_on_interface.hxx compute_v_on_interface/compute_dv_dtt.cxx compute_v_on_interface/compute_v_on_interface.cxx constants.hxx
description:
Fix Cp. This involved getting rid of the eta_jump global, fixing
compute_dv_dvv where I was implicitly assuming norm=(1,0), getting rid
of version of compute_v_on_interface that did not get ddv, and some
bugs in compute_Cp where there should be a += instead of =.
diff -r f376a8388cd2 -r a659d92e837f compute_Cp.cxx
--- a/compute_Cp.cxx Tue Mar 27 16:08:25 2012 -0700
+++ b/compute_Cp.cxx Tue Mar 27 16:14:49 2012 -0700
@@ -3,6 +3,13 @@
#include "compute_v_on_interface.hxx"
#include "FTensor.hpp"
#include "sign.hxx"
+
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
+ const FTensor::Tensor1<double,2> &dv,
+ const FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor1<double,2> &z_jump,
+ FTensor::Tensor2<double,2,2> &dz_jump);
void compute_Cp(const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -26,7 +33,7 @@ void compute_Cp(const double zx[Nx+1][Ny
interface_pos[0](0)=i*h-distx[i][j]*h/(distx[i+1][j]-distx[i][j]);
interface_pos[0](1)=(j+0.5)*h;
}
- else if(intersects[1])
+ if(intersects[1])
{
interface_pos[1](0)=(i+0.5)*h;
interface_pos[1](1)=j*h-disty[i][j]*h/(disty[i][j+1]-disty[i][j]);
@@ -35,26 +42,36 @@ void compute_Cp(const double zx[Nx+1][Ny
for(int dd=0;dd<2;++dd)
if(intersects[dd])
{
- FTensor::Tensor1<double,2> v, dv;
+ FTensor::Tensor1<double,2> v, dv, ddv;
FTensor::Tensor1<int,2> dir(0,0);
dir(dd)=1;
compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
- interface_pos[dd],dd,v,dv);
+ interface_pos[dd],dd,v,dv,ddv);
- double dx=(pos(dd)>interface_pos[dd](dd))
- ? (interface_pos[dd](dd)-(pos(dd)-h/2))
- : (interface_pos[dd](dd)-(pos(dd)+h/2));
+ FTensor::Tensor1<double,2> z_jump;
+ FTensor::Tensor2<double,2,2> dz_jump;
- FTensor::Tensor2_symmetric<double,2> ddel(0,1,0);
+ double eta_jump(sign(pos(dd)-interface_pos[dd](dd)));
+ switch(dd)
+ {
+ case 0:
+ eta_jump*=exp(log_etax[i+1][j])-exp(log_etax[i][j]);
+ break;
+ case 1:
+ eta_jump*=exp(log_etay[i][j+1])-exp(log_etay[i][j]);
+ break;
+ default:
+ abort();
+ break;
+ }
+ compute_jumps(eta_jump,v,dv,ddv,z_jump,dz_jump);
+
+ double dx=std::fabs(pos(dd)-interface_pos[dd](dd)
+ - (h/2)*sign(pos(dd)-interface_pos[dd](dd)));
+
const FTensor::Index<'a',2> a;
const FTensor::Index<'b',2> b;
- FTensor::Number<0> n;
- FTensor::Number<1> t;
- FTensor::Tensor2<double,2,2> dz_jump;
- dz_jump(a,n)=-eta_jump*dv(b)*ddel(a,b);
- dz_jump(a,t)=eta_jump*dv(a);
-
- Cp=(-eta_jump*v(a)*dir(a) + dx*dz_jump(a,b)*dir(a)*dir(b))/h;
+ Cp+=dir(a)*(z_jump(a) + dx*dz_jump(a,b)*dir(b))/h;
}
}
diff -r f376a8388cd2 -r a659d92e837f compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Tue Mar 27 16:08:25 2012 -0700
+++ b/compute_Cxyz.cxx Tue Mar 27 16:14:49 2012 -0700
@@ -4,7 +4,16 @@
#include "FTensor.hpp"
#include "Interface.hxx"
-void compute_jumps(const FTensor::Tensor1<double,2> &v,
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
+ const FTensor::Tensor1<double,2> &dv,
+ const FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor1<double,2> &z_jump,
+ FTensor::Tensor2<double,2,2> &dz_jump,
+ FTensor::Tensor3_christof<double,2,2> &ddz_jump);
+
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
const FTensor::Tensor1<double,2> &dv,
const FTensor::Tensor1<double,2> &ddv,
FTensor::Tensor1<double,2> &z_jump,
@@ -58,14 +67,50 @@ void compute_Cxyz(const double zx[Nx+1][
double p_jump;
FTensor::Tensor1<double,2> dp_jump;
- compute_jumps(v,dv,ddv,z_jump,dz_jump,ddz_jump,p_jump,dp_jump);
+ double eta_jump;
+ switch(d)
+ {
+ case 0:
+ switch(dd)
+ {
+ case 0:
+ eta_jump=exp(log_etax[i+1][j])-exp(log_etax[i-1][j]);
+ break;
+ case 1:
+ eta_jump=exp(log_etax[i][j+1])-exp(log_etax[i][j-1]);
+ break;
+ default:
+ abort();
+ break;
+ }
+ break;
+ case 1:
+ switch(dd)
+ {
+ case 0:
+ eta_jump=exp(log_etay[i+1][j])-exp(log_etay[i-1][j]);
+ break;
+ case 1:
+ eta_jump=exp(log_etay[i][j+1])-exp(log_etay[i][j-1]);
+ break;
+ default:
+ abort();
+ break;
+ }
+ break;
+ default:
+ abort();
+ break;
+ }
+
+ compute_jumps(eta_jump,v,dv,ddv,z_jump,dz_jump,ddz_jump,
+ p_jump,dp_jump);
double dz_dd_correction=
xyz(a)*(z_jump(a) + dz_jump(a,b)*dx(b)
+ ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
- const int dz_dd_factor(sign(pos(dd)-interface.pos[dd](dd))
- *(d==dd ? 2 : 1));
+ const int dz_dd_factor(-sign(dx(dd))*(d==dd ? 2 : 1));
C+=dz_dd_factor*dz_dd_correction;
@@ -113,10 +158,22 @@ void compute_Cxyz(const double zx[Nx+1][
FTensor::Tensor1<double,2> z_jump;
FTensor::Tensor2<double,2,2> dz_jump;
FTensor::Tensor3_christof<double,2,2> ddz_jump;
- double p_jump;
- FTensor::Tensor1<double,2> dp_jump;
- compute_jumps(v,dv,ddv,z_jump,dz_jump,ddz_jump,p_jump,dp_jump);
+ double eta_jump;
+ switch(d)
+ {
+ case 0:
+ eta_jump=exp(log_etax[i][j])-exp(log_etay[i-1+dd][j+ee]);
+ break;
+
+ case 1:
+ eta_jump=exp(log_etay[i][j])-exp(log_etax[i+dd][j-1+ee]);
+ break;
+ default:
+ abort();
+ break;
+ }
+ compute_jumps(eta_jump,v,dv,ddv,z_jump,dz_jump,ddz_jump);
FTensor::Tensor1<double,2> dx;
for(int aa=0;aa<2;++aa)
diff -r f376a8388cd2 -r a659d92e837f compute_jumps.cxx
--- a/compute_jumps.cxx Tue Mar 27 16:08:25 2012 -0700
+++ b/compute_jumps.cxx Tue Mar 27 16:14:49 2012 -0700
@@ -1,7 +1,8 @@
#include "FTensor.hpp"
#include "constants.hxx"
-void compute_jumps(const FTensor::Tensor1<double,2> &v,
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
const FTensor::Tensor1<double,2> &dv,
const FTensor::Tensor1<double,2> &ddv,
FTensor::Tensor1<double,2> &z_jump,
@@ -34,3 +35,28 @@ void compute_jumps(const FTensor::Tensor
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);
}
+
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
+ const FTensor::Tensor1<double,2> &dv,
+ const FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor1<double,2> &z_jump,
+ FTensor::Tensor2<double,2,2> &dz_jump,
+ FTensor::Tensor3_christof<double,2,2> &ddz_jump)
+{
+ double p_jump;
+ FTensor::Tensor1<double,2> dp_jump;
+ compute_jumps(eta_jump,v,dv,ddv,z_jump,dz_jump,ddz_jump,p_jump,dp_jump);
+}
+
+void compute_jumps(const double &eta_jump,
+ const FTensor::Tensor1<double,2> &v,
+ const FTensor::Tensor1<double,2> &dv,
+ const FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor1<double,2> &z_jump,
+ FTensor::Tensor2<double,2,2> &dz_jump)
+{
+ FTensor::Tensor3_christof<double,2,2> ddz_jump;
+ compute_jumps(eta_jump,v,dv,ddv,z_jump,dz_jump,ddz_jump);
+}
+
diff -r f376a8388cd2 -r a659d92e837f compute_v_on_interface.hxx
--- a/compute_v_on_interface.hxx Tue Mar 27 16:08:25 2012 -0700
+++ b/compute_v_on_interface.hxx Tue Mar 27 16:14:49 2012 -0700
@@ -12,17 +12,6 @@ void compute_v_on_interface(const double
const FTensor::Tensor1<double,2> &pos,
const int &xyz,
FTensor::Tensor1<double,2> &v,
- FTensor::Tensor1<double,2> &dv);
-
-void compute_v_on_interface(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 FTensor::Tensor1<double,2> &pos,
- const int &xyz,
- FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &dv,
FTensor::Tensor1<double,2> &ddv);
diff -r f376a8388cd2 -r a659d92e837f compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Tue Mar 27 16:08:25 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Tue Mar 27 16:14:49 2012 -0700
@@ -45,6 +45,7 @@ void compute_dv_dtt(const double zx[Nx+1
double eta_pm[2];
+ /* First find which points can give valid derivatives */
for(int d=0;d<2;++d)
{
FTensor::Tensor2_symmetric<std::vector<Valid>,2> dd_valid;
@@ -366,28 +367,31 @@ void compute_dv_dtt(const double zx[Nx+1
for(int d=0;d<2;++d)
dv_xy(a)+=dv_pm[d](a,b)*tangent(b)*eta_pm[d];
+ dv(n)=dv_xy(a)*norm(a);
+ dv(t)=dv_xy(a)*tangent(a);
+
+
FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
FTensor::Tensor1<double,2> ddz_nt_jump;
double eta_jump(eta_pm[1]-eta_pm[0]);
ddz_nt_jump(a)=-eta_jump*ddv(b)*ddel(a,b);
- dv_xy(a)-=h*h*ddz_nt_jump(a);
- dv_xy(a)/=(eta_pm[0]+eta_pm[1])*h;
- dv(n)=dv_xy(a)*norm(a);
- dv(t)=dv_xy(a)*tangent(a);
+ dv(a)-=h*h*ddz_nt_jump(a);
+ dv(a)/=(eta_pm[0]+eta_pm[1])*h;
/* v */
FTensor::Tensor1<double,2> v_xy(0,0);
for(int d=0;d<2;++d)
v_xy(a)+=v_pm[d](a)*eta_pm[d];
+ v(n)=v_xy(a)*norm(a);
+ v(t)=v_xy(a)*tangent(a);
+
+
FTensor::Tensor1<double,2> dz_n_jump;
dz_n_jump(a)=-eta_jump*dv(b)*ddel(a,b);
- v_xy(a)-=2*h*dz_n_jump(a)/3;
- v_xy(a)/=(eta_pm[0] + eta_pm[1]);
-
- v(n)=v_xy(a)*norm(a);
- v(t)=v_xy(a)*tangent(a);
+ v(a)-=2*h*dz_n_jump(a)/3;
+ v(a)/=(eta_pm[0] + eta_pm[1]);
}
diff -r f376a8388cd2 -r a659d92e837f compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_v_on_interface/compute_v_on_interface.cxx Tue Mar 27 16:08:25 2012 -0700
+++ b/compute_v_on_interface/compute_v_on_interface.cxx Tue Mar 27 16:14:49 2012 -0700
@@ -29,30 +29,12 @@ void compute_v_on_interface(const double
const FTensor::Tensor1<double,2> &pos,
const int &xyz,
FTensor::Tensor1<double,2> &v,
- FTensor::Tensor1<double,2> &dv)
-{
- FTensor::Tensor1<double,2> ddv;
- compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,pos,xyz,v,dv,ddv);
-}
-
-void compute_v_on_interface(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 FTensor::Tensor1<double,2> &pos,
- const int &xyz,
- FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &dv,
FTensor::Tensor1<double,2> &ddv)
{
const FTensor::Index<'a',2> a;
const FTensor::Index<'b',2> b;
int ix(pos(0)/h), iy(pos(0)/h - 0.5), jx(pos(1)/h - 0.5), jy(pos(1)/h);
-
- // FTensor::Tensor2<double,2,2> dx((pos(0)-ix*h)/h,(pos(0)-iy*h)/h-0.5,
- // (pos(1)-jx*h)/h-0.5,(pos(1)-jy*h)/h);
/* TODO: fix for curved interfaces, especially when the interface
point is not along a grid coordinate. */
diff -r f376a8388cd2 -r a659d92e837f constants.hxx
--- a/constants.hxx Tue Mar 27 16:08:25 2012 -0700
+++ b/constants.hxx Tue Mar 27 16:14:49 2012 -0700
@@ -6,7 +6,6 @@ const int Ny(N);
const int Ny(N);
const double min_eta=1;
const double max_eta=1e10;
-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);
More information about the CIG-COMMITS
mailing list