[cig-commits] commit: Fix the way mixed derivatives are calculated.
Mercurial
hg at geodynamics.org
Tue Apr 3 11:23:17 PDT 2012
changeset: 145:67eb40d1a628
user: Walter Landry <wlandry at caltech.edu>
date: Mon Apr 02 20:03:01 2012 -0700
files: compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx
description:
Fix the way mixed derivatives are calculated.
diff -r b7879e4e3cd7 -r 67eb40d1a628 compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx
--- a/compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx Fri Mar 30 12:52:10 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx Mon Apr 02 20:03:01 2012 -0700
@@ -12,18 +12,28 @@ void compute_2nd_derivs(const int &d, co
FTensor::Tensor3_christof<double,2,2> &ddv_pm,
double &eta_pm)
{
+ /* Mixed derivative */
if(d0!=d1)
{
- if(v.i==nx-1 || v.j==ny-1)
+ if(d==0)
{
- ddv_pm(d,d0,d1)=0;
+ if(v.i==0 || v.i==nx)
+ ddv_pm(d,d0,d1)=0;
+ else
+ ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
+ - vel(z,log_eta,v.i,v.j)
+ - vel(z,log_eta,v.i+1,v.j-1)
+ + vel(z,log_eta,v.i,v.j-1);
}
else
{
- ddv_pm(d,d0,d1)=vel(z,log_eta,v.i+1,v.j+1)
- - vel(z,log_eta,v.i,v.j+1)
- - vel(z,log_eta,v.i+1,v.j)
- + vel(z,log_eta,v.i,v.j);
+ if(v.j==0 || v.j==ny)
+ ddv_pm(d,d0,d1)=0;
+ else
+ ddv_pm(d,d0,d1)=vel(z,log_eta,v.i,v.j+1)
+ - vel(z,log_eta,v.i,v.j)
+ - vel(z,log_eta,v.i-1,v.j+1)
+ + vel(z,log_eta,v.i-1,v.j);
}
}
else
diff -r b7879e4e3cd7 -r 67eb40d1a628 compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx Fri Mar 30 12:52:10 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx Mon Apr 02 20:03:01 2012 -0700
@@ -9,7 +9,6 @@
#include "compute_values.hxx"
#include <tuple>
#include <algorithm>
-#include <queue>
void compute_dv_dtt(const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -144,16 +143,16 @@ void compute_dv_dtt(const double zx[Nx+1
if(d==0)
{
/* ddv */
- bool v;
+ bool v_dd;
if(j==0)
- v=sign(distx[i][j])==sign(distx[i][j+1]);
+ v_dd=sign(distx[i][j])==sign(distx[i][j+1]);
else if(j==max_y)
- v=sign(distx[i][j])==sign(distx[i][j-1]);
+ v_dd=sign(distx[i][j])==sign(distx[i][j-1]);
else
- v=sign(distx[i][j])==sign(distx[i][j+1])
+ v_dd=sign(distx[i][j])==sign(distx[i][j+1])
&& sign(distx[i][j])==sign(distx[i][j-1]);
dd_valid(1,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(distx[i][j]),v,i,j,dd_diff);
+ Valid(sign(distx[i][j]),v_dd,i,j,dd_diff);
/* dv */
bool v_d;
@@ -192,39 +191,38 @@ void compute_dv_dtt(const double zx[Nx+1
for(int i=std::max(starti-max_r,0);i<=std::min(starti+max_r,max_x);++i)
for(int j=std::max(startj-max_r,0);j<=std::min(startj+max_r,max_y);++j)
{
- FTensor::Tensor1<double,2> p_dd(i*h+dx+h/2,j*h+dy+h/2), dd_diff;
+ /* Switch dx and dy since the derivative of u_x ends up on
+ the u_y grid and vice versa */
+ FTensor::Tensor1<double,2> p_dd(i*h+dy,j*h+dx), dd_diff;
dd_diff(a)=p_dd(a)-pos(a);
+ bool v_dd;
if(d==0)
{
- if(i<max_x)
- {
- bool v;
- if(j==max_y)
- v=sign(distx[i][j])==sign(distx[i+1][j]);
- else
- v=sign(distx[i][j])==sign(distx[i][j+1])
- && sign(distx[i][j])==sign(distx[i+1][j])
- && sign(distx[i][j])==sign(distx[i+1][j+1]);
-
- dd_valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(distx[i][j]),v,i,j,dd_diff);
- }
+ if(i==0)
+ v_dd=(sign(distx[i+1][j])==sign(distx[i+1][j-1]));
+ else if(i==max_x)
+ v_dd=(sign(distx[i][j])==sign(distx[i][j-1]));
+ else
+ v_dd=(sign(distx[i+1][j])==sign(distx[i][j])
+ && sign(distx[i+1][j])==sign(distx[i+1][j-1])
+ && sign(distx[i+1][j])==sign(distx[i][j-1]));
+
+ dd_valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),v_dd,i,j,dd_diff);
}
else
{
- if(j<max_y)
- {
- bool v;
- if(i==max_x)
- v=sign(distx[i][j])==sign(distx[i][j+1]);
- else
- v=sign(distx[i][j])==sign(distx[i][j+1])
- && sign(distx[i][j])==sign(distx[i+1][j])
- && sign(distx[i][j])==sign(distx[i+1][j+1]);
-
- dd_valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(distx[i][j]),v,i,j,dd_diff);
- }
+ if(j==0)
+ v_dd=sign(disty[i][j+1])==sign(disty[i-1][j+1]);
+ else if(j==max_y)
+ v_dd=sign(disty[i][j])==sign(disty[i-1][j]);
+ else
+ v_dd=sign(disty[i][j+1])==sign(disty[i][j])
+ && sign(disty[i][j+1])==sign(disty[i-1][j+1])
+ && sign(disty[i][j+1])==sign(disty[i-1][j]);
+
+ dd_valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),v_dd,i,j,dd_diff);
}
}
@@ -273,9 +271,11 @@ void compute_dv_dtt(const double zx[Nx+1
++y;
if(values[pm].size()==5)
{
- if(find(x_coords.begin(),x_coords.end(),value.i)==x_coords.end())
+ if(find(x_coords.begin(),x_coords.end(),value.i)
+ ==x_coords.end())
x_coords.push_back(value.i);
- if(find(y_coords.begin(),y_coords.end(),value.j)==y_coords.end())
+ if(find(y_coords.begin(),y_coords.end(),value.j)
+ ==y_coords.end())
y_coords.push_back(value.j);
}
}
@@ -345,7 +345,8 @@ void compute_dv_dtt(const double zx[Nx+1
for(int pm=0;pm<2;++pm)
{
- if(V.sign==(pm==0 ? -1 : 1) && V.valid && !is_set(pm,d0,d1))
+ if(V.sign==(pm==0 ? -1 : 1) && V.valid
+ && !is_set(pm,d0,d1))
{
is_set(pm,d0,d1)=true;
if(d==0)
@@ -381,9 +382,9 @@ void compute_dv_dtt(const double zx[Nx+1
dv(a)=nt_to_xy(b,a)*dv_xy(b);
FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
+ double eta_jump(eta_pm[1]-eta_pm[0]);
+ // /* TODO: fix for curved surfaces */
FTensor::Tensor1<double,2> ddz_nt_jump;
- double eta_jump(eta_pm[1]-eta_pm[0]);
- /* TODO: fix for curved surfaces */
ddz_nt_jump(a)=-eta_jump*ddv(b)*ddel(a,b);
dv(a)-=h*length*ddz_nt_jump(a);
More information about the CIG-COMMITS
mailing list