[cig-commits] commit: move a lot of boilerplate out of residual_2D into computeCompositeResidualOnLevel
Mercurial
hg at geodynamics.org
Tue Mar 15 01:26:15 PDT 2011
changeset: 124:7844b9f83588
user: Walter Landry <wlandry at caltech.edu>
date: Mon Mar 14 16:21:25 2011 -0700
files: StokesFACOps.h StokesFACOps/computeCompositeResidualOnLevel.C StokesFACOps/residual_2D.C wscript
description:
move a lot of boilerplate out of residual_2D into computeCompositeResidualOnLevel
diff -r b4007fa687ae -r 7844b9f83588 StokesFACOps.h
--- a/StokesFACOps.h Mon Mar 14 15:36:38 2011 -0700
+++ b/StokesFACOps.h Mon Mar 14 16:21:25 2011 -0700
@@ -467,18 +467,19 @@ public:
const SAMRAIVectorReal<double>& solution,
const SAMRAIVectorReal<double>& rhs,
int ln,
- bool error_equation_indicator)
- {
- if(d_dim.getValue()==2)
- residual_2D(residual,solution,rhs,ln,error_equation_indicator);
- }
+ bool error_equation_indicator);
- void
- residual_2D(SAMRAIVectorReal<double>& residual,
- const SAMRAIVectorReal<double>& solution,
- const SAMRAIVectorReal<double>& rhs,
- int ln,
- bool error_equation_indicator);
+ void residual_2D
+ (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(
diff -r b4007fa687ae -r 7844b9f83588 StokesFACOps/computeCompositeResidualOnLevel.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Mon Mar 14 16:21:25 2011 -0700
@@ -0,0 +1,151 @@
+#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"
+
+/*
+********************************************************************
+* FACOperatorStrategy virtual *
+* computeCompositeResidualOnLevel function *
+********************************************************************
+*/
+
+void SAMRAI::solv::StokesFACOps::computeCompositeResidualOnLevel
+(SAMRAIVectorReal<double>& residual,
+ const SAMRAIVectorReal<double>& solution,
+ const SAMRAIVectorReal<double>& rhs,
+ int ln,
+ bool error_equation_indicator) {
+
+ t_compute_composite_residual->start();
+
+ checkInputPatchDataIndices();
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (residual.getPatchHierarchy() != d_hierarchy
+ || solution.getPatchHierarchy() != d_hierarchy
+ || rhs.getPatchHierarchy() != d_hierarchy) {
+ TBOX_ERROR(d_object_name << ": Vector hierarchy does not match\n"
+ "internal hierarchy.");
+ }
+#endif
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
+
+ /*
+ * Set up the bc helper so that when we use a refine schedule
+ * to fill ghosts, the correct data is operated on.
+ */
+ const int p_id = solution.getComponentDescriptorIndex(0);
+ const int v_id = solution.getComponentDescriptorIndex(1);
+ p_refine_patch_strategy.setTargetDataId(p_id);
+ v_refine_patch_strategy.setTargetDataId(v_id);
+ // v_refine_patch_strategy.setHomogeneousBc(error_equation_indicator);
+
+ /*
+ * Assumptions:
+ * 1. Data does not yet exist in ghost boundaries.
+ * 2. Residual data on next finer grid (if any)
+ * has been computed already.
+ * 3. Flux data from next finer grid (if any) has
+ * been computed but has not been coarsened to
+ * this level.
+ *
+ * Steps:
+ * S1. Fill solution ghost data by refinement
+ * or setting physical boundary conditions.
+ * This also brings in information from coarser
+ * to form the composite grid flux.
+ * S2. Compute flux on ln.
+ * S3. If next finer is available,
+ * Coarsen flux data on next finer level,
+ * overwriting flux computed from coarse data.
+ * S4. Compute residual data from flux.
+ */
+
+ /* S1. Fill solution ghost data. */
+
+ set_boundaries(v_id,ln,error_equation_indicator);
+ if (ln > d_ln_min) {
+ /* Fill from current, next coarser level and physical boundary */
+ xeqScheduleGhostFill(p_id, v_id, ln);
+ } else {
+ /* Fill from current and physical boundary */
+ xeqScheduleGhostFillNoCoarse(p_id, v_id, ln);
+ }
+
+ /*
+ * S4. Compute residual on patches in level.
+ */
+
+ for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
+ tbox::Pointer<hier::Patch> patch = *pi;
+ tbox::Pointer<pdat::CellData<double> >
+ p_ptr = solution.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::SideData<double> >
+ v_ptr = solution.getComponentPatchData(1, *patch);
+ tbox::Pointer<pdat::CellData<double> >
+ cell_viscosity_ptr = patch->getPatchData(cell_viscosity_id);
+ tbox::Pointer<pdat::CellData<double> >
+ p_rhs_ptr = rhs.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::SideData<double> >
+ v_rhs_ptr = rhs.getComponentPatchData(1, *patch);
+ tbox::Pointer<pdat::CellData<double> >
+ p_resid_ptr = residual.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::SideData<double> >
+ v_resid_ptr = residual.getComponentPatchData(1, *patch);
+
+ hier::Box pbox=patch->getBox();
+ pbox.growUpper(hier::IntVector::getOne(d_dim));
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = patch->getPatchGeometry();
+
+ switch(d_dim.getValue())
+ {
+ case 2:
+ residual_2D(*p_ptr,*v_ptr,*cell_viscosity_ptr,*p_rhs_ptr,*v_rhs_ptr,
+ *p_resid_ptr,*v_resid_ptr,*patch,pbox,*geom);
+ break;
+ case 3:
+ // residual_3D();
+ break;
+ default:
+ abort();
+ }
+ }
+
+ /* We also need to set the boundaries of the rhs so that coarsening
+ works correctly. */
+ const int v_rhs_id = rhs.getComponentDescriptorIndex(1);
+ set_boundaries(v_rhs_id,ln,true);
+ xeqScheduleGhostFillNoCoarse(invalid_id, v_rhs_id, ln);
+
+ t_compute_composite_residual->stop();
+}
+
diff -r b4007fa687ae -r 7844b9f83588 StokesFACOps/residual_2D.C
--- a/StokesFACOps/residual_2D.C Mon Mar 14 15:36:38 2011 -0700
+++ b/StokesFACOps/residual_2D.C Mon Mar 14 16:21:25 2011 -0700
@@ -1,12 +1,3 @@
-/*************************************************************************
- *
- * This file is part of the SAMRAI distribution. For full copyright
- * information, see COPYRIGHT and COPYING.LESSER.
- *
- * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
- * Description: Operator class for cell-centered scalar Stokes using FAC
- *
- ************************************************************************/
#include "StokesFACOps.h"
#include IOMANIP_HEADER_FILE
@@ -42,195 +33,101 @@
#include "Boundary.h"
#include "dRc_dp.h"
-/*
-********************************************************************
-* FACOperatorStrategy virtual *
-* computeCompositeResidualOnLevel function *
-********************************************************************
-*/
void SAMRAI::solv::StokesFACOps::residual_2D
-(SAMRAIVectorReal<double>& residual,
- const SAMRAIVectorReal<double>& solution,
- const SAMRAIVectorReal<double>& rhs,
- int ln,
- bool error_equation_indicator) {
+(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), jp(0,1);
- t_compute_composite_residual->start();
+ tbox::Pointer<pdat::NodeData<double> >
+ edge_viscosity_ptr = patch.getPatchData(edge_viscosity_id);
+ pdat::NodeData<double> &edge_viscosity(*edge_viscosity_ptr);
- checkInputPatchDataIndices();
-#ifdef DEBUG_CHECK_ASSERTIONS
- if (residual.getPatchHierarchy() != d_hierarchy
- || solution.getPatchHierarchy() != d_hierarchy
- || rhs.getPatchHierarchy() != d_hierarchy) {
- TBOX_ERROR(d_object_name << ": Vector hierarchy does not match\n"
- "internal hierarchy.");
- }
-#endif
- tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
+ double dx = geom.getDx()[0];
+ double dy = geom.getDx()[1];
- /*
- * Set up the bc helper so that when we use a refine schedule
- * to fill ghosts, the correct data is operated on.
- */
- const int p_id = solution.getComponentDescriptorIndex(0);
- const int v_id = solution.getComponentDescriptorIndex(1);
- p_refine_patch_strategy.setTargetDataId(p_id);
- v_refine_patch_strategy.setTargetDataId(v_id);
- // v_refine_patch_strategy.setHomogeneousBc(error_equation_indicator);
+ for(pdat::CellIterator ci(pbox); ci; ci++)
+ {
+ pdat::CellIndex center(*ci);
+ pdat::CellIndex up(center), down(center), right(center),
+ left(center);
- /*
- * Assumptions:
- * 1. Data does not yet exist in ghost boundaries.
- * 2. Residual data on next finer grid (if any)
- * has been computed already.
- * 3. Flux data from next finer grid (if any) has
- * been computed but has not been coarsened to
- * this level.
- *
- * Steps:
- * S1. Fill solution ghost data by refinement
- * or setting physical boundary conditions.
- * This also brings in information from coarser
- * to form the composite grid flux.
- * S2. Compute flux on ln.
- * S3. If next finer is available,
- * Coarsen flux data on next finer level,
- * overwriting flux computed from coarse data.
- * S4. Compute residual data from flux.
- */
+ ++up[1];
+ --down[1];
+ ++right[0];
+ --left[0];
- /* S1. Fill solution ghost data. */
+ 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);
+ const pdat::NodeIndex
+ center_e(center,pdat::NodeIndex::LowerLeft),
+ up_e(up,pdat::NodeIndex::LowerLeft),
+ right_e(right,pdat::NodeIndex::LowerLeft);
- set_boundaries(v_id,ln,error_equation_indicator);
- if (ln > d_ln_min) {
- /* Fill from current, next coarser level and physical boundary */
- xeqScheduleGhostFill(p_id, v_id, ln);
- } else {
- /* Fill from current and physical boundary */
- xeqScheduleGhostFillNoCoarse(p_id, v_id, ln);
- }
+ /* 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;
+ p_resid(center)=p_rhs(center) - dvx_dx - dvy_dy;
+ }
- /*
- * S4. Compute residual on patches in level.
- */
+ /* vx */
+ 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))
+
+ {
+ v_resid(center_x)=0;
+ }
+ else
+ {
+ v_resid(center_x)=v_rhs(center_x)
+ - v_operator(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);
+ }
+ }
- const hier::Index ip(1,0), jp(0,1);
- for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
- tbox::Pointer<hier::Patch> patch = *pi;
- tbox::Pointer<pdat::CellData<double> >
- p_ptr = solution.getComponentPatchData(0, *patch);
- pdat::CellData<double> &p(*p_ptr);
- tbox::Pointer<pdat::SideData<double> >
- v_ptr = solution.getComponentPatchData(1, *patch);
- pdat::SideData<double> &v(*v_ptr);
- tbox::Pointer<pdat::CellData<double> >
- cell_viscosity_ptr = patch->getPatchData(cell_viscosity_id);
- pdat::CellData<double> &cell_viscosity(*cell_viscosity_ptr);
- tbox::Pointer<pdat::NodeData<double> >
- edge_viscosity_ptr = patch->getPatchData(edge_viscosity_id);
- pdat::NodeData<double> &edge_viscosity(*edge_viscosity_ptr);
- tbox::Pointer<pdat::CellData<double> >
- p_rhs_ptr = rhs.getComponentPatchData(0, *patch);
- pdat::CellData<double> &p_rhs(*p_rhs_ptr);
- tbox::Pointer<pdat::SideData<double> >
- v_rhs_ptr = rhs.getComponentPatchData(1, *patch);
- pdat::SideData<double> &v_rhs(*v_rhs_ptr);
- tbox::Pointer<pdat::CellData<double> >
- p_resid_ptr = residual.getComponentPatchData(0, *patch);
- pdat::CellData<double> &p_resid(*p_resid_ptr);
- tbox::Pointer<pdat::SideData<double> >
- v_resid_ptr = residual.getComponentPatchData(1, *patch);
- pdat::SideData<double> &v_resid(*v_resid_ptr);
-
- hier::Box pbox=patch->getBox();
- tbox::Pointer<geom::CartesianPatchGeometry> geom = patch->getPatchGeometry();
- double dx = geom->getDx()[0];
- double dy = geom->getDx()[1];
-
- pbox.growUpper(hier::IntVector::getOne(d_dim));
- for(pdat::CellIterator ci(pbox); ci; ci++)
- {
- pdat::CellIndex center(*ci);
- pdat::CellIndex up(center), down(center), right(center),
- left(center);
-
- ++up[1];
- --down[1];
- ++right[0];
- --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);
- const pdat::NodeIndex
- center_e(center,pdat::NodeIndex::LowerLeft),
- up_e(up,pdat::NodeIndex::LowerLeft),
- right_e(right,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;
- p_resid(center)=p_rhs(center) - dvx_dx - dvy_dy;
- }
-
- /* vx */
- 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))
-
- {
- v_resid(center_x)=0;
- }
- else
- {
- v_resid(center_x)=v_rhs(center_x)
- - v_operator(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);
- }
- }
-
- /* vy */
- 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))
- {
- v_resid(center_y)=0;
- }
- else
- {
- v_resid(center_y)=v_rhs(center_y)
- - v_operator(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);
- }
- }
- }
- }
-
- /* We also need to set the boundaries of the rhs so that coarsening
- works correctly. */
- const int v_rhs_id = rhs.getComponentDescriptorIndex(1);
- set_boundaries(v_rhs_id,ln,true);
- xeqScheduleGhostFillNoCoarse(invalid_id, v_rhs_id, ln);
-
- t_compute_composite_residual->stop();
+ /* vy */
+ 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))
+ {
+ v_resid(center_y)=0;
+ }
+ else
+ {
+ v_resid(center_y)=v_rhs(center_y)
+ - v_operator(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);
+ }
+ }
+ }
}
diff -r b4007fa687ae -r 7844b9f83588 wscript
--- a/wscript Mon Mar 14 15:36:38 2011 -0700
+++ b/wscript Mon Mar 14 16:21:25 2011 -0700
@@ -33,6 +33,7 @@ def build(bld):
'set_V_boundary.C',
'StokesFACOps/StokesFACOps.C',
'StokesFACOps/checkInputPatchDataIndices.C',
+ 'StokesFACOps/computeCompositeResidualOnLevel.C',
'StokesFACOps/residual_2D.C',
'StokesFACOps/computeFluxOnPatch.C',
'StokesFACOps/computeResidualNorm.C',
More information about the CIG-COMMITS
mailing list