[cig-commits] commit: Added some support for computing 1st derivatives in a general manner.
Mercurial
hg at geodynamics.org
Wed Mar 21 03:47:56 PDT 2012
changeset: 102:018b068cbb9f
user: Walter Landry <wlandry at caltech.edu>
date: Tue Mar 20 16:31:16 2012 -0700
files: compute_v_on_interface/Valid.hxx compute_v_on_interface/compute_1st_derivs.hxx compute_v_on_interface/compute_2nd_derivs.hxx compute_v_on_interface/compute_dv_dtt.cxx
description:
Added some support for computing 1st derivatives in a general manner.
diff -r 65c1c31f67b1 -r 018b068cbb9f compute_v_on_interface/Valid.hxx
--- a/compute_v_on_interface/Valid.hxx Tue Mar 20 14:53:35 2012 -0700
+++ b/compute_v_on_interface/Valid.hxx Tue Mar 20 16:31:16 2012 -0700
@@ -8,14 +8,17 @@ public:
public:
int sign,i,j;
bool valid;
+ FTensor::Tensor1<double,2> diff;
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) {}
+ Valid(const int &Sign, const bool &val, const int &I, const int &J,
+ const FTensor::Tensor1<double,2> &Diff):
+ sign(Sign), i(I), j(J), valid(val), diff(Diff),
+ r(diff(0)*diff(0)+diff(1)*diff(1)) {}
bool operator<(const Valid &v) const
{
- return r<v.r;
+ return r<v.r;
}
};
diff -r 65c1c31f67b1 -r 018b068cbb9f compute_v_on_interface/compute_1st_derivs.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_v_on_interface/compute_1st_derivs.hxx Tue Mar 20 16:31:16 2012 -0700
@@ -0,0 +1,64 @@
+#ifndef GAMR_COMPUTE_1ST_DERIVS_HXX
+#define GAMR_COMPUTE_1ST_DERIVS_HXX
+
+#include <cassert>
+#include "Valid.hxx"
+#include <gsl/gsl_linalg.h>
+
+template<int ny>
+void compute_1st_derivs(const int &d, const std::vector<Valid> &derivs,
+ const int &d0,
+ const int &nx,
+ const double z[][ny],
+ const double log_eta[][ny],
+ FTensor::Tensor2<double,2,2> &dv_pm)
+{
+ /* Create and invert the Vandermonde matrix */
+ assert(derivs.size()==3);
+ double m[3*3], rhs[3];
+ for(int i=0;i<3;++i)
+ {
+ switch(d0)
+ {
+ case 0:
+ if(derivs[i].i==nx-1)
+ {
+ rhs[i]=0;
+ }
+ else
+ {
+ rhs[i]=vel(z,log_eta,derivs[i].i+1,derivs[i].j)
+ - vel(z,log_eta,derivs[i].i,derivs[i].j);
+ }
+ break;
+ case 1:
+ if(derivs[i].j==ny-1)
+ {
+ rhs[i]=0;
+ }
+ else
+ {
+ rhs[i]=vel(z,log_eta,derivs[i].i,derivs[i].j+1)
+ - vel(z,log_eta,derivs[i].i,derivs[i].j);
+ }
+ break;
+ default:
+ abort();
+ }
+ m[0+3*i]=1;
+ m[1+3*i]=derivs[i].diff(0);
+ m[2+3*i]=derivs[i].diff(1);
+ }
+
+ gsl_matrix_view mv=gsl_matrix_view_array(m,3,3);
+ gsl_vector_view rhsv=gsl_vector_view_array(rhs,3);
+ gsl_vector *x=gsl_vector_alloc(3);
+ int s;
+ gsl_permutation *perm=gsl_permutation_alloc(3);
+ gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
+ gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,x);
+
+ dv_pm(d,d0)=gsl_vector_get(x,0);
+}
+
+#endif
diff -r 65c1c31f67b1 -r 018b068cbb9f compute_v_on_interface/compute_2nd_derivs.hxx
--- a/compute_v_on_interface/compute_2nd_derivs.hxx Tue Mar 20 14:53:35 2012 -0700
+++ b/compute_v_on_interface/compute_2nd_derivs.hxx Tue Mar 20 16:31:16 2012 -0700
@@ -4,23 +4,23 @@
#include "Valid.hxx"
template<int ny>
-void compute_2nd_derivs(const int &pm, const int &d, const Valid &v,
+void compute_2nd_derivs(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])
+ FTensor::Tensor3_christof<double,2,2> &ddv_pm,
+ double &eta_pm)
{
if(d0!=d1)
{
if(v.i==nx-1 || v.j==ny-1)
{
- ddv_pm[pm](d,d0,d1)=0;
+ ddv_pm(d,d0,d1)=0;
}
else
{
- ddv_pm[pm](d,d0,d1)=vel(z,log_eta,v.i+1,v.j+1)
+ 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);
@@ -33,34 +33,34 @@ void compute_2nd_derivs(const int &pm, c
case 0:
if(v.i==0)
{
- ddv_pm[pm](d,d0,d1)=vel(z,log_eta,v.i+1,v.j)
+ ddv_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)=
+ ddv_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)
+ ddv_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)
+ ddv_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)=
+ ddv_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)
+ ddv_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;
@@ -68,7 +68,7 @@ void compute_2nd_derivs(const int &pm, c
abort();
}
}
- eta_pm[pm]=exp(log_eta[v.i][v.j]);
+ eta_pm=exp(log_eta[v.i][v.j]);
}
#endif
diff -r 65c1c31f67b1 -r 018b068cbb9f compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Tue Mar 20 14:53:35 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Tue Mar 20 16:31:16 2012 -0700
@@ -5,6 +5,7 @@
#include "../sign.hxx"
#include "vel.hxx"
#include "compute_2nd_derivs.hxx"
+#include "compute_1st_derivs.hxx"
#include <tuple>
#include <algorithm>
#include <iostream>
@@ -30,6 +31,10 @@ void compute_dv_dtt(const double zx[Nx+1
FTensor::Tensor3_christof<double,2,2> ddv_pm[2];
ddv_pm[0](a,b,c)=0;
ddv_pm[1](a,b,c)=0;
+
+ FTensor::Tensor2<double,2,2> dv_pm[2];
+ dv_pm[0](a,b)=0;
+ dv_pm[1](a,b)=0;
double eta_pm[2];
@@ -66,12 +71,12 @@ void compute_dv_dtt(const double zx[Nx+1
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,dd_diff(a)*dd_diff(a));
+ i,j,dd_diff);
if(i<max_x)
d_valid(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]),
- i,j,d_diff(a)*d_diff(a));
+ i,j,d_diff);
}
else
@@ -86,7 +91,7 @@ void compute_dv_dtt(const double zx[Nx+1
&& 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(a)*dd_diff(a));
+ Valid(sign(disty[i][j]),v_dd,i,j,dd_diff);
bool v_d;
/* This bothers me a little. What if the interface
@@ -98,7 +103,7 @@ void compute_dv_dtt(const double zx[Nx+1
v_d=(sign(disty[i][j])==sign(disty[i+1][j]));
d_valid(0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(disty[i][j]),v_d,i,j,d_diff(a)*d_diff(a));
+ Valid(sign(disty[i][j]),v_d,i,j,d_diff);
}
}
@@ -122,7 +127,7 @@ void compute_dv_dtt(const double zx[Nx+1
v=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(a)*dd_diff(a));
+ Valid(sign(distx[i][j]),v,i,j,dd_diff);
bool v_d;
if(j==max_y)
@@ -130,7 +135,7 @@ void compute_dv_dtt(const double zx[Nx+1
else
v_d=sign(distx[i][j])==sign(distx[i][j+1]);
d_valid(1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(distx[i][j]),v_d,i,j,d_diff(a)*d_diff(a));
+ Valid(sign(distx[i][j]),v_d,i,j,d_diff);
}
else
{
@@ -138,11 +143,11 @@ void compute_dv_dtt(const double zx[Nx+1
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])
&& sign(disty[i][j])==sign(disty[i][j-1]),
- i,j,dd_diff(a)*dd_diff(a));
+ i,j,dd_diff);
if(j<max_y)
d_valid(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]),
- i,j,d_diff(a)*d_diff(a));
+ i,j,d_diff);
}
}
@@ -165,7 +170,7 @@ void compute_dv_dtt(const double zx[Nx+1
&& 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(a)*dd_diff(a));
+ Valid(sign(distx[i][j]),v,i,j,dd_diff);
}
}
else
@@ -181,7 +186,7 @@ void compute_dv_dtt(const double zx[Nx+1
&& 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(a)*dd_diff(a));
+ Valid(sign(distx[i][j]),v,i,j,dd_diff);
}
}
}
@@ -194,26 +199,63 @@ void compute_dv_dtt(const double zx[Nx+1
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: dd_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_2nd_derivs(pm,0,v,d0,d1,Nx+1,zx,log_etax,
- ddv_pm,eta_pm);
- else
- compute_2nd_derivs(pm,1,v,d0,d1,Nx,zy,log_etay,
- ddv_pm,eta_pm);
- }
- }
- }
- }
+ {
+ std::vector<Valid> derivs[2];
+ /* First derivative */
+ for(auto &v: d_valid(d0))
+ {
+ if(v.valid)
+ {
+ const int pm=(v.sign==-1 ? 0 : 1);
+ switch(derivs[pm].size())
+ {
+ default:
+ break;
+ case 0:
+ case 1:
+ derivs[pm].push_back(v);
+ break;
+ case 2:
+ Valid &v0(derivs[pm][0]);
+ Valid &v1(derivs[pm][1]);
+ if(!((v0.i==v1.i && v1.i==v.i)
+ || (v0.j==v1.j && v1.j==v.j)))
+ {
+ derivs[pm].push_back(v);
+ if(d==0)
+ compute_1st_derivs(0,derivs[pm],d0,Nx+1,
+ zx,log_etax,dv_pm[pm]);
+ else
+ compute_1st_derivs(1,derivs[pm],d0,Nx,
+ zy,log_etay,dv_pm[pm]);
+ }
+ break;
+ }
+ }
+ }
+ /* Second derivative */
+ for(int d1=d0;d1<2;++d1)
+ {
+ for(auto &v: dd_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_2nd_derivs(0,v,d0,d1,Nx+1,zx,log_etax,
+ ddv_pm[pm],eta_pm[pm]);
+ else
+ compute_2nd_derivs(1,v,d0,d1,Nx,zy,log_etay,
+ ddv_pm[pm],eta_pm[pm]);
+ }
+ }
+ }
+ }
+ }
}
double temp(0);
FTensor::Tensor1<double,2> ddv_xy(0,0);
More information about the CIG-COMMITS
mailing list