[cig-commits] commit: Add residual_3D and refactor dtau_mixed to accomodate
Mercurial
hg at geodynamics.org
Tue Mar 15 01:26:18 PDT 2011
changeset: 127:3a1236261b98
user: Walter Landry <wlandry at caltech.edu>
date: Tue Mar 15 00:31:30 2011 -0700
files: StokesFACOps.h StokesFACOps/Update_V.C StokesFACOps/residual_2D.C StokesFACOps/residual_3D.C wscript
description:
Add residual_3D and refactor dtau_mixed to accomodate
diff -r 020a6c5ebd2a -r 3a1236261b98 StokesFACOps.h
--- a/StokesFACOps.h Mon Mar 14 17:08:00 2011 -0700
+++ b/StokesFACOps.h Tue Mar 15 00:31:30 2011 -0700
@@ -481,6 +481,18 @@ public:
const hier::Box &pbox,
const geom::CartesianPatchGeometry &geom);
+ void residual_3D
+ (pdat::CellData<double> &p,
+ pdat::SideData<double> &v,
+ pdat::CellData<double> &cell_viscosity,
+ pdat::CellData<double> &p_rhs,
+ pdat::SideData<double> &v_rhs,
+ pdat::CellData<double> &p_resid,
+ pdat::SideData<double> &v_resid,
+ hier::Patch &patch,
+ const hier::Box &pbox,
+ const geom::CartesianPatchGeometry &geom);
+
virtual double
computeResidualNorm(
const SAMRAIVectorReal<double>& residual,
@@ -574,26 +586,24 @@ private:
const double &theta_momentum);
/* The mixed derivative of the stress. We have to use a template
- because 2D uses Node for the edge viscosity, while 3D uses
- Edge. Written as if it is dtau_xy_dy. */
+ because 2D uses Node's for the edge viscosity, while 3D uses
+ Edge's. Written as if it is dtau_xy_dy. */
+
template<class E_data, class E_index>
double dtau_mixed(pdat::SideData<double> &v,
const E_data &edge_viscosity,
- const pdat::SideIndex ¢er_x,
- const pdat::SideIndex &up_x,
- const pdat::SideIndex &down_x,
- const pdat::SideIndex ¢er_y,
- const pdat::SideIndex &up_y,
- const E_index ¢er_e,
- const E_index &up_e,
+ const pdat::SideIndex &x,
+ const pdat::SideIndex &y,
+ const E_index &edge,
const hier::Index &ip,
+ const hier::Index &jp,
const double &dx,
const double &dy)
{
- return edge_viscosity(up_e)*((v(up_x)-v(center_x))/(dy*dy)
- + (v(up_y)-v(up_y-ip))/(dx*dy))
- - edge_viscosity(center_e)*((v(center_x)-v(down_x))/(dy*dy)
- + (v(center_y)-v(center_y-ip))/(dx*dy));
+ return ((v(x+jp)-v(x))/(dy*dy)
+ + (v(y+jp)-v(y+jp-ip))/(dx*dy)) * edge_viscosity(edge+jp)
+ - ((v(x)-v(x-jp))/(dy*dy)
+ + (v(y)-v(y-ip))/(dx*dy)) * edge_viscosity(edge);
}
/* The action of the velocity operator. It is written from the
@@ -605,26 +615,18 @@ private:
pdat::CellData<double> &cell_viscosity,
pdat::NodeData<double> &edge_viscosity,
const pdat::CellIndex ¢er,
- const pdat::CellIndex &left,
- const pdat::SideIndex ¢er_x,
- const pdat::SideIndex &right_x,
- const pdat::SideIndex &left_x,
- const pdat::SideIndex &up_x,
- const pdat::SideIndex &down_x,
- const pdat::SideIndex ¢er_y,
- const pdat::SideIndex &up_y,
- const pdat::NodeIndex ¢er_e,
- const pdat::NodeIndex &up_e,
+ const pdat::NodeIndex &edge,
+ const pdat::SideIndex &x,
+ const pdat::SideIndex &y,
const hier::Index &ip,
+ const hier::Index &jp,
const double &dx,
const double &dy)
{
- double dtau_xx_dx=
- 2*((v(right_x)-v(center_x))*cell_viscosity(center)
- - (v(center_x)-v(left_x))*cell_viscosity(left))/(dx*dx);
- double dtau_xy_dy=dtau_mixed(v,edge_viscosity,center_x,up_x,down_x,
- center_y,up_y,center_e,up_e,ip,dx,dy);
- double dp_dx=(p(center)-p(left))/dx;
+ double dtau_xx_dx=2*((v(x+ip)-v(x))*cell_viscosity(center)
+ - (v(x)-v(x-ip))*cell_viscosity(center-ip))/(dx*dx);
+ double dtau_xy_dy=dtau_mixed(v,edge_viscosity,x,y,edge,ip,jp,dx,dy);
+ double dp_dx=(p(center)-p(center-ip))/dx;
return dtau_xx_dx + dtau_xy_dy - dp_dx;
}
@@ -634,32 +636,23 @@ private:
pdat::CellData<double> &cell_viscosity,
pdat::EdgeData<double> &edge_viscosity,
const pdat::CellIndex ¢er,
- const pdat::CellIndex &left,
- const pdat::SideIndex ¢er_x,
- const pdat::SideIndex &right_x,
- const pdat::SideIndex &left_x,
- const pdat::SideIndex &up_x,
- const pdat::SideIndex &down_x,
- const pdat::SideIndex ¢er_y,
- const pdat::SideIndex &up_y,
- const pdat::SideIndex ¢er_z,
- const pdat::SideIndex &front_z,
- const pdat::EdgeIndex ¢er_e,
- const pdat::EdgeIndex &up_e,
- const pdat::EdgeIndex &front_e,
+ const pdat::EdgeIndex &edge_y,
+ const pdat::EdgeIndex &edge_z,
+ const pdat::SideIndex &x,
+ const pdat::SideIndex &y,
+ const pdat::SideIndex &z,
const hier::Index &ip,
+ const hier::Index &jp,
+ const hier::Index &kp,
const double &dx,
const double &dy,
const double &dz)
{
- double dtau_xx_dx=
- 2*((v(right_x)-v(center_x))*cell_viscosity(center)
- - (v(center_x)-v(left_x))*cell_viscosity(left))/(dx*dx);
- double dtau_xy_dy=dtau_mixed(v,edge_viscosity,center_x,up_x,down_x,
- center_y,up_y,center_e,up_e,ip,dx,dy);
- double dtau_xz_dz=dtau_mixed(v,edge_viscosity,center_x,up_x,down_x,
- center_z,front_z,center_e,front_e,ip,dx,dz);
- double dp_dx=(p(center)-p(left))/dx;
+ double dtau_xx_dx=2*((v(x+ip)-v(x))*cell_viscosity(center)
+ - (v(x)-v(x-ip))*cell_viscosity(center-ip))/(dx*dx);
+ double dtau_xy_dy=dtau_mixed(v,edge_viscosity,x,y,edge_z,ip,jp,dx,dy);
+ double dtau_xz_dz=dtau_mixed(v,edge_viscosity,x,z,edge_y,ip,kp,dx,dz);
+ double dp_dx=(p(center)-p(center-ip))/dx;
return dtau_xx_dx + dtau_xy_dy + dtau_xz_dz - dp_dx;
}
diff -r 020a6c5ebd2a -r 3a1236261b98 StokesFACOps/Update_V.C
--- a/StokesFACOps/Update_V.C Mon Mar 14 17:08:00 2011 -0700
+++ b/StokesFACOps/Update_V.C Tue Mar 15 00:31:30 2011 -0700
@@ -120,10 +120,7 @@ void SAMRAI::solv::StokesFACOps::Update_
double delta_Rx=v_rhs(center_x)
- v_operator_2D(v,p,cell_viscosity,edge_viscosity,center,
- left,center_x,
- right_x,left_x,up_x,down_x,center_y,up_y,
- center_e,up_e,
- ip,dx,dy);
+ center_e,center_x,center_y,ip,jp,dx,dy);
/* No scaling here, though there should be. */
maxres=std::max(maxres,std::fabs(delta_Rx));
diff -r 020a6c5ebd2a -r 3a1236261b98 StokesFACOps/residual_2D.C
--- a/StokesFACOps/residual_2D.C Mon Mar 14 17:08:00 2011 -0700
+++ b/StokesFACOps/residual_2D.C Tue Mar 15 00:31:30 2011 -0700
@@ -67,26 +67,16 @@ void SAMRAI::solv::StokesFACOps::residua
--left[0];
const pdat::SideIndex
- center_x(center,0,pdat::SideIndex::Lower),
- left_x(left,0,pdat::SideIndex::Lower),
- right_x(right,0,pdat::SideIndex::Lower),
- down_x(down,0,pdat::SideIndex::Lower),
- up_x(up,0,pdat::SideIndex::Lower),
- center_y(center,1,pdat::SideIndex::Lower),
- left_y(left,1,pdat::SideIndex::Lower),
- right_y(right,1,pdat::SideIndex::Lower),
- down_y(down,1,pdat::SideIndex::Lower),
- up_y(up,1,pdat::SideIndex::Lower);
+ x(center,0,pdat::SideIndex::Lower),
+ y(center,1,pdat::SideIndex::Lower);
const pdat::NodeIndex
- center_e(center,pdat::NodeIndex::LowerLeft),
- up_e(up,pdat::NodeIndex::LowerLeft),
- right_e(right,pdat::NodeIndex::LowerLeft);
+ edge(center,pdat::NodeIndex::LowerLeft);
/* p */
if(center[0]!=pbox.upper(0) && center[1]!=pbox.upper(1))
{
- double dvx_dx=(v(right_x) - v(center_x))/dx;
- double dvy_dy=(v(up_y) - v(center_y))/dy;
+ double dvx_dx=(v(x+ip) - v(x))/dx;
+ double dvy_dy=(v(y+jp) - v(y))/dy;
p_resid(center)=p_rhs(center) - dvx_dx - dvy_dy;
}
@@ -94,19 +84,16 @@ void SAMRAI::solv::StokesFACOps::residua
if(center[1]!=pbox.upper(1))
{
/* If x==0 */
- if((center[0]==pbox.lower(0) && v(left_x)==boundary_value)
- || (center[0]==pbox.upper(0)
- && v(right_x)==boundary_value))
-
+ if((center[0]==pbox.lower(0) && v(x-ip)==boundary_value)
+ || (center[0]==pbox.upper(0) && v(x+ip)==boundary_value))
{
- v_resid(center_x)=0;
+ v_resid(x)=0;
}
else
{
- v_resid(center_x)=v_rhs(center_x)
+ v_resid(x)=v_rhs(x)
- v_operator_2D(v,p,cell_viscosity,edge_viscosity,center,
- left,center_x,right_x,left_x,up_x,down_x,
- center_y,up_y,center_e,up_e,ip,dx,dy);
+ edge,x,y,ip,jp,dx,dy);
}
}
@@ -114,18 +101,16 @@ void SAMRAI::solv::StokesFACOps::residua
if(center[0]!=pbox.upper(0))
{
/* If y==0 */
- if((center[1]==pbox.lower(1) && v(down_y)==boundary_value)
- || (center[1]==pbox.upper(1) && v(up_y)==boundary_value))
+ if((center[1]==pbox.lower(1) && v(y-jp)==boundary_value)
+ || (center[1]==pbox.upper(1) && v(y+jp)==boundary_value))
{
- v_resid(center_y)=0;
+ v_resid(y)=0;
}
else
{
- v_resid(center_y)=v_rhs(center_y)
+ v_resid(y)=v_rhs(y)
- v_operator_2D(v,p,cell_viscosity,edge_viscosity,center,
- down,center_y,up_y,down_y,right_y,left_y,
- center_x,right_x,center_e,right_e,jp,
- dy,dx);
+ edge,y,x,jp,ip,dy,dx);
}
}
}
diff -r 020a6c5ebd2a -r 3a1236261b98 StokesFACOps/residual_3D.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/residual_3D.C Tue Mar 15 00:31:30 2011 -0700
@@ -0,0 +1,142 @@
+#include "StokesFACOps.h"
+
+#include IOMANIP_HEADER_FILE
+
+#include "SAMRAI/hier/BoundaryBoxUtils.h"
+#include "SAMRAI/geom/CartesianGridGeometry.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/hier/Index.h"
+#include "SAMRAI/hier/Variable.h"
+#include "SAMRAI/hier/VariableDatabase.h"
+#include "SAMRAI/pdat/CellDoubleConstantRefine.h"
+#include "SAMRAI/pdat/CellVariable.h"
+#include "SAMRAI/pdat/OutersideData.h"
+#include "SAMRAI/pdat/OutersideVariable.h"
+#include "SAMRAI/hier/PatchData.h"
+#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/solv/FACPreconditioner.h"
+#include "StokesHypreSolver.h"
+#include "SAMRAI/tbox/Array.h"
+#include "SAMRAI/tbox/MathUtilities.h"
+#include "SAMRAI/tbox/StartupShutdownManager.h"
+#include "SAMRAI/tbox/Timer.h"
+#include "SAMRAI/tbox/TimerManager.h"
+#include "SAMRAI/tbox/Utilities.h"
+#include "SAMRAI/tbox/MathUtilities.h"
+#include "SAMRAI/xfer/CoarsenAlgorithm.h"
+#include "SAMRAI/xfer/CoarsenOperator.h"
+#include "SAMRAI/xfer/CoarsenSchedule.h"
+#include "SAMRAI/xfer/RefineAlgorithm.h"
+#include "SAMRAI/xfer/RefineOperator.h"
+#include "SAMRAI/xfer/RefineSchedule.h"
+#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
+
+#include "Boundary.h"
+
+void SAMRAI::solv::StokesFACOps::residual_3D
+(pdat::CellData<double> &p,
+ pdat::SideData<double> &v,
+ pdat::CellData<double> &cell_viscosity,
+ pdat::CellData<double> &p_rhs,
+ pdat::SideData<double> &v_rhs,
+ pdat::CellData<double> &p_resid,
+ pdat::SideData<double> &v_resid,
+ hier::Patch &patch,
+ const hier::Box &pbox,
+ const geom::CartesianPatchGeometry &geom)
+{
+ const hier::Index ip(1,0,0), jp(0,1,0), kp(0,0,1);
+
+ tbox::Pointer<pdat::EdgeData<double> >
+ edge_viscosity_ptr = patch.getPatchData(edge_viscosity_id);
+ pdat::EdgeData<double> &edge_viscosity(*edge_viscosity_ptr);
+
+ double dx = geom.getDx()[0];
+ double dy = geom.getDx()[1];
+ double dz = geom.getDx()[2];
+
+ for(pdat::CellIterator ci(pbox); ci; ci++)
+ {
+ pdat::CellIndex center(*ci);
+ pdat::CellIndex up(center), down(center), right(center),
+ left(center), front(center), back(center);
+
+ ++right[0];
+ --left[0];
+ ++up[1];
+ --down[1];
+ ++front[2];
+ --back[2];
+
+ const pdat::SideIndex
+ x(center,0,pdat::SideIndex::Lower),
+ y(center,1,pdat::SideIndex::Lower),
+ z(center,2,pdat::SideIndex::Lower);
+ const pdat::EdgeIndex
+ edge_x(center,0,pdat::EdgeIndex::LowerLeft),
+ edge_y(center,1,pdat::EdgeIndex::LowerLeft),
+ edge_z(center,2,pdat::EdgeIndex::LowerLeft);
+
+ /* p */
+ if(center[0]!=pbox.upper(0) && center[1]!=pbox.upper(1)
+ && center[2]!=pbox.upper(2))
+ {
+ double dvx_dx=(v(x+ip) - v(x))/dx;
+ double dvy_dy=(v(y+jp) - v(y))/dy;
+ double dvz_dz=(v(z+kp) - v(z))/dy;
+ p_resid(center)=p_rhs(center) - dvx_dx - dvy_dy - dvz_dz;
+ }
+
+ /* vx */
+ if(center[1]!=pbox.upper(1) && center[2]!=pbox.upper(2))
+ {
+ /* If x==0 */
+ if((center[0]==pbox.lower(0) && v(x-ip)==boundary_value)
+ || (center[0]==pbox.upper(0) && v(x+ip)==boundary_value))
+ {
+ v_resid(x)=0;
+ }
+ else
+ {
+ v_resid(x)=v_rhs(x)
+ - v_operator_3D(v,p,cell_viscosity,edge_viscosity,
+ center,edge_y,edge_z,x,y,z,ip,jp,kp,dx,dy,dz);
+ }
+ }
+
+ /* vy */
+ if(center[0]!=pbox.upper(0) && center[2]!=pbox.upper(2))
+ {
+ /* If y==0 */
+ if((center[1]==pbox.lower(1) && v(y-jp)==boundary_value)
+ || (center[1]==pbox.upper(1) && v(y+jp)==boundary_value))
+ {
+ v_resid(y)=0;
+ }
+ else
+ {
+ v_resid(y)=v_rhs(y)
+ - v_operator_3D(v,p,cell_viscosity,edge_viscosity,
+ center,edge_z,edge_x,y,z,x,jp,kp,ip,dy,dz,dx);
+ }
+ }
+
+ /* vz */
+ if(center[1]!=pbox.upper(0) && center[2]!=pbox.upper(2))
+ {
+ /* If z==0 */
+ if((center[2]==pbox.lower(1) && v(z-kp)==boundary_value)
+ || (center[2]==pbox.upper(1) && v(z+kp)==boundary_value))
+ {
+ v_resid(z)=0;
+ }
+ else
+ {
+ v_resid(z)=v_rhs(z)
+ - v_operator_3D(v,p,cell_viscosity,edge_viscosity,
+ center,edge_x,edge_y,z,x,y,kp,ip,jp,dz,dx,dy);
+ }
+ }
+ }
+}
+
diff -r 020a6c5ebd2a -r 3a1236261b98 wscript
--- a/wscript Mon Mar 14 17:08:00 2011 -0700
+++ b/wscript Tue Mar 15 00:31:30 2011 -0700
@@ -35,6 +35,7 @@ def build(bld):
'StokesFACOps/checkInputPatchDataIndices.C',
'StokesFACOps/computeCompositeResidualOnLevel.C',
'StokesFACOps/residual_2D.C',
+ 'StokesFACOps/residual_3D.C',
'StokesFACOps/computeFluxOnPatch.C',
'StokesFACOps/computeResidualNorm.C',
'StokesFACOps/computeVectorWeights.C',
More information about the CIG-COMMITS
mailing list