[cig-commits] commit: Some fixes for derivatives near the boundary. Also, if multiple 2nd
Mercurial
hg at geodynamics.org
Fri Apr 6 10:29:03 PDT 2012
changeset: 150:597cdfd1cf51
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Apr 06 10:28:49 2012 -0700
files: compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx
description:
Some fixes for derivatives near the boundary. Also, if multiple 2nd
derivatives are the same distance, it uses an average.
diff -r 32f74f34089c -r 597cdfd1cf51 compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx
--- a/compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx Fri Apr 06 10:20:15 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx Fri Apr 06 10:28:49 2012 -0700
@@ -21,7 +21,7 @@ void compute_1st_derivs(const int &d, co
switch(d0)
{
case 0:
- if(derivs[i].i==nx-1)
+ if(derivs[i].i==nx-1 || derivs[i].i==-1)
{
rhs[i]=0;
}
@@ -32,7 +32,7 @@ void compute_1st_derivs(const int &d, co
}
break;
case 1:
- if(derivs[i].j==ny-1)
+ if(derivs[i].j==ny-1 || derivs[i].j==-1)
{
rhs[i]=0;
}
diff -r 32f74f34089c -r 597cdfd1cf51 compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx
--- a/compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx Fri Apr 06 10:20:15 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_2nd_derivs.hxx Fri Apr 06 10:28:49 2012 -0700
@@ -17,7 +17,7 @@ void compute_2nd_derivs(const int &d, co
{
if(d==0)
{
- if(v.i==0 || v.i==nx)
+ 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+1,v.j)
@@ -27,7 +27,7 @@ void compute_2nd_derivs(const int &d, co
}
else
{
- if(v.j==0 || v.j==ny)
+ 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,v.j+1)
diff -r 32f74f34089c -r 597cdfd1cf51 compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx Fri Apr 06 10:20:15 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx Fri Apr 06 10:28:49 2012 -0700
@@ -9,6 +9,11 @@
#include "compute_values.hxx"
#include <tuple>
#include <algorithm>
+
+inline bool about_equal(const double &r1, const double &r2)
+{
+ return std::fabs(r1-r2)<1e-8*std::fabs(r1+r2);
+}
void compute_dv_dtt(const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -57,13 +62,13 @@ void compute_dv_dtt(const double zx[Nx+1
std::vector<Valid> valid((2*max_r+1)*(2*max_r+1));
- int max_x(Nx-d), max_y(Ny+d-1);
+ int max_x(Nx-d), max_y(Ny+d-1), min_x(-d), min_y(d-1);
double dx(d*h/2), dy((1-d)*h/2);
int starti((pos(0)-dx)/h), startj((pos(1)-dy)/h);
/* d/dx^2 */
- for(int i=std::max(starti-max_r,0);i<=std::min(starti+max_r,max_x);++i)
+ for(int i=std::max(starti-max_r,min_x);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(i*h+dx,j*h+dy), p_d(i*h+dx+h/2,j*h+dy),
@@ -73,8 +78,9 @@ void compute_dv_dtt(const double zx[Nx+1
if(d==0)
{
+ dd_diff(a)-=length*norm(a)*sign(distx[i][j]);
/* ddv */
- if(i>0 && i<max_x)
+ if(i>min_x && i<max_x)
dd_valid(0,0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
Valid(sign(distx[i][j]),
sign(distx[i][j])==sign(distx[i+1][j])
@@ -90,32 +96,35 @@ void compute_dv_dtt(const double zx[Nx+1
i,j,d_diff);
}
/* v */
- diff(a)=dd_diff(a)-1.5*length*norm(a)*sign(distx[i][j]);
+ diff(a)=dd_diff(a)-0.5*length*norm(a)*sign(distx[i][j]);
valid[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
Valid(sign(distx[i][j]),true,i,j,diff);
}
else
{
+ dd_diff(a)-=length*norm(a)*sign(disty[i][j]);
/* ddv */
- bool v_dd;
- if(i==0)
- v_dd=(sign(disty[i][j])==sign(disty[i+1][j]));
- else if(i==max_x)
- v_dd=(sign(disty[i][j])==sign(disty[i-1][j]));
- else
- v_dd=(sign(disty[i][j])==sign(disty[i+1][j])
- && sign(disty[i][j])==sign(disty[i-1][j]));
+ if(i>min_x)
+ {
+ bool v_dd;
+ if(i==0)
+ v_dd=(sign(disty[i][j])==sign(disty[i+1][j]));
+ else if(i==max_x)
+ v_dd=(sign(disty[i][j])==sign(disty[i-1][j]));
+ else
+ v_dd=(sign(disty[i][j])==sign(disty[i+1][j])
+ && sign(disty[i][j])==sign(disty[i-1][j]));
- dd_valid(0,0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(disty[i][j]),v_dd,i,j,dd_diff);
-
+ dd_valid(0,0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),v_dd,i,j,dd_diff);
+ }
/* dv */
bool v_d;
/* TODO: This bothers me a little. What if the
interface lies between the point and the boundary?
Shouldn't that make the derivative calculation
invalid? */
- if(i==max_x)
+ if(i==min_x || i==max_x)
v_d=true;
else
v_d=(sign(disty[i][j])==sign(disty[i+1][j]));
@@ -125,15 +134,17 @@ void compute_dv_dtt(const double zx[Nx+1
Valid(sign(disty[i][j]),v_d,i,j,d_diff);
/* v */
- diff(a)=dd_diff(a)-1.5*length*norm(a)*sign(disty[i][j]);
- valid[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(disty[i][j]),true,i,j,diff);
+ diff(a)=dd_diff(a)-0.5*length*norm(a)*sign(disty[i][j]);
+ if(i>min_x)
+ valid[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),true,i,j,diff);
}
}
/* d/dy^2 */
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)
+ for(int j=std::max(startj-max_r,min_y);j<=std::min(startj+max_r,max_y);
+ ++j)
{
FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), p_d(i*h+dx,j*h+dy+h/2),
dd_diff, d_diff;
@@ -143,23 +154,27 @@ void compute_dv_dtt(const double zx[Nx+1
if(d==0)
{
/* ddv */
- bool v_dd;
- if(j==0)
- v_dd=sign(distx[i][j])==sign(distx[i][j+1]);
- else if(j==max_y)
- v_dd=sign(distx[i][j])==sign(distx[i][j-1]);
- else
- 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_dd,i,j,dd_diff);
+ dd_diff(a)-=length*norm(a)*sign(distx[i][j]);
+ if(i>min_y)
+ {
+ bool v_dd;
+ if(j==0)
+ v_dd=(sign(distx[i][j])==sign(distx[i][j+1]));
+ else if(j==max_y)
+ v_dd=(sign(distx[i][j])==sign(distx[i][j-1]));
+ else
+ 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_dd,i,j,dd_diff);
+ }
/* dv */
bool v_d;
- if(j==max_y)
+ if(j==min_y || j==max_y)
v_d=true;
else
- v_d=sign(distx[i][j])==sign(distx[i][j+1]);
+ v_d=(sign(distx[i][j])==sign(distx[i][j+1]));
d_diff(a)-=length*norm(a)*sign(distx[i][j]);
d_valid(1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
@@ -168,7 +183,8 @@ void compute_dv_dtt(const double zx[Nx+1
else
{
/* ddv */
- if(j>0 && j<max_y)
+ dd_diff(a)-=length*norm(a)*sign(disty[i][j]);
+ if(j>min_y && j<max_y)
dd_valid(1,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
Valid(sign(disty[i][j]),
sign(disty[i][j])==sign(disty[i][j+1])
@@ -188,8 +204,8 @@ void compute_dv_dtt(const double zx[Nx+1
}
/* d/dxdy */
- 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)
+ for(int i=std::max(starti-max_r,0);i<=std::min(starti+max_r,Nx+d-1);++i)
+ for(int j=std::max(startj-max_r,0);j<=std::min(startj+max_r,Ny-d);++j)
{
/* Switch dx and dy since the derivative of u_x ends up on
the u_y grid and vice versa */
@@ -198,29 +214,31 @@ void compute_dv_dtt(const double zx[Nx+1
bool v_dd;
if(d==0)
{
- 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]));
+ if(j==0)
+ v_dd=(sign(distx[i+1][j])==sign(distx[i][j]));
+ else if(j==Ny)
+ v_dd=(sign(distx[i+1][j-1])==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_diff(a)-=length*norm(a)*sign(distx[i][j]);
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==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]);
+ if(i==0)
+ v_dd=(sign(disty[i][j+1])==sign(disty[i][j]));
+ else if(i==Nx)
+ v_dd=(sign(disty[i-1][j+1])==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_diff(a)-=length*norm(a)*sign(disty[i][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);
}
@@ -338,16 +356,19 @@ void compute_dv_dtt(const double zx[Nx+1
/* Second derivative ddv */
for(int d1=d0;d1<2;++d1)
{
+ double distances[]={0,0};
+ int num[]={0,0};
for(auto &V: dd_valid(d0,d1))
{
- if(is_set(0,d0,d1) && is_set(1,d0,d1))
- break;
-
for(int pm=0;pm<2;++pm)
{
if(V.sign==(pm==0 ? -1 : 1) && V.valid
- && !is_set(pm,d0,d1))
+ && (!is_set(pm,d0,d1)
+ || about_equal(V.r,distances[pm])))
{
+ distances[pm]=V.r;
+ ++num[pm];
+
is_set(pm,d0,d1)=true;
if(d==0)
compute_2nd_derivs(0,V,d0,d1,Nx+1,zx,log_etax,
@@ -358,13 +379,15 @@ void compute_dv_dtt(const double zx[Nx+1
}
}
}
+ for(int pm=0;pm<2;++pm)
+ ddv_pm[pm](d,d0,d1)/=num[pm];
}
}
}
/* ddv */
+ FTensor::Tensor1<double,2> ddv_xy(0,0);
double temp(0);
- FTensor::Tensor1<double,2> ddv_xy(0,0);
for(int d=0;d<2;++d)
{
temp+=eta_pm[d]*eta_pm[d];
More information about the CIG-COMMITS
mailing list