[cig-commits] commit: Add 3D. 2D still works, but 3D is completely untested.
Mercurial
hg at geodynamics.org
Sat Mar 19 20:16:48 PDT 2011
changeset: 141:66f150fd3e3b
user: Walter Landry <wlandry at caltech.edu>
date: Sat Mar 19 13:39:48 2011 -0700
files: P_MDPI_Refine.C StokesFACOps.h StokesFACOps/smooth_Gerya.C StokesFACOps/smooth_Tackley_2D.C StokesFACOps/smooth_V.C StokesFACOps/smooth_V_2D.C StokesFACOps/smooth_V_3D.C dRc_dp.h dRm_dv.h wscript
description:
Add 3D. 2D still works, but 3D is completely untested.
diff -r 532a610b132a -r 66f150fd3e3b P_MDPI_Refine.C
--- a/P_MDPI_Refine.C Fri Mar 18 20:55:09 2011 -0700
+++ b/P_MDPI_Refine.C Sat Mar 19 13:39:48 2011 -0700
@@ -121,8 +121,8 @@ void SAMRAI::geom::P_MDPI_Refine::refine
y(c_fine,1,pdat::SideIndex::Lower);
double dRc_dp_weight=
- dRc_dp(fine_box,c_fine,x,y,
- cell_viscosity,edge_viscosity,v,dx,dy);
+ dRc_dp_2D(fine_box,c_fine,x,y,
+ cell_viscosity,edge_viscosity,v,dx,dy);
if(c_fine==fine)
diff -r 532a610b132a -r 66f150fd3e3b StokesFACOps.h
--- a/StokesFACOps.h Fri Mar 18 20:55:09 2011 -0700
+++ b/StokesFACOps.h Sat Mar 19 13:39:48 2011 -0700
@@ -558,6 +558,14 @@ private:
double residual_tolerance = -1.0);
void
+ smooth_Tackley_3D(
+ SAMRAIVectorReal<double>& error,
+ const SAMRAIVectorReal<double>& residual,
+ int ln,
+ int num_sweeps,
+ double residual_tolerance = -1.0);
+
+ void
smooth_Gerya(
SAMRAIVectorReal<double>& error,
const SAMRAIVectorReal<double>& residual,
@@ -565,7 +573,7 @@ private:
int num_sweeps,
double residual_tolerance = -1.0);
- void smooth_V
+ void smooth_V_2D
(const int &axis,
const hier::Box &pbox,
tbox::Pointer<geom::CartesianPatchGeometry> &geom,
@@ -581,6 +589,21 @@ private:
pdat::CellData<double> &cell_viscosity,
pdat::NodeData<double> &edge_viscosity,
const double &theta_momentum);
+
+ void smooth_V_3D
+ (const int &ix,
+ const hier::Box &pbox,
+ tbox::Pointer<geom::CartesianPatchGeometry> &geom,
+ pdat::CellData<double> &p,
+ pdat::SideData<double> &v,
+ pdat::SideData<double> &v_rhs,
+ pdat::CellData<double> &cell_viscosity,
+ pdat::EdgeData<double> &edge_viscosity,
+ const pdat::CellIndex ¢er,
+ const double dx[3],
+ const double &theta_momentum,
+ const hier::Index pp[3],
+ double &maxres);
/* The mixed derivative of the stress. We have to use a template
because 2D uses Node's for the edge viscosity, while 3D uses
diff -r 532a610b132a -r 66f150fd3e3b StokesFACOps/smooth_Gerya.C
--- a/StokesFACOps/smooth_Gerya.C Fri Mar 18 20:55:09 2011 -0700
+++ b/StokesFACOps/smooth_Gerya.C Sat Mar 19 13:39:48 2011 -0700
@@ -146,15 +146,15 @@ void SAMRAI::solv::StokesFACOps::smooth_
/* Update v */
if(j<pbox.upper(1)+1)
{
- smooth_V(0,pbox,geom,center,ip,jp,
- p,v,v_rhs,maxres,dx,dy,cell_viscosity,
- edge_viscosity,theta_momentum);
+ smooth_V_2D(0,pbox,geom,center,ip,jp,
+ p,v,v_rhs,maxres,dx,dy,cell_viscosity,
+ edge_viscosity,theta_momentum);
}
if(i<pbox.upper(0)+1)
{
- smooth_V(1,pbox,geom,center,jp,ip,
- p,v,v_rhs,maxres,dy,dx,cell_viscosity,
- edge_viscosity,theta_momentum);
+ smooth_V_2D(1,pbox,geom,center,jp,ip,
+ p,v,v_rhs,maxres,dy,dx,cell_viscosity,
+ edge_viscosity,theta_momentum);
}
}
}
diff -r 532a610b132a -r 66f150fd3e3b StokesFACOps/smooth_Tackley_2D.C
--- a/StokesFACOps/smooth_Tackley_2D.C Fri Mar 18 20:55:09 2011 -0700
+++ b/StokesFACOps/smooth_Tackley_2D.C Sat Mar 19 13:39:48 2011 -0700
@@ -115,9 +115,9 @@ void SAMRAI::solv::StokesFACOps::smooth_
center[1]=j;
/* Update v */
- smooth_V(0,pbox,geom,center,ip,jp,
- p,v,v_rhs,maxres,dx,dy,cell_viscosity,
- edge_viscosity,theta_momentum);
+ smooth_V_2D(0,pbox,geom,center,ip,jp,
+ p,v,v_rhs,maxres,dx,dy,cell_viscosity,
+ edge_viscosity,theta_momentum);
}
}
}
@@ -169,9 +169,9 @@ void SAMRAI::solv::StokesFACOps::smooth_
center[1]=j;
/* Update v */
- smooth_V(1,pbox,geom,center,jp,ip,
- p,v,v_rhs,maxres,dy,dx,cell_viscosity,
- edge_viscosity,theta_momentum);
+ smooth_V_2D(1,pbox,geom,center,jp,ip,
+ p,v,v_rhs,maxres,dy,dx,cell_viscosity,
+ edge_viscosity,theta_momentum);
}
}
}
@@ -234,7 +234,7 @@ void SAMRAI::solv::StokesFACOps::smooth_
maxres=std::max(maxres,std::fabs(delta_R_continuity));
dp(center)=delta_R_continuity*theta_continuity
- /dRc_dp(pbox,center,x,y,cell_viscosity,edge_viscosity,v,dx,dy);
+ /dRc_dp_2D(pbox,center,x,y,cell_viscosity,edge_viscosity,v,dx,dy);
p(center)+=dp(center);
}
}
@@ -285,8 +285,8 @@ void SAMRAI::solv::StokesFACOps::smooth_
|| (center[0]==pbox.upper(0)
&& v(x+ip)==boundary_value)))
v(x)+=(dp(center) - dp(center-ip))
- /(dx*dRm_dv(cell_viscosity,edge_viscosity,center,
- center-ip,edge+jp,edge,dx,dy));
+ /(dx*dRm_dv_2D(cell_viscosity,edge_viscosity,center,
+ center-ip,edge+jp,edge,dx,dy));
}
if(center[0]<pbox.upper(0))
{
@@ -294,8 +294,8 @@ void SAMRAI::solv::StokesFACOps::smooth_
|| (center[1]==pbox.upper(1)
&& v(y+jp)==boundary_value)))
v(y)+=(dp(center) - dp(center-jp))
- /(dy*dRm_dv(cell_viscosity,edge_viscosity,center,
- center-jp,edge+ip,edge,dy,dx));
+ /(dy*dRm_dv_2D(cell_viscosity,edge_viscosity,center,
+ center-jp,edge+ip,edge,dy,dx));
}
}
}
diff -r 532a610b132a -r 66f150fd3e3b StokesFACOps/smooth_V.C
--- a/StokesFACOps/smooth_V.C Fri Mar 18 20:55:09 2011 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,115 +0,0 @@
-#include "StokesFACOps.h"
-#include "Boundary.h"
-#include "dRc_dp.h"
-/*
-********************************************************************
-* Updates one component of the velocity during a red-black *
-* Gauss-Seidel iteration. *
-********************************************************************
-*/
-void SAMRAI::solv::StokesFACOps::smooth_V
-(const int &axis,
- const hier::Box &pbox,
- tbox::Pointer<geom::CartesianPatchGeometry> &geom,
- const pdat::CellIndex ¢er,
- const hier::Index &ip,
- const hier::Index &jp,
- pdat::CellData<double> &p,
- pdat::SideData<double> &v,
- pdat::SideData<double> &v_rhs,
- double &maxres,
- const double &dx,
- const double &dy,
- pdat::CellData<double> &cell_viscosity,
- pdat::NodeData<double> &edge_viscosity,
- const double &theta_momentum)
-{
- const int off_axis=(axis==0) ? 1 : 0;
-
- const pdat::SideIndex x(center,axis,pdat::SideIndex::Lower),
- y(center,off_axis,pdat::SideIndex::Lower);
- const pdat::NodeIndex edge(center,pdat::NodeIndex::LowerLeft);
-
- /* If at the 'x' boundaries, leave vx as is */
- if(!((center[axis]==pbox.lower(axis) && v(x-ip)==boundary_value)
- || (center[axis]==pbox.upper(axis)+1 && v(x+ip)==boundary_value)))
- {
- /* If at the boundary, set things up so that the derivative does
- not change. */
- hier::Index offset(0,0);
- offset[axis]=2;
- bool set_lower_boundary(false), set_upper_boundary(false);
- double dv_lower(0), dv_upper(0);
- if(center[axis]==pbox.lower(axis)+1
- && !geom->getTouchesRegularBoundary(axis,0))
- {
- set_lower_boundary=true;
- dv_lower=v(x-offset) - v(x);
- }
- if(center[axis]==pbox.upper(axis)
- && !geom->getTouchesRegularBoundary(axis,1))
- {
- set_upper_boundary=true;
- dv_upper=v(x+offset) - v(x);
- }
-
- double C_vx=dRm_dv(cell_viscosity,edge_viscosity,center,center-ip,
- edge+jp,edge,dx,dy);
-
- double delta_Rx=v_rhs(x)
- - v_operator_2D(v,p,cell_viscosity,edge_viscosity,center,
- edge,x,y,ip,jp,dx,dy);
-
- /* No scaling here, though there should be. */
- maxres=std::max(maxres,std::fabs(delta_Rx));
-
- v(x)+=delta_Rx*theta_momentum/C_vx;
-
- // tbox::plog << "smooth V "
- // << dx << " "
- // << axis << " "
- // << center[0] << " "
- // << center[1] << " "
- // << v(x) << " "
- // << delta_Rx << " "
- // << v(x+ip) << " "
- // << v(x-ip) << " "
- // << v(x+jp) << " "
- // << v(x-jp) << " "
- // << v(y+jp) << " "
- // << v(y+jp-ip) << " "
- // << v(y) << " "
- // << v(y-ip) << " "
- // << "\n";
-
- /* Set the boundary elements so that the
- derivative is unchanged. */
- if(set_lower_boundary)
- {
- v(x-offset)=v(x) + dv_lower;
- // tbox::plog << "offset minus "
- // << dx << " "
- // << axis << " "
- // << (x-offset)[0] << " "
- // << (x-offset)[1] << " "
- // << v(x) << " "
- // << v(x-offset) << " "
- // << dv_lower << " "
- // << "\n";
- }
- if(set_upper_boundary)
- {
- v(x+offset)=v(x) + dv_upper;
- // tbox::plog << "offset plus "
- // << dx << " "
- // << axis << " "
- // << (x+offset)[0] << " "
- // << (x+offset)[1] << " "
- // << v(x) << " "
- // << v(x+offset) << " "
- // << dv_upper << " "
- // << "\n";
- }
- }
-}
-
diff -r 532a610b132a -r 66f150fd3e3b StokesFACOps/smooth_V_2D.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/smooth_V_2D.C Sat Mar 19 13:39:48 2011 -0700
@@ -0,0 +1,80 @@
+#include "StokesFACOps.h"
+#include "Boundary.h"
+#include "dRm_dv.h"
+/*
+********************************************************************
+* Updates one component of the velocity during a red-black *
+* Gauss-Seidel iteration. *
+********************************************************************
+*/
+void SAMRAI::solv::StokesFACOps::smooth_V_2D
+(const int &axis,
+ const hier::Box &pbox,
+ tbox::Pointer<geom::CartesianPatchGeometry> &geom,
+ const pdat::CellIndex ¢er,
+ const hier::Index &ip,
+ const hier::Index &jp,
+ pdat::CellData<double> &p,
+ pdat::SideData<double> &v,
+ pdat::SideData<double> &v_rhs,
+ double &maxres,
+ const double &dx,
+ const double &dy,
+ pdat::CellData<double> &cell_viscosity,
+ pdat::NodeData<double> &edge_viscosity,
+ const double &theta_momentum)
+{
+ const int off_axis=(axis==0) ? 1 : 0;
+
+ const pdat::SideIndex x(center,axis,pdat::SideIndex::Lower),
+ y(center,off_axis,pdat::SideIndex::Lower);
+ const pdat::NodeIndex edge(center,pdat::NodeIndex::LowerLeft);
+
+ /* If at the 'x' boundaries, leave vx as is */
+ if(!((center[axis]==pbox.lower(axis) && v(x-ip)==boundary_value)
+ || (center[axis]==pbox.upper(axis)+1 && v(x+ip)==boundary_value)))
+ {
+ /* If at the boundary, set things up so that the derivative does
+ not change. */
+ hier::Index offset(0,0);
+ offset[axis]=2;
+ bool set_lower_boundary(false), set_upper_boundary(false);
+ double dv_lower(0), dv_upper(0);
+ if(center[axis]==pbox.lower(axis)+1
+ && !geom->getTouchesRegularBoundary(axis,0))
+ {
+ set_lower_boundary=true;
+ dv_lower=v(x-offset) - v(x);
+ }
+ if(center[axis]==pbox.upper(axis)
+ && !geom->getTouchesRegularBoundary(axis,1))
+ {
+ set_upper_boundary=true;
+ dv_upper=v(x+offset) - v(x);
+ }
+
+ double C_vx=dRm_dv_2D(cell_viscosity,edge_viscosity,center,center-ip,
+ edge+jp,edge,dx,dy);
+
+ double delta_Rx=v_rhs(x)
+ - v_operator_2D(v,p,cell_viscosity,edge_viscosity,center,
+ edge,x,y,ip,jp,dx,dy);
+
+ /* No scaling here, though there should be. */
+ maxres=std::max(maxres,std::fabs(delta_Rx));
+
+ v(x)+=delta_Rx*theta_momentum/C_vx;
+
+ /* Set the boundary elements so that the derivative is
+ unchanged. */
+ if(set_lower_boundary)
+ {
+ v(x-offset)=v(x) + dv_lower;
+ }
+ if(set_upper_boundary)
+ {
+ v(x+offset)=v(x) + dv_upper;
+ }
+ }
+}
+
diff -r 532a610b132a -r 66f150fd3e3b StokesFACOps/smooth_V_3D.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/smooth_V_3D.C Sat Mar 19 13:39:48 2011 -0700
@@ -0,0 +1,79 @@
+#include "StokesFACOps.h"
+#include "Boundary.h"
+#include "dRc_dp.h"
+/*
+********************************************************************
+* Updates one component of the velocity during a red-black *
+* Gauss-Seidel iteration. *
+********************************************************************
+*/
+void SAMRAI::solv::StokesFACOps::smooth_V_3D
+(const int &ix,
+ const hier::Box &pbox,
+ tbox::Pointer<geom::CartesianPatchGeometry> &geom,
+ pdat::CellData<double> &p,
+ pdat::SideData<double> &v,
+ pdat::SideData<double> &v_rhs,
+ pdat::CellData<double> &cell_viscosity,
+ pdat::EdgeData<double> &edge_viscosity,
+ const pdat::CellIndex ¢er,
+ const double Dx[3],
+ const double &theta_momentum,
+ const hier::Index pp[3],
+ double &maxres)
+{
+ const int iy((ix+1)%3), iz((ix+2)%3);
+ const pdat::SideIndex x(center,ix,pdat::SideIndex::Lower),
+ y(center,iy,pdat::SideIndex::Lower),
+ z(center,iz,pdat::SideIndex::Lower);
+ const pdat::EdgeIndex edge_y(center,iy,pdat::EdgeIndex::LowerLeft),
+ edge_z(center,iz,pdat::EdgeIndex::LowerLeft);
+
+ /* If at the 'x' boundaries, leave vx as is */
+ if(!((center[ix]==pbox.lower(ix) && v(x-pp[ix])==boundary_value)
+ || (center[ix]==pbox.upper(ix)+1 && v(x+pp[ix])==boundary_value)))
+ {
+ /* If at the boundary, set things up so that the derivative does
+ not change. */
+ hier::Index offset(0,0,0);
+ offset[ix]=2;
+ bool set_lower_boundary(false), set_upper_boundary(false);
+ double dv_lower(0), dv_upper(0);
+ if(center[ix]==pbox.lower(ix)+1
+ && !geom->getTouchesRegularBoundary(ix,0))
+ {
+ set_lower_boundary=true;
+ dv_lower=v(x-offset) - v(x);
+ }
+ if(center[ix]==pbox.upper(ix) && !geom->getTouchesRegularBoundary(ix,1))
+ {
+ set_upper_boundary=true;
+ dv_upper=v(x+offset) - v(x);
+ }
+
+ double C_vx=dRm_dv_3D(cell_viscosity,edge_viscosity,center,center-pp[ix],
+ edge_y+pp[iy],edge_y,edge_z+pp[iz],edge_z,
+ Dx[ix],Dx[iy],Dx[iz]);
+
+ double delta_Rx=v_rhs(x)
+ - v_operator_3D(v,p,cell_viscosity,edge_viscosity,center,edge_y,edge_z,
+ x,y,z,pp[ix],pp[iy],pp[iz],Dx[ix],Dx[iy],Dx[iz]);
+
+ /* No scaling here, though there should be. */
+ maxres=std::max(maxres,std::fabs(delta_Rx));
+
+ v(x)+=delta_Rx*theta_momentum/C_vx;
+
+ /* Set the boundary elements so that the derivative is
+ unchanged. */
+ if(set_lower_boundary)
+ {
+ v(x-offset)=v(x) + dv_lower;
+ }
+ if(set_upper_boundary)
+ {
+ v(x+offset)=v(x) + dv_upper;
+ }
+ }
+}
+
diff -r 532a610b132a -r 66f150fd3e3b dRc_dp.h
--- a/dRc_dp.h Fri Mar 18 20:55:09 2011 -0700
+++ b/dRc_dp.h Sat Mar 19 13:39:48 2011 -0700
@@ -12,15 +12,15 @@
equation, so we use Tackley's method to chain together
derivatives */
-inline double dRc_dp(const SAMRAI::hier::Box &pbox,
- const SAMRAI::pdat::CellIndex ¢er,
- const SAMRAI::pdat::SideIndex &x,
- const SAMRAI::pdat::SideIndex &y,
- SAMRAI::pdat::CellData<double> &cell_viscosity,
- SAMRAI::pdat::NodeData<double> &edge_viscosity,
- SAMRAI::pdat::SideData<double> &v,
- const double &dx,
- const double &dy)
+inline double dRc_dp_2D(const SAMRAI::hier::Box &pbox,
+ const SAMRAI::pdat::CellIndex ¢er,
+ const SAMRAI::pdat::SideIndex &x,
+ const SAMRAI::pdat::SideIndex &y,
+ SAMRAI::pdat::CellData<double> &cell_viscosity,
+ SAMRAI::pdat::NodeData<double> &edge_viscosity,
+ SAMRAI::pdat::SideData<double> &v,
+ const double &dx,
+ const double &dy)
{
const SAMRAI::hier::Index ip(1,0), jp(0,1);
const SAMRAI::pdat::NodeIndex
@@ -35,24 +35,59 @@ inline double dRc_dp(const SAMRAI::hier:
double result(0);
if(!(center[0]==pbox.lower(0) && v(x-ip)==boundary_value))
- result+=dRc_dvx_p * dRm_dp_xp/dRm_dv(cell_viscosity,edge_viscosity,
- center+ip,
- center,up_e+ip,center_e+ip,dx,dy);
+ result+=dRc_dvx_p * dRm_dp_xp/dRm_dv_2D(cell_viscosity,edge_viscosity,
+ center+ip,
+ center,up_e+ip,center_e+ip,dx,dy);
if(!(center[0]==pbox.upper(0)+1 && v(x+ip)==boundary_value))
- result+=dRc_dvx_m * dRm_dp_xm/dRm_dv(cell_viscosity,edge_viscosity,
- center,center-ip,up_e,center_e,dx,dy);
+ result+=dRc_dvx_m * dRm_dp_xm/dRm_dv_2D(cell_viscosity,edge_viscosity,
+ center,center-ip,up_e,center_e,dx,dy);
if(!(center[1]==pbox.lower(1) && v(y-jp)==boundary_value))
- result+=dRc_dvy_p * dRm_dp_yp/dRm_dv(cell_viscosity,edge_viscosity,center+jp,
- center,right_e+jp,center_e+jp,dy,dx);
+ result+=dRc_dvy_p * dRm_dp_yp/dRm_dv_2D(cell_viscosity,edge_viscosity,center+jp,
+ center,right_e+jp,center_e+jp,dy,dx);
if(!(center[1]==pbox.upper(1)+1 && v(y+jp)==boundary_value))
- result+=dRc_dvy_m * dRm_dp_ym/dRm_dv(cell_viscosity,edge_viscosity,center,
- center-jp,right_e,center_e,dy,dx);
+ result+=dRc_dvy_m * dRm_dp_ym/dRm_dv_2D(cell_viscosity,edge_viscosity,center,
+ center-jp,right_e,center_e,dy,dx);
return result;
}
+inline double dRc_dp_3D(const SAMRAI::hier::Box &pbox,
+ const SAMRAI::pdat::CellIndex ¢er,
+ SAMRAI::pdat::CellData<double> &cell_viscosity,
+ SAMRAI::pdat::EdgeData<double> &edge_viscosity,
+ SAMRAI::pdat::SideData<double> &v,
+ const double Dx[3],
+ const SAMRAI::hier::Index pp[3])
+{
+ double result(0);
+ for(int ix=0;ix<3;++ix)
+ {
+ const int iy((ix+1)%3), iz((ix+2)%3);
+ const SAMRAI::pdat::SideIndex x(center,ix,SAMRAI::pdat::SideIndex::Lower);
+ const SAMRAI::pdat::EdgeIndex
+ center_y(center,iy,SAMRAI::pdat::EdgeIndex::LowerLeft),
+ up_y(center_y+pp[iy]),
+ center_z(center,iz,SAMRAI::pdat::EdgeIndex::LowerLeft),
+ front_z(center_z+pp[iz]);
+
+ if(!(center[ix]==pbox.lower(ix) && v(x-pp[ix])==boundary_value))
+ result+=-1.0/(dRm_dv_3D(cell_viscosity,edge_viscosity,
+ center+pp[ix],center,up_y+pp[ix],center_y+pp[ix],
+ front_z+pp[ix],center_z+pp[ix],
+ Dx[ix],Dx[iy],Dx[iz])
+ * Dx[ix]*Dx[ix]);
+
+
+ if(!(center[ix]==pbox.upper(ix)+1 && v(x+pp[ix])==boundary_value))
+ result+=-1.0/(dRm_dv_3D(cell_viscosity,edge_viscosity,
+ center,center-pp[ix],up_y,center_y,
+ front_z,center_z,Dx[ix],Dx[iy],Dx[iz])
+ * Dx[ix]*Dx[ix]);
+ }
+ return result;
+}
#endif
diff -r 532a610b132a -r 66f150fd3e3b dRm_dv.h
--- a/dRm_dv.h Fri Mar 18 20:55:09 2011 -0700
+++ b/dRm_dv.h Sat Mar 19 13:39:48 2011 -0700
@@ -2,25 +2,42 @@
#define GAMR_DRM_DV
#include "SAMRAI/pdat/CellData.h"
-#include "SAMRAI/pdat/NodeData.h"
+#include "SAMRAI/pdat/EdgeData.h"
/* The derivative of the momentum equation w/respect to velocity. It
is written from the perspective of vx(center_x), but pass in
different values for center etc. to get vy or vx(!center_x). */
-inline double dRm_dv(SAMRAI::pdat::CellData<double> &cell_viscosity,
- SAMRAI::pdat::NodeData<double> &edge_viscosity,
- const SAMRAI::pdat::CellIndex ¢er,
- const SAMRAI::pdat::CellIndex &left,
- const SAMRAI::pdat::NodeIndex &up_e,
- const SAMRAI::pdat::NodeIndex ¢er_e,
- const double &dx,
- const double &dy)
+template<class E_data, class E_index>
+double dRm_dv_2D(SAMRAI::pdat::CellData<double> &cell_viscosity,
+ E_data &edge_viscosity,
+ const SAMRAI::pdat::CellIndex ¢er,
+ const SAMRAI::pdat::CellIndex &left,
+ const E_index &up_e,
+ const E_index ¢er_e,
+ const double &dx,
+ const double &dy)
{
return -2*(cell_viscosity(center) + cell_viscosity(left))/(dx*dx)
- (edge_viscosity(up_e) + edge_viscosity(center_e))/(dy*dy);
}
+inline double dRm_dv_3D(SAMRAI::pdat::CellData<double> &cell_viscosity,
+ SAMRAI::pdat::EdgeData<double> &edge_viscosity,
+ const SAMRAI::pdat::CellIndex ¢er,
+ const SAMRAI::pdat::CellIndex &left,
+ const SAMRAI::pdat::EdgeIndex &up_y,
+ const SAMRAI::pdat::EdgeIndex ¢er_y,
+ const SAMRAI::pdat::EdgeIndex &front_z,
+ const SAMRAI::pdat::EdgeIndex ¢er_z,
+ const double &dx,
+ const double &dy,
+ const double &dz)
+{
+ return dRm_dv_2D(cell_viscosity,edge_viscosity,center,left,up_y,center_y,dx,dy)
+ - (edge_viscosity(front_z) + edge_viscosity(center_z))/(dz*dz);
+}
+
diff -r 532a610b132a -r 66f150fd3e3b wscript
--- a/wscript Fri Mar 18 20:55:09 2011 -0700
+++ b/wscript Sat Mar 19 13:39:48 2011 -0700
@@ -49,11 +49,13 @@ def build(bld):
'StokesFACOps/restrictSolution.C',
'StokesFACOps/smoothError.C',
'StokesFACOps/smooth_Tackley_2D.C',
+ 'StokesFACOps/smooth_Tackley_3D.C',
'StokesFACOps/smooth_Gerya.C',
'StokesFACOps/set_boundaries.C',
'StokesFACOps/solveCoarsestLevel.C',
'StokesFACOps/solveCoarsestLevel_HYPRE.C',
- 'StokesFACOps/smooth_V.C',
+ 'StokesFACOps/smooth_V_2D.C',
+ 'StokesFACOps/smooth_V_3D.C',
'StokesFACOps/xeqScheduleFluxCoarsen.C',
'StokesFACOps/xeqScheduleGhostFill.C',
'StokesFACOps/xeqScheduleGhostFillNoCoarse.C',
More information about the CIG-COMMITS
mailing list