[cig-commits] commit: Make the interface interpolation stencil one-sided, since it was not really using the other side anyway.
Mercurial
hg at geodynamics.org
Wed Feb 15 06:37:21 PST 2012
changeset: 48:85bab1d05473
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 15 06:34:25 2012 -0800
files: compute_v_on_interface.cxx
description:
Make the interface interpolation stencil one-sided, since it was not really using the other side anyway.
diff -r 06750548dbd3 -r 85bab1d05473 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Tue Feb 14 16:23:18 2012 -0800
+++ b/compute_v_on_interface.cxx Wed Feb 15 06:34:25 2012 -0800
@@ -10,7 +10,7 @@ namespace {
template<int Ny>
void add_to_matrix(const int &n_terms, const int n, double z[][Ny],
double log_eta[][Ny], const int i, const int j,
- const double &u0, const double &dxm, const double &dxp,
+ const double &u0, const double &dx,
const double &dy,
double m[], double rhs[], const int sign=1)
{
@@ -20,12 +20,9 @@ namespace {
m[0+n*n_terms]+=sign*u0;
m[1+n*n_terms]+=sign*dy;
m[2+n*n_terms]+=sign*dy*dy/2;
- m[3+n*n_terms]+=sign*dxm;
- m[4+n*n_terms]+=sign*dxm*dy;
- m[5+n*n_terms]+=sign*dxm*dxm/2;
- m[6+n*n_terms]+=sign*dxp;
- m[7+n*n_terms]+=sign*dxp*dy;
- m[8+n*n_terms]+=sign*dxp*dxp/2;
+ m[3+n*n_terms]+=sign*dx;
+ m[4+n*n_terms]+=sign*dx*dy;
+ m[5+n*n_terms]+=sign*dx*dx/2;
}
}
@@ -45,7 +42,7 @@ void interpolate_to_surface(double z[][N
const double &delta_x, const bool &aligned,
double &v0, double &dv_y, double &dv_yy)
{
- const int n_terms(9);
+ const int n_terms(6);
/* Format of m is
u0, u,y u,yy um,x um,xy um,xx up,x up,xy up,xx
@@ -83,23 +80,15 @@ void interpolate_to_surface(double z[][N
i_center+1,j+1 - i_center-1,j-1
*/
- add_to_matrix(n_terms,0,z,log_eta,i_center,j,1,0,di*dx,0,m,rhs);
- add_to_matrix(n_terms,1,z,log_eta,i_center,j+1,1,0,di*dx,h,m,rhs);
- add_to_matrix(n_terms,2,z,log_eta,i_center,j-1,1,0,di*dx,-h,m,rhs);
- add_to_matrix(n_terms,3,z,log_eta,i_center-di,j,1,0,di*(dx+h),0,m,rhs);
- add_to_matrix(n_terms,4,z,log_eta,i_center-2*di,j,1,0,di*(dx+2*h),
+ add_to_matrix(n_terms,0,z,log_eta,i_center,j,1,di*dx,0,m,rhs);
+ add_to_matrix(n_terms,1,z,log_eta,i_center,j+1,1,di*dx,h,m,rhs);
+ add_to_matrix(n_terms,2,z,log_eta,i_center,j-1,1,di*dx,-h,m,rhs);
+ add_to_matrix(n_terms,3,z,log_eta,i_center-di,j,1,di*(dx+h),0,m,rhs);
+ add_to_matrix(n_terms,4,z,log_eta,i_center-2*di,j,1,di*(dx+2*h),
0,m,rhs);
- add_to_matrix(n_terms,5,z,log_eta,i_center-di,j+1,1,0,di*(dx+h),h,m,rhs);
- add_to_matrix(n_terms,5,z,log_eta,i_center-di,j-1,1,0,di*(dx+h),
+ add_to_matrix(n_terms,5,z,log_eta,i_center-di,j+1,1,di*(dx+h),h,m,rhs);
+ add_to_matrix(n_terms,5,z,log_eta,i_center-di,j-1,1,di*(dx+h),
-h,m,rhs,-1);
-
- add_to_matrix(n_terms,6,z,log_eta,i_center+di,j,1,di*(dx-h),0,0,m,rhs);
- add_to_matrix(n_terms,7,z,log_eta,i_center+2*di,j,1,di*(dx-2*h),0,
- 0,m,rhs);
- add_to_matrix(n_terms,8,z,log_eta,i_center+di,j+1,1,di*(dx-h),0,h,m,rhs);
- add_to_matrix(n_terms,8,z,log_eta,i_center+di,j-1,1,di*(dx-h),0,
- -h,m,rhs,-1);
-
// if(j>0)
// {
@@ -118,9 +107,6 @@ void interpolate_to_surface(double z[][N
// rhs[4]=analytic(di*(dx+2*h),0);
// rhs[5]=analytic(di*(dx+h),h)-analytic(di*(dx+h),-h);
- // rhs[6]=analytic(di*(dx-h),0);
- // rhs[7]=analytic(di*(dx-2*h),0);
- // rhs[8]=analytic(di*(dx-h),h)-analytic(di*(dx-h),-h);
// }
/* At the j==0 boundary, all of the j-1 terms are not applied.
@@ -151,41 +137,33 @@ void interpolate_to_surface(double z[][N
i_center+2,j + i_center+2,j+1
*/
- add_to_matrix(n_terms,0,z,log_eta,i_center,j,1,0,-di*dx,-h/2,m,rhs);
- add_to_matrix(n_terms,1,z,log_eta,i_center,j+1,1,0,-di*dx,h/2,m,rhs);
- add_to_matrix(n_terms,2,z,log_eta,i_center-di,j,1,0,-di*(dx+h),
+ add_to_matrix(n_terms,0,z,log_eta,i_center,j,1,-di*dx,-h/2,m,rhs);
+ add_to_matrix(n_terms,1,z,log_eta,i_center,j+1,1,-di*dx,h/2,m,rhs);
+ add_to_matrix(n_terms,2,z,log_eta,i_center-di,j,1,-di*(dx+h),
-h/2,m,rhs);
- add_to_matrix(n_terms,3,z,log_eta,i_center-di,j+1,1,0,-di*(dx+h),
+ add_to_matrix(n_terms,3,z,log_eta,i_center-di,j+1,1,-di*(dx+h),
h/2,m,rhs);
- add_to_matrix(n_terms,4,z,log_eta,i_center,j+2,1,0,-di*dx,3*h/2,m,rhs);
- add_to_matrix(n_terms,4,z,log_eta,i_center,j-1,1,0,-di*dx,-3*h/2,m,rhs);
- add_to_matrix(n_terms,5,z,log_eta,i_center-2*di,j,1,0,-di*(dx+2*h),
+ add_to_matrix(n_terms,4,z,log_eta,i_center,j+2,1,-di*dx,3*h/2,m,rhs);
+ add_to_matrix(n_terms,4,z,log_eta,i_center,j-1,1,-di*dx,-3*h/2,m,rhs);
+ add_to_matrix(n_terms,5,z,log_eta,i_center-2*di,j,1,-di*(dx+2*h),
-h/2,m,rhs);
- add_to_matrix(n_terms,5,z,log_eta,i_center-2*di,j+1,1,0,-di*(dx+2*h),
- h/2,m,rhs);
-
- add_to_matrix(n_terms,6,z,log_eta,i_center+di,j,1,-di*(dx-h),0,-h/2,m,rhs);
- add_to_matrix(n_terms,7,z,log_eta,i_center+di,j+1,1,-di*(dx-h),0,
- h/2,m,rhs);
- add_to_matrix(n_terms,8,z,log_eta,i_center+2*di,j,1,-di*(dx-2*h),0,
- -h/2,m,rhs);
- add_to_matrix(n_terms,8,z,log_eta,i_center+2*di,j+1,1,-di*(dx-2*h),0,
+ add_to_matrix(n_terms,5,z,log_eta,i_center-2*di,j+1,1,-di*(dx+2*h),
h/2,m,rhs);
}
// if(j>0)
- // for(int a=0;a<9;++a)
- // {
- // std::cout << a << " "
- // << std::setw(15)
- // << std::setprecision(10)
- // << rhs[a] << " ";
- // for(int b=0;b<9;++b)
- // std::cout << std::setw(15)
- // << std::setprecision(10)
- // << m[b+9*a] << " ";
- // std::cout << "\n";
- // }
+ // for(int a=0;a<n_terms;++a)
+ // {
+ // std::cout << a << " "
+ // << std::setw(15)
+ // << std::setprecision(10)
+ // << rhs[a] << " ";
+ // for(int b=0;b<n_terms;++b)
+ // std::cout << std::setw(15)
+ // << std::setprecision(10)
+ // << m[b+n_terms*a] << " ";
+ // std::cout << "\n";
+ // }
gsl_matrix_view mv=gsl_matrix_view_array(m,n_terms,n_terms);
gsl_vector_view rhsv=gsl_vector_view_array(rhs,n_terms);
@@ -198,7 +176,7 @@ void interpolate_to_surface(double z[][N
// if(j>0)
// {
// std::cout << "vector ";
- // for(int ii=0;ii<9;++ii)
+ // for(int ii=0;ii<n_terms;++ii)
// std::cout << gsl_vector_get(result,ii) << " ";
// std::cout << "\n";
// exit(0);
More information about the CIG-COMMITS
mailing list