[cig-commits] commit: Compute d/dtt using a general method
Mercurial
hg at geodynamics.org
Thu Mar 15 06:53:08 PDT 2012
changeset: 94:d97a502abc10
user: Walter Landry <wlandry at caltech.edu>
date: Wed Mar 14 20:00:23 2012 -0700
files: compute_Cp.cxx compute_Cxyz.cxx compute_v_on_interface.cxx compute_v_on_interface.hxx compute_v_on_interface/compute_dv_dtt.cxx compute_v_on_interface/compute_v_on_interface.cxx wscript
description:
Compute d/dtt using a general method
diff -r 94b794cb763c -r d97a502abc10 compute_Cp.cxx
--- a/compute_Cp.cxx Tue Mar 13 07:48:23 2012 -0700
+++ b/compute_Cp.cxx Wed Mar 14 20:00:23 2012 -0700
@@ -38,9 +38,8 @@ void compute_Cp(const Model &model, cons
FTensor::Tensor1<double,2> v, dv;
FTensor::Tensor1<int,2> dir(0,0);
dir(dd)=1;
- compute_v_on_interface(zx,zy,log_etax,log_etay,
- interface_pos[dd],
- i,j,dd,v,dv);
+ compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
+ interface_pos[dd],i,j,dd,v,dv);
double dx=(pos(dd)>interface_pos[dd](dd))
? (interface_pos[dd](dd)-(pos(dd)-h/2))
diff -r 94b794cb763c -r d97a502abc10 compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Tue Mar 13 07:48:23 2012 -0700
+++ b/compute_Cxyz.cxx Wed Mar 14 20:00:23 2012 -0700
@@ -39,7 +39,8 @@ void compute_Cxyz(const Model &model, co
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[dd],i,j,d,
+ compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
+ interface.pos[dd],i,j,d,
v,dv,ddv,v_dv,dv_dv,ddv_dv);
FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
@@ -104,6 +105,7 @@ void compute_Cxyz(const Model &model, co
dC_dv+=(p_dv + dx*dp_dv(a)*dir(a))/h;
}
+ /* TODO handle corner cutting */
if(interface.intersect_xy[dd])
{
FTensor::Tensor1<int,2> dir2(0,0), dir3(0,0);
diff -r 94b794cb763c -r d97a502abc10 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Tue Mar 13 07:48:23 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,292 +0,0 @@
-#include <algorithm>
-#include <cmath>
-#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)
-{
- double ux(1.1), ux_xp(1.2), ux_xxp(1.3), ux_xyp(1.4),
- ux_xm(1.5), ux_xxm(1.6), ux_xym(1.7),
- uy(2.1), uy_xp(2.2), uy_xxp(2.3), uy_xyp(2.4),
- uy_xm(2.5), uy_xxm(2.6), uy_xym(2.7);
- double ux_y(-(max_eta*uy_xp - min_eta*uy_xm)/eta_jump),
- ux_yy(-(max_eta*uy_xyp - min_eta*uy_xym)/eta_jump),
- uy_y(-(max_eta*ux_xp - min_eta*ux_xm)/eta_jump),
- uy_yy(-(max_eta*ux_xyp - min_eta*ux_xym)/eta_jump);
-
- // std::cout << "u y "
- // << ux << " "
- // << ux_y << " "
- // << ux_yy << " "
- // << uy << " "
- // << uy_y << " "
- // << uy_yy << " "
- // << "\n";
- // exit(0);
-
- if(return_ux)
- {
- if(x>0)
- return ux + x*ux_xp + x*x*ux_xxp/2 + x*y*ux_xyp
- + y*ux_y + y*y*ux_yy/2;
- else
- return ux + x*ux_xm + x*x*ux_xxm/2 + x*y*ux_xym
- + y*ux_y + y*y*ux_yy/2;
- }
- else
- {
- if(x>0)
- return uy + x*uy_xp + x*x*uy_xxp/2 + x*y*uy_xyp
- + y*uy_y + y*y*uy_yy/2;
- else
- return uy + x*uy_xm + x*x*uy_xxm/2 + x*y*uy_xym
- + y*uy_y + y*y*uy_yy/2;
- }
-
-}
-
-
-template<int ny>
-double vel(const double z[][ny], const double log_eta[][ny],
- const int i, const int j)
-{
- // if(j>ny/8-4 && j<ny/8+4)
- // {
- // double x,y;
- // if(ny%2==0)
- // {
- // x=i*h-middle;
- // y=(j+0.5)*h-1.0/8;
- // // y=j*h-1.0/8;
- // return analytic(x,y,true);
- // }
- // else
- // {
- // x=(i+0.5)*h-middle;
- // y=j*h-1.0/8;
- // // y=(j-0.5)*h-1.0/8;
- // return analytic(x,y,false);
- // }
- // }
-
- return z[i][j]*std::exp(-log_eta[i][j]);
-}
-
-template<int ny>
-double compute_dv_yy(const double z[][ny],
- const double log_eta[][ny], const double &dx,
- const int i, const int j)
-{
- double dv_yy_p;
- if(j==0)
- {
- dv_yy_p=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/(h*h);
- }
- else if(j==ny-1)
- {
- dv_yy_p=(-vel(z,log_eta,i+1,j) + vel(z,log_eta,i+1,j-1))/(h*h);
- }
- else
- {
- dv_yy_p=(vel(z,log_eta,i+1,j+1) - 2*vel(z,log_eta,i+1,j)
- + vel(z,log_eta,i+1,j-1))/(h*h);
- }
-
- double dv_yy_m;
- if(j==0)
- {
- dv_yy_m=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/(h*h);
- }
- else if(j==ny-1)
- {
- dv_yy_m=(-vel(z,log_eta,i,j) + vel(z,log_eta,i,j-1))/(h*h);
- }
- else
- {
- dv_yy_m=(vel(z,log_eta,i,j+1) - 2*vel(z,log_eta,i,j)
- + vel(z,log_eta,i,j-1))/(h*h);
- }
-
- double eta_p(exp(log_eta[i+1][j])), eta_m(exp(log_eta[i][j]));
-
- return (dv_yy_p*eta_p*eta_p + dv_yy_m*eta_m*eta_m)/(eta_m*eta_m + eta_p*eta_p);
-}
-
-template<int ny>
-double compute_dv_y(const double z[][ny],
- const double log_eta[][ny],
- const double &jump,
- const double &dx,
- const int i, const int j)
-{
- double dv_y_p1=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/h;
- double dv_y_p2=(vel(z,log_eta,i+2,j+1) - vel(z,log_eta,i+2,j))/h;
- double dv_y_p=(1-dx)*dv_y_p1 + dx*dv_y_p2;
-
- double dv_y_m0=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/h;
- double dv_y_m1=(vel(z,log_eta,i-1,j+1) - vel(z,log_eta,i-1,j))/h;
- double dv_y_m=(1-dx)*dv_y_m1 + dx*dv_y_m0;
-
- return (max_eta*dv_y_p + min_eta*dv_y_m - h*jump)
- /(max_eta + min_eta);
-}
-
-template<int ny>
-double compute_v(const double z[][ny],
- const double log_eta[][ny],
- const double &jump,
- const double &dx,
- const int i, const int j)
-{
- /* v=v0 + dv*x + ddv*h*h/2, where x is centered on v_p2. The
- distance between points is 1 here, which simplifies the
- expressions. */
-
- double v_p1(vel(z,log_eta,i+1,j)), v_p2(vel(z,log_eta,i+2,j)),
- v_p3(vel(z,log_eta,i+3,j));
- double vp0(v_p2), dvp((v_p3-v_p1)/2), ddvp(v_p3 - 2*v_p2 + v_p1);
- double v_p(vp0 - (1-dx)*dvp + (1-dx)*(1-dx)*ddvp/2),
- v_pp(vp0 + dx*dvp + dx*dx*ddvp/2);
-
- double v_m0(vel(z,log_eta,i,j)), v_m1(vel(z,log_eta,i-1,j)),
- v_m2(vel(z,log_eta,i-2,j));
- double vm0(v_m1), dvm((v_m0-v_m2)/2), ddvm(v_m2 - 2*v_m1 + v_m0);
- double v_m(vm0 + dx*dvm + dx*dx*ddvm/2),
- v_mm(vm0 - (1-dx)*dvm + (1-dx)*(1-dx)*ddvm/2);
-
- return (4*(max_eta*v_p + min_eta*v_m) - (max_eta*v_pp + min_eta*v_mm)
- - 2*h*jump)
- /(3*(max_eta+min_eta));
-}
-
-/* For now, hard code the solCx answer */
-double dzy_yx_jump(const double &dvx_yy)
-{
- return -eta_jump*dvx_yy;
-}
-
-double dzx_x_jump(const double &dvy_y)
-{
- return -eta_jump*dvy_y;
-}
-
-double dzx_yx_jump(const double &dvy_yy)
-{
- return -eta_jump*dvy_yy;
-}
-
-double dzy_x_jump(const double &dvx_y)
-{
- return -eta_jump*dvx_y;
-}
-
-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 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,pos,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 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,
- FTensor::Tensor1<double,2> &v_dv,
- FTensor::Tensor1<double,2> &dv_dv,
- FTensor::Tensor1<double,2> &ddv_dv)
-{
- 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);
- 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);
-
- 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
- {
- 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);
-
- v(0)=0;
- dv(1)=0;
- ddv(0)=0;
-
- double eta=exp(log_etay[i][j]);
- double x_i((i+0.5)*h);
- double delta_x=(x_i>middle) ? (middle-(x_i-h))
- : (middle-(x_i+h));
- double delta=(h-std::fabs(delta_x))/h;
- double dvy_yy_dv=-2*eta*eta/((h*h)*(min_eta*min_eta + max_eta*max_eta));
- double dzx_yx_dv=-eta_jump*dvy_yy_dv;
- double dvx_y_dv=(-h/(min_eta+max_eta))*dzx_yx_dv;
- double dzy_x_dv=-eta_jump*dvx_y_dv;
- double dvp_dv=delta/2 + delta*delta/2;
- double dvpp_dv=-(1-delta)/2 + (1-delta)*(1-delta)/2;
- double dvy_dv=
- (eta*(4*dvp_dv - dvpp_dv) - 2*h*dzy_x_dv)/(3*(min_eta+max_eta));
-
- ddv_dv(0)=0;
- ddv_dv(1)=dvy_yy_dv;
- dv_dv(0)=dvx_y_dv;
- dv_dv(1)=0;
- v_dv(0)=0;
- v_dv(1)=dvy_dv;
-
-
- // if(j==1)
- // std::cout << "v "
- // << ddv_dv(1) << " "
- // << dv_dv(0) << " "
- // << v_dv(1) << " "
- // << delta << " "
- // << dy << " "
- // << "\n";
- }
-}
diff -r 94b794cb763c -r d97a502abc10 compute_v_on_interface.hxx
--- a/compute_v_on_interface.hxx Tue Mar 13 07:48:23 2012 -0700
+++ b/compute_v_on_interface.hxx Wed Mar 14 20:00:23 2012 -0700
@@ -3,28 +3,32 @@
#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 FTensor::Tensor1<double,2> &pos,
- const int &i, const int &j,
- 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 &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 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,
- FTensor::Tensor1<double,2> &v_dv,
- FTensor::Tensor1<double,2> &dv_dv,
- FTensor::Tensor1<double,2> &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 distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ 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,
+ FTensor::Tensor1<double,2> &v_dv,
+ FTensor::Tensor1<double,2> &dv_dv,
+ FTensor::Tensor1<double,2> &ddv_dv);
#endif
diff -r 94b794cb763c -r d97a502abc10 compute_v_on_interface/compute_dv_dtt.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Wed Mar 14 20:00:23 2012 -0700
@@ -0,0 +1,272 @@
+#include "../constants.hxx"
+#include <cmath>
+#include "FTensor.hpp"
+#include <vector>
+#include "../sign.hxx"
+#include <tuple>
+#include <algorithm>
+#include <iostream>
+#include <limits>
+
+class Valid
+{
+public:
+ int sign,i,j;
+ bool valid;
+ double r;
+
+ Valid(): sign(0), valid(false), r(std::numeric_limits<double>::infinity()) {}
+ Valid(const int &Sign, const bool &val, const int &I, const int &J, const double &R):
+ sign(Sign), i(I), j(J), valid(val), r(R) {}
+ bool operator<(const Valid &v) const
+ {
+ return r<v.r;
+ }
+};
+
+template<int ny>
+double vel(const double z[][ny], const double log_eta[][ny],
+ const int i, const int j)
+{
+ return z[i][j]*std::exp(-log_eta[i][j]);
+}
+
+template<int ny>
+void compute_derivatives(const int &pm, const int &d, const Valid &v,
+ const int &d0, const int &d1,
+ const int &nx,
+ const double z[][ny],
+ const double log_eta[][ny],
+ FTensor::Tensor3_christof<double,2,2> ddv_pm[2],
+ double eta_pm[2])
+{
+ if(d0!=d1)
+ {
+ if(v.i==nx-1 || v.j==ny-1)
+ {
+ ddv_pm[pm](d,d0,d1)=0;
+ }
+ else
+ {
+ ddv_pm[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);
+ }
+ }
+ else
+ {
+ switch(d0)
+ {
+ case 0:
+ if(v.i==0)
+ {
+ ddv_pm[pm](d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
+ - vel(z,log_eta,v.i,v.j);
+ }
+ else if(v.i==nx-1)
+ {
+ ddv_pm[pm](d,d0,d1)=
+ - vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i-1,v.j);
+ }
+ else
+ {
+ ddv_pm[pm](d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
+ - 2*vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i-1,v.j);
+ }
+ break;
+ case 1:
+ if(v.j==0)
+ {
+ ddv_pm[pm](d,d0,d1)=vel(z,log_eta,v.i,v.j+1)
+ - vel(z,log_eta,v.i,v.j);
+ }
+ else if(v.j==ny-1)
+ {
+ ddv_pm[pm](d,d0,d1)=
+ - vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i,v.j-1);
+ }
+ else
+ {
+ ddv_pm[pm](d,d0,d1)=vel(z,log_eta,v.i,v.j+1)
+ - 2*vel(z,log_eta,v.i,v.j) + vel(z,log_eta,v.i,v.j-1);
+ }
+ break;
+ default:
+ abort();
+ }
+ }
+ eta_pm[pm]=exp(log_eta[v.i][v.j]);
+}
+
+
+double compute_dv_dtt(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,
+ FTensor::Tensor1<double,2> &tangent)
+{
+ /* vx */
+
+ const int max_r=4;
+ FTensor::Tensor2_symmetric<std::vector<Valid>,2> valid;
+ valid(0,0).resize((2*max_r+1)*(2*max_r+1));
+ valid(0,1).resize((2*max_r+1)*(2*max_r+1));
+ valid(1,1).resize((2*max_r+1)*(2*max_r+1));
+
+ const FTensor::Index<'a',2> a;
+ const FTensor::Index<'b',2> b;
+ const FTensor::Index<'c',2> c;
+ FTensor::Tensor3_christof<double,2,2> ddv_pm[2];
+ ddv_pm[0](a,b,c)=0;
+ ddv_pm[1](a,b,c)=0;
+
+ double eta_pm[2];
+
+ for(int d=0;d<2;++d)
+ {
+ int max_x(Nx-d), max_y(Ny+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 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), diff;
+ diff(a)=p(a)-pos(a);
+
+ if(d==0)
+ {
+ if(i>0 && i<max_x)
+ 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])
+ && sign(distx[i][j])==sign(distx[i-1][j]),
+ i,j,diff(a)*diff(a));
+ }
+ else
+ {
+ bool v;
+ if(i==0)
+ v=sign(disty[i][j])==sign(disty[i+1][j]);
+ else if(i==max_x)
+ v=sign(disty[i][j])==sign(disty[i-1][j]);
+ else
+ v=sign(disty[i][j])==sign(disty[i+1][j])
+ && sign(disty[i][j])==sign(disty[i-1][j]);
+
+ valid(0,0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),v,i,j,diff(a)*diff(a));
+ }
+ }
+
+ /* 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)
+ {
+ FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), diff;
+ diff(a)=p(a)-pos(a);
+ if(d==0)
+ {
+ bool v;
+ if(j==0)
+ v=sign(distx[i][j])==sign(distx[i][j+1]);
+ else if(j==max_y)
+ 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][j-1]);
+ valid(1,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),v,i,j,diff(a)*diff(a));
+ }
+ else if(j>0 && j<max_y)
+ 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])
+ && sign(disty[i][j])==sign(disty[i][j-1]),
+ i,j,diff(a)*diff(a));
+ }
+
+ /* 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)
+ {
+ FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), diff;
+ diff(a)=p(a)-pos(a);
+ if(d==0)
+ {
+ if(i>0 && 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]);
+
+ valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),v,i,j,diff(a)*diff(a));
+ }
+ }
+ else
+ {
+ if(j>0 && 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]);
+
+ valid(0,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),v,i,j,diff(a)*diff(a));
+ }
+ }
+ }
+
+ std::sort(valid(0,0).begin(),valid(0,0).end());
+ std::sort(valid(0,1).begin(),valid(0,1).end());
+ std::sort(valid(1,1).begin(),valid(1,1).end());
+
+ FTensor::Tensor3_christof<bool,2,2> is_set(false,false,false,false,
+ false,false);
+ for(int d0=0;d0<2;++d0)
+ for(int d1=d0;d1<2;++d1)
+ {
+ for(auto &v: valid(d0,d1))
+ {
+ 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)=true;
+ if(d==0)
+ compute_derivatives(pm,0,v,d0,d1,Nx+1,zx,log_etax,
+ ddv_pm,eta_pm);
+ else
+ compute_derivatives(pm,1,v,d0,d1,Nx,zy,log_etay,
+ ddv_pm,eta_pm);
+
+ }
+ }
+ }
+ }
+ }
+ double temp(0);
+ FTensor::Tensor1<double,2> ddv(0,0);
+ for(int d=0;d<2;++d)
+ {
+ temp+=eta_pm[d]*eta_pm[d];
+ ddv(a)+=ddv_pm[d](a,b,c)*tangent(b)*tangent(c)*eta_pm[d]*eta_pm[d];
+ }
+
+ return ddv(0)/(temp*h*h);
+}
+
diff -r 94b794cb763c -r d97a502abc10 compute_v_on_interface/compute_v_on_interface.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_v_on_interface/compute_v_on_interface.cxx Wed Mar 14 20:00:23 2012 -0700
@@ -0,0 +1,325 @@
+#include <algorithm>
+#include <cmath>
+#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)
+{
+ double ux(1.1), ux_xp(1.2), ux_xxp(1.3), ux_xyp(1.4),
+ ux_xm(1.5), ux_xxm(1.6), ux_xym(1.7),
+ uy(2.1), uy_xp(2.2), uy_xxp(2.3), uy_xyp(2.4),
+ uy_xm(2.5), uy_xxm(2.6), uy_xym(2.7);
+ double ux_y(-(max_eta*uy_xp - min_eta*uy_xm)/eta_jump),
+ ux_yy(-(max_eta*uy_xyp - min_eta*uy_xym)/eta_jump),
+ uy_y(-(max_eta*ux_xp - min_eta*ux_xm)/eta_jump),
+ uy_yy(-(max_eta*ux_xyp - min_eta*ux_xym)/eta_jump);
+
+ // std::cout << "u y "
+ // << ux << " "
+ // << ux_y << " "
+ // << ux_yy << " "
+ // << uy << " "
+ // << uy_y << " "
+ // << uy_yy << " "
+ // << "\n";
+ // exit(0);
+
+ if(return_ux)
+ {
+ if(x>0)
+ return ux + x*ux_xp + x*x*ux_xxp/2 + x*y*ux_xyp
+ + y*ux_y + y*y*ux_yy/2;
+ else
+ return ux + x*ux_xm + x*x*ux_xxm/2 + x*y*ux_xym
+ + y*ux_y + y*y*ux_yy/2;
+ }
+ else
+ {
+ if(x>0)
+ return uy + x*uy_xp + x*x*uy_xxp/2 + x*y*uy_xyp
+ + y*uy_y + y*y*uy_yy/2;
+ else
+ return uy + x*uy_xm + x*x*uy_xxm/2 + x*y*uy_xym
+ + y*uy_y + y*y*uy_yy/2;
+ }
+
+}
+
+
+template<int ny>
+double vel(const double z[][ny], const double log_eta[][ny],
+ const int i, const int j)
+{
+ // if(j>ny/8-4 && j<ny/8+4)
+ // {
+ // double x,y;
+ // if(ny%2==0)
+ // {
+ // x=i*h-middle;
+ // y=(j+0.5)*h-1.0/8;
+ // // y=j*h-1.0/8;
+ // return analytic(x,y,true);
+ // }
+ // else
+ // {
+ // x=(i+0.5)*h-middle;
+ // y=j*h-1.0/8;
+ // // y=(j-0.5)*h-1.0/8;
+ // return analytic(x,y,false);
+ // }
+ // }
+
+ return z[i][j]*std::exp(-log_eta[i][j]);
+}
+
+template<int ny>
+double compute_dv_yy(const double z[][ny],
+ const double log_eta[][ny], const double &dx,
+ const int i, const int j)
+{
+ double dv_yy_p;
+ if(j==0)
+ {
+ dv_yy_p=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/(h*h);
+ }
+ else if(j==ny-1)
+ {
+ dv_yy_p=(-vel(z,log_eta,i+1,j) + vel(z,log_eta,i+1,j-1))/(h*h);
+ }
+ else
+ {
+ dv_yy_p=(vel(z,log_eta,i+1,j+1) - 2*vel(z,log_eta,i+1,j)
+ + vel(z,log_eta,i+1,j-1))/(h*h);
+ }
+
+ double dv_yy_m;
+ if(j==0)
+ {
+ dv_yy_m=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/(h*h);
+ }
+ else if(j==ny-1)
+ {
+ dv_yy_m=(-vel(z,log_eta,i,j) + vel(z,log_eta,i,j-1))/(h*h);
+ }
+ else
+ {
+ dv_yy_m=(vel(z,log_eta,i,j+1) - 2*vel(z,log_eta,i,j)
+ + vel(z,log_eta,i,j-1))/(h*h);
+ }
+
+ double eta_p(exp(log_eta[i+1][j])), eta_m(exp(log_eta[i][j]));
+
+ return (dv_yy_p*eta_p*eta_p + dv_yy_m*eta_m*eta_m)/(eta_m*eta_m + eta_p*eta_p);
+}
+
+template<int ny>
+double compute_dv_y(const double z[][ny],
+ const double log_eta[][ny],
+ const double &jump,
+ const double &dx,
+ const int i, const int j)
+{
+ double dv_y_p1=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/h;
+ double dv_y_p2=(vel(z,log_eta,i+2,j+1) - vel(z,log_eta,i+2,j))/h;
+ double dv_y_p=(1-dx)*dv_y_p1 + dx*dv_y_p2;
+
+ double dv_y_m0=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/h;
+ double dv_y_m1=(vel(z,log_eta,i-1,j+1) - vel(z,log_eta,i-1,j))/h;
+ double dv_y_m=(1-dx)*dv_y_m1 + dx*dv_y_m0;
+
+ return (max_eta*dv_y_p + min_eta*dv_y_m - h*jump)
+ /(max_eta + min_eta);
+}
+
+template<int ny>
+double compute_v(const double z[][ny],
+ const double log_eta[][ny],
+ const double &jump,
+ const double &dx,
+ const int i, const int j)
+{
+ /* v=v0 + dv*x + ddv*h*h/2, where x is centered on v_p2. The
+ distance between points is 1 here, which simplifies the
+ expressions. */
+
+ double v_p1(vel(z,log_eta,i+1,j)), v_p2(vel(z,log_eta,i+2,j)),
+ v_p3(vel(z,log_eta,i+3,j));
+ double vp0(v_p2), dvp((v_p3-v_p1)/2), ddvp(v_p3 - 2*v_p2 + v_p1);
+ double v_p(vp0 - (1-dx)*dvp + (1-dx)*(1-dx)*ddvp/2),
+ v_pp(vp0 + dx*dvp + dx*dx*ddvp/2);
+
+ double v_m0(vel(z,log_eta,i,j)), v_m1(vel(z,log_eta,i-1,j)),
+ v_m2(vel(z,log_eta,i-2,j));
+ double vm0(v_m1), dvm((v_m0-v_m2)/2), ddvm(v_m2 - 2*v_m1 + v_m0);
+ double v_m(vm0 + dx*dvm + dx*dx*ddvm/2),
+ v_mm(vm0 - (1-dx)*dvm + (1-dx)*(1-dx)*ddvm/2);
+
+ return (4*(max_eta*v_p + min_eta*v_m) - (max_eta*v_pp + min_eta*v_mm)
+ - 2*h*jump)
+ /(3*(max_eta+min_eta));
+}
+
+/* For now, hard code the solCx answer */
+double dzy_yx_jump(const double &dvx_yy)
+{
+ return -eta_jump*dvx_yy;
+}
+
+double dzx_x_jump(const double &dvy_y)
+{
+ return -eta_jump*dvy_y;
+}
+
+double dzx_yx_jump(const double &dvy_yy)
+{
+ return -eta_jump*dvy_yy;
+}
+
+double dzy_x_jump(const double &dvx_y)
+{
+ return -eta_jump*dvx_y;
+}
+
+double compute_dv_dtt(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,
+ FTensor::Tensor1<double,2> &tangent);
+
+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 &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,distx,disty,pos,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 distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ 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,
+ FTensor::Tensor1<double,2> &v_dv,
+ FTensor::Tensor1<double,2> &dv_dv,
+ FTensor::Tensor1<double,2> &ddv_dv)
+{
+ 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);
+ double dx_i((pos(0)-ix*h)/h), dx_j((pos(0)-iy*h)/h-0.5),
+ dy_i((pos(1)-jx*h)/h-0.5), dy_j((pos(1)-jy*h)/h);
+
+ if(xyz==0)
+ {
+ FTensor::Tensor1<double,2> norm, tangent;
+ norm(0)=((distx[iy+1][jy] - disty[iy][jy])*(h-dx_j)
+ + (disty[iy+2][jy] - disty[iy+1][jy])*dx_j)/(h*h);
+ norm(1)=((disty[iy][jy+1] - disty[iy][jy])*(h-dx_j)
+ + (disty[iy+1][jy+1] - disty[iy+1][jy])*dx_j)/(h*h);
+
+ norm(a)/=sqrt(norm(b)*norm(b));
+
+ tangent(0)=-norm(1);
+ tangent(1)=norm(0);
+
+ double length=1/std::max(std::abs(norm(0)),std::abs(norm(1)));
+
+ FTensor::Tensor2<double,2,2> ddv_xy;
+ ddv(0)=compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,tangent);
+
+ // FTensor::Tensor1<double,2> pos_p, pos_pp, pos_m, pos_mm;
+ // pos_p(a)=pos(a)+length*norm(a);
+ // pos_pp(a)=pos(a)+2*norm(a)*length;
+ // pos_m(a)=pos(a)-length*norm(a);
+ // pos_mm(a)=pos(a)-2*norm(a)*length;
+
+
+ // ddv(0)=compute_dv_yy(zx,log_etax,dx_i,ix,j);
+
+
+ dv(1)=compute_dv_y(zy,log_etay,dzy_yx_jump(ddv(0)),dx_j,iy,j);
+ v(0)=compute_v(zx,log_etax,dzx_x_jump(dv(1)),dx_i,ix,j);
+
+ 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
+ {
+ ddv(1)=compute_dv_yy(zy,log_etay,dx_j,iy,j);
+ dv(0)=compute_dv_y(zx,log_etax,dzx_yx_jump(ddv(1)),dx_i,ix,j-1);
+ v(1)=compute_v(zy,log_etay,dzy_x_jump(dv(0)),dx_j,iy,j);
+
+ v(0)=0;
+ dv(1)=0;
+ ddv(0)=0;
+
+ double eta=exp(log_etay[i][j]);
+ double x_i((i+0.5)*h);
+ double delta_x=(x_i>middle) ? (middle-(x_i-h))
+ : (middle-(x_i+h));
+ double delta=(h-std::fabs(delta_x))/h;
+ double dvy_yy_dv=-2*eta*eta/((h*h)*(min_eta*min_eta + max_eta*max_eta));
+ double dzx_yx_dv=-eta_jump*dvy_yy_dv;
+ double dvx_y_dv=(-h/(min_eta+max_eta))*dzx_yx_dv;
+ double dzy_x_dv=-eta_jump*dvx_y_dv;
+ double dvp_dv=delta/2 + delta*delta/2;
+ double dvpp_dv=-(1-delta)/2 + (1-delta)*(1-delta)/2;
+ double dvy_dv=
+ (eta*(4*dvp_dv - dvpp_dv) - 2*h*dzy_x_dv)/(3*(min_eta+max_eta));
+
+ ddv_dv(0)=0;
+ ddv_dv(1)=dvy_yy_dv;
+ dv_dv(0)=dvx_y_dv;
+ dv_dv(1)=0;
+ v_dv(0)=0;
+ v_dv(1)=dvy_dv;
+
+ }
+}
diff -r 94b794cb763c -r d97a502abc10 wscript
--- a/wscript Tue Mar 13 07:48:23 2012 -0700
+++ b/wscript Wed Mar 14 20:00:23 2012 -0700
@@ -8,12 +8,13 @@ def build(bld):
features = ['cxx','cprogram'],
source = ['main.cxx',
'initial.cxx',
- 'compute_v_on_interface.cxx',
+ 'compute_v_on_interface/compute_v_on_interface.cxx',
+ 'compute_v_on_interface/compute_dv_dtt.cxx',
'compute_Cxyz.cxx',
'compute_Cp.cxx'],
target = 'Gamr_disc',
- # cxxflags = ['-std=c++0x','-g','-D_GLIBCXX_DEBUG','-Wall'],
+ # cxxflags = ['-std=c++0x','-g','-D_GLIBCXX_DEBUG','-Wall','-Drestrict='],
cxxflags = ['-std=c++0x','-O3','-Wall','-Drestrict='],
lib = ['boost_filesystem', 'boost_system', 'gsl', 'gslcblas'],
includes = ['FTensor'],
More information about the CIG-COMMITS
mailing list