[cig-commits] commit: Add an Interface class and make compute_Cxyz more able to compute Cy as well as Cx
Mercurial
hg at geodynamics.org
Mon Mar 12 18:01:09 PDT 2012
changeset: 84:14f98b7e3c8c
user: Walter Landry <wlandry at caltech.edu>
date: Mon Mar 12 14:32:21 2012 -0700
files: Interface.hxx compute_Cp.cxx compute_Cxyz.cxx compute_v_on_interface.cxx compute_v_on_interface.hxx main.cxx
description:
Add an Interface class and make compute_Cxyz more able to compute Cy as well as Cx
diff -r 5fe6aa6142b8 -r 14f98b7e3c8c Interface.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Interface.hxx Mon Mar 12 14:32:21 2012 -0700
@@ -0,0 +1,48 @@
+#ifndef GAMR_INTERFACE_HXX
+#define GAMR_INTERFACE_HXX
+
+class Interface
+{
+public:
+ FTensor::Tensor1<double,2> grid_pos, pos;
+ FTensor::Tensor1<bool,2> intersect_dd, intersect_xy;
+ bool intersect_dp;
+
+ Interface(const int &d, const int &i, const int &j)
+ {
+ switch(d)
+ {
+ case 0:
+ grid_pos(0)=i*h;
+ grid_pos(1)=(j+0.5)*h;
+ break;
+ case 1:
+ grid_pos(0)=(i+0.5)*h;
+ grid_pos(1)=j*h;
+ break;
+ default:
+ abort();
+ }
+
+ pos(0)=middle;
+ pos(1)=grid_pos(1);
+
+ intersect_dd(0)=(grid_pos(0) + h > pos(0) && grid_pos(0)-h < pos(0));
+ intersect_dd(1)=false;
+
+ intersect_xy(0)=(grid_pos(0) + h/2 > pos(0) && grid_pos(0)-h/2 < pos(0));
+ intersect_xy(1)=false;
+
+ intersect_dp=(d==1 ? false : intersect_xy(0));
+ }
+
+ bool intersects() const
+ {
+ return intersect_dd(0) || intersect_dd(1)
+ || intersect_xy(0) || intersect_xy(1) || intersect_dp;
+ }
+
+
+};
+
+#endif
diff -r 5fe6aa6142b8 -r 14f98b7e3c8c compute_Cp.cxx
--- a/compute_Cp.cxx Sun Mar 11 04:57:54 2012 -0700
+++ b/compute_Cp.cxx Mon Mar 12 14:32:21 2012 -0700
@@ -18,7 +18,9 @@ void compute_Cp(const Model &model, cons
if(i<Nx && j<Ny && x+h/2>middle && x-h/2<middle)
{
FTensor::Tensor1<double,2> v, dv, dir(1,0);
- compute_v_on_interface(zx,zy,log_etax,log_etay,middle,i,j,0,v,dv);
+ compute_v_on_interface(zx,zy,log_etax,log_etay,
+ FTensor::Tensor1<double,2>(middle,0),
+ i,j,0,v,dv);
double dx=(x>middle) ? (middle-(x-h/2)) : (middle-(x+h/2));
diff -r 5fe6aa6142b8 -r 14f98b7e3c8c compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Sun Mar 11 04:57:54 2012 -0700
+++ b/compute_Cxyz.cxx Mon Mar 12 14:32:21 2012 -0700
@@ -3,38 +3,48 @@
#include <iostream>
#include "compute_v_on_interface.hxx"
#include "FTensor.hpp"
+#include "Interface.hxx"
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 &d,
const int &i, const int &j,
double &C, double &dC_dz)
{
C=0;
dC_dz=0;
- if(model==Model::solCx)
+
+ const Interface interface(d,i,j);
+ if(!interface.intersects())
+ return;
+ FTensor::Tensor1<double,2> pos;
+ const FTensor::Index<'a',2> a;
+ const FTensor::Index<'b',2> b;
+ const FTensor::Index<'c',2> c;
+ pos(a)=interface.grid_pos(a);
+
+ double dC_dv(0);
+ for(int dd=0;dd<2;++dd)
{
- const double x(i*h);
- if(j<Ny && i!=0 && i!=Nx && x-h<middle && x+h>middle)
+ if(interface.intersect_dd(dd))
{
- double dC_dv(0);
/* xyz= the component we are correcting, Cx, Cy, or Cz.
- dir=direction of derivative
+ dir=direction of 2nd 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,
+ FTensor::Tensor1<double,2> v, dv, ddv, v_dv, dv_dv, ddv_dv;
+ FTensor::Tensor1<int,2> xyz(0,0), dir(0,0);
+ xyz(d)=1;
+ dir(dd)=1;
+ compute_v_on_interface(zx,zy,log_etax,log_etay,interface.pos,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));
+ FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
+
+ double dx=(pos(dd)>interface.pos(dd)) ? (interface.pos(dd)-(pos(dd)-h))
+ : (interface.pos(dd)-(pos(dd)+h));
double z_jump=eta_jump*v(a)*xyz(a);
FTensor::Tensor2<double,2,2> dz_jump;
@@ -54,14 +64,13 @@ void compute_Cxyz(const Model &model, co
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;
+ -(z_jump - dx*dz_d_jump + dx*dx*ddz_dd_jump/2)/(h*h);
+ const int dz_dd_factor((pos(dd)>interface.pos(dd) ? -1 : 1)*(d==dd ? 2 : 1));
+ C=dz_dd_factor*dz_dd_correction;
FTensor::Tensor1<double,2> dp_dv(2*eta_jump*ddv_dv(n),
-2*eta_jump*ddv_dv(t));
@@ -76,25 +85,46 @@ void compute_Cxyz(const Model &model, co
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);
+ dC_dv=-dz_dd_factor
+ *(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)
+ if(d==dd && interface.intersect_dp)
{
- dx=(x>middle) ? ((x-h/2)-middle)
- : ((x+h/2)-middle);
+ dx=(pos(dd)>interface.pos(dd)) ? ((pos(dd)-h/2)-interface.pos(dd))
+ : ((pos(dd)+h/2)-interface.pos(dd));
- 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;
+ C+=(p_jump + dx*dp_jump(a)*dir(a))/h;
+ dC_dv+=(p_dv + dx*dp_dv(a)*dir(a))/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;
+ if(interface.intersect_xy(dd))
+ {
+ FTensor::Tensor1<int,2> dir2(0,0), dir3(0,0);
+ dir2((d+1)%2)=1;
+ dir3((dd+1)%2)=1;
+ dx=(pos(dd)>interface.pos(dd)) ? ((pos(dd)-h/2)-interface.pos(dd))
+ : ((pos(dd)+h/2)-interface.pos(dd));
+
+ C+=-(dz_jump(a,b)*dir2(a)*dir3(b)
+ + dx*ddz_jump(a,b,c)*dir2(a)*dir3(b)*dir(c))/h;
+
+ dC_dv+=-(dz_dv(a,b)*dir2(a)*dir3(b)
+ + dx*ddz_dv(a,b,c)*dir2(a)*dir3(b)*dir(c))/h;
}
- dC_dz=dC_dv*exp(-log_etax[i][j]);
}
}
+
+ switch(d)
+ {
+ case 0:
+ dC_dz=dC_dv*exp(-log_etax[i][j]);
+ break;
+ case 1:
+ dC_dz=dC_dv*exp(-log_etay[i][j]);
+ break;
+ default:
+ abort();
+ }
}
diff -r 5fe6aa6142b8 -r 14f98b7e3c8c compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Sun Mar 11 04:57:54 2012 -0700
+++ b/compute_v_on_interface.cxx Mon Mar 12 14:32:21 2012 -0700
@@ -202,7 +202,8 @@ void compute_v_on_interface(const double
double &dvx_yy, double &dvy_yy)
{
FTensor::Tensor1<double,2> v, dv, ddv, v_dv, dv_dv, ddv_dv;
- compute_v_on_interface(zx,zy,log_etax,log_etay,x,i,j,xyz,v,dv,ddv,
+ compute_v_on_interface(zx,zy,log_etax,log_etay,
+ FTensor::Tensor1<double,2>(x,0),i,j,xyz,v,dv,ddv,
v_dv,dv_dv,ddv_dv);
vx=v(0);
vy=v(1);
@@ -216,13 +217,14 @@ void compute_v_on_interface(const double
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
- const double &x, const int &i, const int &j,
+ const FTensor::Tensor1<double,2> &pos,
+ const int &i, const int &j,
const int &xyz,
FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &dv)
{
FTensor::Tensor1<double,2> ddv, v_dv, dv_dv, ddv_dv;
- compute_v_on_interface(zx,zy,log_etax,log_etay,x,i,j,xyz,v,dv,ddv,
+ compute_v_on_interface(zx,zy,log_etax,log_etay,pos,i,j,xyz,v,dv,ddv,
v_dv,dv_dv,ddv_dv);
}
@@ -230,7 +232,8 @@ void compute_v_on_interface(const double
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
- const double &x, const int &i, const int &j,
+ const FTensor::Tensor1<double,2> &pos,
+ const int &i, const int &j,
const int &xyz,
FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &dv,
@@ -239,8 +242,8 @@ void compute_v_on_interface(const double
FTensor::Tensor1<double,2> &dv_dv,
FTensor::Tensor1<double,2> &ddv_dv)
{
- int ix(x/h), iy(x/h - 0.5);
- double dx((x-ix*h)/h), dy((x-iy*h)/h-0.5);
+ int ix(pos(0)/h), iy(pos(0)/h - 0.5);
+ double dx((pos(0)-ix*h)/h), dy((pos(0)-iy*h)/h-0.5);
if(xyz==0)
{
ddv(0)=compute_dv_yy(zx,log_etax,dx,ix,j);
diff -r 5fe6aa6142b8 -r 14f98b7e3c8c compute_v_on_interface.hxx
--- a/compute_v_on_interface.hxx Sun Mar 11 04:57:54 2012 -0700
+++ b/compute_v_on_interface.hxx Mon Mar 12 14:32:21 2012 -0700
@@ -17,7 +17,8 @@ extern void compute_v_on_interface(const
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
- const double &x, const int &i, const int &j,
+ const FTensor::Tensor1<double,2> &pos,
+ const int &i, const int &j,
const int &xyz,
FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &dv);
@@ -26,7 +27,8 @@ extern void compute_v_on_interface(const
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
- const double &x, const int &i, const int &j,
+ const FTensor::Tensor1<double,2> &pos,
+ const int &i, const int &j,
const int &xyz,
FTensor::Tensor1<double,2> &v,
FTensor::Tensor1<double,2> &dv,
diff -r 5fe6aa6142b8 -r 14f98b7e3c8c main.cxx
--- a/main.cxx Sun Mar 11 04:57:54 2012 -0700
+++ b/main.cxx Mon Mar 12 14:32:21 2012 -0700
@@ -14,6 +14,7 @@ extern void compute_Cxyz(const Model &mo
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const int &d,
const int &i, const int &j,
double &C, double &dCx_dz);
@@ -160,7 +161,7 @@ int main()
/* Compute the residual and update zx */
double Cx, dCx_dzx;
- compute_Cxyz(model,zx,zy,log_etax,log_etay,i,j,Cx,
+ compute_Cxyz(model,zx,zy,log_etax,log_etay,0,i,j,Cx,
dCx_dzx);
double Rx=2*dzx_xx + dzx_yy + dzy_xy
@@ -391,8 +392,7 @@ int main()
}
double Cx, dCx_dzx;
- compute_Cxyz(model,zx,zy,log_etax,log_etay,i,j,Cx,
- dCx_dzx);
+ compute_Cxyz(model,zx,zy,log_etax,log_etay,0,i,j,Cx,dCx_dzx);
double dRx_dzx=-6/(h*h);
if(j==0 || j==Ny-1)
More information about the CIG-COMMITS
mailing list