[cig-commits] commit: Split StokesFACOps into individual files
Mercurial
hg at geodynamics.org
Fri Feb 25 14:13:10 PST 2011
changeset: 19:04ca745a8fed
user: Walter Landry <wlandry at caltech.edu>
date: Sun Jan 02 00:12:38 2011 -0800
files: Makefile StokesFACOps.C StokesFACOps.h StokesFACOps/StokesFACOps.C StokesFACOps/checkInputPatchDataIndices.C StokesFACOps/computeCompositeResidualOnLevel.C StokesFACOps/computeFluxOnPatch.C StokesFACOps/computeResidualNorm.C StokesFACOps/computeResidualOnPatch.C StokesFACOps/computeVectorWeights.C StokesFACOps/deallocateOperatorState.C StokesFACOps/ewingFixFlux.C StokesFACOps/finalizeCallback.C StokesFACOps/initializeOperatorState.C StokesFACOps/postprocessOneCycle.C StokesFACOps/prolongErrorAndCorrect.C StokesFACOps/redOrBlackSmoothingOnPatch.C StokesFACOps/restrictResidual.C StokesFACOps/restrictSolution.C StokesFACOps/smoothError.C StokesFACOps/smoothErrorByRedBlack.C StokesFACOps/solveCoarsestLevel.C StokesFACOps/solveCoarsestLevel_HYPRE.C StokesFACOps/xeqScheduleFluxCoarsen.C StokesFACOps/xeqScheduleGhostFill.C StokesFACOps/xeqScheduleGhostFillNoCoarse.C StokesFACOps/xeqScheduleProlongation.C StokesFACOps/xeqScheduleRRestriction.C StokesFACOps/xeqScheduleURestriction.C
description:
Split StokesFACOps into individual files
diff -r 0e7b866a1bb0 -r 04ca745a8fed Makefile
--- a/Makefile Sat Jan 01 23:30:12 2011 -0800
+++ b/Makefile Sun Jan 02 00:12:38 2011 -0800
@@ -29,7 +29,33 @@ CXX_OBJS = main.o FACStokes/FACStok
FACStokes/packDerivedDataIntoDoubleBuffer.o \
FACStokes/resetHierarchyConfiguration.o \
FACStokes/setupPlotter.o \
- FACStokes/solveStokes.o StokesFACOps.o \
+ FACStokes/solveStokes.o \
+ StokesFACOps/StokesFACOps.o \
+ StokesFACOps/checkInputPatchDataIndices.o \
+ StokesFACOps/computeCompositeResidualOnLevel.o \
+ StokesFACOps/computeFluxOnPatch.o \
+ StokesFACOps/computeResidualNorm.o \
+ StokesFACOps/computeResidualOnPatch.o \
+ StokesFACOps/computeVectorWeights.o \
+ StokesFACOps/deallocateOperatorState.o \
+ StokesFACOps/ewingFixFlux.o \
+ StokesFACOps/finalizeCallback.o \
+ StokesFACOps/initializeOperatorState.o \
+ StokesFACOps/postprocessOneCycle.o \
+ StokesFACOps/prolongErrorAndCorrect.o \
+ StokesFACOps/redOrBlackSmoothingOnPatch.o \
+ StokesFACOps/restrictResidual.o \
+ StokesFACOps/restrictSolution.o \
+ StokesFACOps/smoothError.o \
+ StokesFACOps/smoothErrorByRedBlack.o \
+ StokesFACOps/solveCoarsestLevel.o \
+ StokesFACOps/solveCoarsestLevel_HYPRE.o \
+ StokesFACOps/xeqScheduleFluxCoarsen.o \
+ StokesFACOps/xeqScheduleGhostFill.o \
+ StokesFACOps/xeqScheduleGhostFillNoCoarse.o \
+ StokesFACOps/xeqScheduleProlongation.o \
+ StokesFACOps/xeqScheduleRRestriction.o \
+ StokesFACOps/xeqScheduleURestriction.o \
StokesFACSolver/StokesFACSolver.o \
StokesFACSolver/StokesFACSolver_Destructor.o \
StokesFACSolver/createVectorWrappers.o \
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps.C
--- a/StokesFACOps.C Sat Jan 01 23:30:12 2011 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,2830 +0,0 @@
-/*************************************************************************
- *
- * 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
- *
- ************************************************************************/
-#ifndef included_solv_StokesFACOps_C
-#define included_solv_StokesFACOps_C
-
-#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"
-
-#ifndef SAMRAI_INLINE
-#include "StokesFACOps.I"
-#endif
-
-namespace SAMRAI {
-namespace solv {
-
-tbox::Pointer<pdat::CellVariable<double> >
-StokesFACOps::s_cell_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
-
-tbox::Pointer<pdat::SideVariable<double> >
-StokesFACOps::s_flux_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
-
-tbox::Pointer<pdat::OutersideVariable<double> >
-StokesFACOps::s_oflux_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
-
-tbox::StartupShutdownManager::Handler
-StokesFACOps::s_finalize_handler(
- 0,
- 0,
- 0,
- StokesFACOps::finalizeCallback,
- tbox::StartupShutdownManager::priorityVariables);
-
-extern "C" {
-
-#ifdef __INTEL_COMPILER
-#pragma warning (disable:1419)
-#endif
-
-void F77_FUNC(compfluxvardc2d, COMPFLUXVARDC2D) (
- double* xflux,
- double* yflux,
- const int* fluxgi,
- const int* fluxgj,
- const double* xdiff_coef,
- const double* ydiff_coef,
- const int* dcgi,
- const int* dcgj,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const double* dx);
-void F77_FUNC(compfluxcondc2d, COMPFLUXCONDC2D) (
- double* xflux,
- double* yflux,
- const int* fluxgi,
- const int* fluxgj,
- const double & diff_coef,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const double* dx);
-void F77_FUNC(rbgswithfluxmaxvardcvarsf2d, RBGSWITHFLUXMAXVARDCVARSF2D) (
- const double* xflux,
- const double* yflux,
- const int* fluxgi,
- const int* fluxgj,
- const double* xdiff_coef,
- const double* ydiff_coef,
- const int* dcgi,
- const int* dcgj,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- const double* scalar_field,
- const int* scalar_field_gi,
- const int* scalar_field_gj,
- double* soln,
- const int* solngi,
- const int* solngj,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const double* dx,
- const int* offset,
- const double* maxres);
-void F77_FUNC(rbgswithfluxmaxcondcvarsf2d, RBGSWITHFLUXMAXCONDCVARSF2D) (
- const double* xflux,
- const double* yflux,
- const int* fluxgi,
- const int* fluxgj,
- const double & dc,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- const double* scalar_field,
- const int* scalar_field_gi,
- const int* scalar_field_gj,
- double* soln,
- const int* solngi,
- const int* solngj,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const double* dx,
- const int* offset,
- const double* maxres);
-void F77_FUNC(rbgswithfluxmaxvardcconsf2d, RBGSWITHFLUXMAXVARDCCONSF2D) (
- const double* xflux,
- const double* yflux,
- const int* fluxgi,
- const int* fluxgj,
- const double* xdiff_coef,
- const double* ydiff_coef,
- const int* dcgi,
- const int* dcgj,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- const double & scalar_field,
- double* soln,
- const int* solngi,
- const int* solngj,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const double* dx,
- const int* offset,
- const double* maxres);
-void F77_FUNC(rbgswithfluxmaxcondcconsf2d, RBGSWITHFLUXMAXCONDCCONSF2D) (
- const double* xflux,
- const double* yflux,
- const int* fluxgi,
- const int* fluxgj,
- const double & dc,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- const double & scalar_field,
- double* soln,
- const int* solngi,
- const int* solngj,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const double* dx,
- const int* offset,
- const double* maxres);
-void F77_FUNC(compresvarsca2d, COMPRESVARSCA2D) (
- const double* xflux,
- const double* yflux,
- const int* fluxgi,
- const int* fluxgj,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- double* residual,
- const int* residualgi,
- const int* residualgj,
- const double* scalar_field,
- const int* scalar_field_gi,
- const int* scalar_field_gj,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const double* dx);
-void F77_FUNC(compresconsca2d, COMPRESCONSCA2D) (
- const double* xflux,
- const double* yflux,
- const int* fluxgi,
- const int* fluxgj,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- double* residual,
- const int* residualgi,
- const int* residualgj,
- const double & scalar_field,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const double* dx);
-void F77_FUNC(ewingfixfluxvardc2d, EWINGFIXFLUXVARDC2D) (
- const double* xflux,
- const double* yflux,
- const int* fluxgi,
- const int* fluxgj,
- const double* xdiff_coef,
- const double* ydiff_coef,
- const int* dcgi,
- const int* dcgj,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* location_index,
- const int* ratio_to_coarser,
- const int* blower,
- const int* bupper,
- const double* dx);
-void F77_FUNC(ewingfixfluxcondc2d, EWINGFIXFLUXCONDC2D) (
- const double* xflux,
- const double* yflux,
- const int* fluxgi,
- const int* fluxgj,
- const double & diff_coef,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* location_index,
- const int* ratio_to_coarser,
- const int* blower,
- const int* bupper,
- const double* dx);
-
-void F77_FUNC(compfluxvardc3d, COMPFLUXVARDC3D) (
- double* xflux,
- double* yflux,
- double* zflux,
- const int* fluxgi,
- const int* fluxgj,
- const int* fluxgk,
- const double* xdiff_coef,
- const double* ydiff_coef,
- const double* zdiff_coef,
- const int* dcgi,
- const int* dcgj,
- const int* dcgk,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* solngk,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* kfirst,
- const int* klast,
- const double* dx);
-void F77_FUNC(compfluxcondc3d, COMPFLUXCONDC3D) (
- double* xflux,
- double* yflux,
- double* zflux,
- const int* fluxgi,
- const int* fluxgj,
- const int* fluxgk,
- const double & diff_coef,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* solngk,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* kfirst,
- const int* klast,
- const double* dx);
-void F77_FUNC(rbgswithfluxmaxvardcvarsf3d, RBGSWITHFLUXMAXVARDCVARSF3D) (
- const double* xflux,
- const double* yflux,
- const double* zflux,
- const int* fluxgi,
- const int* fluxgj,
- const int* fluxgk,
- const double* xdiff_coef,
- const double* ydiff_coef,
- const double* zdiff_coef,
- const int* dcgi,
- const int* dcgj,
- const int* dcgk,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- const int* rhsgk,
- const double* scalar_field,
- const int* scalar_field_gi,
- const int* scalar_field_gj,
- const int* scalar_field_gk,
- double* soln,
- const int* solngi,
- const int* solngj,
- const int* solngk,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* kfirst,
- const int* klast,
- const double* dx,
- const int* offset,
- const double* maxres);
-void F77_FUNC(rbgswithfluxmaxcondcvarsf3d, RBGSWITHFLUXMAXCONDCVARSF3D) (
- const double* xflux,
- const double* yflux,
- const double* zflux,
- const int* fluxgi,
- const int* fluxgj,
- const int* fluxgk,
- const double & dc,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- const int* rhsgk,
- const double* scalar_field,
- const int* scalar_field_gi,
- const int* scalar_field_gj,
- const int* scalar_field_gk,
- double* soln,
- const int* solngi,
- const int* solngj,
- const int* solngk,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* kfirst,
- const int* klast,
- const double* dx,
- const int* offset,
- const double* maxres);
-void F77_FUNC(rbgswithfluxmaxvardcconsf3d, RBGSWITHFLUXMAXVARDCCONSF3D) (
- const double* xflux,
- const double* yflux,
- const double* zflux,
- const int* fluxgi,
- const int* fluxgj,
- const int* fluxgk,
- const double* xdiff_coef,
- const double* ydiff_coef,
- const double* zdiff_coef,
- const int* dcgi,
- const int* dcgj,
- const int* dcgk,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- const int* rhsgk,
- const double & scalar_field,
- double* soln,
- const int* solngi,
- const int* solngj,
- const int* solngk,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* kfirst,
- const int* klast,
- const double* dx,
- const int* offset,
- const double* maxres);
-void F77_FUNC(rbgswithfluxmaxcondcconsf3d, RBGSWITHFLUXMAXCONDCCONSF3D) (
- const double* xflux,
- const double* yflux,
- const double* zflux,
- const int* fluxgi,
- const int* fluxgj,
- const int* fluxgk,
- const double & dc,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- const int* rhsgk,
- const double & scalar_field,
- double* soln,
- const int* solngi,
- const int* solngj,
- const int* solngk,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* kfirst,
- const int* klast,
- const double* dx,
- const int* offset,
- const double* maxres);
-void F77_FUNC(compresvarsca3d, COMPRESVARSCA3D) (
- const double* xflux,
- const double* yflux,
- const double* zflux,
- const int* fluxgi,
- const int* fluxgj,
- const int* fluxgk,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- const int* rhsgk,
- double* residual,
- const int* residualgi,
- const int* residualgj,
- const int* residualgk,
- const double* scalar_field,
- const int* scalar_field_gi,
- const int* scalar_field_gj,
- const int* scalar_field_gk,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* solngk,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* kfirst,
- const int* klast,
- const double* dx);
-void F77_FUNC(compresconsca3d, COMPRESCONSCA3D) (
- const double* xflux,
- const double* yflux,
- const double* zflux,
- const int* fluxgi,
- const int* fluxgj,
- const int* fluxgk,
- const double* rhs,
- const int* rhsgi,
- const int* rhsgj,
- const int* rhsgk,
- double* residual,
- const int* residualgi,
- const int* residualgj,
- const int* residualgk,
- const double & scalar_field,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* solngk,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* kfirst,
- const int* klast,
- const double* dx);
-void F77_FUNC(ewingfixfluxvardc3d, EWINGFIXFLUXVARDC3D) (
- const double* xflux,
- const double* yflux,
- const double* zflux,
- const int* fluxgi,
- const int* fluxgj,
- const int* fluxgk,
- const double* xdiff_coef,
- const double* ydiff_coef,
- const double* zdiff_coef,
- const int* dcgi,
- const int* dcgj,
- const int* dcgk,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* solngk,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* kfirst,
- const int* klast,
- const int* location_index,
- const int* ratio_to_coarser,
- const int* blower,
- const int* bupper,
- const double* dx);
-void F77_FUNC(ewingfixfluxcondc3d, EWINGFIXFLUXCONDC3D) (
- const double* xflux,
- const double* yflux,
- const double* zflux,
- const int* fluxgi,
- const int* fluxgj,
- const int* fluxgk,
- const double & diff_coef,
- const double* soln,
- const int* solngi,
- const int* solngj,
- const int* solngk,
- const int* ifirst,
- const int* ilast,
- const int* jfirst,
- const int* jlast,
- const int* kfirst,
- const int* klast,
- const int* location_index,
- const int* ratio_to_coarser,
- const int* blower,
- const int* bupper,
- const double* dx);
-
-}
-
-/*
- ********************************************************************
- * Constructor. *
- ********************************************************************
- */
-StokesFACOps::StokesFACOps(
- const tbox::Dimension& dim,
- const std::string& object_name,
- tbox::Pointer<tbox::Database> database):
- d_dim(dim),
- d_object_name(object_name),
- d_hierarchy(),
- d_ln_min(-1),
- d_ln_max(-1),
- d_cf_boundary(),
- d_stokes_spec(object_name + "::Stokes specs"),
- d_smoothing_choice("redblack"),
- d_coarse_solver_choice(
-#ifdef HAVE_HYPRE
- "hypre"
-#else
- "redblack"
-#endif
-
- ),
- d_cf_discretization("Ewing"),
- d_prolongation_method("CONSTANT_REFINE"),
- d_coarse_solver_tolerance(1.e-8),
- d_coarse_solver_max_iterations(10),
- d_residual_tolerance_during_smoothing(-1.0),
- d_flux_id(-1),
-#ifdef HAVE_HYPRE
- d_hypre_solver(dim,
- object_name + "::hypre_solver",
- database && database->isDatabase("hypre_solver") ?
- database->getDatabase("hypre_solver"):
- tbox::Pointer<tbox::Database>(NULL)),
-#endif
- d_physical_bc_coef(NULL),
- d_context(hier::VariableDatabase::getDatabase()
- ->getContext(object_name + "::PRIVATE_CONTEXT")),
- d_cell_scratch_id(-1),
- d_flux_scratch_id(-1),
- d_oflux_scratch_id(-1),
- d_prolongation_refine_operator(),
- d_prolongation_refine_algorithm(),
- d_prolongation_refine_schedules(),
- d_urestriction_coarsen_operator(),
- d_urestriction_coarsen_algorithm(),
- d_urestriction_coarsen_schedules(),
- d_rrestriction_coarsen_operator(),
- d_rrestriction_coarsen_algorithm(),
- d_rrestriction_coarsen_schedules(),
- d_flux_coarsen_operator(),
- d_flux_coarsen_algorithm(),
- d_flux_coarsen_schedules(),
- d_ghostfill_refine_operator(),
- d_ghostfill_refine_algorithm(),
- d_ghostfill_refine_schedules(),
- d_ghostfill_nocoarse_refine_operator(),
- d_ghostfill_nocoarse_refine_algorithm(),
- d_ghostfill_nocoarse_refine_schedules(),
- d_bc_helper(dim,
- d_object_name + "::bc helper"),
- d_enable_logging(false),
- d_preconditioner(NULL),
- d_hopscell(),
- d_hopsside()
-{
-
- t_restrict_solution = tbox::TimerManager::getManager()->
- getTimer("solv::StokesFACOps::restrictSolution()");
- t_restrict_residual = tbox::TimerManager::getManager()->
- getTimer("solv::StokesFACOps::restrictResidual()");
- t_prolong = tbox::TimerManager::getManager()->
- getTimer("solv::StokesFACOps::prolongErrorAndCorrect()");
- t_smooth_error = tbox::TimerManager::getManager()->
- getTimer("solv::StokesFACOps::smoothError()");
- t_solve_coarsest = tbox::TimerManager::getManager()->
- getTimer("solv::StokesFACOps::solveCoarsestLevel()");
- t_compute_composite_residual = tbox::TimerManager::getManager()->
- getTimer("solv::StokesFACOps::computeCompositeResidualOnLevel()");
- t_compute_residual_norm = tbox::TimerManager::getManager()->
- getTimer("solv::StokesFACOps::computeResidualNorm()");
-
- if (d_dim == tbox::Dimension(1) || d_dim > tbox::Dimension(3)) {
- TBOX_ERROR("StokesFACOps : DIM == 1 or > 3 not implemented yet.\n");
- }
-
- if (s_cell_scratch_var[dim.getValue() - 1].isNull()) {
- TBOX_ASSERT(s_cell_scratch_var[dim.getValue() - 1].isNull());
- TBOX_ASSERT(s_cell_scratch_var[dim.getValue() - 1].isNull());
-
- std::ostringstream ss;
- ss << "StokesFACOps::private_cell_scratch" << dim.getValue();
- s_cell_scratch_var[dim.getValue() - 1] = new pdat::CellVariable<double>
- (dim, ss.str());
- ss.str("");
- ss << "StokesFACOps::private_flux_scratch" << dim.getValue();
- s_flux_scratch_var[dim.getValue() - 1] = new pdat::SideVariable<double>
- (dim, ss.str());
- ss.str("");
- ss << "StokesFACOps::private_oflux_scratch" << dim.getValue();
- s_oflux_scratch_var[dim.getValue() - 1] = new pdat::OutersideVariable<double>
- (dim, ss.str());
- }
-
- hier::VariableDatabase* vdb = hier::VariableDatabase::getDatabase();
- d_cell_scratch_id = vdb->
- registerVariableAndContext(s_cell_scratch_var[dim.getValue() - 1],
- d_context,
- hier::IntVector::getOne(dim));
- d_flux_scratch_id = vdb->
- registerVariableAndContext(s_flux_scratch_var[dim.getValue() - 1],
- d_context,
- hier::IntVector::getZero(d_dim));
- d_oflux_scratch_id = vdb->
- registerVariableAndContext(s_oflux_scratch_var[dim.getValue() - 1],
- d_context,
- hier::IntVector::getZero(d_dim));
-
- /*
- * Some variables initialized by default are overriden by input.
- */
- if (database) {
-
- d_coarse_solver_choice =
- database->getStringWithDefault("coarse_solver_choice",
- d_coarse_solver_choice);
- d_coarse_solver_tolerance =
- database->getDoubleWithDefault("coarse_solver_tolerance",
- d_coarse_solver_tolerance);
- d_coarse_solver_max_iterations =
- database->getIntegerWithDefault("coarse_solver_max_iterations",
- d_coarse_solver_max_iterations);
- d_smoothing_choice =
- database->getStringWithDefault("smoothing_choice",
- d_smoothing_choice);
-
- d_cf_discretization =
- database->getStringWithDefault("cf_discretization",
- d_cf_discretization);
-
- d_prolongation_method =
- database->getStringWithDefault("prolongation_method",
- d_prolongation_method);
-
- d_enable_logging =
- database->getBoolWithDefault("enable_logging",
- d_enable_logging);
-
- }
-
- /*
- * Check input validity and correctness.
- */
- checkInputPatchDataIndices();
-
-}
-
-StokesFACOps::~StokesFACOps(
- void)
-{
-}
-
-/*
- ************************************************************************
- * FACOperatorStrategy virtual initializeOperatorState function. *
- * *
- * Set internal variables to correspond to the solution passed in. *
- * Look up transfer operators. *
- ************************************************************************
- */
-
-void StokesFACOps::initializeOperatorState(
- const SAMRAIVectorReal<double>& solution,
- const SAMRAIVectorReal<double>& rhs)
-{
- deallocateOperatorState();
- int ln;
- hier::VariableDatabase* vdb = hier::VariableDatabase::getDatabase();
-
- d_hierarchy = solution.getPatchHierarchy();
- d_ln_min = solution.getCoarsestLevelNumber();
- d_ln_max = solution.getFinestLevelNumber();
- d_hopscell = new math::HierarchyCellDataOpsReal<double>(d_hierarchy,
- d_ln_min,
- d_ln_max);
- d_hopsside = new math::HierarchySideDataOpsReal<double>(d_hierarchy,
- d_ln_min,
- d_ln_max);
-
-#ifdef DEBUG_CHECK_ASSERTIONS
-
- if (d_physical_bc_coef == NULL) {
- /*
- * It's an error not to have bc object set.
- * Note that the bc object cannot be passed in through
- * the argument because the interface is inherited.
- */
- TBOX_ERROR(
- d_object_name << ": No physical bc object in\n"
- <<
- "StokesFACOps::initializeOperatorState\n"
- << "You must use "
- <<
- "StokesFACOps::setPhysicalBcCoefObject\n"
- <<
- "to set one before calling initializeOperatorState\n");
- }
-
- if (solution.getNumberOfComponents() != 1) {
- TBOX_WARNING(d_object_name
- << ": Solution vector has multiple components.\n"
- << "Solver is for component 0 only.\n");
- }
- if (rhs.getNumberOfComponents() != 1) {
- TBOX_WARNING(d_object_name
- << ": RHS vector has multiple components.\n"
- << "Solver is for component 0 only.\n");
- }
-
- /*
- * Make sure that solution and rhs data
- * are of correct type
- * are allocated
- * has sufficient ghost width
- */
- tbox::Pointer<hier::Variable> var;
- {
- vdb->mapIndexToVariable(rhs.getComponentDescriptorIndex(0),
- var);
- if (!var) {
- TBOX_ERROR(d_object_name << ": RHS component does not\n"
- << "correspond to a variable.\n");
- }
- tbox::Pointer<pdat::CellVariable<double> > cell_var = var;
- if (!cell_var) {
- TBOX_ERROR(d_object_name
- << ": RHS variable is not cell-centered double\n");
- }
- }
- {
- vdb->mapIndexToVariable(solution.getComponentDescriptorIndex(0),
- var);
- if (!var) {
- TBOX_ERROR(d_object_name << ": Solution component does not\n"
- << "correspond to a variable.\n");
- }
- tbox::Pointer<pdat::CellVariable<double> > cell_var = var;
- if (!cell_var) {
- TBOX_ERROR(d_object_name
- << ": Solution variable is not cell-centered double\n");
- }
- }
- for (ln = d_ln_min; ln <= d_ln_max; ++ln) {
- tbox::Pointer<hier::PatchLevel> level_ptr =
- d_hierarchy->getPatchLevel(ln);
- hier::PatchLevel& level = *level_ptr;
- for (hier::PatchLevel::Iterator pi(level); pi; pi++) {
- hier::Patch& patch = **pi;
- tbox::Pointer<hier::PatchData> fd =
- patch.getPatchData(rhs.getComponentDescriptorIndex(0));
- if (fd) {
- /*
- * Some data checks can only be done if the data already exists.
- */
- tbox::Pointer<pdat::CellData<double> > cd = fd;
- if (!cd) {
- TBOX_ERROR(d_object_name
- << ": RHS data is not cell-centered double\n");
- }
- if (cd->getDepth() > 1) {
- TBOX_WARNING(d_object_name
- << ": RHS data has multiple depths.\n"
- << "Solver is for depth 0 only.\n");
- }
- }
- tbox::Pointer<hier::PatchData> ud =
- patch.getPatchData(solution.getComponentDescriptorIndex(0));
- if (ud) {
- /*
- * Some data checks can only be done if the data already exists.
- */
- tbox::Pointer<pdat::CellData<double> > cd = ud;
- if (!cd) {
- TBOX_ERROR(d_object_name
- << ": Solution data is not cell-centered double\n");
- }
- if (cd->getDepth() > 1) {
- TBOX_WARNING(d_object_name
- << ": Solution data has multiple depths.\n"
- << "Solver is for depth 0 only.\n");
- }
- if (cd->getGhostCellWidth() < hier::IntVector::getOne(d_dim)) {
- TBOX_ERROR(d_object_name
- << ": Solution data has insufficient ghost width\n");
- }
- }
- }
- }
-
- /*
- * Solution and rhs must have some similar properties.
- */
- if (rhs.getPatchHierarchy() != d_hierarchy
- || rhs.getCoarsestLevelNumber() != d_ln_min
- || rhs.getFinestLevelNumber() != d_ln_max) {
- TBOX_ERROR(d_object_name << ": solution and rhs do not have\n"
- << "the same set of patch levels.\n");
- }
-
-#endif
-
- /*
- * Initialize the coarse-fine boundary description for the
- * hierarchy.
- */
- d_cf_boundary.resizeArray(d_hierarchy->getNumberOfLevels());
-
- hier::IntVector max_gcw(d_dim, 1);
- for (ln = d_ln_min; ln <= d_ln_max; ++ln) {
- d_cf_boundary[ln] = new hier::CoarseFineBoundary(*d_hierarchy,
- ln,
- max_gcw);
- }
-#ifdef HAVE_HYPRE
- if (d_coarse_solver_choice == "hypre") {
- d_hypre_solver.initializeSolverState(d_hierarchy, d_ln_min);
- /*
- * Share the boundary condition object with the hypre solver
- * to make sure that boundary condition settings are consistent
- * between the two objects.
- */
- d_hypre_solver.setPhysicalBcCoefObject(d_physical_bc_coef);
- d_hypre_solver.setMatrixCoefficients(d_stokes_spec);
- }
-#endif
-
- /*
- * Get the transfer operators.
- * Flux coarsening is conservative.
- * Cell (solution, error, etc) coarsening is conservative.
- * Cell refinement from same level is constant refinement.
- * Cell refinement from coarser level is chosen by the
- * choice of coarse-fine discretization, d_cf_discretization,
- * which should be set to either "Ewing" or one of the
- * acceptable strings for looking up the refine operator.
- */
- tbox::Pointer<geom::CartesianGridGeometry> geometry =
- d_hierarchy->getGridGeometry();
- tbox::Pointer<hier::Variable> variable;
-
- vdb->mapIndexToVariable(d_cell_scratch_id, variable);
- d_prolongation_refine_operator =
- geometry->lookupRefineOperator(variable,
- d_prolongation_method);
-
- vdb->mapIndexToVariable(d_cell_scratch_id, variable);
- d_urestriction_coarsen_operator =
- d_rrestriction_coarsen_operator =
- geometry->lookupCoarsenOperator(variable,
- "CONSERVATIVE_COARSEN");
-
- vdb->mapIndexToVariable(d_oflux_scratch_id, variable);
- d_flux_coarsen_operator =
- geometry->lookupCoarsenOperator(variable,
- "CONSERVATIVE_COARSEN");
-
- vdb->mapIndexToVariable(d_cell_scratch_id, variable);
- d_ghostfill_refine_operator =
- geometry->lookupRefineOperator(variable,
- d_cf_discretization == "Ewing" ?
- "CONSTANT_REFINE" : d_cf_discretization);
-
- vdb->mapIndexToVariable(d_cell_scratch_id, variable);
- d_ghostfill_nocoarse_refine_operator =
- geometry->lookupRefineOperator(variable,
- "CONSTANT_REFINE");
-
-#ifdef DEBUG_CHECK_ASSERTIONS
- if (!d_prolongation_refine_operator) {
- TBOX_ERROR(d_object_name
- << ": Cannot find prolongation refine operator");
- }
- if (!d_urestriction_coarsen_operator) {
- TBOX_ERROR(d_object_name
- << ": Cannot find restriction coarsening operator");
- }
- if (!d_rrestriction_coarsen_operator) {
- TBOX_ERROR(d_object_name
- << ": Cannot find restriction coarsening operator");
- }
- if (!d_flux_coarsen_operator) {
- TBOX_ERROR(d_object_name
- << ": Cannot find flux coarsening operator");
- }
- if (!d_ghostfill_refine_operator) {
- TBOX_ERROR(d_object_name
- << ": Cannot find ghost filling refinement operator");
- }
- if (!d_ghostfill_nocoarse_refine_operator) {
- TBOX_ERROR(d_object_name
- << ": Cannot find ghost filling refinement operator");
- }
-#endif
-
- for (ln = d_ln_min + 1; ln <= d_ln_max; ++ln) {
- d_hierarchy->getPatchLevel(ln)->
- allocatePatchData(d_oflux_scratch_id);
- }
-
- /*
- * Make space for saving communication schedules.
- * There is no need to delete the old schedules first
- * because we have deallocated the solver state above.
- */
- d_prolongation_refine_schedules.resizeArray(d_ln_max + 1);
- d_ghostfill_refine_schedules.resizeArray(d_ln_max + 1);
- d_ghostfill_nocoarse_refine_schedules.resizeArray(d_ln_max + 1);
- d_urestriction_coarsen_schedules.resizeArray(d_ln_max + 1);
- d_rrestriction_coarsen_schedules.resizeArray(d_ln_max + 1);
- d_flux_coarsen_schedules.resizeArray(d_ln_max + 1);
-
- d_prolongation_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
- d_urestriction_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
- d_rrestriction_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
- d_flux_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
- d_ghostfill_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
- d_ghostfill_nocoarse_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
-
- d_prolongation_refine_algorithm->
- registerRefine(d_cell_scratch_id,
- solution.getComponentDescriptorIndex(0),
- d_cell_scratch_id,
- d_prolongation_refine_operator);
- d_urestriction_coarsen_algorithm->
- registerCoarsen(solution.getComponentDescriptorIndex(0),
- solution.getComponentDescriptorIndex(0),
- d_urestriction_coarsen_operator);
- d_rrestriction_coarsen_algorithm->
- registerCoarsen(rhs.getComponentDescriptorIndex(0),
- rhs.getComponentDescriptorIndex(0),
- d_rrestriction_coarsen_operator);
- d_ghostfill_refine_algorithm->
- registerRefine(solution.getComponentDescriptorIndex(0),
- solution.getComponentDescriptorIndex(0),
- solution.getComponentDescriptorIndex(0),
- d_ghostfill_refine_operator);
- d_flux_coarsen_algorithm->
- registerCoarsen(((d_flux_id != -1) ? d_flux_id : d_flux_scratch_id),
- d_oflux_scratch_id,
- d_flux_coarsen_operator);
- d_ghostfill_nocoarse_refine_algorithm->
- registerRefine(solution.getComponentDescriptorIndex(0),
- solution.getComponentDescriptorIndex(0),
- solution.getComponentDescriptorIndex(0),
- d_ghostfill_nocoarse_refine_operator);
-
- for (int dest_ln = d_ln_min + 1; dest_ln <= d_ln_max; ++dest_ln) {
-
- tbox::Pointer<xfer::PatchLevelFullFillPattern> fill_pattern(
- new xfer::PatchLevelFullFillPattern());
- d_prolongation_refine_schedules[dest_ln] =
- d_prolongation_refine_algorithm->
- createSchedule(fill_pattern,
- d_hierarchy->getPatchLevel(dest_ln),
- tbox::Pointer<hier::PatchLevel>(),
- dest_ln - 1,
- d_hierarchy,
- &d_bc_helper);
- if (!d_prolongation_refine_schedules[dest_ln]) {
- TBOX_ERROR(d_object_name
- << ": Cannot create a refine schedule for prolongation!\n");
- }
- d_ghostfill_refine_schedules[dest_ln] =
- d_ghostfill_refine_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- dest_ln - 1,
- d_hierarchy,
- &d_bc_helper);
- if (!d_ghostfill_refine_schedules[dest_ln]) {
- TBOX_ERROR(d_object_name
- << ": Cannot create a refine schedule for ghost filling!\n");
- }
- d_ghostfill_nocoarse_refine_schedules[dest_ln] =
- d_ghostfill_nocoarse_refine_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- &d_bc_helper);
- if (!d_ghostfill_nocoarse_refine_schedules[dest_ln]) {
- TBOX_ERROR(
- d_object_name
- <<
- ": Cannot create a refine schedule for ghost filling on bottom level!\n");
- }
- }
- for (int dest_ln = d_ln_min; dest_ln < d_ln_max; ++dest_ln) {
- d_urestriction_coarsen_schedules[dest_ln] =
- d_urestriction_coarsen_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- d_hierarchy->getPatchLevel(dest_ln + 1));
- if (!d_urestriction_coarsen_schedules[dest_ln]) {
- TBOX_ERROR(d_object_name
- << ": Cannot create a coarsen schedule for U restriction!\n");
- }
- d_rrestriction_coarsen_schedules[dest_ln] =
- d_rrestriction_coarsen_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- d_hierarchy->getPatchLevel(dest_ln + 1));
- if (!d_rrestriction_coarsen_schedules[dest_ln]) {
- TBOX_ERROR(d_object_name
- << ": Cannot create a coarsen schedule for R restriction!\n");
- }
- d_flux_coarsen_schedules[dest_ln] =
- d_flux_coarsen_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- d_hierarchy->getPatchLevel(dest_ln + 1));
- if (!d_flux_coarsen_schedules[dest_ln]) {
- TBOX_ERROR(d_object_name
- << ": Cannot create a coarsen schedule for flux transfer!\n");
- }
- }
- d_ghostfill_nocoarse_refine_schedules[d_ln_min] =
- d_ghostfill_nocoarse_refine_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(d_ln_min),
- &d_bc_helper);
- if (!d_ghostfill_nocoarse_refine_schedules[d_ln_min]) {
- TBOX_ERROR(
- d_object_name
- <<
- ": Cannot create a refine schedule for ghost filling on bottom level!\n");
- }
-}
-
-/*
- ********************************************************************
- * FACOperatorStrategy virtual deallocateOperatorState *
- * function. Deallocate internal hierarchy-dependent data. *
- * State is allocated iff hierarchy is set. *
- ********************************************************************
- */
-
-void StokesFACOps::deallocateOperatorState()
-{
- if (d_hierarchy) {
- int ln;
- for (ln = d_ln_min + 1; ln <= d_ln_max; ++ln) {
- d_hierarchy->getPatchLevel(ln)->
- deallocatePatchData(d_oflux_scratch_id);
- }
- d_cf_boundary.resizeArray(0);
-#ifdef HAVE_HYPRE
- d_hypre_solver.deallocateSolverState();
-#endif
- d_hierarchy.setNull();
- d_ln_min = -1;
- d_ln_max = -1;
-
- d_prolongation_refine_algorithm.setNull();
- d_prolongation_refine_schedules.setNull();
-
- d_urestriction_coarsen_algorithm.setNull();
- d_urestriction_coarsen_schedules.setNull();
-
- d_rrestriction_coarsen_algorithm.setNull();
- d_rrestriction_coarsen_schedules.setNull();
-
- d_flux_coarsen_algorithm.setNull();
- d_flux_coarsen_schedules.setNull();
-
- d_ghostfill_refine_algorithm.setNull();
- d_ghostfill_refine_schedules.setNull();
-
- d_ghostfill_nocoarse_refine_algorithm.setNull();
- d_ghostfill_nocoarse_refine_schedules.setNull();
-
- }
-}
-
-/*
- ********************************************************************
- * FACOperatorStrategy virtual postprocessOneCycle function. *
- ********************************************************************
- */
-
-void StokesFACOps::postprocessOneCycle(
- int fac_cycle_num,
- const SAMRAIVectorReal<double>& current_soln,
- const SAMRAIVectorReal<double>& residual)
-{
- NULL_USE(current_soln);
- NULL_USE(residual);
-
- if (d_enable_logging) {
- if (d_preconditioner) {
- /*
- * Output convergence progress. This is probably only appropriate
- * if the solver is NOT being used as a preconditioner.
- */
- double avg_factor, final_factor;
- d_preconditioner->getConvergenceFactors(avg_factor, final_factor);
- tbox::plog
- << "iter=" << std::setw(4) << fac_cycle_num
- << " resid=" << d_preconditioner->getResidualNorm()
- << " net conv=" << d_preconditioner->getNetConvergenceFactor()
- << " final conv=" << d_preconditioner->getNetConvergenceFactor()
- << " avg conv=" << d_preconditioner->getAvgConvergenceFactor()
- << std::endl;
- }
- }
-}
-
-/*
- ********************************************************************
- * FACOperatorStrategy virtual restrictSolution function. *
- * After restricting solution, update ghost cells of the affected *
- * level. *
- ********************************************************************
- */
-
-void StokesFACOps::restrictSolution(
- const SAMRAIVectorReal<double>& s,
- SAMRAIVectorReal<double>& d,
- int dest_ln) {
-
- t_restrict_solution->start();
-
- xeqScheduleURestriction(d.getComponentDescriptorIndex(0),
- s.getComponentDescriptorIndex(0),
- dest_ln);
-
- d_bc_helper.setHomogeneousBc(false);
- d_bc_helper.setTargetDataId(d.getComponentDescriptorIndex(0));
-
- if (dest_ln == d_ln_min) {
- xeqScheduleGhostFillNoCoarse(d.getComponentDescriptorIndex(0),
- dest_ln);
- } else {
- xeqScheduleGhostFill(d.getComponentDescriptorIndex(0),
- dest_ln);
- }
-
- t_restrict_solution->stop();
-}
-
-/*
- ********************************************************************
- * FACOperatorStrategy virtual restrictresidual function. *
- ********************************************************************
- */
-
-void StokesFACOps::restrictResidual(
- const SAMRAIVectorReal<double>& s,
- SAMRAIVectorReal<double>& d,
- int dest_ln) {
-
- t_restrict_residual->start();
-
- xeqScheduleRRestriction(d.getComponentDescriptorIndex(0),
- s.getComponentDescriptorIndex(0),
- dest_ln);
-
- t_restrict_residual->stop();
-}
-
-/*
- ***********************************************************************
- * FACOperatorStrategy virtual prolongErrorAndCorrect function. *
- * After the prolongation, we set the physical boundary condition *
- * for the correction, which is zero. Other ghost cell values, *
- * which are preset to zero, need not be set. *
- ***********************************************************************
- */
-
-void StokesFACOps::prolongErrorAndCorrect(
- const SAMRAIVectorReal<double>& s,
- SAMRAIVectorReal<double>& d,
- int dest_ln) {
-
- t_prolong->start();
-
-#ifdef DEBUG_CHECK_ASSERTIONS
- if (s.getPatchHierarchy() != d_hierarchy
- || d.getPatchHierarchy() != d_hierarchy) {
- TBOX_ERROR(d_object_name << ": Vector hierarchy does not match\n"
- "internal state hierarchy.");
- }
-#endif
-
- tbox::Pointer<hier::PatchLevel> coarse_level =
- d_hierarchy->getPatchLevel(dest_ln - 1);
- tbox::Pointer<hier::PatchLevel> fine_level =
- d_hierarchy->getPatchLevel(dest_ln);
-
- /*
- * Data is prolonged into the scratch space corresponding
- * to index d_cell_scratch_id and allocated here.
- */
- fine_level->allocatePatchData(d_cell_scratch_id);
-
- /*
- * Refine solution into scratch space to fill the fine level
- * interior in the scratch space, then use that refined data
- * to correct the fine level error.
- */
- d_bc_helper.setTargetDataId(d_cell_scratch_id);
- d_bc_helper.setHomogeneousBc(true);
- const int src_index = s.getComponentDescriptorIndex(0);
- xeqScheduleProlongation(d_cell_scratch_id,
- src_index,
- d_cell_scratch_id,
- dest_ln);
-
- /*
- * Add the refined error in the scratch space
- * to the error currently residing in the destination level.
- */
- math::HierarchyCellDataOpsReal<double>
- hierarchy_math_ops(d_hierarchy, dest_ln, dest_ln);
- const int dst_index = d.getComponentDescriptorIndex(0);
- hierarchy_math_ops.add(dst_index, dst_index, d_cell_scratch_id);
-
- fine_level->deallocatePatchData(d_cell_scratch_id);
-
- t_prolong->stop();
-
-}
-
-/*
- ********************************************************************
- ********************************************************************
- */
-
-void StokesFACOps::smoothError(
- SAMRAIVectorReal<double>& data,
- const SAMRAIVectorReal<double>& residual,
- int ln,
- int num_sweeps)
-{
-
- t_smooth_error->start();
-
- checkInputPatchDataIndices();
- if (d_smoothing_choice == "redblack") {
- smoothErrorByRedBlack(data,
- residual,
- ln,
- num_sweeps,
- d_residual_tolerance_during_smoothing);
- } else {
- TBOX_ERROR(d_object_name << ": Bad smoothing choice '"
- << d_smoothing_choice
- << "' in StokesFACOps.");
- }
-
- t_smooth_error->stop();
-}
-
-/*
- ********************************************************************
- * Workhorse function to smooth error using red-black *
- * Gauss-Seidel iterations. *
- ********************************************************************
- */
-
-void StokesFACOps::smoothErrorByRedBlack(
- SAMRAIVectorReal<double>& data,
- const SAMRAIVectorReal<double>& residual,
- int ln,
- int num_sweeps,
- double residual_tolerance)
-{
-
- checkInputPatchDataIndices();
-
-#ifdef DEBUG_CHECK_ASSERTIONS
- if (data.getPatchHierarchy() != d_hierarchy
- || residual.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);
-
- const int data_id = data.getComponentDescriptorIndex(0);
-
- const int flux_id = (d_flux_id != -1) ? d_flux_id : d_flux_scratch_id;
-
- d_bc_helper.setTargetDataId(data_id);
- d_bc_helper.setHomogeneousBc(true);
- xeqScheduleGhostFillNoCoarse(data_id, ln);
-
- if (ln > d_ln_min) {
- /*
- * Perform a one-time transfer of data from coarser level,
- * to fill ghost boundaries that will not change through
- * the smoothing loop.
- */
- xeqScheduleGhostFill(data_id, ln);
- }
-
- /*
- * Smooth the number of sweeps specified or until
- * the convergence is satisfactory.
- */
- int isweep;
- double red_maxres, blk_maxres, maxres = 0;
- red_maxres = blk_maxres = residual_tolerance + 1;
- /*
- * Instead of checking residual convergence globally,
- * we check the not_converged flag. This avoids possible
- * round-off errors affecting different processes differently,
- * leading to disagreement on whether to continue smoothing.
- */
- int not_converged = 1;
- for (isweep = 0; isweep < num_sweeps && not_converged; ++isweep) {
- red_maxres = blk_maxres = 0;
-
- // Red sweep.
- xeqScheduleGhostFillNoCoarse(data_id, ln);
- for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
- tbox::Pointer<hier::Patch> patch = *pi;
-
- bool deallocate_flux_data_when_done = false;
- if (flux_id == d_flux_scratch_id) {
- /*
- * Using internal temporary storage for flux.
- * For each patch, make sure the internal
- * side-centered data is allocated and note
- * whether that data should be deallocated when done.
- */
- if (!patch->checkAllocated(flux_id)) {
- patch->allocatePatchData(flux_id);
- deallocate_flux_data_when_done = true;
- }
- }
-
- tbox::Pointer<pdat::CellData<double> >
- scalar_field_data = d_stokes_spec.cIsVariable() ?
- patch->getPatchData(d_stokes_spec.getCPatchDataId()) :
- tbox::Pointer<hier::PatchData>(NULL);
- tbox::Pointer<pdat::CellData<double> >
- err_data = data.getComponentPatchData(0, *patch);
- tbox::Pointer<pdat::CellData<double> >
- residual_data = residual.getComponentPatchData(0, *patch);
- tbox::Pointer<pdat::SideData<double> >
- flux_data = patch->getPatchData(flux_id);
-
- computeFluxOnPatch(
- *patch,
- level->getRatioToCoarserLevel(),
- *err_data,
- *flux_data);
-
- redOrBlackSmoothingOnPatch(*patch,
- *flux_data,
- *residual_data,
- *err_data,
- 'r',
- &red_maxres);
-
- if (deallocate_flux_data_when_done) {
- patch->deallocatePatchData(flux_id);
- }
- } // End patch number *pi
- xeqScheduleGhostFillNoCoarse(data_id, ln);
-
- // Black sweep.
- for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
- tbox::Pointer<hier::Patch> patch = *pi;
-
- bool deallocate_flux_data_when_done = false;
- if (flux_id == d_flux_scratch_id) {
- /*
- * Using internal temporary storage for flux.
- * For each patch, make sure the internal
- * side-centered data is allocated and note
- * whether that data should be deallocated when done.
- */
- if (!patch->checkAllocated(flux_id)) {
- patch->allocatePatchData(flux_id);
- deallocate_flux_data_when_done = true;
- }
- }
-
- tbox::Pointer<pdat::CellData<double> >
- scalar_field_data = d_stokes_spec.cIsVariable() ?
- patch->getPatchData(d_stokes_spec.getCPatchDataId()) :
- tbox::Pointer<hier::PatchData>(NULL);
- tbox::Pointer<pdat::CellData<double> >
- err_data = data.getComponentPatchData(0, *patch);
- tbox::Pointer<pdat::CellData<double> >
- residual_data = residual.getComponentPatchData(0, *patch);
- tbox::Pointer<pdat::SideData<double> >
- flux_data = patch->getPatchData(flux_id);
-
- computeFluxOnPatch(
- *patch,
- level->getRatioToCoarserLevel(),
- *err_data,
- *flux_data);
-
- redOrBlackSmoothingOnPatch(*patch,
- *flux_data,
- *residual_data,
- *err_data,
- 'b',
- &blk_maxres);
-
- if (deallocate_flux_data_when_done) {
- patch->deallocatePatchData(flux_id);
- }
- } // End patch number *pi
- xeqScheduleGhostFillNoCoarse(data_id, ln);
- if (residual_tolerance >= 0.0) {
- /*
- * Check for early end of sweeps due to convergence
- * only if it is numerically possible (user gave a
- * non negative value for residual tolerance).
- */
- maxres = tbox::MathUtilities<double>::Max(red_maxres, blk_maxres);
- not_converged = maxres > residual_tolerance;
- const tbox::SAMRAI_MPI& mpi(d_hierarchy->getDomainMappedBoxLevel().getMPI());
- if (mpi.getSize() > 1) {
- mpi.AllReduce(¬_converged, 1, MPI_MAX);
- }
- }
- } // End sweep number isweep
- if (d_enable_logging) tbox::plog
- << d_object_name << " RBGS smoothing maxres = " << maxres << "\n"
- << " after " << isweep << " sweeps.\n";
-
-}
-
-/*
- ********************************************************************
- * Fix flux on coarse-fine boundaries computed from a *
- * constant-refine interpolation of coarse level data. *
- ********************************************************************
- */
-
-void StokesFACOps::ewingFixFlux(
- const hier::Patch& patch,
- const pdat::CellData<double>& soln_data,
- pdat::SideData<double>& flux_data,
- const hier::IntVector& ratio_to_coarser) const
-{
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(d_dim, patch, soln_data, flux_data,
- ratio_to_coarser);
-
- const int patch_ln = patch.getPatchLevelNumber();
- const hier::GlobalId id = patch.getGlobalId();
- tbox::Pointer<geom::CartesianPatchGeometry> patch_geom =
- patch.getPatchGeometry();
- const double* dx = patch_geom->getDx();
- const hier::Box& patch_box(patch.getBox());
- const hier::Index& plower = patch_box.lower();
- const hier::Index& pupper = patch_box.upper();
-
- const tbox::Array<hier::BoundaryBox>& bboxes =
- d_cf_boundary[patch_ln]->getBoundaries(id, 1);
- int bn, nboxes = bboxes.getSize();
-
- if (d_stokes_spec.dIsVariable()) {
-
- tbox::Pointer<pdat::SideData<double> > diffcoef_data;
- diffcoef_data = patch.getPatchData(d_stokes_spec.getDPatchDataId());
-
- for (bn = 0; bn < nboxes; ++bn) {
- const hier::BoundaryBox& boundary_box = bboxes[bn];
-
- TBOX_ASSERT(boundary_box.getBoundaryType() == 1);
-
- const hier::Box& bdry_box = boundary_box.getBox();
- const hier::Index& blower = bdry_box.lower();
- const hier::Index& bupper = bdry_box.upper();
- const int location_index = boundary_box.getLocationIndex();
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(ewingfixfluxvardc2d, EWINGFIXFLUXVARDC2D) (
- flux_data.getPointer(0), flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- diffcoef_data->getPointer(0), diffcoef_data->getPointer(1),
- &diffcoef_data->getGhostCellWidth()[0],
- &diffcoef_data->getGhostCellWidth()[1],
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &plower[0], &pupper[0], &plower[1], &pupper[1],
- &location_index,
- &ratio_to_coarser[0],
- &blower[0], &bupper[0],
- dx);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(ewingfixfluxvardc3d, EWINGFIXFLUXVARDC3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- diffcoef_data->getPointer(0),
- diffcoef_data->getPointer(1),
- diffcoef_data->getPointer(2),
- &diffcoef_data->getGhostCellWidth()[0],
- &diffcoef_data->getGhostCellWidth()[1],
- &diffcoef_data->getGhostCellWidth()[2],
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &plower[0], &pupper[0],
- &plower[1], &pupper[1],
- &plower[2], &pupper[2],
- &location_index,
- &ratio_to_coarser[0],
- &blower[0], &bupper[0],
- dx);
- } else {
- TBOX_ERROR("StokesFACOps : DIM > 3 not supported" << std::endl);
- }
-
- }
- } else {
-
- const double diffcoef_constant = d_stokes_spec.getDConstant();
-
- for (bn = 0; bn < nboxes; ++bn) {
- const hier::BoundaryBox& boundary_box = bboxes[bn];
-
- TBOX_ASSERT(boundary_box.getBoundaryType() == 1);
-
- const hier::Box& bdry_box = boundary_box.getBox();
- const hier::Index& blower = bdry_box.lower();
- const hier::Index& bupper = bdry_box.upper();
- const int location_index = boundary_box.getLocationIndex();
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(ewingfixfluxcondc2d, EWINGFIXFLUXCONDC2D) (
- flux_data.getPointer(0), flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- diffcoef_constant,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &plower[0], &pupper[0],
- &plower[1], &pupper[1],
- &location_index,
- &ratio_to_coarser[0],
- &blower[0], &bupper[0],
- dx);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(ewingfixfluxcondc3d, EWINGFIXFLUXCONDC3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- diffcoef_constant,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &plower[0], &pupper[0],
- &plower[1], &pupper[1],
- &plower[2], &pupper[2],
- &location_index,
- &ratio_to_coarser[0],
- &blower[0], &bupper[0],
- dx);
- }
- }
- }
-}
-
-/*
- ********************************************************************
- * FACOperatorStrategy virtual solveCoarsestLevel *
- * function *
- ********************************************************************
- */
-
-int StokesFACOps::solveCoarsestLevel(
- SAMRAIVectorReal<double>& data,
- const SAMRAIVectorReal<double>& residual,
- int coarsest_ln) {
-
- t_solve_coarsest->start();
-
- checkInputPatchDataIndices();
-
- int return_value = 0;
-
- if (d_coarse_solver_choice == "jacobi") {
- d_residual_tolerance_during_smoothing = d_coarse_solver_tolerance;
- smoothError(data,
- residual,
- coarsest_ln,
- d_coarse_solver_max_iterations);
- d_residual_tolerance_during_smoothing = -1.0;
- } else if (d_coarse_solver_choice == "redblack") {
- d_residual_tolerance_during_smoothing = d_coarse_solver_tolerance;
- smoothError(data,
- residual,
- coarsest_ln,
- d_coarse_solver_max_iterations);
- d_residual_tolerance_during_smoothing = -1.0;
- } else if (d_coarse_solver_choice == "hypre") {
-#ifndef HAVE_HYPRE
- TBOX_ERROR(d_object_name << ": Coarse level solver choice '"
- << d_coarse_solver_choice
- << "' unavailable in "
- << "scapStokesOps::solveCoarsestLevel.");
-#else
- return_value = solveCoarsestLevel_HYPRE(data, residual, coarsest_ln);
-#endif
- } else {
- TBOX_ERROR(
- d_object_name << ": Bad coarse level solver choice '"
- << d_coarse_solver_choice
- <<
- "' in scapStokesOps::solveCoarsestLevel.");
- }
-
- xeqScheduleGhostFillNoCoarse(data.getComponentDescriptorIndex(0),
- coarsest_ln);
-
- t_solve_coarsest->stop();
-
- return return_value;
-}
-
-#ifdef HAVE_HYPRE
-/*
- ********************************************************************
- * Solve coarsest level using Hypre *
- * We only solve for the error, so we always use homogeneous bc. *
- ********************************************************************
- */
-
-int StokesFACOps::solveCoarsestLevel_HYPRE(
- SAMRAIVectorReal<double>& data,
- const SAMRAIVectorReal<double>& residual,
- int coarsest_ln) {
-
- NULL_USE(coarsest_ln);
-
-#ifndef HAVE_HYPRE
- TBOX_ERROR(d_object_name << ": Coarse level solver choice '"
- << d_coarse_solver_choice
- << "' unavailable in "
- << "StokesFACOps::solveCoarsestLevel.");
-
- return 0;
-
-#else
-
- checkInputPatchDataIndices();
- d_hypre_solver.setStoppingCriteria(d_coarse_solver_max_iterations,
- d_coarse_solver_tolerance);
- const int solver_ret =
- d_hypre_solver.solveSystem(
- data.getComponentDescriptorIndex(0),
- residual.getComponentDescriptorIndex(0),
- true);
- /*
- * Present data on the solve.
- * The Hypre solver returns 0 if converged.
- */
- if (d_enable_logging) tbox::plog
- << d_object_name << " Hypre solve " << (solver_ret ? "" : "NOT ")
- << "converged\n"
- << "\titerations: " << d_hypre_solver.getNumberOfIterations() << "\n"
- << "\tresidual: " << d_hypre_solver.getRelativeResidualNorm() << "\n";
-
- return !solver_ret;
-
-#endif
-
-}
-#endif
-
-/*
- ********************************************************************
- * FACOperatorStrategy virtual *
- * computeCompositeResidualOnLevel function *
- ********************************************************************
- */
-
-void 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 soln_id = solution.getComponentDescriptorIndex(0);
- d_bc_helper.setTargetDataId(soln_id);
- d_bc_helper.setHomogeneousBc(error_equation_indicator);
-
- const int flux_id = (d_flux_id != -1) ? d_flux_id : d_flux_scratch_id;
-
- /*
- * 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. */
- {
- tbox::Pointer<xfer::RefineSchedule> ln_refine_schedule;
- if (ln > d_ln_min) {
- /* Fill from current, next coarser level and physical boundary */
- xeqScheduleGhostFill(soln_id, ln);
- } else {
- /* Fill from current and physical boundary */
- xeqScheduleGhostFillNoCoarse(soln_id, ln);
- }
- }
-
- /*
- * For the whole level, make sure the internal
- * side-centered data is allocated and note
- * whether that data should be deallocated when done.
- * We do this for the whole level because the data
- * undergoes transfer operations which require the
- * whole level data.
- */
- bool deallocate_flux_data_when_done = false;
- if (flux_id == d_flux_scratch_id) {
- if (!level->checkAllocated(flux_id)) {
- level->allocatePatchData(flux_id);
- deallocate_flux_data_when_done = true;
- }
- }
-
- /*
- * S2. Compute flux on patches in level.
- */
- for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
- tbox::Pointer<hier::Patch> patch = *pi;
-
- tbox::Pointer<pdat::CellData<double> >
- soln_data = solution.getComponentPatchData(0, *patch);
- tbox::Pointer<pdat::CellData<double> >
- rhs_data = rhs.getComponentPatchData(0, *patch);
- tbox::Pointer<pdat::CellData<double> >
- residual_data = residual.getComponentPatchData(0, *patch);
- tbox::Pointer<pdat::SideData<double> >
- flux_data = patch->getPatchData(flux_id);
- computeFluxOnPatch(
- *patch,
- level->getRatioToCoarserLevel(),
- *soln_data,
- *flux_data);
-
- }
-
- /*
- * S3. Coarsen oflux data from next finer level so that
- * the computed flux becomes the composite grid flux.
- */
- if (ln < d_ln_max) {
- xeqScheduleFluxCoarsen(flux_id, d_oflux_scratch_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> >
- soln_data = solution.getComponentPatchData(0, *patch);
- tbox::Pointer<pdat::CellData<double> >
- rhs_data = rhs.getComponentPatchData(0, *patch);
- tbox::Pointer<pdat::CellData<double> >
- residual_data = residual.getComponentPatchData(0, *patch);
- tbox::Pointer<pdat::SideData<double> >
- flux_data = patch->getPatchData(flux_id);
- computeResidualOnPatch(*patch,
- *flux_data,
- *soln_data,
- *rhs_data,
- *residual_data);
-
- if (ln > d_ln_min) {
- /*
- * Save outerflux data so that next coarser level
- * can compute its coarse-fine composite flux.
- * This is not strictly needed in this "compute residual"
- * loop through the patches, but we put it here to
- * avoid writing another loop for it.
- */
- tbox::Pointer<pdat::OutersideData<double> >
- oflux_data = patch->getPatchData(d_oflux_scratch_id);
-
- TBOX_ASSERT(oflux_data);
-
- oflux_data->copy(*flux_data);
- }
- }
-
- if (deallocate_flux_data_when_done) {
- level->deallocatePatchData(flux_id);
- }
-
- t_compute_composite_residual->stop();
-}
-
-/*
- ********************************************************************
- * FACOperatorStrategy virtual computeResidualNorm *
- * function *
- ********************************************************************
- */
-
-double StokesFACOps::computeResidualNorm(
- const SAMRAIVectorReal<double>& residual,
- int fine_ln,
- int coarse_ln)
-{
-
- if (coarse_ln != residual.getCoarsestLevelNumber() ||
- fine_ln != residual.getFinestLevelNumber()) {
- TBOX_ERROR("StokesFACOps::computeResidualNorm() is not\n"
- << "set up to compute residual except on the range of\n"
- << "levels defining the vector.\n");
- }
- t_compute_residual_norm->start();
- /*
- * The residual vector was cloned from vectors that has
- * the proper weights associated with them, so we do not
- * have to explicitly weight the residuals.
- *
- * maxNorm: not good to use because Hypre's norm does not
- * correspond to it. Also maybe too sensitive to spikes.
- * L2Norm: maybe good. But does not correspond to the
- * scale of the quantity.
- * L1Norm: maybe good. Correspond to scale of quantity,
- * but may be too insensitive to spikes.
- * RMSNorm: maybe good.
- */
- double norm = residual.RMSNorm();
- t_compute_residual_norm->stop();
- return norm;
-}
-
-/*
- ********************************************************************
- * Compute the vector weight and put it at a specified patch data *
- * index. *
- ********************************************************************
- */
-
-void StokesFACOps::computeVectorWeights(
- tbox::Pointer<hier::PatchHierarchy> hierarchy,
- int weight_id,
- int coarsest_ln,
- int finest_ln) const
-{
- TBOX_ASSERT(!hierarchy.isNull());
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS1(d_dim, *hierarchy);
-
- if (coarsest_ln == -1) coarsest_ln = 0;
- if (finest_ln == -1) finest_ln = hierarchy->getFinestLevelNumber();
- if (finest_ln < coarsest_ln) {
- TBOX_ERROR(d_object_name
- << ": Illegal level number range. finest_ln < coarsest_ln.");
- }
-
- int ln;
- for (ln = finest_ln; ln >= coarsest_ln; --ln) {
-
- /*
- * On every level, first assign cell volume to vector weight.
- */
-
- tbox::Pointer<hier::PatchLevel> level =
- hierarchy->getPatchLevel(ln);
- for (hier::PatchLevel::Iterator p(level); p; p++) {
- tbox::Pointer<hier::Patch> patch = *p;
- tbox::Pointer<geom::CartesianPatchGeometry> patch_geometry =
- patch->getPatchGeometry();
- const double* dx = patch_geometry->getDx();
- double cell_vol = dx[0];
- if (d_dim > tbox::Dimension(1)) {
- cell_vol *= dx[1];
- }
-
- if (d_dim > tbox::Dimension(2)) {
- cell_vol *= dx[2];
- }
-
- tbox::Pointer<pdat::CellData<double> > w =
- patch->getPatchData(weight_id);
- if (!w) {
- TBOX_ERROR(d_object_name
- << ": weight id must refer to a pdat::CellVariable");
- }
- w->fillAll(cell_vol);
- }
-
- /*
- * On all but the finest level, assign 0 to vector
- * weight to cells covered by finer cells.
- */
-
- if (ln < finest_ln) {
-
- /*
- * First get the boxes that describe index space of the next finer
- * level and coarsen them to describe corresponding index space
- * at this level.
- */
-
- tbox::Pointer<hier::PatchLevel> next_finer_level =
- hierarchy->getPatchLevel(ln + 1);
- hier::BoxArray coarsened_boxes = next_finer_level->getBoxes();
- hier::IntVector coarsen_ratio(next_finer_level->getRatioToLevelZero());
- coarsen_ratio /= level->getRatioToLevelZero();
- coarsened_boxes.coarsen(coarsen_ratio);
-
- /*
- * Then set vector weight to 0 wherever there is
- * a nonempty intersection with the next finer level.
- * Note that all assignments are local.
- */
-
- for (hier::PatchLevel::Iterator p(level); p; p++) {
-
- tbox::Pointer<hier::Patch> patch = *p;
- for (int i = 0; i < coarsened_boxes.getNumberOfBoxes(); i++) {
-
- hier::Box coarse_box = coarsened_boxes[i];
- hier::Box intersection = coarse_box * (patch->getBox());
- if (!intersection.empty()) {
- tbox::Pointer<pdat::CellData<double> > w =
- patch->getPatchData(weight_id);
- w->fillAll(0.0, intersection);
-
- } // assignment only in non-empty intersection
- } // loop over coarsened boxes from finer level
- } // loop over patches in level
- } // all levels except finest
- } // loop over levels
-}
-
-/*
- ********************************************************************
- * Check the validity and correctness of input data for this class. *
- ********************************************************************
- */
-
-void StokesFACOps::checkInputPatchDataIndices() const {
- /*
- * Check input validity and correctness.
- */
- hier::VariableDatabase& vdb(*hier::VariableDatabase::getDatabase());
-
- if (!d_stokes_spec.dIsConstant()
- && d_stokes_spec.getDPatchDataId() != -1) {
- tbox::Pointer<hier::Variable> var;
- vdb.mapIndexToVariable(d_stokes_spec.getDPatchDataId(), var);
- tbox::Pointer<pdat::SideVariable<double> > diffcoef_var = var;
- if (!diffcoef_var) {
- TBOX_ERROR(d_object_name
- << ": Bad diffusion coefficient patch data index.");
- }
- }
-
- if (!d_stokes_spec.cIsConstant() && !d_stokes_spec.cIsZero()) {
- tbox::Pointer<hier::Variable> var;
- vdb.mapIndexToVariable(d_stokes_spec.getCPatchDataId(), var);
- tbox::Pointer<pdat::CellVariable<double> > scalar_field_var = var;
- if (!scalar_field_var) {
- TBOX_ERROR(d_object_name << ": Bad linear term patch data index.");
- }
- }
-
- if (d_flux_id != -1) {
- tbox::Pointer<hier::Variable> var;
- vdb.mapIndexToVariable(d_flux_id, var);
- tbox::Pointer<pdat::SideVariable<double> > flux_var = var;
-
- TBOX_ASSERT(flux_var);
- }
-
-}
-
-/*
- *******************************************************************
- * *
- * AMR-unaware patch-centered computational kernels. *
- * *
- *******************************************************************
- */
-
-void StokesFACOps::computeFluxOnPatch(
- const hier::Patch& patch,
- const hier::IntVector& ratio_to_coarser_level,
- const pdat::CellData<double>& w_data,
- pdat::SideData<double>& Dgradw_data) const
-{
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(d_dim, patch, ratio_to_coarser_level, w_data,
- Dgradw_data);
- TBOX_ASSERT(patch.inHierarchy());
- TBOX_ASSERT(w_data.getGhostCellWidth() >=
- hier::IntVector::getOne(ratio_to_coarser_level.getDim()));
-
- tbox::Pointer<geom::CartesianPatchGeometry> patch_geom =
- patch.getPatchGeometry();
- const hier::Box& box = patch.getBox();
- const int* lower = &box.lower()[0];
- const int* upper = &box.upper()[0];
- const double* dx = patch_geom->getDx();
-
- double D_value;
- tbox::Pointer<pdat::SideData<double> > D_data;
- if (d_stokes_spec.dIsConstant()) {
- D_value = d_stokes_spec.getDConstant();
- } else {
- D_data = patch.getPatchData(d_stokes_spec.getDPatchDataId());
- }
-
- if (d_stokes_spec.dIsConstant()) {
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(compfluxcondc2d, COMPFLUXCONDC2D) (
- Dgradw_data.getPointer(0),
- Dgradw_data.getPointer(1),
- &Dgradw_data.getGhostCellWidth()[0],
- &Dgradw_data.getGhostCellWidth()[1],
- D_value,
- w_data.getPointer(),
- &w_data.getGhostCellWidth()[0],
- &w_data.getGhostCellWidth()[1],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- dx);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(compfluxcondc3d, COMPFLUXCONDC3D) (
- Dgradw_data.getPointer(0),
- Dgradw_data.getPointer(1),
- Dgradw_data.getPointer(2),
- &Dgradw_data.getGhostCellWidth()[0],
- &Dgradw_data.getGhostCellWidth()[1],
- &Dgradw_data.getGhostCellWidth()[2],
- D_value,
- w_data.getPointer(),
- &w_data.getGhostCellWidth()[0],
- &w_data.getGhostCellWidth()[1],
- &w_data.getGhostCellWidth()[2],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- &lower[2], &upper[2],
- dx);
- }
- } else {
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(compfluxvardc2d, COMPFLUXVARDC2D) (
- Dgradw_data.getPointer(0),
- Dgradw_data.getPointer(1),
- &Dgradw_data.getGhostCellWidth()[0],
- &Dgradw_data.getGhostCellWidth()[1],
- D_data->getPointer(0),
- D_data->getPointer(1),
- &D_data->getGhostCellWidth()[0],
- &D_data->getGhostCellWidth()[1],
- w_data.getPointer(),
- &w_data.getGhostCellWidth()[0],
- &w_data.getGhostCellWidth()[1],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- dx);
- }
- if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(compfluxvardc3d, COMPFLUXVARDC3D) (
- Dgradw_data.getPointer(0),
- Dgradw_data.getPointer(1),
- Dgradw_data.getPointer(2),
- &Dgradw_data.getGhostCellWidth()[0],
- &Dgradw_data.getGhostCellWidth()[1],
- &Dgradw_data.getGhostCellWidth()[2],
- D_data->getPointer(0),
- D_data->getPointer(1),
- D_data->getPointer(2),
- &D_data->getGhostCellWidth()[0],
- &D_data->getGhostCellWidth()[1],
- &D_data->getGhostCellWidth()[2],
- w_data.getPointer(),
- &w_data.getGhostCellWidth()[0],
- &w_data.getGhostCellWidth()[1],
- &w_data.getGhostCellWidth()[2],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- &lower[2], &upper[2],
- dx);
- }
- }
-
- const int patch_ln = patch.getPatchLevelNumber();
-
- if (d_cf_discretization == "Ewing" && patch_ln > d_ln_min) {
- ewingFixFlux(patch,
- w_data,
- Dgradw_data,
- ratio_to_coarser_level);
- }
-
-}
-
-void StokesFACOps::computeResidualOnPatch(
- const hier::Patch& patch,
- const pdat::SideData<double>& flux_data,
- const pdat::CellData<double>& soln_data,
- const pdat::CellData<double>& rhs_data,
- pdat::CellData<double>& residual_data) const
-{
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS5(d_dim, patch, flux_data, soln_data, rhs_data,
- residual_data);
-
- tbox::Pointer<geom::CartesianPatchGeometry> patch_geom =
- patch.getPatchGeometry();
- const hier::Box& box = patch.getBox();
- const int* lower = &box.lower()[0];
- const int* upper = &box.upper()[0];
- const double* dx = patch_geom->getDx();
-
- tbox::Pointer<pdat::CellData<double> > scalar_field_data;
- double scalar_field_constant;
- if (d_stokes_spec.cIsVariable()) {
- scalar_field_data =
- patch.getPatchData(d_stokes_spec.getCPatchDataId());
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(compresvarsca2d, COMPRESVARSCA2D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- residual_data.getPointer(),
- &residual_data.getGhostCellWidth()[0],
- &residual_data.getGhostCellWidth()[1],
- scalar_field_data->getPointer(),
- &scalar_field_data->getGhostCellWidth()[0],
- &scalar_field_data->getGhostCellWidth()[1],
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &lower[0], &upper[0], &lower[1], &upper[1],
- dx);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(compresvarsca3d, COMPRESVARSCA3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- &rhs_data.getGhostCellWidth()[2],
- residual_data.getPointer(),
- &residual_data.getGhostCellWidth()[0],
- &residual_data.getGhostCellWidth()[1],
- &residual_data.getGhostCellWidth()[2],
- scalar_field_data->getPointer(),
- &scalar_field_data->getGhostCellWidth()[0],
- &scalar_field_data->getGhostCellWidth()[1],
- &scalar_field_data->getGhostCellWidth()[2],
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &lower[0], &upper[0], &lower[1], &upper[1], &lower[2], &upper[2],
- dx);
- }
- } else if (d_stokes_spec.cIsConstant()) {
- scalar_field_constant = d_stokes_spec.getCConstant();
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(compresconsca2d, COMPRESCONSCA2D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- residual_data.getPointer(),
- &residual_data.getGhostCellWidth()[0],
- &residual_data.getGhostCellWidth()[1],
- scalar_field_constant,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &lower[0], &upper[0], &lower[1], &upper[1],
- dx);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(compresconsca3d, COMPRESCONSCA3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- &rhs_data.getGhostCellWidth()[2],
- residual_data.getPointer(),
- &residual_data.getGhostCellWidth()[0],
- &residual_data.getGhostCellWidth()[1],
- &residual_data.getGhostCellWidth()[2],
- scalar_field_constant,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &lower[0], &upper[0], &lower[1], &upper[1], &lower[2], &upper[2],
- dx);
- }
- } else {
- scalar_field_constant = 0.0;
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(compresconsca2d, COMPRESCONSCA2D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- residual_data.getPointer(),
- &residual_data.getGhostCellWidth()[0],
- &residual_data.getGhostCellWidth()[1],
- 0.0,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &lower[0], &upper[0], &lower[1], &upper[1],
- dx);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(compresconsca3d, COMPRESCONSCA3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- &rhs_data.getGhostCellWidth()[2],
- residual_data.getPointer(),
- &residual_data.getGhostCellWidth()[0],
- &residual_data.getGhostCellWidth()[1],
- &residual_data.getGhostCellWidth()[2],
- 0.0,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &lower[0], &upper[0], &lower[1], &upper[1], &lower[2], &upper[2],
- dx);
- }
- }
-}
-
-void StokesFACOps::redOrBlackSmoothingOnPatch(
- const hier::Patch& patch,
- const pdat::SideData<double>& flux_data,
- const pdat::CellData<double>& rhs_data,
- pdat::CellData<double>& soln_data,
- char red_or_black,
- double* p_maxres) const
-{
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(d_dim, patch, flux_data, soln_data, rhs_data);
- TBOX_ASSERT(red_or_black == 'r' || red_or_black == 'b');
-
- const int offset = red_or_black == 'r' ? 0 : 1;
- tbox::Pointer<geom::CartesianPatchGeometry> patch_geom =
- patch.getPatchGeometry();
- const hier::Box& box = patch.getBox();
- const int* lower = &box.lower()[0];
- const int* upper = &box.upper()[0];
- const double* dx = patch_geom->getDx();
-
- tbox::Pointer<pdat::CellData<double> > scalar_field_data;
- double scalar_field_constant;
- tbox::Pointer<pdat::SideData<double> > diffcoef_data;
- double diffcoef_constant;
-
- if (d_stokes_spec.cIsVariable()) {
- scalar_field_data =
- patch.getPatchData(d_stokes_spec.getCPatchDataId());
- } else if (d_stokes_spec.cIsConstant()) {
- scalar_field_constant = d_stokes_spec.getCConstant();
- } else {
- scalar_field_constant = 0.0;
- }
- if (d_stokes_spec.dIsVariable()) {
- diffcoef_data = patch.getPatchData(d_stokes_spec.getDPatchDataId());
- } else {
- diffcoef_constant = d_stokes_spec.getDConstant();
- }
-
- double maxres = 0.0;
- if (d_stokes_spec.dIsVariable() && d_stokes_spec.cIsVariable()) {
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(rbgswithfluxmaxvardcvarsf2d, RBGSWITHFLUXMAXVARDCVARSF2D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- diffcoef_data->getPointer(0),
- diffcoef_data->getPointer(1),
- &diffcoef_data->getGhostCellWidth()[0],
- &diffcoef_data->getGhostCellWidth()[1],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- scalar_field_data->getPointer(),
- &scalar_field_data->getGhostCellWidth()[0],
- &scalar_field_data->getGhostCellWidth()[1],
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- dx,
- &offset, &maxres);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(rbgswithfluxmaxvardcvarsf3d, RBGSWITHFLUXMAXVARDCVARSF3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- diffcoef_data->getPointer(0),
- diffcoef_data->getPointer(1),
- diffcoef_data->getPointer(2),
- &diffcoef_data->getGhostCellWidth()[0],
- &diffcoef_data->getGhostCellWidth()[1],
- &diffcoef_data->getGhostCellWidth()[2],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- &rhs_data.getGhostCellWidth()[2],
- scalar_field_data->getPointer(),
- &scalar_field_data->getGhostCellWidth()[0],
- &scalar_field_data->getGhostCellWidth()[1],
- &scalar_field_data->getGhostCellWidth()[2],
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- &lower[2], &upper[2],
- dx,
- &offset, &maxres);
- }
- } else if (d_stokes_spec.dIsVariable() && d_stokes_spec.cIsConstant()) {
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(rbgswithfluxmaxvardcconsf2d, RBGSWITHFLUXMAXVARDCCONSF2D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- diffcoef_data->getPointer(0),
- diffcoef_data->getPointer(1),
- &diffcoef_data->getGhostCellWidth()[0],
- &diffcoef_data->getGhostCellWidth()[1],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- scalar_field_constant,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- dx,
- &offset, &maxres);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(rbgswithfluxmaxvardcconsf3d, RBGSWITHFLUXMAXVARDCCONSF3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- diffcoef_data->getPointer(0),
- diffcoef_data->getPointer(1),
- diffcoef_data->getPointer(2),
- &diffcoef_data->getGhostCellWidth()[0],
- &diffcoef_data->getGhostCellWidth()[1],
- &diffcoef_data->getGhostCellWidth()[2],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- &rhs_data.getGhostCellWidth()[2],
- scalar_field_constant,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- &lower[2], &upper[2],
- dx,
- &offset, &maxres);
- }
- } else if (d_stokes_spec.dIsVariable() && d_stokes_spec.cIsZero()) {
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(rbgswithfluxmaxvardcconsf2d, RBGSWITHFLUXMAXVARDCCONSF2D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- diffcoef_data->getPointer(0),
- diffcoef_data->getPointer(1),
- &diffcoef_data->getGhostCellWidth()[0],
- &diffcoef_data->getGhostCellWidth()[1],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- 0.0,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- dx,
- &offset, &maxres);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(rbgswithfluxmaxvardcconsf3d, RBGSWITHFLUXMAXVARDCCONSF3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- diffcoef_data->getPointer(0),
- diffcoef_data->getPointer(1),
- diffcoef_data->getPointer(2),
- &diffcoef_data->getGhostCellWidth()[0],
- &diffcoef_data->getGhostCellWidth()[1],
- &diffcoef_data->getGhostCellWidth()[2],
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- &rhs_data.getGhostCellWidth()[2],
- 0.0,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- &lower[2], &upper[2],
- dx,
- &offset, &maxres);
- }
- } else if (!d_stokes_spec.dIsVariable() && d_stokes_spec.cIsVariable()) {
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(rbgswithfluxmaxcondcvarsf2d, RBGSWITHFLUXMAXCONDCVARSF2D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- diffcoef_constant,
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- scalar_field_data->getPointer(),
- &scalar_field_data->getGhostCellWidth()[0],
- &scalar_field_data->getGhostCellWidth()[1],
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- dx,
- &offset, &maxres);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(rbgswithfluxmaxcondcvarsf3d, RBGSWITHFLUXMAXCONDCVARSF3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- diffcoef_constant,
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- &rhs_data.getGhostCellWidth()[2],
- scalar_field_data->getPointer(),
- &scalar_field_data->getGhostCellWidth()[0],
- &scalar_field_data->getGhostCellWidth()[1],
- &scalar_field_data->getGhostCellWidth()[2],
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- &lower[2], &upper[2],
- dx,
- &offset, &maxres);
- }
- } else if (!d_stokes_spec.dIsVariable() && d_stokes_spec.cIsConstant()) {
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(rbgswithfluxmaxcondcconsf2d, RBGSWITHFLUXMAXCONDCCONSF2D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- diffcoef_constant,
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- scalar_field_constant,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- dx,
- &offset, &maxres);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(rbgswithfluxmaxcondcconsf3d, RBGSWITHFLUXMAXCONDCCONSF3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- diffcoef_constant,
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- &rhs_data.getGhostCellWidth()[2],
- scalar_field_constant,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- &lower[2], &upper[2],
- dx,
- &offset, &maxres);
- }
- } else if (!d_stokes_spec.dIsVariable() && d_stokes_spec.cIsZero()) {
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(rbgswithfluxmaxcondcconsf2d, RBGSWITHFLUXMAXCONDCCONSF2D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- diffcoef_constant,
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- 0.0,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- dx,
- &offset, &maxres);
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(rbgswithfluxmaxcondcconsf3d, RBGSWITHFLUXMAXCONDCCONSF3D) (
- flux_data.getPointer(0),
- flux_data.getPointer(1),
- flux_data.getPointer(2),
- &flux_data.getGhostCellWidth()[0],
- &flux_data.getGhostCellWidth()[1],
- &flux_data.getGhostCellWidth()[2],
- diffcoef_constant,
- rhs_data.getPointer(),
- &rhs_data.getGhostCellWidth()[0],
- &rhs_data.getGhostCellWidth()[1],
- &rhs_data.getGhostCellWidth()[2],
- 0.0,
- soln_data.getPointer(),
- &soln_data.getGhostCellWidth()[0],
- &soln_data.getGhostCellWidth()[1],
- &soln_data.getGhostCellWidth()[2],
- &lower[0], &upper[0],
- &lower[1], &upper[1],
- &lower[2], &upper[2],
- dx,
- &offset, &maxres);
- }
- }
-
- *p_maxres = maxres;
-}
-
-void
-StokesFACOps::xeqScheduleProlongation(
- int dst_id,
- int src_id,
- int scr_id,
- int dest_ln)
-{
- if (!d_prolongation_refine_schedules[dest_ln]) {
- TBOX_ERROR("Expected schedule not found.");
- }
- xfer::RefineAlgorithm refiner(d_dim);
- refiner.
- registerRefine(dst_id,
- src_id,
- scr_id,
- d_prolongation_refine_operator);
- refiner.
- resetSchedule(d_prolongation_refine_schedules[dest_ln]);
- d_prolongation_refine_schedules[dest_ln]->fillData(0.0);
- d_prolongation_refine_algorithm->
- resetSchedule(d_prolongation_refine_schedules[dest_ln]);
-}
-
-void
-StokesFACOps::xeqScheduleURestriction(
- int dst_id,
- int src_id,
- int dest_ln)
-{
- if (!d_urestriction_coarsen_schedules[dest_ln]) {
- TBOX_ERROR("Expected schedule not found.");
- }
-
- xfer::CoarsenAlgorithm coarsener(d_dim);
- coarsener.registerCoarsen(dst_id,
- src_id,
- d_urestriction_coarsen_operator);
- coarsener.resetSchedule(d_urestriction_coarsen_schedules[dest_ln]);
- d_urestriction_coarsen_schedules[dest_ln]->coarsenData();
- d_urestriction_coarsen_algorithm->
- resetSchedule(d_urestriction_coarsen_schedules[dest_ln]);
-}
-
-void
-StokesFACOps::xeqScheduleRRestriction(
- int dst_id,
- int src_id,
- int dest_ln)
-{
- if (!d_rrestriction_coarsen_schedules[dest_ln]) {
- TBOX_ERROR("Expected schedule not found.");
- }
-
- xfer::CoarsenAlgorithm coarsener(d_dim);
- coarsener.registerCoarsen(dst_id,
- src_id,
- d_rrestriction_coarsen_operator);
- coarsener.resetSchedule(d_rrestriction_coarsen_schedules[dest_ln]);
- d_rrestriction_coarsen_schedules[dest_ln]->coarsenData();
- d_rrestriction_coarsen_algorithm->
- resetSchedule(d_rrestriction_coarsen_schedules[dest_ln]);
-}
-
-void
-StokesFACOps::xeqScheduleFluxCoarsen(
- int dst_id,
- int src_id,
- int dest_ln)
-{
- if (!d_flux_coarsen_schedules[dest_ln]) {
- TBOX_ERROR("Expected schedule not found.");
- }
-
- xfer::CoarsenAlgorithm coarsener(d_dim);
- coarsener.registerCoarsen(dst_id,
- src_id,
- d_flux_coarsen_operator);
-
- coarsener.resetSchedule(d_flux_coarsen_schedules[dest_ln]);
- d_flux_coarsen_schedules[dest_ln]->coarsenData();
- d_flux_coarsen_algorithm->
- resetSchedule(d_flux_coarsen_schedules[dest_ln]);
-}
-
-void
-StokesFACOps::xeqScheduleGhostFill(
- int dst_id,
- int dest_ln)
-{
- if (!d_ghostfill_refine_schedules[dest_ln]) {
- TBOX_ERROR("Expected schedule not found.");
- }
- xfer::RefineAlgorithm refiner(d_dim);
- refiner.
- registerRefine(dst_id,
- dst_id,
- dst_id,
- d_ghostfill_refine_operator);
- refiner.
- resetSchedule(d_ghostfill_refine_schedules[dest_ln]);
- d_ghostfill_refine_schedules[dest_ln]->fillData(0.0);
- d_ghostfill_refine_algorithm->
- resetSchedule(d_ghostfill_refine_schedules[dest_ln]);
-}
-
-void
-StokesFACOps::xeqScheduleGhostFillNoCoarse(
- int dst_id,
- int dest_ln)
-{
- if (!d_ghostfill_nocoarse_refine_schedules[dest_ln]) {
- TBOX_ERROR("Expected schedule not found.");
- }
- xfer::RefineAlgorithm refiner(d_dim);
- refiner.
- registerRefine(dst_id,
- dst_id,
- dst_id,
- d_ghostfill_nocoarse_refine_operator);
- refiner.
- resetSchedule(d_ghostfill_nocoarse_refine_schedules[dest_ln]);
- d_ghostfill_nocoarse_refine_schedules[dest_ln]->fillData(0.0);
- d_ghostfill_nocoarse_refine_algorithm->
- resetSchedule(d_ghostfill_nocoarse_refine_schedules[dest_ln]);
-}
-
-void
-StokesFACOps::finalizeCallback()
-{
- for (int d = 0; d < tbox::Dimension::MAXIMUM_DIMENSION_VALUE; ++d) {
- s_cell_scratch_var[d].setNull();
- s_flux_scratch_var[d].setNull();
- s_oflux_scratch_var[d].setNull();
- }
-}
-
-}
-}
-#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps.h
--- a/StokesFACOps.h Sat Jan 01 23:30:12 2011 -0800
+++ b/StokesFACOps.h Sun Jan 02 00:12:38 2011 -0800
@@ -137,8 +137,7 @@ public:
*
* Deallocate internal data.
*/
- ~StokesFACOps(
- void);
+ ~StokesFACOps(void) {}
/*!
* @brief Set the scalar Stokes equation specifications.
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/StokesFACOps.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/StokesFACOps.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,232 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+#ifndef SAMRAI_INLINE
+#include "StokesFACOps.I"
+#endif
+
+namespace SAMRAI {
+ namespace solv {
+
+ tbox::Pointer<pdat::CellVariable<double> >
+ StokesFACOps::s_cell_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+
+ tbox::Pointer<pdat::SideVariable<double> >
+ StokesFACOps::s_flux_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+
+ tbox::Pointer<pdat::OutersideVariable<double> >
+ StokesFACOps::s_oflux_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+
+ tbox::StartupShutdownManager::Handler
+ StokesFACOps::s_finalize_handler(
+ 0,
+ 0,
+ 0,
+ StokesFACOps::finalizeCallback,
+ tbox::StartupShutdownManager::priorityVariables);
+
+ /*
+********************************************************************
+* Constructor. *
+********************************************************************
+*/
+ StokesFACOps::StokesFACOps(
+ const tbox::Dimension& dim,
+ const std::string& object_name,
+ tbox::Pointer<tbox::Database> database):
+ d_dim(dim),
+ d_object_name(object_name),
+ d_hierarchy(),
+ d_ln_min(-1),
+ d_ln_max(-1),
+ d_cf_boundary(),
+ d_stokes_spec(object_name + "::Stokes specs"),
+ d_smoothing_choice("redblack"),
+ d_coarse_solver_choice(
+#ifdef HAVE_HYPRE
+ "hypre"
+#else
+ "redblack"
+#endif
+
+ ),
+ d_cf_discretization("Ewing"),
+ d_prolongation_method("CONSTANT_REFINE"),
+ d_coarse_solver_tolerance(1.e-8),
+ d_coarse_solver_max_iterations(10),
+ d_residual_tolerance_during_smoothing(-1.0),
+ d_flux_id(-1),
+#ifdef HAVE_HYPRE
+ d_hypre_solver(dim,
+ object_name + "::hypre_solver",
+ database && database->isDatabase("hypre_solver") ?
+ database->getDatabase("hypre_solver"):
+ tbox::Pointer<tbox::Database>(NULL)),
+#endif
+ d_physical_bc_coef(NULL),
+ d_context(hier::VariableDatabase::getDatabase()
+ ->getContext(object_name + "::PRIVATE_CONTEXT")),
+ d_cell_scratch_id(-1),
+ d_flux_scratch_id(-1),
+ d_oflux_scratch_id(-1),
+ d_prolongation_refine_operator(),
+ d_prolongation_refine_algorithm(),
+ d_prolongation_refine_schedules(),
+ d_urestriction_coarsen_operator(),
+ d_urestriction_coarsen_algorithm(),
+ d_urestriction_coarsen_schedules(),
+ d_rrestriction_coarsen_operator(),
+ d_rrestriction_coarsen_algorithm(),
+ d_rrestriction_coarsen_schedules(),
+ d_flux_coarsen_operator(),
+ d_flux_coarsen_algorithm(),
+ d_flux_coarsen_schedules(),
+ d_ghostfill_refine_operator(),
+ d_ghostfill_refine_algorithm(),
+ d_ghostfill_refine_schedules(),
+ d_ghostfill_nocoarse_refine_operator(),
+ d_ghostfill_nocoarse_refine_algorithm(),
+ d_ghostfill_nocoarse_refine_schedules(),
+ d_bc_helper(dim,
+ d_object_name + "::bc helper"),
+ d_enable_logging(false),
+ d_preconditioner(NULL),
+ d_hopscell(),
+ d_hopsside()
+ {
+
+ t_restrict_solution = tbox::TimerManager::getManager()->
+ getTimer("solv::StokesFACOps::restrictSolution()");
+ t_restrict_residual = tbox::TimerManager::getManager()->
+ getTimer("solv::StokesFACOps::restrictResidual()");
+ t_prolong = tbox::TimerManager::getManager()->
+ getTimer("solv::StokesFACOps::prolongErrorAndCorrect()");
+ t_smooth_error = tbox::TimerManager::getManager()->
+ getTimer("solv::StokesFACOps::smoothError()");
+ t_solve_coarsest = tbox::TimerManager::getManager()->
+ getTimer("solv::StokesFACOps::solveCoarsestLevel()");
+ t_compute_composite_residual = tbox::TimerManager::getManager()->
+ getTimer("solv::StokesFACOps::computeCompositeResidualOnLevel()");
+ t_compute_residual_norm = tbox::TimerManager::getManager()->
+ getTimer("solv::StokesFACOps::computeResidualNorm()");
+
+ if (d_dim == tbox::Dimension(1) || d_dim > tbox::Dimension(3)) {
+ TBOX_ERROR("StokesFACOps : DIM == 1 or > 3 not implemented yet.\n");
+ }
+
+ if (s_cell_scratch_var[dim.getValue() - 1].isNull()) {
+ TBOX_ASSERT(s_cell_scratch_var[dim.getValue() - 1].isNull());
+ TBOX_ASSERT(s_cell_scratch_var[dim.getValue() - 1].isNull());
+
+ std::ostringstream ss;
+ ss << "StokesFACOps::private_cell_scratch" << dim.getValue();
+ s_cell_scratch_var[dim.getValue() - 1] = new pdat::CellVariable<double>
+ (dim, ss.str());
+ ss.str("");
+ ss << "StokesFACOps::private_flux_scratch" << dim.getValue();
+ s_flux_scratch_var[dim.getValue() - 1] = new pdat::SideVariable<double>
+ (dim, ss.str());
+ ss.str("");
+ ss << "StokesFACOps::private_oflux_scratch" << dim.getValue();
+ s_oflux_scratch_var[dim.getValue() - 1] = new pdat::OutersideVariable<double>
+ (dim, ss.str());
+ }
+
+ hier::VariableDatabase* vdb = hier::VariableDatabase::getDatabase();
+ d_cell_scratch_id = vdb->
+ registerVariableAndContext(s_cell_scratch_var[dim.getValue() - 1],
+ d_context,
+ hier::IntVector::getOne(dim));
+ d_flux_scratch_id = vdb->
+ registerVariableAndContext(s_flux_scratch_var[dim.getValue() - 1],
+ d_context,
+ hier::IntVector::getZero(d_dim));
+ d_oflux_scratch_id = vdb->
+ registerVariableAndContext(s_oflux_scratch_var[dim.getValue() - 1],
+ d_context,
+ hier::IntVector::getZero(d_dim));
+
+ /*
+ * Some variables initialized by default are overriden by input.
+ */
+ if (database) {
+
+ d_coarse_solver_choice =
+ database->getStringWithDefault("coarse_solver_choice",
+ d_coarse_solver_choice);
+ d_coarse_solver_tolerance =
+ database->getDoubleWithDefault("coarse_solver_tolerance",
+ d_coarse_solver_tolerance);
+ d_coarse_solver_max_iterations =
+ database->getIntegerWithDefault("coarse_solver_max_iterations",
+ d_coarse_solver_max_iterations);
+ d_smoothing_choice =
+ database->getStringWithDefault("smoothing_choice",
+ d_smoothing_choice);
+
+ d_cf_discretization =
+ database->getStringWithDefault("cf_discretization",
+ d_cf_discretization);
+
+ d_prolongation_method =
+ database->getStringWithDefault("prolongation_method",
+ d_prolongation_method);
+
+ d_enable_logging =
+ database->getBoolWithDefault("enable_logging",
+ d_enable_logging);
+
+ }
+
+ /*
+ * Check input validity and correctness.
+ */
+ checkInputPatchDataIndices();
+
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/checkInputPatchDataIndices.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/checkInputPatchDataIndices.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,93 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+* Check the validity and correctness of input data for this class. *
+********************************************************************
+*/
+
+ void StokesFACOps::checkInputPatchDataIndices() const {
+ /*
+ * Check input validity and correctness.
+ */
+ hier::VariableDatabase& vdb(*hier::VariableDatabase::getDatabase());
+
+ if (!d_stokes_spec.dIsConstant()
+ && d_stokes_spec.getDPatchDataId() != -1) {
+ tbox::Pointer<hier::Variable> var;
+ vdb.mapIndexToVariable(d_stokes_spec.getDPatchDataId(), var);
+ tbox::Pointer<pdat::SideVariable<double> > diffcoef_var = var;
+ if (!diffcoef_var) {
+ TBOX_ERROR(d_object_name
+ << ": Bad diffusion coefficient patch data index.");
+ }
+ }
+
+ if (!d_stokes_spec.cIsConstant() && !d_stokes_spec.cIsZero()) {
+ tbox::Pointer<hier::Variable> var;
+ vdb.mapIndexToVariable(d_stokes_spec.getCPatchDataId(), var);
+ tbox::Pointer<pdat::CellVariable<double> > scalar_field_var = var;
+ if (!scalar_field_var) {
+ TBOX_ERROR(d_object_name << ": Bad linear term patch data index.");
+ }
+ }
+
+ if (d_flux_id != -1) {
+ tbox::Pointer<hier::Variable> var;
+ vdb.mapIndexToVariable(d_flux_id, var);
+ tbox::Pointer<pdat::SideVariable<double> > flux_var = var;
+
+ TBOX_ASSERT(flux_var);
+ }
+
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/computeCompositeResidualOnLevel.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,210 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+* FACOperatorStrategy virtual *
+* computeCompositeResidualOnLevel function *
+********************************************************************
+*/
+
+ void 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 soln_id = solution.getComponentDescriptorIndex(0);
+ d_bc_helper.setTargetDataId(soln_id);
+ d_bc_helper.setHomogeneousBc(error_equation_indicator);
+
+ const int flux_id = (d_flux_id != -1) ? d_flux_id : d_flux_scratch_id;
+
+ /*
+ * 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. */
+ {
+ tbox::Pointer<xfer::RefineSchedule> ln_refine_schedule;
+ if (ln > d_ln_min) {
+ /* Fill from current, next coarser level and physical boundary */
+ xeqScheduleGhostFill(soln_id, ln);
+ } else {
+ /* Fill from current and physical boundary */
+ xeqScheduleGhostFillNoCoarse(soln_id, ln);
+ }
+ }
+
+ /*
+ * For the whole level, make sure the internal
+ * side-centered data is allocated and note
+ * whether that data should be deallocated when done.
+ * We do this for the whole level because the data
+ * undergoes transfer operations which require the
+ * whole level data.
+ */
+ bool deallocate_flux_data_when_done = false;
+ if (flux_id == d_flux_scratch_id) {
+ if (!level->checkAllocated(flux_id)) {
+ level->allocatePatchData(flux_id);
+ deallocate_flux_data_when_done = true;
+ }
+ }
+
+ /*
+ * S2. Compute flux on patches in level.
+ */
+ for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
+ tbox::Pointer<hier::Patch> patch = *pi;
+
+ tbox::Pointer<pdat::CellData<double> >
+ soln_data = solution.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::CellData<double> >
+ rhs_data = rhs.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::CellData<double> >
+ residual_data = residual.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::SideData<double> >
+ flux_data = patch->getPatchData(flux_id);
+ computeFluxOnPatch(
+ *patch,
+ level->getRatioToCoarserLevel(),
+ *soln_data,
+ *flux_data);
+
+ }
+
+ /*
+ * S3. Coarsen oflux data from next finer level so that
+ * the computed flux becomes the composite grid flux.
+ */
+ if (ln < d_ln_max) {
+ xeqScheduleFluxCoarsen(flux_id, d_oflux_scratch_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> >
+ soln_data = solution.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::CellData<double> >
+ rhs_data = rhs.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::CellData<double> >
+ residual_data = residual.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::SideData<double> >
+ flux_data = patch->getPatchData(flux_id);
+ computeResidualOnPatch(*patch,
+ *flux_data,
+ *soln_data,
+ *rhs_data,
+ *residual_data);
+
+ if (ln > d_ln_min) {
+ /*
+ * Save outerflux data so that next coarser level
+ * can compute its coarse-fine composite flux.
+ * This is not strictly needed in this "compute residual"
+ * loop through the patches, but we put it here to
+ * avoid writing another loop for it.
+ */
+ tbox::Pointer<pdat::OutersideData<double> >
+ oflux_data = patch->getPatchData(d_oflux_scratch_id);
+
+ TBOX_ASSERT(oflux_data);
+
+ oflux_data->copy(*flux_data);
+ }
+ }
+
+ if (deallocate_flux_data_when_done) {
+ level->deallocatePatchData(flux_id);
+ }
+
+ t_compute_composite_residual->stop();
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/computeFluxOnPatch.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/computeFluxOnPatch.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,251 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+extern "C" {
+
+void F77_FUNC(compfluxvardc2d, COMPFLUXVARDC2D) (
+ double* xflux,
+ double* yflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const double* xdiff_coef,
+ const double* ydiff_coef,
+ const int* dcgi,
+ const int* dcgj,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* dx);
+void F77_FUNC(compfluxcondc2d, COMPFLUXCONDC2D) (
+ double* xflux,
+ double* yflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const double & diff_coef,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* dx);
+
+void F77_FUNC(compfluxvardc3d, COMPFLUXVARDC3D) (
+ double* xflux,
+ double* yflux,
+ double* zflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const int* fluxgk,
+ const double* xdiff_coef,
+ const double* ydiff_coef,
+ const double* zdiff_coef,
+ const int* dcgi,
+ const int* dcgj,
+ const int* dcgk,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* solngk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* dx);
+void F77_FUNC(compfluxcondc3d, COMPFLUXCONDC3D) (
+ double* xflux,
+ double* yflux,
+ double* zflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const int* fluxgk,
+ const double & diff_coef,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* solngk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* dx);
+
+}
+ /*
+*******************************************************************
+* *
+* AMR-unaware patch-centered computational kernels. *
+* *
+*******************************************************************
+*/
+
+ void StokesFACOps::computeFluxOnPatch(
+ const hier::Patch& patch,
+ const hier::IntVector& ratio_to_coarser_level,
+ const pdat::CellData<double>& w_data,
+ pdat::SideData<double>& Dgradw_data) const
+ {
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(d_dim, patch, ratio_to_coarser_level, w_data,
+ Dgradw_data);
+ TBOX_ASSERT(patch.inHierarchy());
+ TBOX_ASSERT(w_data.getGhostCellWidth() >=
+ hier::IntVector::getOne(ratio_to_coarser_level.getDim()));
+
+ tbox::Pointer<geom::CartesianPatchGeometry> patch_geom =
+ patch.getPatchGeometry();
+ const hier::Box& box = patch.getBox();
+ const int* lower = &box.lower()[0];
+ const int* upper = &box.upper()[0];
+ const double* dx = patch_geom->getDx();
+
+ double D_value;
+ tbox::Pointer<pdat::SideData<double> > D_data;
+ if (d_stokes_spec.dIsConstant()) {
+ D_value = d_stokes_spec.getDConstant();
+ } else {
+ D_data = patch.getPatchData(d_stokes_spec.getDPatchDataId());
+ }
+
+ if (d_stokes_spec.dIsConstant()) {
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(compfluxcondc2d, COMPFLUXCONDC2D) (
+ Dgradw_data.getPointer(0),
+ Dgradw_data.getPointer(1),
+ &Dgradw_data.getGhostCellWidth()[0],
+ &Dgradw_data.getGhostCellWidth()[1],
+ D_value,
+ w_data.getPointer(),
+ &w_data.getGhostCellWidth()[0],
+ &w_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ dx);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(compfluxcondc3d, COMPFLUXCONDC3D) (
+ Dgradw_data.getPointer(0),
+ Dgradw_data.getPointer(1),
+ Dgradw_data.getPointer(2),
+ &Dgradw_data.getGhostCellWidth()[0],
+ &Dgradw_data.getGhostCellWidth()[1],
+ &Dgradw_data.getGhostCellWidth()[2],
+ D_value,
+ w_data.getPointer(),
+ &w_data.getGhostCellWidth()[0],
+ &w_data.getGhostCellWidth()[1],
+ &w_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ &lower[2], &upper[2],
+ dx);
+ }
+ } else {
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(compfluxvardc2d, COMPFLUXVARDC2D) (
+ Dgradw_data.getPointer(0),
+ Dgradw_data.getPointer(1),
+ &Dgradw_data.getGhostCellWidth()[0],
+ &Dgradw_data.getGhostCellWidth()[1],
+ D_data->getPointer(0),
+ D_data->getPointer(1),
+ &D_data->getGhostCellWidth()[0],
+ &D_data->getGhostCellWidth()[1],
+ w_data.getPointer(),
+ &w_data.getGhostCellWidth()[0],
+ &w_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ dx);
+ }
+ if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(compfluxvardc3d, COMPFLUXVARDC3D) (
+ Dgradw_data.getPointer(0),
+ Dgradw_data.getPointer(1),
+ Dgradw_data.getPointer(2),
+ &Dgradw_data.getGhostCellWidth()[0],
+ &Dgradw_data.getGhostCellWidth()[1],
+ &Dgradw_data.getGhostCellWidth()[2],
+ D_data->getPointer(0),
+ D_data->getPointer(1),
+ D_data->getPointer(2),
+ &D_data->getGhostCellWidth()[0],
+ &D_data->getGhostCellWidth()[1],
+ &D_data->getGhostCellWidth()[2],
+ w_data.getPointer(),
+ &w_data.getGhostCellWidth()[0],
+ &w_data.getGhostCellWidth()[1],
+ &w_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ &lower[2], &upper[2],
+ dx);
+ }
+ }
+
+ const int patch_ln = patch.getPatchLevelNumber();
+
+ if (d_cf_discretization == "Ewing" && patch_ln > d_ln_min) {
+ ewingFixFlux(patch,
+ w_data,
+ Dgradw_data,
+ ratio_to_coarser_level);
+ }
+
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/computeResidualNorm.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/computeResidualNorm.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,88 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+* FACOperatorStrategy virtual computeResidualNorm *
+* function *
+********************************************************************
+*/
+
+ double StokesFACOps::computeResidualNorm(const SAMRAIVectorReal<double>& residual,
+ int fine_ln,
+ int coarse_ln)
+ {
+
+ if (coarse_ln != residual.getCoarsestLevelNumber() ||
+ fine_ln != residual.getFinestLevelNumber()) {
+ TBOX_ERROR("StokesFACOps::computeResidualNorm() is not\n"
+ << "set up to compute residual except on the range of\n"
+ << "levels defining the vector.\n");
+ }
+ t_compute_residual_norm->start();
+ /*
+ * The residual vector was cloned from vectors that has
+ * the proper weights associated with them, so we do not
+ * have to explicitly weight the residuals.
+ *
+ * maxNorm: not good to use because Hypre's norm does not
+ * correspond to it. Also maybe too sensitive to spikes.
+ * L2Norm: maybe good. But does not correspond to the
+ * scale of the quantity.
+ * L1Norm: maybe good. Correspond to scale of quantity,
+ * but may be too insensitive to spikes.
+ * RMSNorm: maybe good.
+ */
+ double norm = residual.RMSNorm();
+ t_compute_residual_norm->stop();
+ return norm;
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/computeResidualOnPatch.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/computeResidualOnPatch.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,313 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+extern "C" {
+void F77_FUNC(compresvarsca2d, COMPRESVARSCA2D) (
+ const double* xflux,
+ const double* yflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ double* residual,
+ const int* residualgi,
+ const int* residualgj,
+ const double* scalar_field,
+ const int* scalar_field_gi,
+ const int* scalar_field_gj,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* dx);
+void F77_FUNC(compresconsca2d, COMPRESCONSCA2D) (
+ const double* xflux,
+ const double* yflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ double* residual,
+ const int* residualgi,
+ const int* residualgj,
+ const double & scalar_field,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* dx);
+void F77_FUNC(compresvarsca3d, COMPRESVARSCA3D) (
+ const double* xflux,
+ const double* yflux,
+ const double* zflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const int* fluxgk,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ const int* rhsgk,
+ double* residual,
+ const int* residualgi,
+ const int* residualgj,
+ const int* residualgk,
+ const double* scalar_field,
+ const int* scalar_field_gi,
+ const int* scalar_field_gj,
+ const int* scalar_field_gk,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* solngk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* dx);
+void F77_FUNC(compresconsca3d, COMPRESCONSCA3D) (
+ const double* xflux,
+ const double* yflux,
+ const double* zflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const int* fluxgk,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ const int* rhsgk,
+ double* residual,
+ const int* residualgi,
+ const int* residualgj,
+ const int* residualgk,
+ const double & scalar_field,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* solngk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* dx);
+}
+
+ void StokesFACOps::computeResidualOnPatch(
+ const hier::Patch& patch,
+ const pdat::SideData<double>& flux_data,
+ const pdat::CellData<double>& soln_data,
+ const pdat::CellData<double>& rhs_data,
+ pdat::CellData<double>& residual_data) const
+ {
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS5(d_dim, patch, flux_data, soln_data, rhs_data,
+ residual_data);
+
+ tbox::Pointer<geom::CartesianPatchGeometry> patch_geom =
+ patch.getPatchGeometry();
+ const hier::Box& box = patch.getBox();
+ const int* lower = &box.lower()[0];
+ const int* upper = &box.upper()[0];
+ const double* dx = patch_geom->getDx();
+
+ tbox::Pointer<pdat::CellData<double> > scalar_field_data;
+ double scalar_field_constant;
+ if (d_stokes_spec.cIsVariable()) {
+ scalar_field_data =
+ patch.getPatchData(d_stokes_spec.getCPatchDataId());
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(compresvarsca2d, COMPRESVARSCA2D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ residual_data.getPointer(),
+ &residual_data.getGhostCellWidth()[0],
+ &residual_data.getGhostCellWidth()[1],
+ scalar_field_data->getPointer(),
+ &scalar_field_data->getGhostCellWidth()[0],
+ &scalar_field_data->getGhostCellWidth()[1],
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0], &lower[1], &upper[1],
+ dx);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(compresvarsca3d, COMPRESVARSCA3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ &rhs_data.getGhostCellWidth()[2],
+ residual_data.getPointer(),
+ &residual_data.getGhostCellWidth()[0],
+ &residual_data.getGhostCellWidth()[1],
+ &residual_data.getGhostCellWidth()[2],
+ scalar_field_data->getPointer(),
+ &scalar_field_data->getGhostCellWidth()[0],
+ &scalar_field_data->getGhostCellWidth()[1],
+ &scalar_field_data->getGhostCellWidth()[2],
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0], &lower[1], &upper[1], &lower[2], &upper[2],
+ dx);
+ }
+ } else if (d_stokes_spec.cIsConstant()) {
+ scalar_field_constant = d_stokes_spec.getCConstant();
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(compresconsca2d, COMPRESCONSCA2D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ residual_data.getPointer(),
+ &residual_data.getGhostCellWidth()[0],
+ &residual_data.getGhostCellWidth()[1],
+ scalar_field_constant,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0], &lower[1], &upper[1],
+ dx);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(compresconsca3d, COMPRESCONSCA3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ &rhs_data.getGhostCellWidth()[2],
+ residual_data.getPointer(),
+ &residual_data.getGhostCellWidth()[0],
+ &residual_data.getGhostCellWidth()[1],
+ &residual_data.getGhostCellWidth()[2],
+ scalar_field_constant,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0], &lower[1], &upper[1], &lower[2], &upper[2],
+ dx);
+ }
+ } else {
+ scalar_field_constant = 0.0;
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(compresconsca2d, COMPRESCONSCA2D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ residual_data.getPointer(),
+ &residual_data.getGhostCellWidth()[0],
+ &residual_data.getGhostCellWidth()[1],
+ 0.0,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0], &lower[1], &upper[1],
+ dx);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(compresconsca3d, COMPRESCONSCA3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ &rhs_data.getGhostCellWidth()[2],
+ residual_data.getPointer(),
+ &residual_data.getGhostCellWidth()[0],
+ &residual_data.getGhostCellWidth()[1],
+ &residual_data.getGhostCellWidth()[2],
+ 0.0,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0], &lower[1], &upper[1], &lower[2], &upper[2],
+ dx);
+ }
+ }
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/computeVectorWeights.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/computeVectorWeights.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,151 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+* Compute the vector weight and put it at a specified patch data *
+* index. *
+********************************************************************
+*/
+
+ void StokesFACOps::computeVectorWeights(
+ tbox::Pointer<hier::PatchHierarchy> hierarchy,
+ int weight_id,
+ int coarsest_ln,
+ int finest_ln) const
+ {
+ TBOX_ASSERT(!hierarchy.isNull());
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS1(d_dim, *hierarchy);
+
+ if (coarsest_ln == -1) coarsest_ln = 0;
+ if (finest_ln == -1) finest_ln = hierarchy->getFinestLevelNumber();
+ if (finest_ln < coarsest_ln) {
+ TBOX_ERROR(d_object_name
+ << ": Illegal level number range. finest_ln < coarsest_ln.");
+ }
+
+ int ln;
+ for (ln = finest_ln; ln >= coarsest_ln; --ln) {
+
+ /*
+ * On every level, first assign cell volume to vector weight.
+ */
+
+ tbox::Pointer<hier::PatchLevel> level =
+ hierarchy->getPatchLevel(ln);
+ for (hier::PatchLevel::Iterator p(level); p; p++) {
+ tbox::Pointer<hier::Patch> patch = *p;
+ tbox::Pointer<geom::CartesianPatchGeometry> patch_geometry =
+ patch->getPatchGeometry();
+ const double* dx = patch_geometry->getDx();
+ double cell_vol = dx[0];
+ if (d_dim > tbox::Dimension(1)) {
+ cell_vol *= dx[1];
+ }
+
+ if (d_dim > tbox::Dimension(2)) {
+ cell_vol *= dx[2];
+ }
+
+ tbox::Pointer<pdat::CellData<double> > w =
+ patch->getPatchData(weight_id);
+ if (!w) {
+ TBOX_ERROR(d_object_name
+ << ": weight id must refer to a pdat::CellVariable");
+ }
+ w->fillAll(cell_vol);
+ }
+
+ /*
+ * On all but the finest level, assign 0 to vector
+ * weight to cells covered by finer cells.
+ */
+
+ if (ln < finest_ln) {
+
+ /*
+ * First get the boxes that describe index space of the next finer
+ * level and coarsen them to describe corresponding index space
+ * at this level.
+ */
+
+ tbox::Pointer<hier::PatchLevel> next_finer_level =
+ hierarchy->getPatchLevel(ln + 1);
+ hier::BoxArray coarsened_boxes = next_finer_level->getBoxes();
+ hier::IntVector coarsen_ratio(next_finer_level->getRatioToLevelZero());
+ coarsen_ratio /= level->getRatioToLevelZero();
+ coarsened_boxes.coarsen(coarsen_ratio);
+
+ /*
+ * Then set vector weight to 0 wherever there is
+ * a nonempty intersection with the next finer level.
+ * Note that all assignments are local.
+ */
+
+ for (hier::PatchLevel::Iterator p(level); p; p++) {
+
+ tbox::Pointer<hier::Patch> patch = *p;
+ for (int i = 0; i < coarsened_boxes.getNumberOfBoxes(); i++) {
+
+ hier::Box coarse_box = coarsened_boxes[i];
+ hier::Box intersection = coarse_box * (patch->getBox());
+ if (!intersection.empty()) {
+ tbox::Pointer<pdat::CellData<double> > w =
+ patch->getPatchData(weight_id);
+ w->fillAll(0.0, intersection);
+
+ } // assignment only in non-empty intersection
+ } // loop over coarsened boxes from finer level
+ } // loop over patches in level
+ } // all levels except finest
+ } // loop over levels
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/deallocateOperatorState.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/deallocateOperatorState.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,96 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+* FACOperatorStrategy virtual deallocateOperatorState *
+* function. Deallocate internal hierarchy-dependent data. *
+* State is allocated iff hierarchy is set. *
+********************************************************************
+*/
+
+ void StokesFACOps::deallocateOperatorState()
+ {
+ if (d_hierarchy) {
+ int ln;
+ for (ln = d_ln_min + 1; ln <= d_ln_max; ++ln) {
+ d_hierarchy->getPatchLevel(ln)->
+ deallocatePatchData(d_oflux_scratch_id);
+ }
+ d_cf_boundary.resizeArray(0);
+#ifdef HAVE_HYPRE
+ d_hypre_solver.deallocateSolverState();
+#endif
+ d_hierarchy.setNull();
+ d_ln_min = -1;
+ d_ln_max = -1;
+
+ d_prolongation_refine_algorithm.setNull();
+ d_prolongation_refine_schedules.setNull();
+
+ d_urestriction_coarsen_algorithm.setNull();
+ d_urestriction_coarsen_schedules.setNull();
+
+ d_rrestriction_coarsen_algorithm.setNull();
+ d_rrestriction_coarsen_schedules.setNull();
+
+ d_flux_coarsen_algorithm.setNull();
+ d_flux_coarsen_schedules.setNull();
+
+ d_ghostfill_refine_algorithm.setNull();
+ d_ghostfill_refine_schedules.setNull();
+
+ d_ghostfill_nocoarse_refine_algorithm.setNull();
+ d_ghostfill_nocoarse_refine_schedules.setNull();
+
+ }
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/ewingFixFlux.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/ewingFixFlux.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,286 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+extern "C" {
+void F77_FUNC(ewingfixfluxvardc2d, EWINGFIXFLUXVARDC2D) (
+ const double* xflux,
+ const double* yflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const double* xdiff_coef,
+ const double* ydiff_coef,
+ const int* dcgi,
+ const int* dcgj,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* location_index,
+ const int* ratio_to_coarser,
+ const int* blower,
+ const int* bupper,
+ const double* dx);
+void F77_FUNC(ewingfixfluxcondc2d, EWINGFIXFLUXCONDC2D) (
+ const double* xflux,
+ const double* yflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const double & diff_coef,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* location_index,
+ const int* ratio_to_coarser,
+ const int* blower,
+ const int* bupper,
+ const double* dx);
+
+void F77_FUNC(ewingfixfluxvardc3d, EWINGFIXFLUXVARDC3D) (
+ const double* xflux,
+ const double* yflux,
+ const double* zflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const int* fluxgk,
+ const double* xdiff_coef,
+ const double* ydiff_coef,
+ const double* zdiff_coef,
+ const int* dcgi,
+ const int* dcgj,
+ const int* dcgk,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* solngk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const int* location_index,
+ const int* ratio_to_coarser,
+ const int* blower,
+ const int* bupper,
+ const double* dx);
+void F77_FUNC(ewingfixfluxcondc3d, EWINGFIXFLUXCONDC3D) (
+ const double* xflux,
+ const double* yflux,
+ const double* zflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const int* fluxgk,
+ const double & diff_coef,
+ const double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* solngk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const int* location_index,
+ const int* ratio_to_coarser,
+ const int* blower,
+ const int* bupper,
+ const double* dx);
+}
+ /*
+********************************************************************
+* Fix flux on coarse-fine boundaries computed from a *
+* constant-refine interpolation of coarse level data. *
+********************************************************************
+*/
+
+ void StokesFACOps::ewingFixFlux(
+ const hier::Patch& patch,
+ const pdat::CellData<double>& soln_data,
+ pdat::SideData<double>& flux_data,
+ const hier::IntVector& ratio_to_coarser) const
+ {
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(d_dim, patch, soln_data, flux_data,
+ ratio_to_coarser);
+
+ const int patch_ln = patch.getPatchLevelNumber();
+ const hier::GlobalId id = patch.getGlobalId();
+ tbox::Pointer<geom::CartesianPatchGeometry> patch_geom =
+ patch.getPatchGeometry();
+ const double* dx = patch_geom->getDx();
+ const hier::Box& patch_box(patch.getBox());
+ const hier::Index& plower = patch_box.lower();
+ const hier::Index& pupper = patch_box.upper();
+
+ const tbox::Array<hier::BoundaryBox>& bboxes =
+ d_cf_boundary[patch_ln]->getBoundaries(id, 1);
+ int bn, nboxes = bboxes.getSize();
+
+ if (d_stokes_spec.dIsVariable()) {
+
+ tbox::Pointer<pdat::SideData<double> > diffcoef_data;
+ diffcoef_data = patch.getPatchData(d_stokes_spec.getDPatchDataId());
+
+ for (bn = 0; bn < nboxes; ++bn) {
+ const hier::BoundaryBox& boundary_box = bboxes[bn];
+
+ TBOX_ASSERT(boundary_box.getBoundaryType() == 1);
+
+ const hier::Box& bdry_box = boundary_box.getBox();
+ const hier::Index& blower = bdry_box.lower();
+ const hier::Index& bupper = bdry_box.upper();
+ const int location_index = boundary_box.getLocationIndex();
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(ewingfixfluxvardc2d, EWINGFIXFLUXVARDC2D) (
+ flux_data.getPointer(0), flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ diffcoef_data->getPointer(0), diffcoef_data->getPointer(1),
+ &diffcoef_data->getGhostCellWidth()[0],
+ &diffcoef_data->getGhostCellWidth()[1],
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &plower[0], &pupper[0], &plower[1], &pupper[1],
+ &location_index,
+ &ratio_to_coarser[0],
+ &blower[0], &bupper[0],
+ dx);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(ewingfixfluxvardc3d, EWINGFIXFLUXVARDC3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ diffcoef_data->getPointer(0),
+ diffcoef_data->getPointer(1),
+ diffcoef_data->getPointer(2),
+ &diffcoef_data->getGhostCellWidth()[0],
+ &diffcoef_data->getGhostCellWidth()[1],
+ &diffcoef_data->getGhostCellWidth()[2],
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &plower[0], &pupper[0],
+ &plower[1], &pupper[1],
+ &plower[2], &pupper[2],
+ &location_index,
+ &ratio_to_coarser[0],
+ &blower[0], &bupper[0],
+ dx);
+ } else {
+ TBOX_ERROR("StokesFACOps : DIM > 3 not supported" << std::endl);
+ }
+
+ }
+ } else {
+
+ const double diffcoef_constant = d_stokes_spec.getDConstant();
+
+ for (bn = 0; bn < nboxes; ++bn) {
+ const hier::BoundaryBox& boundary_box = bboxes[bn];
+
+ TBOX_ASSERT(boundary_box.getBoundaryType() == 1);
+
+ const hier::Box& bdry_box = boundary_box.getBox();
+ const hier::Index& blower = bdry_box.lower();
+ const hier::Index& bupper = bdry_box.upper();
+ const int location_index = boundary_box.getLocationIndex();
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(ewingfixfluxcondc2d, EWINGFIXFLUXCONDC2D) (
+ flux_data.getPointer(0), flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ diffcoef_constant,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &plower[0], &pupper[0],
+ &plower[1], &pupper[1],
+ &location_index,
+ &ratio_to_coarser[0],
+ &blower[0], &bupper[0],
+ dx);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(ewingfixfluxcondc3d, EWINGFIXFLUXCONDC3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ diffcoef_constant,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &plower[0], &pupper[0],
+ &plower[1], &pupper[1],
+ &plower[2], &pupper[2],
+ &location_index,
+ &ratio_to_coarser[0],
+ &blower[0], &bupper[0],
+ dx);
+ }
+ }
+ }
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/finalizeCallback.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/finalizeCallback.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,61 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ void
+ StokesFACOps::finalizeCallback()
+ {
+ for (int d = 0; d < tbox::Dimension::MAXIMUM_DIMENSION_VALUE; ++d) {
+ s_cell_scratch_var[d].setNull();
+ s_flux_scratch_var[d].setNull();
+ s_oflux_scratch_var[d].setNull();
+ }
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/initializeOperatorState.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/initializeOperatorState.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,420 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+************************************************************************
+* FACOperatorStrategy virtual initializeOperatorState function. *
+* *
+* Set internal variables to correspond to the solution passed in. *
+* Look up transfer operators. *
+************************************************************************
+*/
+
+ void StokesFACOps::initializeOperatorState(
+ const SAMRAIVectorReal<double>& solution,
+ const SAMRAIVectorReal<double>& rhs)
+ {
+ deallocateOperatorState();
+ int ln;
+ hier::VariableDatabase* vdb = hier::VariableDatabase::getDatabase();
+
+ d_hierarchy = solution.getPatchHierarchy();
+ d_ln_min = solution.getCoarsestLevelNumber();
+ d_ln_max = solution.getFinestLevelNumber();
+ d_hopscell = new math::HierarchyCellDataOpsReal<double>(d_hierarchy,
+ d_ln_min,
+ d_ln_max);
+ d_hopsside = new math::HierarchySideDataOpsReal<double>(d_hierarchy,
+ d_ln_min,
+ d_ln_max);
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+
+ if (d_physical_bc_coef == NULL) {
+ /*
+ * It's an error not to have bc object set.
+ * Note that the bc object cannot be passed in through
+ * the argument because the interface is inherited.
+ */
+ TBOX_ERROR(
+ d_object_name << ": No physical bc object in\n"
+ <<
+ "StokesFACOps::initializeOperatorState\n"
+ << "You must use "
+ <<
+ "StokesFACOps::setPhysicalBcCoefObject\n"
+ <<
+ "to set one before calling initializeOperatorState\n");
+ }
+
+ if (solution.getNumberOfComponents() != 1) {
+ TBOX_WARNING(d_object_name
+ << ": Solution vector has multiple components.\n"
+ << "Solver is for component 0 only.\n");
+ }
+ if (rhs.getNumberOfComponents() != 1) {
+ TBOX_WARNING(d_object_name
+ << ": RHS vector has multiple components.\n"
+ << "Solver is for component 0 only.\n");
+ }
+
+ /*
+ * Make sure that solution and rhs data
+ * are of correct type
+ * are allocated
+ * has sufficient ghost width
+ */
+ tbox::Pointer<hier::Variable> var;
+ {
+ vdb->mapIndexToVariable(rhs.getComponentDescriptorIndex(0),
+ var);
+ if (!var) {
+ TBOX_ERROR(d_object_name << ": RHS component does not\n"
+ << "correspond to a variable.\n");
+ }
+ tbox::Pointer<pdat::CellVariable<double> > cell_var = var;
+ if (!cell_var) {
+ TBOX_ERROR(d_object_name
+ << ": RHS variable is not cell-centered double\n");
+ }
+ }
+ {
+ vdb->mapIndexToVariable(solution.getComponentDescriptorIndex(0),
+ var);
+ if (!var) {
+ TBOX_ERROR(d_object_name << ": Solution component does not\n"
+ << "correspond to a variable.\n");
+ }
+ tbox::Pointer<pdat::CellVariable<double> > cell_var = var;
+ if (!cell_var) {
+ TBOX_ERROR(d_object_name
+ << ": Solution variable is not cell-centered double\n");
+ }
+ }
+ for (ln = d_ln_min; ln <= d_ln_max; ++ln) {
+ tbox::Pointer<hier::PatchLevel> level_ptr =
+ d_hierarchy->getPatchLevel(ln);
+ hier::PatchLevel& level = *level_ptr;
+ for (hier::PatchLevel::Iterator pi(level); pi; pi++) {
+ hier::Patch& patch = **pi;
+ tbox::Pointer<hier::PatchData> fd =
+ patch.getPatchData(rhs.getComponentDescriptorIndex(0));
+ if (fd) {
+ /*
+ * Some data checks can only be done if the data already exists.
+ */
+ tbox::Pointer<pdat::CellData<double> > cd = fd;
+ if (!cd) {
+ TBOX_ERROR(d_object_name
+ << ": RHS data is not cell-centered double\n");
+ }
+ if (cd->getDepth() > 1) {
+ TBOX_WARNING(d_object_name
+ << ": RHS data has multiple depths.\n"
+ << "Solver is for depth 0 only.\n");
+ }
+ }
+ tbox::Pointer<hier::PatchData> ud =
+ patch.getPatchData(solution.getComponentDescriptorIndex(0));
+ if (ud) {
+ /*
+ * Some data checks can only be done if the data already exists.
+ */
+ tbox::Pointer<pdat::CellData<double> > cd = ud;
+ if (!cd) {
+ TBOX_ERROR(d_object_name
+ << ": Solution data is not cell-centered double\n");
+ }
+ if (cd->getDepth() > 1) {
+ TBOX_WARNING(d_object_name
+ << ": Solution data has multiple depths.\n"
+ << "Solver is for depth 0 only.\n");
+ }
+ if (cd->getGhostCellWidth() < hier::IntVector::getOne(d_dim)) {
+ TBOX_ERROR(d_object_name
+ << ": Solution data has insufficient ghost width\n");
+ }
+ }
+ }
+ }
+
+ /*
+ * Solution and rhs must have some similar properties.
+ */
+ if (rhs.getPatchHierarchy() != d_hierarchy
+ || rhs.getCoarsestLevelNumber() != d_ln_min
+ || rhs.getFinestLevelNumber() != d_ln_max) {
+ TBOX_ERROR(d_object_name << ": solution and rhs do not have\n"
+ << "the same set of patch levels.\n");
+ }
+
+#endif
+
+ /*
+ * Initialize the coarse-fine boundary description for the
+ * hierarchy.
+ */
+ d_cf_boundary.resizeArray(d_hierarchy->getNumberOfLevels());
+
+ hier::IntVector max_gcw(d_dim, 1);
+ for (ln = d_ln_min; ln <= d_ln_max; ++ln) {
+ d_cf_boundary[ln] = new hier::CoarseFineBoundary(*d_hierarchy,
+ ln,
+ max_gcw);
+ }
+#ifdef HAVE_HYPRE
+ if (d_coarse_solver_choice == "hypre") {
+ d_hypre_solver.initializeSolverState(d_hierarchy, d_ln_min);
+ /*
+ * Share the boundary condition object with the hypre solver
+ * to make sure that boundary condition settings are consistent
+ * between the two objects.
+ */
+ d_hypre_solver.setPhysicalBcCoefObject(d_physical_bc_coef);
+ d_hypre_solver.setMatrixCoefficients(d_stokes_spec);
+ }
+#endif
+
+ /*
+ * Get the transfer operators.
+ * Flux coarsening is conservative.
+ * Cell (solution, error, etc) coarsening is conservative.
+ * Cell refinement from same level is constant refinement.
+ * Cell refinement from coarser level is chosen by the
+ * choice of coarse-fine discretization, d_cf_discretization,
+ * which should be set to either "Ewing" or one of the
+ * acceptable strings for looking up the refine operator.
+ */
+ tbox::Pointer<geom::CartesianGridGeometry> geometry =
+ d_hierarchy->getGridGeometry();
+ tbox::Pointer<hier::Variable> variable;
+
+ vdb->mapIndexToVariable(d_cell_scratch_id, variable);
+ d_prolongation_refine_operator =
+ geometry->lookupRefineOperator(variable,
+ d_prolongation_method);
+
+ vdb->mapIndexToVariable(d_cell_scratch_id, variable);
+ d_urestriction_coarsen_operator =
+ d_rrestriction_coarsen_operator =
+ geometry->lookupCoarsenOperator(variable,
+ "CONSERVATIVE_COARSEN");
+
+ vdb->mapIndexToVariable(d_oflux_scratch_id, variable);
+ d_flux_coarsen_operator =
+ geometry->lookupCoarsenOperator(variable,
+ "CONSERVATIVE_COARSEN");
+
+ vdb->mapIndexToVariable(d_cell_scratch_id, variable);
+ d_ghostfill_refine_operator =
+ geometry->lookupRefineOperator(variable,
+ d_cf_discretization == "Ewing" ?
+ "CONSTANT_REFINE" : d_cf_discretization);
+
+ vdb->mapIndexToVariable(d_cell_scratch_id, variable);
+ d_ghostfill_nocoarse_refine_operator =
+ geometry->lookupRefineOperator(variable,
+ "CONSTANT_REFINE");
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (!d_prolongation_refine_operator) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot find prolongation refine operator");
+ }
+ if (!d_urestriction_coarsen_operator) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot find restriction coarsening operator");
+ }
+ if (!d_rrestriction_coarsen_operator) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot find restriction coarsening operator");
+ }
+ if (!d_flux_coarsen_operator) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot find flux coarsening operator");
+ }
+ if (!d_ghostfill_refine_operator) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot find ghost filling refinement operator");
+ }
+ if (!d_ghostfill_nocoarse_refine_operator) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot find ghost filling refinement operator");
+ }
+#endif
+
+ for (ln = d_ln_min + 1; ln <= d_ln_max; ++ln) {
+ d_hierarchy->getPatchLevel(ln)->
+ allocatePatchData(d_oflux_scratch_id);
+ }
+
+ /*
+ * Make space for saving communication schedules.
+ * There is no need to delete the old schedules first
+ * because we have deallocated the solver state above.
+ */
+ d_prolongation_refine_schedules.resizeArray(d_ln_max + 1);
+ d_ghostfill_refine_schedules.resizeArray(d_ln_max + 1);
+ d_ghostfill_nocoarse_refine_schedules.resizeArray(d_ln_max + 1);
+ d_urestriction_coarsen_schedules.resizeArray(d_ln_max + 1);
+ d_rrestriction_coarsen_schedules.resizeArray(d_ln_max + 1);
+ d_flux_coarsen_schedules.resizeArray(d_ln_max + 1);
+
+ d_prolongation_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
+ d_urestriction_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
+ d_rrestriction_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
+ d_flux_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
+ d_ghostfill_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
+ d_ghostfill_nocoarse_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
+
+ d_prolongation_refine_algorithm->
+ registerRefine(d_cell_scratch_id,
+ solution.getComponentDescriptorIndex(0),
+ d_cell_scratch_id,
+ d_prolongation_refine_operator);
+ d_urestriction_coarsen_algorithm->
+ registerCoarsen(solution.getComponentDescriptorIndex(0),
+ solution.getComponentDescriptorIndex(0),
+ d_urestriction_coarsen_operator);
+ d_rrestriction_coarsen_algorithm->
+ registerCoarsen(rhs.getComponentDescriptorIndex(0),
+ rhs.getComponentDescriptorIndex(0),
+ d_rrestriction_coarsen_operator);
+ d_ghostfill_refine_algorithm->
+ registerRefine(solution.getComponentDescriptorIndex(0),
+ solution.getComponentDescriptorIndex(0),
+ solution.getComponentDescriptorIndex(0),
+ d_ghostfill_refine_operator);
+ d_flux_coarsen_algorithm->
+ registerCoarsen(((d_flux_id != -1) ? d_flux_id : d_flux_scratch_id),
+ d_oflux_scratch_id,
+ d_flux_coarsen_operator);
+ d_ghostfill_nocoarse_refine_algorithm->
+ registerRefine(solution.getComponentDescriptorIndex(0),
+ solution.getComponentDescriptorIndex(0),
+ solution.getComponentDescriptorIndex(0),
+ d_ghostfill_nocoarse_refine_operator);
+
+ for (int dest_ln = d_ln_min + 1; dest_ln <= d_ln_max; ++dest_ln) {
+
+ tbox::Pointer<xfer::PatchLevelFullFillPattern> fill_pattern(
+ new xfer::PatchLevelFullFillPattern());
+ d_prolongation_refine_schedules[dest_ln] =
+ d_prolongation_refine_algorithm->
+ createSchedule(fill_pattern,
+ d_hierarchy->getPatchLevel(dest_ln),
+ tbox::Pointer<hier::PatchLevel>(),
+ dest_ln - 1,
+ d_hierarchy,
+ &d_bc_helper);
+ if (!d_prolongation_refine_schedules[dest_ln]) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot create a refine schedule for prolongation!\n");
+ }
+ d_ghostfill_refine_schedules[dest_ln] =
+ d_ghostfill_refine_algorithm->
+ createSchedule(d_hierarchy->getPatchLevel(dest_ln),
+ dest_ln - 1,
+ d_hierarchy,
+ &d_bc_helper);
+ if (!d_ghostfill_refine_schedules[dest_ln]) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot create a refine schedule for ghost filling!\n");
+ }
+ d_ghostfill_nocoarse_refine_schedules[dest_ln] =
+ d_ghostfill_nocoarse_refine_algorithm->
+ createSchedule(d_hierarchy->getPatchLevel(dest_ln),
+ &d_bc_helper);
+ if (!d_ghostfill_nocoarse_refine_schedules[dest_ln]) {
+ TBOX_ERROR(
+ d_object_name
+ <<
+ ": Cannot create a refine schedule for ghost filling on bottom level!\n");
+ }
+ }
+ for (int dest_ln = d_ln_min; dest_ln < d_ln_max; ++dest_ln) {
+ d_urestriction_coarsen_schedules[dest_ln] =
+ d_urestriction_coarsen_algorithm->
+ createSchedule(d_hierarchy->getPatchLevel(dest_ln),
+ d_hierarchy->getPatchLevel(dest_ln + 1));
+ if (!d_urestriction_coarsen_schedules[dest_ln]) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot create a coarsen schedule for U restriction!\n");
+ }
+ d_rrestriction_coarsen_schedules[dest_ln] =
+ d_rrestriction_coarsen_algorithm->
+ createSchedule(d_hierarchy->getPatchLevel(dest_ln),
+ d_hierarchy->getPatchLevel(dest_ln + 1));
+ if (!d_rrestriction_coarsen_schedules[dest_ln]) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot create a coarsen schedule for R restriction!\n");
+ }
+ d_flux_coarsen_schedules[dest_ln] =
+ d_flux_coarsen_algorithm->
+ createSchedule(d_hierarchy->getPatchLevel(dest_ln),
+ d_hierarchy->getPatchLevel(dest_ln + 1));
+ if (!d_flux_coarsen_schedules[dest_ln]) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot create a coarsen schedule for flux transfer!\n");
+ }
+ }
+ d_ghostfill_nocoarse_refine_schedules[d_ln_min] =
+ d_ghostfill_nocoarse_refine_algorithm->
+ createSchedule(d_hierarchy->getPatchLevel(d_ln_min),
+ &d_bc_helper);
+ if (!d_ghostfill_nocoarse_refine_schedules[d_ln_min]) {
+ TBOX_ERROR(
+ d_object_name
+ <<
+ ": Cannot create a refine schedule for ghost filling on bottom level!\n");
+ }
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/postprocessOneCycle.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/postprocessOneCycle.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,84 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+* FACOperatorStrategy virtual postprocessOneCycle function. *
+********************************************************************
+*/
+
+ void StokesFACOps::postprocessOneCycle(
+ int fac_cycle_num,
+ const SAMRAIVectorReal<double>& current_soln,
+ const SAMRAIVectorReal<double>& residual)
+ {
+ NULL_USE(current_soln);
+ NULL_USE(residual);
+
+ if (d_enable_logging) {
+ if (d_preconditioner) {
+ /*
+ * Output convergence progress. This is probably only appropriate
+ * if the solver is NOT being used as a preconditioner.
+ */
+ double avg_factor, final_factor;
+ d_preconditioner->getConvergenceFactors(avg_factor, final_factor);
+ tbox::plog
+ << "iter=" << std::setw(4) << fac_cycle_num
+ << " resid=" << d_preconditioner->getResidualNorm()
+ << " net conv=" << d_preconditioner->getNetConvergenceFactor()
+ << " final conv=" << d_preconditioner->getNetConvergenceFactor()
+ << " avg conv=" << d_preconditioner->getAvgConvergenceFactor()
+ << std::endl;
+ }
+ }
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/prolongErrorAndCorrect.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/prolongErrorAndCorrect.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,114 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+***********************************************************************
+* FACOperatorStrategy virtual prolongErrorAndCorrect function. *
+* After the prolongation, we set the physical boundary condition *
+* for the correction, which is zero. Other ghost cell values, *
+* which are preset to zero, need not be set. *
+***********************************************************************
+*/
+
+ void StokesFACOps::prolongErrorAndCorrect(
+ const SAMRAIVectorReal<double>& s,
+ SAMRAIVectorReal<double>& d,
+ int dest_ln) {
+
+ t_prolong->start();
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (s.getPatchHierarchy() != d_hierarchy
+ || d.getPatchHierarchy() != d_hierarchy) {
+ TBOX_ERROR(d_object_name << ": Vector hierarchy does not match\n"
+ "internal state hierarchy.");
+ }
+#endif
+
+ tbox::Pointer<hier::PatchLevel> coarse_level =
+ d_hierarchy->getPatchLevel(dest_ln - 1);
+ tbox::Pointer<hier::PatchLevel> fine_level =
+ d_hierarchy->getPatchLevel(dest_ln);
+
+ /*
+ * Data is prolonged into the scratch space corresponding
+ * to index d_cell_scratch_id and allocated here.
+ */
+ fine_level->allocatePatchData(d_cell_scratch_id);
+
+ /*
+ * Refine solution into scratch space to fill the fine level
+ * interior in the scratch space, then use that refined data
+ * to correct the fine level error.
+ */
+ d_bc_helper.setTargetDataId(d_cell_scratch_id);
+ d_bc_helper.setHomogeneousBc(true);
+ const int src_index = s.getComponentDescriptorIndex(0);
+ xeqScheduleProlongation(d_cell_scratch_id,
+ src_index,
+ d_cell_scratch_id,
+ dest_ln);
+
+ /*
+ * Add the refined error in the scratch space
+ * to the error currently residing in the destination level.
+ */
+ math::HierarchyCellDataOpsReal<double>
+ hierarchy_math_ops(d_hierarchy, dest_ln, dest_ln);
+ const int dst_index = d.getComponentDescriptorIndex(0);
+ hierarchy_math_ops.add(dst_index, dst_index, d_cell_scratch_id);
+
+ fine_level->deallocatePatchData(d_cell_scratch_id);
+
+ t_prolong->stop();
+
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/redOrBlackSmoothingOnPatch.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/redOrBlackSmoothingOnPatch.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,600 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+extern "C" {
+void F77_FUNC(rbgswithfluxmaxvardcvarsf2d, RBGSWITHFLUXMAXVARDCVARSF2D) (
+ const double* xflux,
+ const double* yflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const double* xdiff_coef,
+ const double* ydiff_coef,
+ const int* dcgi,
+ const int* dcgj,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ const double* scalar_field,
+ const int* scalar_field_gi,
+ const int* scalar_field_gj,
+ double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* dx,
+ const int* offset,
+ const double* maxres);
+void F77_FUNC(rbgswithfluxmaxcondcvarsf2d, RBGSWITHFLUXMAXCONDCVARSF2D) (
+ const double* xflux,
+ const double* yflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const double & dc,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ const double* scalar_field,
+ const int* scalar_field_gi,
+ const int* scalar_field_gj,
+ double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* dx,
+ const int* offset,
+ const double* maxres);
+void F77_FUNC(rbgswithfluxmaxvardcconsf2d, RBGSWITHFLUXMAXVARDCCONSF2D) (
+ const double* xflux,
+ const double* yflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const double* xdiff_coef,
+ const double* ydiff_coef,
+ const int* dcgi,
+ const int* dcgj,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ const double & scalar_field,
+ double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* dx,
+ const int* offset,
+ const double* maxres);
+void F77_FUNC(rbgswithfluxmaxcondcconsf2d, RBGSWITHFLUXMAXCONDCCONSF2D) (
+ const double* xflux,
+ const double* yflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const double & dc,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ const double & scalar_field,
+ double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* dx,
+ const int* offset,
+ const double* maxres);
+void F77_FUNC(rbgswithfluxmaxvardcvarsf3d, RBGSWITHFLUXMAXVARDCVARSF3D) (
+ const double* xflux,
+ const double* yflux,
+ const double* zflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const int* fluxgk,
+ const double* xdiff_coef,
+ const double* ydiff_coef,
+ const double* zdiff_coef,
+ const int* dcgi,
+ const int* dcgj,
+ const int* dcgk,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ const int* rhsgk,
+ const double* scalar_field,
+ const int* scalar_field_gi,
+ const int* scalar_field_gj,
+ const int* scalar_field_gk,
+ double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* solngk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* dx,
+ const int* offset,
+ const double* maxres);
+void F77_FUNC(rbgswithfluxmaxcondcvarsf3d, RBGSWITHFLUXMAXCONDCVARSF3D) (
+ const double* xflux,
+ const double* yflux,
+ const double* zflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const int* fluxgk,
+ const double & dc,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ const int* rhsgk,
+ const double* scalar_field,
+ const int* scalar_field_gi,
+ const int* scalar_field_gj,
+ const int* scalar_field_gk,
+ double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* solngk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* dx,
+ const int* offset,
+ const double* maxres);
+void F77_FUNC(rbgswithfluxmaxvardcconsf3d, RBGSWITHFLUXMAXVARDCCONSF3D) (
+ const double* xflux,
+ const double* yflux,
+ const double* zflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const int* fluxgk,
+ const double* xdiff_coef,
+ const double* ydiff_coef,
+ const double* zdiff_coef,
+ const int* dcgi,
+ const int* dcgj,
+ const int* dcgk,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ const int* rhsgk,
+ const double & scalar_field,
+ double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* solngk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* dx,
+ const int* offset,
+ const double* maxres);
+void F77_FUNC(rbgswithfluxmaxcondcconsf3d, RBGSWITHFLUXMAXCONDCCONSF3D) (
+ const double* xflux,
+ const double* yflux,
+ const double* zflux,
+ const int* fluxgi,
+ const int* fluxgj,
+ const int* fluxgk,
+ const double & dc,
+ const double* rhs,
+ const int* rhsgi,
+ const int* rhsgj,
+ const int* rhsgk,
+ const double & scalar_field,
+ double* soln,
+ const int* solngi,
+ const int* solngj,
+ const int* solngk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* dx,
+ const int* offset,
+ const double* maxres);
+
+}
+ void StokesFACOps::redOrBlackSmoothingOnPatch(
+ const hier::Patch& patch,
+ const pdat::SideData<double>& flux_data,
+ const pdat::CellData<double>& rhs_data,
+ pdat::CellData<double>& soln_data,
+ char red_or_black,
+ double* p_maxres) const
+ {
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(d_dim, patch, flux_data, soln_data, rhs_data);
+ TBOX_ASSERT(red_or_black == 'r' || red_or_black == 'b');
+
+ const int offset = red_or_black == 'r' ? 0 : 1;
+ tbox::Pointer<geom::CartesianPatchGeometry> patch_geom =
+ patch.getPatchGeometry();
+ const hier::Box& box = patch.getBox();
+ const int* lower = &box.lower()[0];
+ const int* upper = &box.upper()[0];
+ const double* dx = patch_geom->getDx();
+
+ tbox::Pointer<pdat::CellData<double> > scalar_field_data;
+ double scalar_field_constant;
+ tbox::Pointer<pdat::SideData<double> > diffcoef_data;
+ double diffcoef_constant;
+
+ if (d_stokes_spec.cIsVariable()) {
+ scalar_field_data =
+ patch.getPatchData(d_stokes_spec.getCPatchDataId());
+ } else if (d_stokes_spec.cIsConstant()) {
+ scalar_field_constant = d_stokes_spec.getCConstant();
+ } else {
+ scalar_field_constant = 0.0;
+ }
+ if (d_stokes_spec.dIsVariable()) {
+ diffcoef_data = patch.getPatchData(d_stokes_spec.getDPatchDataId());
+ } else {
+ diffcoef_constant = d_stokes_spec.getDConstant();
+ }
+
+ double maxres = 0.0;
+ if (d_stokes_spec.dIsVariable() && d_stokes_spec.cIsVariable()) {
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(rbgswithfluxmaxvardcvarsf2d, RBGSWITHFLUXMAXVARDCVARSF2D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ diffcoef_data->getPointer(0),
+ diffcoef_data->getPointer(1),
+ &diffcoef_data->getGhostCellWidth()[0],
+ &diffcoef_data->getGhostCellWidth()[1],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ scalar_field_data->getPointer(),
+ &scalar_field_data->getGhostCellWidth()[0],
+ &scalar_field_data->getGhostCellWidth()[1],
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ dx,
+ &offset, &maxres);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(rbgswithfluxmaxvardcvarsf3d, RBGSWITHFLUXMAXVARDCVARSF3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ diffcoef_data->getPointer(0),
+ diffcoef_data->getPointer(1),
+ diffcoef_data->getPointer(2),
+ &diffcoef_data->getGhostCellWidth()[0],
+ &diffcoef_data->getGhostCellWidth()[1],
+ &diffcoef_data->getGhostCellWidth()[2],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ &rhs_data.getGhostCellWidth()[2],
+ scalar_field_data->getPointer(),
+ &scalar_field_data->getGhostCellWidth()[0],
+ &scalar_field_data->getGhostCellWidth()[1],
+ &scalar_field_data->getGhostCellWidth()[2],
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ &lower[2], &upper[2],
+ dx,
+ &offset, &maxres);
+ }
+ } else if (d_stokes_spec.dIsVariable() && d_stokes_spec.cIsConstant()) {
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(rbgswithfluxmaxvardcconsf2d, RBGSWITHFLUXMAXVARDCCONSF2D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ diffcoef_data->getPointer(0),
+ diffcoef_data->getPointer(1),
+ &diffcoef_data->getGhostCellWidth()[0],
+ &diffcoef_data->getGhostCellWidth()[1],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ scalar_field_constant,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ dx,
+ &offset, &maxres);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(rbgswithfluxmaxvardcconsf3d, RBGSWITHFLUXMAXVARDCCONSF3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ diffcoef_data->getPointer(0),
+ diffcoef_data->getPointer(1),
+ diffcoef_data->getPointer(2),
+ &diffcoef_data->getGhostCellWidth()[0],
+ &diffcoef_data->getGhostCellWidth()[1],
+ &diffcoef_data->getGhostCellWidth()[2],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ &rhs_data.getGhostCellWidth()[2],
+ scalar_field_constant,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ &lower[2], &upper[2],
+ dx,
+ &offset, &maxres);
+ }
+ } else if (d_stokes_spec.dIsVariable() && d_stokes_spec.cIsZero()) {
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(rbgswithfluxmaxvardcconsf2d, RBGSWITHFLUXMAXVARDCCONSF2D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ diffcoef_data->getPointer(0),
+ diffcoef_data->getPointer(1),
+ &diffcoef_data->getGhostCellWidth()[0],
+ &diffcoef_data->getGhostCellWidth()[1],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ 0.0,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ dx,
+ &offset, &maxres);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(rbgswithfluxmaxvardcconsf3d, RBGSWITHFLUXMAXVARDCCONSF3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ diffcoef_data->getPointer(0),
+ diffcoef_data->getPointer(1),
+ diffcoef_data->getPointer(2),
+ &diffcoef_data->getGhostCellWidth()[0],
+ &diffcoef_data->getGhostCellWidth()[1],
+ &diffcoef_data->getGhostCellWidth()[2],
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ &rhs_data.getGhostCellWidth()[2],
+ 0.0,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ &lower[2], &upper[2],
+ dx,
+ &offset, &maxres);
+ }
+ } else if (!d_stokes_spec.dIsVariable() && d_stokes_spec.cIsVariable()) {
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(rbgswithfluxmaxcondcvarsf2d, RBGSWITHFLUXMAXCONDCVARSF2D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ diffcoef_constant,
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ scalar_field_data->getPointer(),
+ &scalar_field_data->getGhostCellWidth()[0],
+ &scalar_field_data->getGhostCellWidth()[1],
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ dx,
+ &offset, &maxres);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(rbgswithfluxmaxcondcvarsf3d, RBGSWITHFLUXMAXCONDCVARSF3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ diffcoef_constant,
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ &rhs_data.getGhostCellWidth()[2],
+ scalar_field_data->getPointer(),
+ &scalar_field_data->getGhostCellWidth()[0],
+ &scalar_field_data->getGhostCellWidth()[1],
+ &scalar_field_data->getGhostCellWidth()[2],
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ &lower[2], &upper[2],
+ dx,
+ &offset, &maxres);
+ }
+ } else if (!d_stokes_spec.dIsVariable() && d_stokes_spec.cIsConstant()) {
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(rbgswithfluxmaxcondcconsf2d, RBGSWITHFLUXMAXCONDCCONSF2D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ diffcoef_constant,
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ scalar_field_constant,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ dx,
+ &offset, &maxres);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(rbgswithfluxmaxcondcconsf3d, RBGSWITHFLUXMAXCONDCCONSF3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ diffcoef_constant,
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ &rhs_data.getGhostCellWidth()[2],
+ scalar_field_constant,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ &lower[2], &upper[2],
+ dx,
+ &offset, &maxres);
+ }
+ } else if (!d_stokes_spec.dIsVariable() && d_stokes_spec.cIsZero()) {
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(rbgswithfluxmaxcondcconsf2d, RBGSWITHFLUXMAXCONDCCONSF2D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ diffcoef_constant,
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ 0.0,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ dx,
+ &offset, &maxres);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(rbgswithfluxmaxcondcconsf3d, RBGSWITHFLUXMAXCONDCCONSF3D) (
+ flux_data.getPointer(0),
+ flux_data.getPointer(1),
+ flux_data.getPointer(2),
+ &flux_data.getGhostCellWidth()[0],
+ &flux_data.getGhostCellWidth()[1],
+ &flux_data.getGhostCellWidth()[2],
+ diffcoef_constant,
+ rhs_data.getPointer(),
+ &rhs_data.getGhostCellWidth()[0],
+ &rhs_data.getGhostCellWidth()[1],
+ &rhs_data.getGhostCellWidth()[2],
+ 0.0,
+ soln_data.getPointer(),
+ &soln_data.getGhostCellWidth()[0],
+ &soln_data.getGhostCellWidth()[1],
+ &soln_data.getGhostCellWidth()[2],
+ &lower[0], &upper[0],
+ &lower[1], &upper[1],
+ &lower[2], &upper[2],
+ dx,
+ &offset, &maxres);
+ }
+ }
+
+ *p_maxres = maxres;
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/restrictResidual.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/restrictResidual.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,71 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+* FACOperatorStrategy virtual restrictresidual function. *
+********************************************************************
+*/
+
+ void StokesFACOps::restrictResidual(
+ const SAMRAIVectorReal<double>& s,
+ SAMRAIVectorReal<double>& d,
+ int dest_ln) {
+
+ t_restrict_residual->start();
+
+ xeqScheduleRRestriction(d.getComponentDescriptorIndex(0),
+ s.getComponentDescriptorIndex(0),
+ dest_ln);
+
+ t_restrict_residual->stop();
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/restrictSolution.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/restrictSolution.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,84 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+* FACOperatorStrategy virtual restrictSolution function. *
+* After restricting solution, update ghost cells of the affected *
+* level. *
+********************************************************************
+*/
+
+ void StokesFACOps::restrictSolution(
+ const SAMRAIVectorReal<double>& s,
+ SAMRAIVectorReal<double>& d,
+ int dest_ln) {
+
+ t_restrict_solution->start();
+
+ xeqScheduleURestriction(d.getComponentDescriptorIndex(0),
+ s.getComponentDescriptorIndex(0),
+ dest_ln);
+
+ d_bc_helper.setHomogeneousBc(false);
+ d_bc_helper.setTargetDataId(d.getComponentDescriptorIndex(0));
+
+ if (dest_ln == d_ln_min) {
+ xeqScheduleGhostFillNoCoarse(d.getComponentDescriptorIndex(0),
+ dest_ln);
+ } else {
+ xeqScheduleGhostFill(d.getComponentDescriptorIndex(0),
+ dest_ln);
+ }
+
+ t_restrict_solution->stop();
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/smoothError.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/smoothError.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,81 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+********************************************************************
+*/
+
+ void StokesFACOps::smoothError(
+ SAMRAIVectorReal<double>& data,
+ const SAMRAIVectorReal<double>& residual,
+ int ln,
+ int num_sweeps)
+ {
+
+ t_smooth_error->start();
+
+ checkInputPatchDataIndices();
+ if (d_smoothing_choice == "redblack") {
+ smoothErrorByRedBlack(data,
+ residual,
+ ln,
+ num_sweeps,
+ d_residual_tolerance_during_smoothing);
+ } else {
+ TBOX_ERROR(d_object_name << ": Bad smoothing choice '"
+ << d_smoothing_choice
+ << "' in StokesFACOps.");
+ }
+
+ t_smooth_error->stop();
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/smoothErrorByRedBlack.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/smoothErrorByRedBlack.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,227 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+* Workhorse function to smooth error using red-black *
+* Gauss-Seidel iterations. *
+********************************************************************
+*/
+
+ void StokesFACOps::smoothErrorByRedBlack(
+ SAMRAIVectorReal<double>& data,
+ const SAMRAIVectorReal<double>& residual,
+ int ln,
+ int num_sweeps,
+ double residual_tolerance)
+ {
+
+ checkInputPatchDataIndices();
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (data.getPatchHierarchy() != d_hierarchy
+ || residual.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);
+
+ const int data_id = data.getComponentDescriptorIndex(0);
+
+ const int flux_id = (d_flux_id != -1) ? d_flux_id : d_flux_scratch_id;
+
+ d_bc_helper.setTargetDataId(data_id);
+ d_bc_helper.setHomogeneousBc(true);
+ xeqScheduleGhostFillNoCoarse(data_id, ln);
+
+ if (ln > d_ln_min) {
+ /*
+ * Perform a one-time transfer of data from coarser level,
+ * to fill ghost boundaries that will not change through
+ * the smoothing loop.
+ */
+ xeqScheduleGhostFill(data_id, ln);
+ }
+
+ /*
+ * Smooth the number of sweeps specified or until
+ * the convergence is satisfactory.
+ */
+ int isweep;
+ double red_maxres, blk_maxres, maxres = 0;
+ red_maxres = blk_maxres = residual_tolerance + 1;
+ /*
+ * Instead of checking residual convergence globally,
+ * we check the not_converged flag. This avoids possible
+ * round-off errors affecting different processes differently,
+ * leading to disagreement on whether to continue smoothing.
+ */
+ int not_converged = 1;
+ for (isweep = 0; isweep < num_sweeps && not_converged; ++isweep) {
+ red_maxres = blk_maxres = 0;
+
+ // Red sweep.
+ xeqScheduleGhostFillNoCoarse(data_id, ln);
+ for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
+ tbox::Pointer<hier::Patch> patch = *pi;
+
+ bool deallocate_flux_data_when_done = false;
+ if (flux_id == d_flux_scratch_id) {
+ /*
+ * Using internal temporary storage for flux.
+ * For each patch, make sure the internal
+ * side-centered data is allocated and note
+ * whether that data should be deallocated when done.
+ */
+ if (!patch->checkAllocated(flux_id)) {
+ patch->allocatePatchData(flux_id);
+ deallocate_flux_data_when_done = true;
+ }
+ }
+
+ tbox::Pointer<pdat::CellData<double> >
+ scalar_field_data = d_stokes_spec.cIsVariable() ?
+ patch->getPatchData(d_stokes_spec.getCPatchDataId()) :
+ tbox::Pointer<hier::PatchData>(NULL);
+ tbox::Pointer<pdat::CellData<double> >
+ err_data = data.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::CellData<double> >
+ residual_data = residual.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::SideData<double> >
+ flux_data = patch->getPatchData(flux_id);
+
+ computeFluxOnPatch(
+ *patch,
+ level->getRatioToCoarserLevel(),
+ *err_data,
+ *flux_data);
+
+ redOrBlackSmoothingOnPatch(*patch,
+ *flux_data,
+ *residual_data,
+ *err_data,
+ 'r',
+ &red_maxres);
+
+ if (deallocate_flux_data_when_done) {
+ patch->deallocatePatchData(flux_id);
+ }
+ } // End patch number *pi
+ xeqScheduleGhostFillNoCoarse(data_id, ln);
+
+ // Black sweep.
+ for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
+ tbox::Pointer<hier::Patch> patch = *pi;
+
+ bool deallocate_flux_data_when_done = false;
+ if (flux_id == d_flux_scratch_id) {
+ /*
+ * Using internal temporary storage for flux.
+ * For each patch, make sure the internal
+ * side-centered data is allocated and note
+ * whether that data should be deallocated when done.
+ */
+ if (!patch->checkAllocated(flux_id)) {
+ patch->allocatePatchData(flux_id);
+ deallocate_flux_data_when_done = true;
+ }
+ }
+
+ tbox::Pointer<pdat::CellData<double> >
+ scalar_field_data = d_stokes_spec.cIsVariable() ?
+ patch->getPatchData(d_stokes_spec.getCPatchDataId()) :
+ tbox::Pointer<hier::PatchData>(NULL);
+ tbox::Pointer<pdat::CellData<double> >
+ err_data = data.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::CellData<double> >
+ residual_data = residual.getComponentPatchData(0, *patch);
+ tbox::Pointer<pdat::SideData<double> >
+ flux_data = patch->getPatchData(flux_id);
+
+ computeFluxOnPatch(
+ *patch,
+ level->getRatioToCoarserLevel(),
+ *err_data,
+ *flux_data);
+
+ redOrBlackSmoothingOnPatch(*patch,
+ *flux_data,
+ *residual_data,
+ *err_data,
+ 'b',
+ &blk_maxres);
+
+ if (deallocate_flux_data_when_done) {
+ patch->deallocatePatchData(flux_id);
+ }
+ } // End patch number *pi
+ xeqScheduleGhostFillNoCoarse(data_id, ln);
+ if (residual_tolerance >= 0.0) {
+ /*
+ * Check for early end of sweeps due to convergence
+ * only if it is numerically possible (user gave a
+ * non negative value for residual tolerance).
+ */
+ maxres = tbox::MathUtilities<double>::Max(red_maxres, blk_maxres);
+ not_converged = maxres > residual_tolerance;
+ const tbox::SAMRAI_MPI& mpi(d_hierarchy->getDomainMappedBoxLevel().getMPI());
+ if (mpi.getSize() > 1) {
+ mpi.AllReduce(¬_converged, 1, MPI_MAX);
+ }
+ }
+ } // End sweep number isweep
+ if (d_enable_logging) tbox::plog
+ << d_object_name << " RBGS smoothing maxres = " << maxres << "\n"
+ << " after " << isweep << " sweeps.\n";
+
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/solveCoarsestLevel.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/solveCoarsestLevel.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,108 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ /*
+********************************************************************
+* FACOperatorStrategy virtual solveCoarsestLevel *
+* function *
+********************************************************************
+*/
+
+ int StokesFACOps::solveCoarsestLevel(
+ SAMRAIVectorReal<double>& data,
+ const SAMRAIVectorReal<double>& residual,
+ int coarsest_ln) {
+
+ t_solve_coarsest->start();
+
+ checkInputPatchDataIndices();
+
+ int return_value = 0;
+
+ if (d_coarse_solver_choice == "jacobi") {
+ d_residual_tolerance_during_smoothing = d_coarse_solver_tolerance;
+ smoothError(data,
+ residual,
+ coarsest_ln,
+ d_coarse_solver_max_iterations);
+ d_residual_tolerance_during_smoothing = -1.0;
+ } else if (d_coarse_solver_choice == "redblack") {
+ d_residual_tolerance_during_smoothing = d_coarse_solver_tolerance;
+ smoothError(data,
+ residual,
+ coarsest_ln,
+ d_coarse_solver_max_iterations);
+ d_residual_tolerance_during_smoothing = -1.0;
+ } else if (d_coarse_solver_choice == "hypre") {
+#ifndef HAVE_HYPRE
+ TBOX_ERROR(d_object_name << ": Coarse level solver choice '"
+ << d_coarse_solver_choice
+ << "' unavailable in "
+ << "scapStokesOps::solveCoarsestLevel.");
+#else
+ return_value = solveCoarsestLevel_HYPRE(data, residual, coarsest_ln);
+#endif
+ } else {
+ TBOX_ERROR(
+ d_object_name << ": Bad coarse level solver choice '"
+ << d_coarse_solver_choice
+ <<
+ "' in scapStokesOps::solveCoarsestLevel.");
+ }
+
+ xeqScheduleGhostFillNoCoarse(data.getComponentDescriptorIndex(0),
+ coarsest_ln);
+
+ t_solve_coarsest->stop();
+
+ return return_value;
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/solveCoarsestLevel_HYPRE.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/solveCoarsestLevel_HYPRE.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,101 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+#ifdef HAVE_HYPRE
+ /*
+********************************************************************
+* Solve coarsest level using Hypre *
+* We only solve for the error, so we always use homogeneous bc. *
+********************************************************************
+*/
+
+ int StokesFACOps::solveCoarsestLevel_HYPRE(
+ SAMRAIVectorReal<double>& data,
+ const SAMRAIVectorReal<double>& residual,
+ int coarsest_ln) {
+
+ NULL_USE(coarsest_ln);
+
+#ifndef HAVE_HYPRE
+ TBOX_ERROR(d_object_name << ": Coarse level solver choice '"
+ << d_coarse_solver_choice
+ << "' unavailable in "
+ << "StokesFACOps::solveCoarsestLevel.");
+
+ return 0;
+
+#else
+
+ checkInputPatchDataIndices();
+ d_hypre_solver.setStoppingCriteria(d_coarse_solver_max_iterations,
+ d_coarse_solver_tolerance);
+ const int solver_ret =
+ d_hypre_solver.solveSystem(
+ data.getComponentDescriptorIndex(0),
+ residual.getComponentDescriptorIndex(0),
+ true);
+ /*
+ * Present data on the solve.
+ * The Hypre solver returns 0 if converged.
+ */
+ if (d_enable_logging) tbox::plog
+ << d_object_name << " Hypre solve " << (solver_ret ? "" : "NOT ")
+ << "converged\n"
+ << "\titerations: " << d_hypre_solver.getNumberOfIterations() << "\n"
+ << "\tresidual: " << d_hypre_solver.getRelativeResidualNorm() << "\n";
+
+ return !solver_ret;
+
+#endif
+
+ }
+#endif
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/xeqScheduleFluxCoarsen.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/xeqScheduleFluxCoarsen.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,72 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ void
+ StokesFACOps::xeqScheduleFluxCoarsen(
+ int dst_id,
+ int src_id,
+ int dest_ln)
+ {
+ if (!d_flux_coarsen_schedules[dest_ln]) {
+ TBOX_ERROR("Expected schedule not found.");
+ }
+
+ xfer::CoarsenAlgorithm coarsener(d_dim);
+ coarsener.registerCoarsen(dst_id,
+ src_id,
+ d_flux_coarsen_operator);
+
+ coarsener.resetSchedule(d_flux_coarsen_schedules[dest_ln]);
+ d_flux_coarsen_schedules[dest_ln]->coarsenData();
+ d_flux_coarsen_algorithm->
+ resetSchedule(d_flux_coarsen_schedules[dest_ln]);
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/xeqScheduleGhostFill.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/xeqScheduleGhostFill.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,72 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ void
+ StokesFACOps::xeqScheduleGhostFill(
+ int dst_id,
+ int dest_ln)
+ {
+ if (!d_ghostfill_refine_schedules[dest_ln]) {
+ TBOX_ERROR("Expected schedule not found.");
+ }
+ xfer::RefineAlgorithm refiner(d_dim);
+ refiner.
+ registerRefine(dst_id,
+ dst_id,
+ dst_id,
+ d_ghostfill_refine_operator);
+ refiner.
+ resetSchedule(d_ghostfill_refine_schedules[dest_ln]);
+ d_ghostfill_refine_schedules[dest_ln]->fillData(0.0);
+ d_ghostfill_refine_algorithm->
+ resetSchedule(d_ghostfill_refine_schedules[dest_ln]);
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/xeqScheduleGhostFillNoCoarse.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/xeqScheduleGhostFillNoCoarse.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,72 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ void
+ StokesFACOps::xeqScheduleGhostFillNoCoarse(
+ int dst_id,
+ int dest_ln)
+ {
+ if (!d_ghostfill_nocoarse_refine_schedules[dest_ln]) {
+ TBOX_ERROR("Expected schedule not found.");
+ }
+ xfer::RefineAlgorithm refiner(d_dim);
+ refiner.
+ registerRefine(dst_id,
+ dst_id,
+ dst_id,
+ d_ghostfill_nocoarse_refine_operator);
+ refiner.
+ resetSchedule(d_ghostfill_nocoarse_refine_schedules[dest_ln]);
+ d_ghostfill_nocoarse_refine_schedules[dest_ln]->fillData(0.0);
+ d_ghostfill_nocoarse_refine_algorithm->
+ resetSchedule(d_ghostfill_nocoarse_refine_schedules[dest_ln]);
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/xeqScheduleProlongation.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/xeqScheduleProlongation.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,74 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ void
+ StokesFACOps::xeqScheduleProlongation(
+ int dst_id,
+ int src_id,
+ int scr_id,
+ int dest_ln)
+ {
+ if (!d_prolongation_refine_schedules[dest_ln]) {
+ TBOX_ERROR("Expected schedule not found.");
+ }
+ xfer::RefineAlgorithm refiner(d_dim);
+ refiner.
+ registerRefine(dst_id,
+ src_id,
+ scr_id,
+ d_prolongation_refine_operator);
+ refiner.
+ resetSchedule(d_prolongation_refine_schedules[dest_ln]);
+ d_prolongation_refine_schedules[dest_ln]->fillData(0.0);
+ d_prolongation_refine_algorithm->
+ resetSchedule(d_prolongation_refine_schedules[dest_ln]);
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/xeqScheduleRRestriction.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/xeqScheduleRRestriction.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,71 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ void
+ StokesFACOps::xeqScheduleRRestriction(
+ int dst_id,
+ int src_id,
+ int dest_ln)
+ {
+ if (!d_rrestriction_coarsen_schedules[dest_ln]) {
+ TBOX_ERROR("Expected schedule not found.");
+ }
+
+ xfer::CoarsenAlgorithm coarsener(d_dim);
+ coarsener.registerCoarsen(dst_id,
+ src_id,
+ d_rrestriction_coarsen_operator);
+ coarsener.resetSchedule(d_rrestriction_coarsen_schedules[dest_ln]);
+ d_rrestriction_coarsen_schedules[dest_ln]->coarsenData();
+ d_rrestriction_coarsen_algorithm->
+ resetSchedule(d_rrestriction_coarsen_schedules[dest_ln]);
+ }
+
+ }
+}
+#endif
diff -r 0e7b866a1bb0 -r 04ca745a8fed StokesFACOps/xeqScheduleURestriction.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/xeqScheduleURestriction.C Sun Jan 02 00:12:38 2011 -0800
@@ -0,0 +1,71 @@
+/*************************************************************************
+ *
+ * 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
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesFACOps_C
+#define included_solv_StokesFACOps_C
+
+#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"
+
+namespace SAMRAI {
+ namespace solv {
+
+ void
+ StokesFACOps::xeqScheduleURestriction(
+ int dst_id,
+ int src_id,
+ int dest_ln)
+ {
+ if (!d_urestriction_coarsen_schedules[dest_ln]) {
+ TBOX_ERROR("Expected schedule not found.");
+ }
+
+ xfer::CoarsenAlgorithm coarsener(d_dim);
+ coarsener.registerCoarsen(dst_id,
+ src_id,
+ d_urestriction_coarsen_operator);
+ coarsener.resetSchedule(d_urestriction_coarsen_schedules[dest_ln]);
+ d_urestriction_coarsen_schedules[dest_ln]->coarsenData();
+ d_urestriction_coarsen_algorithm->
+ resetSchedule(d_urestriction_coarsen_schedules[dest_ln]);
+ }
+
+ }
+}
+#endif
More information about the CIG-COMMITS
mailing list