[cig-commits] commit: Split Valid, compute_2nd_derivs, and vel into separate files.
Mercurial
hg at geodynamics.org
Tue Mar 20 14:53:50 PDT 2012
changeset: 96:5adaedce245c
user: Walter Landry <wlandry at caltech.edu>
date: Mon Mar 19 14:25:43 2012 -0700
files: compute_v_on_interface/Valid.hxx compute_v_on_interface/compute_2nd_derivs.hxx compute_v_on_interface/compute_dv_dtt.cxx compute_v_on_interface/compute_v_on_interface.cxx compute_v_on_interface/vel.hxx
description:
Split Valid, compute_2nd_derivs, and vel into separate files.
diff -r a451764d762d -r 5adaedce245c compute_v_on_interface/Valid.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_v_on_interface/Valid.hxx Mon Mar 19 14:25:43 2012 -0700
@@ -0,0 +1,23 @@
+#ifndef GAMR_VALID_HXX
+#define GAMR_VALID_HXX
+
+#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;
+ }
+};
+
+
+#endif
diff -r a451764d762d -r 5adaedce245c compute_v_on_interface/compute_2nd_derivs.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_v_on_interface/compute_2nd_derivs.hxx Mon Mar 19 14:25:43 2012 -0700
@@ -0,0 +1,74 @@
+#ifndef GAMR_COMPUTE_2ND_DERIVS_HXX
+#define GAMR_COMPUTE_2ND_DERIVS_HXX
+
+#include "Valid.hxx"
+
+template<int ny>
+void compute_2nd_derivs(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]);
+}
+
+#endif
diff -r a451764d762d -r 5adaedce245c compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Thu Mar 15 06:52:55 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Mon Mar 19 14:25:43 2012 -0700
@@ -3,112 +3,23 @@
#include "FTensor.hpp"
#include <vector>
#include "../sign.hxx"
+#include "vel.hxx"
+#include "compute_2nd_derivs.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]);
-}
-
-
-FTensor::Tensor1<double,2> 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> &norm,
- FTensor::Tensor1<double,2> &tangent)
+void 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,
+ const FTensor::Tensor1<double,2> &norm,
+ const FTensor::Tensor1<double,2> &tangent,
+ FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor1<double,2> &ddv_dv)
{
/* vx */
@@ -249,11 +160,11 @@ FTensor::Tensor1<double,2> compute_dv_dt
{
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);
+ compute_2nd_derivs(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);
+ compute_2nd_derivs(pm,1,v,d0,d1,Nx,zy,log_etay,
+ ddv_pm,eta_pm);
}
}
@@ -261,14 +172,15 @@ FTensor::Tensor1<double,2> compute_dv_dt
}
}
double temp(0);
- FTensor::Tensor1<double,2> ddv(0,0);
+ FTensor::Tensor1<double,2> ddv_xy(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];
+ ddv_xy(a)+=ddv_pm[d](a,b,c)*tangent(b)*tangent(c)*eta_pm[d]*eta_pm[d];
}
- ddv(a)/=temp*h*h;
+ ddv_xy(a)/=temp*h*h;
- return FTensor::Tensor1<double,2>(ddv(a)*norm(a),ddv(a)*tangent(a));
+ ddv(0)=ddv_xy(a)*norm(a);
+ ddv(1)=ddv_xy(a)*tangent(a);
}
diff -r a451764d762d -r 5adaedce245c compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_v_on_interface/compute_v_on_interface.cxx Thu Mar 15 06:52:55 2012 -0700
+++ b/compute_v_on_interface/compute_v_on_interface.cxx Mon Mar 19 14:25:43 2012 -0700
@@ -5,75 +5,7 @@
#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]);
-}
+#include "vel.hxx"
template<int ny>
double compute_dv_yy(const double z[][ny],
@@ -183,15 +115,17 @@ double dzy_x_jump(const double &dvx_y)
return -eta_jump*dvx_y;
}
-FTensor::Tensor1<double,2> 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> &norm,
- FTensor::Tensor1<double,2> &tangent);
+void 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,
+ const FTensor::Tensor1<double,2> &norm,
+ const FTensor::Tensor1<double,2> &tangent,
+ FTensor::Tensor1<double,2> &ddv,
+ FTensor::Tensor1<double,2> &ddv_dv);
void compute_v_on_interface(const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -245,10 +179,11 @@ void compute_v_on_interface(const double
tangent(0)=-norm(1);
tangent(1)=norm(0);
- double length=1/std::max(std::abs(norm(0)),std::abs(norm(1)));
+ // double length=1/std::max(std::abs(norm(0)),std::abs(norm(1)));
- ddv(a)=compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,
- norm,tangent)(a);
+ FTensor::Tensor1<double,2> ddv_dv_temp;
+ compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,
+ norm,tangent,ddv,ddv_dv_temp);
// FTensor::Tensor1<double,2> pos_p, pos_pp, pos_m, pos_mm;
// pos_p(a)=pos(a)+length*norm(a);
diff -r a451764d762d -r 5adaedce245c compute_v_on_interface/vel.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_v_on_interface/vel.hxx Mon Mar 19 14:25:43 2012 -0700
@@ -0,0 +1,73 @@
+#ifndef GAMR_VEL_HXX
+#define GAMR_VEL_HXX
+
+inline 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]);
+}
+
+#endif
More information about the CIG-COMMITS
mailing list