[cig-commits] commit: Make the stencil smaller since it is only first order.
Mercurial
hg at geodynamics.org
Tue May 8 05:50:57 PDT 2012
changeset: 159:7eed54710466
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Tue May 08 05:50:38 2012 -0700
files: compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx compute_coefficients/compute_v_on_interface/compute_values.hxx
description:
Make the stencil smaller since it is only first order.
diff -r ea067e692ffb -r 7eed54710466 compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx
--- a/compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx Mon May 07 08:50:48 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_1st_derivs.hxx Tue May 08 05:50:38 2012 -0700
@@ -6,61 +6,39 @@
#include <gsl/gsl_linalg.h>
template<int ny>
-void compute_1st_derivs(const int &d, const std::vector<Valid> &derivs,
+void compute_1st_derivs(const int &d, const Valid &V,
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)
+ /* TODO: this does not seem right at the boundary. */
+ switch(d0)
{
- switch(d0)
+ case 0:
+ if(V.i==nx-1 || V.i==-1)
{
- case 0:
- if(derivs[i].i==nx-1 || derivs[i].i==-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 || derivs[i].j==-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();
+ dv_pm(d,d0)=0;
}
- m[0+3*i]=1;
- m[1+3*i]=derivs[i].diff(0);
- m[2+3*i]=derivs[i].diff(1);
+ else
+ {
+ dv_pm(d,d0)=vel(z,log_eta,V.i+1,V.j) - vel(z,log_eta,V.i,V.j);
+ }
+ break;
+ case 1:
+ if(V.j==ny-1 || V.j==-1)
+ {
+ dv_pm(d,d0)=0;
+ }
+ else
+ {
+ dv_pm(d,d0)=vel(z,log_eta,V.i,V.j+1) - vel(z,log_eta,V.i,V.j);
+ }
+ break;
+ default:
+ abort();
}
-
- 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);
- gsl_vector_free(x);
- gsl_permutation_free(perm);
}
#endif
diff -r ea067e692ffb -r 7eed54710466 compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx Mon May 07 08:50:48 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_dv_dtt.cxx Tue May 08 05:50:38 2012 -0700
@@ -56,13 +56,35 @@ void compute_dv_dtt(const double zx[Nx+1
int starti((pos(0)-dx)/h), startj((pos(1)-dy)/h);
- /* d/dx^2 */
+ /* v */
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),
- d_diff, diff;
+ FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), diff;
diff(a)=p(a)-pos(a);
+
+ if(d==0)
+ {
+ diff(a)-=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
+ {
+ diff(a)-=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/dx */
+ 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_d(i*h+dx+h/2,j*h+dy),
+ d_diff;
d_diff(a)=p_d(a)-pos(a);
if(d==0)
@@ -76,10 +98,6 @@ void compute_dv_dtt(const double zx[Nx+1
sign(distx[i][j])==sign(distx[i+1][j]),
i,j,d_diff);
}
- /* v */
- diff(a)-=1.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
{
@@ -97,27 +115,19 @@ void compute_dv_dtt(const double zx[Nx+1
d_diff(a)-=length*norm(a)*sign(disty[i][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);
-
- /* v */
- diff(a)-=1.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 */
+ /* d/dy */
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,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),
- d_diff;
+ FTensor::Tensor1<double,2> p_d(i*h+dx,j*h+dy+h/2), d_diff;
d_diff(a)=p_d(a)-pos(a);
if(d==0)
{
- /* dv */
bool v_d;
if(j==min_y || j==max_y)
v_d=true;
@@ -130,7 +140,6 @@ void compute_dv_dtt(const double zx[Nx+1
}
else
{
- /* dv */
if(j<max_y)
{
d_diff(a)-=length*norm(a)*sign(disty[i][j]);
@@ -147,9 +156,6 @@ void compute_dv_dtt(const double zx[Nx+1
sort(valid.begin(), valid.end());
- FTensor::Tensor3_christof<bool,2,2> is_set(false,false,false,false,
- false,false);
-
/* Compute everything */
/* v */
@@ -164,49 +170,25 @@ void compute_dv_dtt(const double zx[Nx+1
const int pm=(V.sign==-1 ? 0 : 1);
switch(values[pm].size())
{
- case 0:
- case 1:
- case 2:
+ default:
values[pm].push_back(V);
break;
- case 6:
+ case 3:
break;
- default:
- int x(0),y(0);
- std::vector<int> x_coords={V.i};
- std::vector<int> y_coords={V.j};
- for(auto &value: values[pm])
+ case 2:
+ if((V.i!=values[pm][0].i || V.i !=values[pm][1].i)
+ && (V.j!=values[pm][0].j || V.j !=values[pm][1].j))
{
- if(V.i==value.i)
- ++x;
- if(V.j==value.j)
- ++y;
- if(values[pm].size()==5)
- {
- 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())
- y_coords.push_back(value.j);
- }
- }
- if(x<3 && y<3 && !(values[pm].size()==5 &&
- (x_coords.size()<3 || y_coords.size()<3)))
- values[pm].push_back(V);
- if(values[pm].size()==6)
- {
+ values[pm].push_back(V);
if(d==0)
{
- eta_pm[pm]=exp(log_etax[values[pm].begin()->i]
- [values[pm].begin()->j]);
+ eta_pm[pm]=exp(log_etax[values[pm][0].i][values[pm][0].j]);
compute_values(0,values[pm],Nx+1,zx,log_etax,
norm,pm==0 ? -length : length,v_pm[pm]);
}
else
{
- eta_pm[pm]=exp(log_etay[values[pm].begin()->i]
- [values[pm].begin()->j]);
+ eta_pm[pm]=exp(log_etay[values[pm][0].i][values[pm][0].j]);
compute_values(1,values[pm],Nx,zy,log_etay,
norm,pm==0 ? -length : length,v_pm[pm]);
}
@@ -218,68 +200,41 @@ void compute_dv_dtt(const double zx[Nx+1
for(int d0=0;d0<2;++d0)
{
/* First derivative dv */
- std::vector<Valid> derivs[2];
- for(auto &V: d_valid(d0))
- {
- if(derivs[0].size()==3 && derivs[1].size()==3)
- break;
-
- if(V.sign==0)
- continue;
- 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;
- }
- }
- }
+ for(int pm=0; pm<2; ++pm)
+ for(auto &V: d_valid(d0))
+ {
+ if(V.sign==2*pm-1)
+ {
+ switch(d)
+ {
+ case 0:
+ compute_1st_derivs(0,V,d0,Nx+1,zx,log_etax,dv_pm[pm]);
+ break;
+ case 1:
+ compute_1st_derivs(1,V,d0,Nx,zy,log_etay,dv_pm[pm]);
+ break;
+ default:
+ abort();
+ break;
+ }
+ break;
+ }
+ }
}
}
/* dv */
FTensor::Tensor1<double,2> dv_xy(0,0);
- for(int d=0;d<2;++d)
- dv_xy(a)+=dv_pm[d](a,b)*tangent(b)*eta_pm[d];
+ for(int pm=0;pm<2;++pm)
+ dv_xy(a)+=dv_pm[pm](a,b)*tangent(b)*eta_pm[pm];
- dv(a)=nt_to_xy(b,a)*dv_xy(b);
-
- dv(a)/=(eta_pm[0]+eta_pm[1])*h;
+ dv(a)=nt_to_xy(b,a)*dv_xy(b)/(h*(eta_pm[0]+eta_pm[1]));
/* v */
FTensor::Tensor1<double,2> v_xy(0,0);
for(int d=0;d<2;++d)
v_xy(a)+=v_pm[d](a)*eta_pm[d];
- v(a)=nt_to_xy(b,a)*v_xy(b);
-
- FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
- double eta_jump(eta_pm[1]-eta_pm[0]);
- FTensor::Tensor1<double,2> dz_n_jump;
- dz_n_jump(a)=-eta_jump*dv(b)*ddel(a,b);
-
- v(a)-=2*length*dz_n_jump(a)/3;
- v(a)/=(eta_pm[0] + eta_pm[1]);
+ v(a)=nt_to_xy(b,a)*v_xy(b)/(eta_pm[0] + eta_pm[1]);
}
diff -r ea067e692ffb -r 7eed54710466 compute_coefficients/compute_v_on_interface/compute_values.hxx
--- a/compute_coefficients/compute_v_on_interface/compute_values.hxx Mon May 07 08:50:48 2012 -0700
+++ b/compute_coefficients/compute_v_on_interface/compute_values.hxx Tue May 08 05:50:38 2012 -0700
@@ -15,25 +15,23 @@ void compute_values(const int &d, const
FTensor::Tensor1<double,2> &v_pm)
{
/* Create and invert the Vandermonde matrix */
- assert(values.size()==6);
- double m[6*6], rhs[6];
- for(int i=0;i<6;++i)
+ const unsigned int degree(3);
+ assert(values.size()==3);
+ double m[degree*degree], rhs[degree];
+ for(unsigned int i=0;i<degree;++i)
{
rhs[i]=vel(z,log_eta,values[i].i,values[i].j);
- m[0+6*i]=1;
- m[1+6*i]=values[i].diff(0);
- m[2+6*i]=values[i].diff(1);
- m[3+6*i]=values[i].diff(0)*values[i].diff(0)/2;
- m[4+6*i]=values[i].diff(1)*values[i].diff(1)/2;
- m[5+6*i]=values[i].diff(0)*values[i].diff(1);
+ m[0+degree*i]=1;
+ m[1+degree*i]=values[i].diff(0);
+ m[2+degree*i]=values[i].diff(1);
}
- gsl_matrix_view mv=gsl_matrix_view_array(m,6,6);
- gsl_vector_view rhsv=gsl_vector_view_array(rhs,6);
- gsl_vector *x=gsl_vector_alloc(6);
+ gsl_matrix_view mv=gsl_matrix_view_array(m,degree,degree);
+ gsl_vector_view rhsv=gsl_vector_view_array(rhs,degree);
+ gsl_vector *x=gsl_vector_alloc(degree);
int s;
- gsl_permutation *perm=gsl_permutation_alloc(6);
+ gsl_permutation *perm=gsl_permutation_alloc(degree);
gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,x);
@@ -44,15 +42,9 @@ void compute_values(const int &d, const
const FTensor::Index<'b',2> b;
dv_xy(0)=gsl_vector_get(x,1);
dv_xy(1)=gsl_vector_get(x,2);
- double dv=dv_xy(a)*norm(a)*signed_length/2;
+ double dv=dv_xy(a)*norm(a);
- FTensor::Tensor2_symmetric<double,2> ddv_xy;
- ddv_xy(0,0)=gsl_vector_get(x,3);
- ddv_xy(1,1)=gsl_vector_get(x,4);
- ddv_xy(0,1)=gsl_vector_get(x,5);
- double ddv=ddv_xy(a,b)*norm(a)*norm(b)*signed_length*signed_length/4;
-
- v_pm(d)=(3*v - 5*dv + 1.5*ddv)/3;
+ v_pm(d)=v - dv*signed_length;
gsl_vector_free(x);
gsl_permutation_free(perm);
More information about the CIG-COMMITS
mailing list