[cig-commits] commit: Move computation of v_dv, dv_dv, and ddv_dv out of compute_Cx into compute_v_on_interface
Mercurial
hg at geodynamics.org
Sun Mar 11 04:07:10 PDT 2012
changeset: 81:15eb88729afc
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Sun Mar 11 04:06:50 2012 -0700
files: compute_Cp.cxx compute_Cx.cxx compute_Cy.cxx compute_v_on_interface.cxx compute_v_on_interface.hxx
description:
Move computation of v_dv, dv_dv, and ddv_dv out of compute_Cx into compute_v_on_interface
diff -r 85f5a13867e2 -r 15eb88729afc compute_Cp.cxx
--- a/compute_Cp.cxx Sun Mar 11 03:16:41 2012 -0700
+++ b/compute_Cp.cxx Sun Mar 11 04:06:50 2012 -0700
@@ -17,14 +17,8 @@ void compute_Cp(const Model &model, cons
const double x=(i+0.5)*h;
if(i<Nx && j<Ny && x+h/2>middle && x-h/2<middle)
{
- double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
-
- compute_v_on_interface(zx,zy,log_etax,log_etay,
- middle,j,0,
- vx,vy,dvx_y,dvy_y,
- dvx_yy,dvy_yy);
-
- FTensor::Tensor1<double,2> v(vx,0), dv(0,dvy_y), dir(1,0);
+ 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);
double dx=(x>middle) ? (middle-(x-h/2)) : (middle-(x+h/2));
diff -r 85f5a13867e2 -r 15eb88729afc compute_Cx.cxx
--- a/compute_Cx.cxx Sun Mar 11 03:16:41 2012 -0700
+++ b/compute_Cx.cxx Sun Mar 11 04:06:50 2012 -0700
@@ -19,20 +19,15 @@ void compute_Cx(const Model &model, cons
if(j<Ny && i!=0 && i!=Nx && x-h<middle && x+h>middle)
{
double dCx_dvx(0);
- double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
-
- compute_v_on_interface(zx,zy,log_etax,
- log_etay,
- middle,j,0,
- vx,vy,dvx_y,dvy_y,
- dvx_yy,dvy_yy);
-
/* 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(vx,0), dv(0,dvy_y), ddv(dvx_yy,0),
+ 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;
@@ -68,26 +63,6 @@ void compute_Cx(const Model &model, cons
const int dz_dd_sign(x>middle ? -1 : 1);
Cx=dz_dd_sign*2*dz_dd_correction;
-
- double eta=exp(log_etax[i][j]);
- double delta=(h-std::fabs(dx))/h;
- double dvx_yy_dv=-2*eta*eta/(h*h*(min_eta*min_eta + max_eta*max_eta));
- if(j==0 || j==Ny-1)
- dvx_yy_dv/=2;
-
- double dzy_yx_dv=-eta_jump*dvx_yy_dv;
- double dvy_y_dv=-dzy_yx_dv*h/(min_eta+max_eta);
- double dzx_x_dv=eta_jump*dvy_y_dv;
- double dvp_dv=delta/2 + delta*delta/2;
- double dvpp_dv=-(1-delta)/2 + (1-delta)*(1-delta)/2;
- double dvx_dv=(eta*(4*dvp_dv - dvpp_dv) + 2*h*dzx_x_dv)
- /(3*(min_eta+max_eta));
-
-
- FTensor::Tensor1<double,2> ddv_dv(dvx_yy_dv,0);
- FTensor::Tensor1<double,2> dv_dv(0,dvy_y_dv);
- FTensor::Tensor1<double,2> v_dv(dvx_dv,0);
-
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);
@@ -119,7 +94,7 @@ void compute_Cx(const Model &model, cons
- (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/eta;
+ dCx_dzx=dCx_dvx*exp(-log_etax[i][j]);
}
}
}
diff -r 85f5a13867e2 -r 15eb88729afc compute_Cy.cxx
--- a/compute_Cy.cxx Sun Mar 11 03:16:41 2012 -0700
+++ b/compute_Cy.cxx Sun Mar 11 04:06:50 2012 -0700
@@ -22,7 +22,7 @@ void compute_Cy(const Model &model, cons
compute_v_on_interface(zx,zy,log_etax,
log_etay,
- middle,j,1,
+ middle,i,j,1,
vx,vy,dvx_y,dvy_y,
dvx_yy,dvy_yy);
diff -r 85f5a13867e2 -r 15eb88729afc compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Sun Mar 11 03:16:41 2012 -0700
+++ b/compute_v_on_interface.cxx Sun Mar 11 04:06:50 2012 -0700
@@ -3,6 +3,8 @@
#include "constants.hxx"
#include <limits>
#include <iostream>
+#include "FTensor.hpp"
+#include "compute_v_on_interface.hxx"
double analytic(const double x, const double y, const bool return_ux)
{
@@ -193,32 +195,94 @@ 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 &j,
+ const double &x, const int &i, const int &j,
const int &xyz,
double &vx, double &vy,
double &dvx_y, double &dvy_y,
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,
+ v_dv,dv_dv,ddv_dv);
+ vx=v(0);
+ vy=v(1);
+ dvx_y=dv(0);
+ dvy_y=dv(1);
+ dvx_yy=ddv(0);
+ dvy_yy=ddv(1);
+}
+
+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 &x, 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,
+ v_dv,dv_dv,ddv_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 &x, const int &i, const int &j,
+ const int &xyz,
+ FTensor::Tensor1<double,2> &v,
+ FTensor::Tensor1<double,2> &dv,
+ FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor1<double,2> &v_dv,
+ 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);
if(xyz==0)
{
- dvx_yy=compute_dv_yy(zx,log_etax,dx,ix,j);
- dvy_y=compute_dv_y(zy,log_etay,dzy_yx_jump(dvx_yy),dy,iy,j);
- vx=compute_v(zx,log_etax,dzx_x_jump(dvy_y),dx,ix,j);
+ ddv(0)=compute_dv_yy(zx,log_etax,dx,ix,j);
+ dv(1)=compute_dv_y(zy,log_etay,dzy_yx_jump(ddv(0)),dy,iy,j);
+ v(0)=compute_v(zx,log_etax,dzx_x_jump(dv(1)),dx,ix,j);
- vy=std::numeric_limits<double>::infinity();
- dvx_y=std::numeric_limits<double>::infinity();
- dvy_yy=std::numeric_limits<double>::infinity();
+ v(1)=0;
+ dv(0)=0;
+ ddv(1)=0;
+
+ double eta=exp(log_etax[i][j]);
+ double x_i(i*h);
+ double delta_x=(x_i>middle) ? (middle-(x_i-h))
+ : (middle-(x_i+h));
+ double delta=(h-std::fabs(delta_x))/h;
+ double dvx_yy_dv=-2*eta*eta/(h*h*(min_eta*min_eta + max_eta*max_eta));
+ if(j==0 || j==Ny-1)
+ dvx_yy_dv/=2;
+
+ double dzy_yx_dv=-eta_jump*dvx_yy_dv;
+ double dvy_y_dv=-dzy_yx_dv*h/(min_eta+max_eta);
+ double dzx_x_dv=eta_jump*dvy_y_dv;
+ double dvp_dv=delta/2 + delta*delta/2;
+ double dvpp_dv=-(1-delta)/2 + (1-delta)*(1-delta)/2;
+ double dvx_dv=(eta*(4*dvp_dv - dvpp_dv) + 2*h*dzx_x_dv)
+ /(3*(min_eta+max_eta));
+
+ ddv_dv(0)=dvx_yy_dv;
+ ddv_dv(1)=0;
+ dv_dv(0)=0;
+ dv_dv(1)=dvy_y_dv;
+ v_dv(0)=dvx_dv;
+ v_dv(1)=0;
}
else
{
- dvy_yy=compute_dv_yy(zy,log_etay,dy,iy,j);
- dvx_y=compute_dv_y(zx,log_etax,dzx_yx_jump(dvy_yy),dx,ix,j-1);
- vy=compute_v(zy,log_etay,dzy_x_jump(dvx_y),dy,iy,j);
+ ddv(1)=compute_dv_yy(zy,log_etay,dy,iy,j);
+ dv(0)=compute_dv_y(zx,log_etax,dzx_yx_jump(ddv(1)),dx,ix,j-1);
+ v(1)=compute_v(zy,log_etay,dzy_x_jump(dv(0)),dy,iy,j);
- vx=std::numeric_limits<double>::infinity();
- dvy_y=std::numeric_limits<double>::infinity();
- dvx_yy=std::numeric_limits<double>::infinity();
+ v(0)=std::numeric_limits<double>::infinity();
+ dv(1)=std::numeric_limits<double>::infinity();
+ ddv(0)=std::numeric_limits<double>::infinity();
}
}
diff -r 85f5a13867e2 -r 15eb88729afc compute_v_on_interface.hxx
--- a/compute_v_on_interface.hxx Sun Mar 11 03:16:41 2012 -0700
+++ b/compute_v_on_interface.hxx Sun Mar 11 04:06:50 2012 -0700
@@ -1,14 +1,38 @@
#ifndef GAMR_COMPUTE_V_ON_INTERFACE_HXX
#define GAMR_COMPUTE_V_ON_INTERFACE_HXX
+
+#include "FTensor.hpp"
extern 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 &x, const int &j,
+ const double &x, const int &i, const int &j,
const int &xyz,
double &vx, double &vy,
double &dvx_y, double &dvy_y,
double &dvx_yy, double &dvy_yy);
+extern 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 &x, const int &i, const int &j,
+ const int &xyz,
+ FTensor::Tensor1<double,2> &v,
+ FTensor::Tensor1<double,2> &dv);
+
+extern 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 &x, const int &i, const int &j,
+ const int &xyz,
+ FTensor::Tensor1<double,2> &v,
+ FTensor::Tensor1<double,2> &dv,
+ FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor1<double,2> &v_dv,
+ FTensor::Tensor1<double,2> &dv_dv,
+ FTensor::Tensor1<double,2> &ddv_dv);
+
#endif
More information about the CIG-COMMITS
mailing list