[cig-commits] commit: Copy over CellPoisson* stuff from SAMRAI internals and rename it to Stokes
Mercurial
hg at geodynamics.org
Fri Feb 25 14:12:21 PST 2011
changeset: 5:fe9d63509c19
user: Walter Landry <wlandry at caltech.edu>
date: Fri Dec 31 07:12:28 2010 -0800
files: CellStokesFACOps.C CellStokesFACOps.I CellStokesFACOps.h CellStokesFACSolver.C CellStokesFACSolver.I CellStokesFACSolver.h CellStokesHypreSolver.C CellStokesHypreSolver.I CellStokesHypreSolver.h FACStokes.C FACStokes.h Makefile StokesSpecifications.C StokesSpecifications.I StokesSpecifications.h
description:
Copy over CellPoisson* stuff from SAMRAI internals and rename it to Stokes
diff -r cd1d77fa073e -r fe9d63509c19 CellStokesFACOps.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/CellStokesFACOps.C Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,2830 @@
+/*************************************************************************
+ *
+ * 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_CellStokesFACOps_C
+#define included_solv_CellStokesFACOps_C
+
+#include "CellStokesFACOps.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 "CellStokesHypreSolver.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 "CellStokesFACOps.I"
+#endif
+
+namespace SAMRAI {
+namespace solv {
+
+tbox::Pointer<pdat::CellVariable<double> >
+CellStokesFACOps::s_cell_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+
+tbox::Pointer<pdat::SideVariable<double> >
+CellStokesFACOps::s_flux_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+
+tbox::Pointer<pdat::OutersideVariable<double> >
+CellStokesFACOps::s_oflux_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+
+tbox::StartupShutdownManager::Handler
+CellStokesFACOps::s_finalize_handler(
+ 0,
+ 0,
+ 0,
+ CellStokesFACOps::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. *
+ ********************************************************************
+ */
+CellStokesFACOps::CellStokesFACOps(
+ 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::CellStokesFACOps::restrictSolution()");
+ t_restrict_residual = tbox::TimerManager::getManager()->
+ getTimer("solv::CellStokesFACOps::restrictResidual()");
+ t_prolong = tbox::TimerManager::getManager()->
+ getTimer("solv::CellStokesFACOps::prolongErrorAndCorrect()");
+ t_smooth_error = tbox::TimerManager::getManager()->
+ getTimer("solv::CellStokesFACOps::smoothError()");
+ t_solve_coarsest = tbox::TimerManager::getManager()->
+ getTimer("solv::CellStokesFACOps::solveCoarsestLevel()");
+ t_compute_composite_residual = tbox::TimerManager::getManager()->
+ getTimer("solv::CellStokesFACOps::computeCompositeResidualOnLevel()");
+ t_compute_residual_norm = tbox::TimerManager::getManager()->
+ getTimer("solv::CellStokesFACOps::computeResidualNorm()");
+
+ if (d_dim == tbox::Dimension(1) || d_dim > tbox::Dimension(3)) {
+ TBOX_ERROR("CellStokesFACOps : 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 << "CellStokesFACOps::private_cell_scratch" << dim.getValue();
+ s_cell_scratch_var[dim.getValue() - 1] = new pdat::CellVariable<double>
+ (dim, ss.str());
+ ss.str("");
+ ss << "CellStokesFACOps::private_flux_scratch" << dim.getValue();
+ s_flux_scratch_var[dim.getValue() - 1] = new pdat::SideVariable<double>
+ (dim, ss.str());
+ ss.str("");
+ ss << "CellStokesFACOps::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();
+
+}
+
+CellStokesFACOps::~CellStokesFACOps(
+ void)
+{
+}
+
+/*
+ ************************************************************************
+ * FACOperatorStrategy virtual initializeOperatorState function. *
+ * *
+ * Set internal variables to correspond to the solution passed in. *
+ * Look up transfer operators. *
+ ************************************************************************
+ */
+
+void CellStokesFACOps::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"
+ <<
+ "CellStokesFACOps::initializeOperatorState\n"
+ << "You must use "
+ <<
+ "CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps.");
+ }
+
+ t_smooth_error->stop();
+}
+
+/*
+ ********************************************************************
+ * Workhorse function to smooth error using red-black *
+ * Gauss-Seidel iterations. *
+ ********************************************************************
+ */
+
+void CellStokesFACOps::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 CellStokesFACOps::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("CellStokesFACOps : 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 CellStokesFACOps::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 "
+ << "scapCellStokesOps::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 scapCellStokesOps::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 CellStokesFACOps::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 "
+ << "CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps::computeResidualNorm(
+ const SAMRAIVectorReal<double>& residual,
+ int fine_ln,
+ int coarse_ln)
+{
+
+ if (coarse_ln != residual.getCoarsestLevelNumber() ||
+ fine_ln != residual.getFinestLevelNumber()) {
+ TBOX_ERROR("CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps::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 CellStokesFACOps::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
+CellStokesFACOps::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
+CellStokesFACOps::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
+CellStokesFACOps::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
+CellStokesFACOps::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
+CellStokesFACOps::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
+CellStokesFACOps::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
+CellStokesFACOps::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 cd1d77fa073e -r fe9d63509c19 CellStokesFACOps.I
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/CellStokesFACOps.I Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,213 @@
+/*************************************************************************
+ *
+ * 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 solving scalar Stokes using FAC
+ *
+ ************************************************************************/
+namespace SAMRAI {
+namespace solv {
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setPreconditioner(
+ const FACPreconditioner* preconditioner) {
+ d_preconditioner = preconditioner;
+}
+
+#ifdef HAVE_HYPRE
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setUseSMG(
+ bool use_smg)
+{
+ if (d_hierarchy) {
+ TBOX_ERROR(
+ d_object_name << ": setUseSMG(bool) may NOT be called\n"
+ <<
+ "while the solver state is initialized, as that\n"
+ << "would lead to a corrupted solver state.\n");
+ }
+ d_hypre_solver.setUseSMG(use_smg);
+}
+#endif
+
+/*
+ ********************************************************************
+ * Set the physical boundary condition object. *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setPhysicalBcCoefObject(
+ const RobinBcCoefStrategy* physical_bc_coef)
+{
+ d_physical_bc_coef = physical_bc_coef;
+ d_bc_helper.setCoefImplementation(physical_bc_coef);
+#ifdef HAVE_HYPRE
+ d_hypre_solver.setPhysicalBcCoefObject(d_physical_bc_coef);
+#endif
+}
+
+/*
+ ********************************************************************
+ * Set the object specifying the parameters of the Stokes equation *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setStokesSpecifications(
+ const StokesSpecifications& spec)
+{
+ d_stokes_spec = spec;
+}
+
+/*
+ ********************************************************************
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::enableLogging(
+ bool enable_logging)
+{
+ d_enable_logging = enable_logging;
+}
+
+/*
+ ********************************************************************
+ * Set the patch data id for the flux. *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setFluxId(
+ int flux_id) {
+ d_flux_id = flux_id;
+#ifdef DEBUG_CHECK_ASSERTIONS
+ checkInputPatchDataIndices();
+#endif
+}
+
+/*
+ ********************************************************************
+ * Set the choice for smoothing algorithm. *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setSmoothingChoice(
+ const std::string& smoothing_choice)
+{
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (smoothing_choice != "redblack") {
+ TBOX_ERROR(d_object_name << ": Bad smoothing choice '"
+ << smoothing_choice
+ << "' in CellStokesFACOps::setSmoothingChoice.");
+ }
+#endif
+ d_smoothing_choice = smoothing_choice;
+}
+
+/*
+ ********************************************************************
+ * Set the choice for the coarse level solver. *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setCoarsestLevelSolverChoice(
+ const std::string& choice) {
+#ifdef DEBUG_CHECK_ASSERTIONS
+#ifndef HAVE_HYPRE
+ if (choice == "hypre") {
+ TBOX_ERROR(d_object_name << ": HYPRe library is not available.\n");
+ }
+#endif
+#endif
+ if (choice == "redblack"
+ || choice == "hypre") {
+ d_coarse_solver_choice = choice;
+ } else {
+ TBOX_ERROR(
+ d_object_name << ": Bad coarse level solver choice '"
+ << choice
+ <<
+ "' in scapCellStokesOpsX::setCoarseLevelSolver.");
+ }
+}
+
+/*
+ ********************************************************************
+ * Set the tolerance for the coarse level solver. *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setCoarsestLevelSolverTolerance(
+ double tol) {
+ d_coarse_solver_tolerance = tol;
+}
+
+/*
+ ********************************************************************
+ * Set the tolerance for the coarse level solver. *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setCoarsestLevelSolverMaxIterations(
+ int max_iterations) {
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (max_iterations < 0) {
+ TBOX_ERROR(d_object_name << ": Invalid number of max iterations\n");
+ }
+#endif
+ d_coarse_solver_max_iterations = max_iterations;
+}
+
+/*
+ ********************************************************************
+ * Set the coarse-fine discretization method. *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setCoarseFineDiscretization(
+ const std::string& coarsefine_method) {
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (d_hierarchy) {
+ TBOX_ERROR(
+ d_object_name << ": Cannot change coarse-fine\n"
+ <<
+ "discretization method while operator state\n"
+ << "is initialized because that causes a\n"
+ << "corruption in the state.\n");
+ }
+#endif
+ d_cf_discretization = coarsefine_method;
+}
+
+/*
+ ********************************************************************
+ * Set the prolongation method *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACOps::setProlongationMethod(
+ const std::string& prolongation_method) {
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (d_hierarchy) {
+ TBOX_ERROR(
+ d_object_name << ": Cannot change prolongation method\n"
+ <<
+ "while operator state is initialized because that\n"
+ << "causes a corruption in the state.\n");
+ }
+#endif
+ d_prolongation_method = prolongation_method;
+}
+
+}
+}
diff -r cd1d77fa073e -r fe9d63509c19 CellStokesFACOps.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/CellStokesFACOps.h Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,1044 @@
+/*************************************************************************
+ *
+ * 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_CellStokesFACOps
+#define included_solv_CellStokesFACOps
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#include "SAMRAI/solv/CartesianRobinBcHelper.h"
+#include "SAMRAI/solv/FACPreconditioner.h"
+#include "SAMRAI/solv/FACOperatorStrategy.h"
+#include "SAMRAI/solv/RobinBcCoefStrategy.h"
+#include "CellStokesHypreSolver.h"
+#include "SAMRAI/solv/SAMRAIVectorReal.h"
+#include "StokesSpecifications.h"
+#include "SAMRAI/math/HierarchyCellDataOpsReal.h"
+#include "SAMRAI/math/HierarchySideDataOpsReal.h"
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/pdat/CellVariable.h"
+#include "SAMRAI/pdat/CellDoubleConstantRefine.h"
+#include "SAMRAI/pdat/OutersideData.h"
+#include "SAMRAI/pdat/OutersideVariable.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/xfer/CoarsenSchedule.h"
+#include "SAMRAI/xfer/RefineSchedule.h"
+#include "SAMRAI/xfer/CoarsenAlgorithm.h"
+#include "SAMRAI/xfer/CoarsenOperator.h"
+#include "SAMRAI/xfer/RefineAlgorithm.h"
+#include "SAMRAI/xfer/RefineOperator.h"
+#include "SAMRAI/hier/CoarseFineBoundary.h"
+#include "SAMRAI/hier/Patch.h"
+#include "SAMRAI/hier/PatchHierarchy.h"
+#include "SAMRAI/hier/PatchLevel.h"
+#include "SAMRAI/hier/IntVector.h"
+#include "SAMRAI/hier/Box.h"
+#include "SAMRAI/hier/VariableContext.h"
+#include "SAMRAI/tbox/Database.h"
+#include "SAMRAI/tbox/Pointer.h"
+#include "SAMRAI/tbox/Timer.h"
+
+#include <string>
+
+namespace SAMRAI {
+namespace solv {
+
+/*!
+ * @brief FAC operator class to solve Stokes's equation on a SAMR grid,
+ * using cell-centered, second-order finite-volume method, with Robin
+ * boundary conditions.
+ *
+ * This class provides operators that are used by the FAC
+ * preconditioner FACPreconditioner.
+ * It is used to solve the scalar Stokes's equation using a cell-centered
+ * second-order finite-volume discretization.
+ * It is designed to provide all operations specific to
+ * the scalar Stokes's equation,
+ * @f[ \nabla \cdot D \nabla u + C u = f @f]
+ * (see StokesSpecifications) where
+ * - C, D and f are indpendent of u
+ * - C is a cell-centered scalar field
+ * - D is the @em diffusion @em coefficients, stored on faces
+ * - f is a cell-centered scalar function
+ *
+ * You are left to provide the source function, initial guess, etc.,
+ * by specifying them in specific forms.
+ *
+ * This class provides:
+ * -# 5-point (second order), cell-centered stencil operations
+ * for the discrete Laplacian.
+ * -# Red-black Gauss-Seidel smoothing.
+ * -# Provisions for working Robin boundary conditions
+ * (see RobinBcCoefStrategy).
+ *
+ * This class is meant to provide the Stokes-specific operator
+ * used by the FAC preconditioner, FACPreconditioner.
+ * To use the preconditioner with this class, you will have to provide:
+ * -# The solution vector SAMRAIVectorReal,
+ * with appropriate norm weighting for the cell-centered AMR mesh.
+ * This class provides the function computeVectorWeights()
+ * to help with computing the appropriate weights.
+ * Since this is for a scalar equation, only the first depth
+ * of the first component of the vectors are used.
+ * All other parts are ignored.
+ * -# The source vector SAMRAIVectorReal for f.
+ * -# A StokesSpecifications objects to specify
+ * the cell-centered scalar field C and the side-centered
+ * diffusion coefficients D
+ * -# The boundary condition specifications in terms of the coefficients
+ * @f$ \alpha @f$, @f$ \beta @f$ and @f$ \gamma @f$ in the
+ * Robin formula @f$ \alpha u + \beta u_n = \gamma @f$ applied on the
+ * boundary faces. See RobinBcCoefStrategy.
+ *
+ * This class allocates and deallocates only its own scratch data.
+ * Other data that it manipuates are passed in as function arguments.
+ * Hence, it owns none of the solution vectors, error vectors,
+ * diffusion coefficient data, or any such things.
+ *
+ * Input Examples
+ * @verbatim
+ * coarse_solver_choice = "hypre" // see setCoarsestLevelSolverChoice()
+ * coarse_solver_tolerance = 1e-14 // see setCoarsestLevelSolverTolerance()
+ * coarse_solver_max_iterations = 10 // see setCoarsestLevelSolverMaxIterations()
+ * smoothing_choice = "redblack" // see setSmoothingChoice()
+ * cf_discretization = "Ewing" // see setCoarseFineDiscretization()
+ * prolongation_method = "LINEAR_REFINE" // see setProlongationMethod()
+ * hypre_solver = { ... } // tbox::Database for initializing Hypre solver
+ * @endverbatim
+ */
+class CellStokesFACOps:
+ public FACOperatorStrategy
+{
+
+public:
+ /*!
+ * @brief Constructor.
+ *
+ * If you want standard output and logging,
+ * pass in valid pointers for those streams.
+ * @param object_name Ojbect name
+ * @param database Input database
+ */
+ CellStokesFACOps(
+ const tbox::Dimension& dim,
+ const std::string& object_name = std::string(),
+ tbox::Pointer<tbox::Database> database =
+ tbox::Pointer<tbox::Database>(NULL));
+
+ /*!
+ * @brief Destructor.
+ *
+ * Deallocate internal data.
+ */
+ ~CellStokesFACOps(
+ void);
+
+ /*!
+ * @brief Set the scalar Stokes equation specifications.
+ */
+ void
+ setStokesSpecifications(
+ const StokesSpecifications& spec);
+
+ /*!
+ * @brief Enable logging.
+ *
+ * By default, logging is disabled. The logging flag is
+ * propagated to the major components used by this class.
+ */
+ void
+ enableLogging(
+ bool enable_logging);
+
+ //@{
+ /*!
+ * @name Functions for setting solver mathematic algorithm controls
+ */
+
+ /*!
+ * @brief Set the choice of smoothing algorithms.
+ *
+ * Current smoothing choices are:
+ * - "redblack": Red-black Gauss-Seidel smoothing.
+ */
+ void
+ setSmoothingChoice(
+ const std::string& smoothing_choice);
+
+ /*!
+ * @brief Set coarse level solver.
+ *
+ * Select from these:
+ * - @c "redblack" (red-black smoothing until convergence--very slow!)
+ * - @c "hypre" (only if the HYPRE library is available).
+ */
+ void
+ setCoarsestLevelSolverChoice(
+ const std::string& choice);
+
+ /*!
+ * @brief Set tolerance for coarse level solve.
+ *
+ * If the coarse level solver requires a tolerance (currently, they all do),
+ * the specified value is used.
+ */
+ void
+ setCoarsestLevelSolverTolerance(
+ double tol);
+
+ /*!
+ * @brief Set max iterations for coarse level solve.
+ *
+ * If the coarse level solver requires a max iteration limit
+ * (currently, they all do), the specified value is used.
+ */
+ void
+ setCoarsestLevelSolverMaxIterations(
+ int max_iterations);
+
+ /*!
+ * @brief Set the coarse-fine boundary discretization method.
+ *
+ * Specify the @c op_name std::string which will be passed to
+ * xfer::Geometry::lookupRefineOperator() to get the operator
+ * for setting fine grid ghost cells from the coarse grid.
+ * Note that chosing this operator implicitly choses the
+ * discretization method at the coarse-fine boundary.
+ *
+ * There is one important instance where this std::string is
+ * @em not passed to xfer::Geometry::lookupRefineOperator.
+ * If this variable is set to "Ewing", Ewing's coarse-fine
+ * discretization is used (a constant refinement is performed,
+ * and the flux is later corrected to result in Ewing's scheme).
+ * For a reference to Ewing's discretization method, see
+ * "Local Refinement Techniques for Elliptic Problems on Cell-Centered
+ * Grids, I. Error Analysis", Mathematics of Computation, Vol. 56, No. 194,
+ * April 1991, pp. 437-461.
+ *
+ * @param coarsefine_method String selecting the coarse-fine discretization method.
+ */
+ void
+ setCoarseFineDiscretization(
+ const std::string& coarsefine_method);
+
+ /*!
+ * @brief Set the name of the prolongation method.
+ *
+ * Specify the @c op_name std::string which will be passed to
+ * xfer::Geometry::lookupRefineOperator() to get the operator
+ * for prolonging the coarse-grid correction.
+ *
+ * By default, "CONSTANT_REFINE" is used. "LINEAR_REFINE" seems to
+ * to lead to faster convergence, but it does NOT satisfy the Galerkin
+ * condition.
+ *
+ * Prolonging using linear refinement requires a Robin bc
+ * coefficient implementation that is capable of delivering
+ * coefficients for non-hierarchy data, because linear refinement
+ * requires boundary conditions to be set on temporary levels.
+ *
+ * @param prolongation_method String selecting the coarse-fine
+ * discretization method.
+ */
+ void
+ setProlongationMethod(
+ const std::string& prolongation_method);
+
+#ifdef HAVE_HYPRE
+ /*!
+ * @brief Set whether to use Hypre's PFMG algorithm instead of the
+ * SMG algorithm.
+ *
+ * This flag affects the Hypre solver (used to solve the coarsest level).
+ * The flag is used to select which of HYPRE's linear solver algorithms
+ * to use if true, the semicoarsening multigrid algorithm is used, and if
+ * false, the ``PF'' multigrid algorithm is used.
+ * By default, the SMG algorithm is used.
+ *
+ * This setting has effect only when Hypre is chosen for the coarsest
+ * level solver. See setCoarsestLevelSolverChoice().
+ *
+ * Changing the algorithm must be done before initializing the solver
+ * state and must NOT be done while the state is initialized
+ * (the program will exit), as that would corrupt the state.
+ */
+ void
+ setUseSMG(
+ bool use_smg);
+#endif
+
+ //@}
+
+ //@{
+ /*!
+ * @name Functions for setting patch data indices and coefficients
+ */
+
+ /*!
+ * @brief Set the scratch patch data index for the flux.
+ *
+ * The use of this function is optional.
+ * The patch data index should be a pdat::SideData<DIM> type of variable.
+ * If the flux id is -1 (the default initial value), scratch space
+ * for the flux is allocated as needed and immediately deallocated
+ * afterward, level by level. If you have space preallocated for
+ * flux and you would like that to be used, set flux id to the
+ * patch data index of that space.
+ */
+ void
+ setFluxId(
+ int flux_id);
+
+ //@}
+
+ /*!
+ * @brief Provide an implementation for getting the
+ * physical bc coefficients
+ *
+ * If your solution is fixed at the physical boundary
+ * ghost cell centers AND those cells have the correct
+ * values before entering solveSystem(), you may use a
+ * GhostCellRobinBcCoefs object.
+ *
+ * If your solution is @b not fixed at the ghost cell centers,
+ * the ghost cell values will change as the interior
+ * cell values change. In those cases, the flexible
+ * Robin boundary conditions are applied. You must
+ * call this function to provide the implementation for
+ * determining the boundary condition coefficients.
+ *
+ * @param physical_bc_coef tbox::Pointer to an object that can
+ * set the Robin bc coefficients.
+ */
+ void
+ setPhysicalBcCoefObject(
+ const RobinBcCoefStrategy* physical_bc_coef);
+
+ //@{
+
+ /*!
+ * @name Functions for checking validity and correctness of state.
+ */
+
+ /*!
+ * @brief Check validity and correctness of input patch data indices.
+ *
+ * Descriptors checked:
+ * -# Diffusion coefficient (see setDiffcoefId())
+ * -# Flux (see setFluxId())
+ * -# Source (see setScalarFieldId())
+ */
+ void
+ checkInputPatchDataIndices() const;
+
+ //@}
+
+ /*!
+ * @brief Set weight appropriate for computing vector norms.
+ *
+ * If you this function to set the weights used when you
+ * SAMRAIVectorReal::addComponent, you can use the
+ * vector norm functions of SAMRAIVectorReal, and
+ * the weights will be used to blank out coarse grid
+ * regions under fine grids.
+ *
+ * The weights computed are specific to the cell-centered
+ * discretization used by this class. The weight is equal
+ * to the cell volume if the cell has not been refined,
+ * and zero if it has.
+ *
+ * This function is state-independent. All inputs are in
+ * the argument list.
+ *
+ * @param hierarchy Hierarchy configuration to compute weights for
+ * @param weight_id hier::Patch data index of the weight
+ * @param coarsest_ln Coarsest level number. Must be included
+ * in hierarchy. Must not be greater than @c finest_ln.
+ * Default to 0.
+ * @param finest_ln Finest level number. Must be included
+ * in hierarchy. Must not be less than @c coarsest_ln.
+ * Default to finest level in @c hierarchy.
+ */
+ void
+ computeVectorWeights(
+ tbox::Pointer<hier::PatchHierarchy> hierarchy,
+ int weight_id,
+ int coarsest_ln = -1,
+ int finest_ln = -1) const;
+
+ /*!
+ * @brief Set the FAC preconditioner that will be using this object.
+ *
+ * The FAC preconditioner is accessed to get convergence data during
+ * the cycle postprocessing step. It is optional.
+ */
+ void
+ setPreconditioner(
+ const FACPreconditioner* preconditioner);
+
+ /*!
+ * @brief function to compute flux, using general diffusion
+ * coefficient data.
+ *
+ * Recall that this solver class discretizes the PDE
+ * @f[ \nabla \cdot D \nabla u + C u = f @f] on an AMR grid. This member
+ * function allows users of this solver class to compute gradient
+ * terms, @f[ D \nabla w @f], in their code in a manner consistent with the
+ * solver discretization. In particular, when solving PDE systems, it may
+ * be necessary to discretize the gradient operator appearing in equations
+ * not treated by the solver class in the same way as those treated by this
+ * class. These funtions allow users to do this easily. The divergence
+ * operator used in this solver is the standard sum of centered differences
+ * involving flux terms on the cell sides computed by these routines.
+ *
+ * Note that the patch must exist on a level in an AMR hierarchy so that
+ * the discretization can be computed properly at the coarse-fine interface.
+ * Stokes coefficients C and D must exist on the patch, if they are variable.
+ * Also, calling this function does not affect the internal solver state in any
+ * way. However, the solver must be fully initialized before it is called and care
+ * should be exercised to pass arguments so that the solver solution quantity and
+ * other internal solver quantities are not adversely affected.
+ *
+ * @param patch patch on which computation will take place
+ * @param ratio_to_coarser_level refinement ratio from coarser level to level
+ * on which patch lives; if current patch level
+ * is level zero, this is ignored
+ * @param w_data cell-centered data
+ * @param Dgradw_data side-centered flux data (i.e., D (grad w))
+ */
+ void
+ computeFluxOnPatch(
+ const hier::Patch& patch,
+ const hier::IntVector& ratio_to_coarser_level,
+ const pdat::CellData<double>& w_data,
+ pdat::SideData<double>& Dgradw_data) const;
+
+ //@{ @name FACOperatorStrategy virtuals
+
+ virtual void
+ restrictSolution(
+ const SAMRAIVectorReal<double>& source,
+ SAMRAIVectorReal<double>& dest,
+ int dest_ln);
+ virtual void
+ restrictResidual(
+ const SAMRAIVectorReal<double>& source,
+ SAMRAIVectorReal<double>& dest,
+ int dest_ln);
+
+ virtual void
+ prolongErrorAndCorrect(
+ const SAMRAIVectorReal<double>& source,
+ SAMRAIVectorReal<double>& dest,
+ int dest_ln);
+
+ virtual void
+ smoothError(
+ SAMRAIVectorReal<double>& error,
+ const SAMRAIVectorReal<double>& residual,
+ int ln,
+ int num_sweeps);
+
+ virtual int
+ solveCoarsestLevel(
+ SAMRAIVectorReal<double>& error,
+ const SAMRAIVectorReal<double>& residual,
+ int coarsest_ln);
+
+ virtual void
+ computeCompositeResidualOnLevel(
+ SAMRAIVectorReal<double>& residual,
+ const SAMRAIVectorReal<double>& solution,
+ const SAMRAIVectorReal<double>& rhs,
+ int ln,
+ bool error_equation_indicator);
+
+ virtual double
+ computeResidualNorm(
+ const SAMRAIVectorReal<double>& residual,
+ int fine_ln,
+ int coarse_ln);
+
+ virtual void
+ initializeOperatorState(
+ const SAMRAIVectorReal<double>& solution,
+ const SAMRAIVectorReal<double>& rhs);
+
+ virtual void
+ deallocateOperatorState();
+
+ virtual void
+ postprocessOneCycle(
+ int fac_cycle_num,
+ const SAMRAIVectorReal<double>& current_soln,
+ const SAMRAIVectorReal<double>& residual);
+
+ //@}
+
+private:
+ //@{
+ /*!
+ * @name Private workhorse functions.
+ */
+
+ /*!
+ * @brief Red-black Gauss-Seidel error smoothing on a level.
+ *
+ * Smoothes on the residual equation @f$ Ae=r @f$ on a level.
+ *
+ * @param error error vector
+ * @param residual residual vector
+ * @param ln level number
+ * @param num_sweeps number of sweeps
+ * @param residual_tolerance the maximum residual considered to be
+ * converged
+ */
+ void
+ smoothErrorByRedBlack(
+ SAMRAIVectorReal<double>& error,
+ const SAMRAIVectorReal<double>& residual,
+ int ln,
+ int num_sweeps,
+ double residual_tolerance = -1.0);
+
+ /*!
+ * @brief Solve the coarsest level using HYPRE
+ */
+ int
+ solveCoarsestLevel_HYPRE(
+ SAMRAIVectorReal<double>& error,
+ const SAMRAIVectorReal<double>& residual,
+ int ln);
+
+ /*!
+ * @brief Fix flux per Ewing's coarse-fine boundary treatment.
+ *
+ * Ewing's coarse-fine boundary treatment can be implemented
+ * using a constant refinement into the fine-grid ghost boundary,
+ * naively computing the flux using the constant-refined data then
+ * fixing up the flux to correct the error.
+ *
+ * To use this function
+ * -# you must use constant refinement to fill the fine level ghost cells
+ * -# the flux must first be computed and stored
+ *
+ * @param patch patch
+ * @param soln_data cell-centered solution data
+ * @param flux_data side-centered flux data
+ * @param diffcoef_data side-centered diffusion coefficient data
+ * @param cfb coarse-fine boundary object for the level
+ * in which patch resides
+ * @param ratio_to_coarser Refinement ratio to the next coarser level.
+ */
+ void
+ ewingFixFlux(
+ const hier::Patch& patch,
+ const pdat::CellData<double>& soln_data,
+ pdat::SideData<double>& flux_data,
+ const hier::IntVector& ratio_to_coarser) const;
+
+ /*!
+ * @brief AMR-unaware function to compute residual on a single patch,
+ * with variable scalar field.
+ *
+ * @param patch patch
+ * @param flux_data side-centered flux data
+ * @param soln_data cell-centered solution data
+ * @param rhs_data cell-centered rhs data
+ * @param residual_data cell-centered residual data
+ */
+ void
+ 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;
+
+ /*!
+ * @brief AMR-unaware function to red or black smoothing on a single patch,
+ * for variable diffusion coefficient and variable scalar field.
+ *
+ * @param patch patch
+ * @param flux_data side-centered flux data
+ * @param rhs_data cell-centered rhs data
+ * @param scalar_field_data
+ * cell-centered scalar field data
+ * @param soln_data cell-centered solution data
+ * @param red_or_black red-black switch. Set to 'r' or 'b'.
+ * @param p_maxres max residual output. Set to NULL to avoid computing.
+ */
+ void
+ 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 = NULL) const;
+
+ //@}
+
+ //@{ @name For executing, caching and resetting communication schedules.
+
+ /*!
+ * @brief Execute a refinement schedule
+ * for prolonging cell data.
+ *
+ * General notes regarding internal objects for communication:
+ * We maintain objects to support caching schedules to improve
+ * efficiency. Communication is needed in 5 distinct tasks.
+ * -# Prolongation
+ * -# Restriction
+ * -# Flux coarsening. Changing the coarse grid flux to the
+ * composite grid flux by coarsening the fine grid flux
+ * at the coarse-fine boundaries.
+ * -# Fill boundary data from other patches in the same level
+ * and physical boundary condition.
+ * -# Fill boundary data from same level, coarser levels
+ * and physical boundary condition.
+ *
+ * For each task, we maintain a refine or coarsen operator,
+ * and a array of communication schedules (one for each
+ * destination level).
+ *
+ * The 5 member functions named @c xeqSchedule... execute
+ * communication schedules appropriate for five specific tasks.
+ * They use a cached schedule if possible or create and cache
+ * a new schedule if needed. These functions and the data
+ * they manipulate are as follows:
+ * <ol>
+ * <li> xeqScheduleProlongation():
+ * d_prolongation_refine_operator
+ * d_prolongation_refine_schedules
+ * <li> xeqScheduleURestriction():
+ * d_restriction_coarsen_operator,
+ * d_urestriction_coarsen_schedules.
+ * <li> xeqScheduleRRestriction():
+ * d_restriction_coarsen_operator,
+ * d_rrestriction_coarsen_schedules.
+ * <li> xeqScheduleFluxCoarsen():
+ * d_flux_coarsen_operator,
+ * d_flux_coarsen_schedules.
+ * <li> xeqScheduleGhostFill():
+ * d_ghostfill_refine_operator,
+ * d_ghostfill_refine_schedules.
+ * <li> xeqScheduleGhostFillNoCoarse():
+ * d_ghostfill_nocoarse_refine_operator,
+ * d_ghostfill_nocoarse_refine_schedules.
+ * </ol>
+ *
+ * @return refinement schedule for prolongation
+ */
+ void
+ xeqScheduleProlongation(
+ int dst_id,
+ int src_id,
+ int scr_id,
+ int dest_ln);
+
+ /*!
+ * @brief Execute schedule for restricting solution to the specified
+ * level or reregister an existing one.
+ *
+ * See general notes for xeqScheduleProlongation().
+ *
+ * @return coarsening schedule for restriction
+ */
+ void
+ xeqScheduleURestriction(
+ int dst_id,
+ int src_id,
+ int dest_ln);
+
+ /*!
+ * @brief Execute schedule for restricting residual to the specified
+ * level or reregister an existing one.
+ *
+ * See general notes for xeqScheduleProlongation().
+ *
+ * @return coarsening schedule for restriction
+ */
+ void
+ xeqScheduleRRestriction(
+ int dst_id,
+ int src_id,
+ int dest_ln);
+
+ /*!
+ * @brief Execute schedule for coarsening flux to the specified
+ * level or reregister an existing one.
+ *
+ * See general notes for xeqScheduleProlongation().
+ *
+ * @return coarsening schedule for setting composite grid flux at
+ * coarse-fine boundaries.
+ */
+ void
+ xeqScheduleFluxCoarsen(
+ int dst_id,
+ int src_id,
+ int dest_ln);
+
+ /*!
+ * @brief Execute schedule for filling ghosts on the specified
+ * level or reregister an existing one.
+ *
+ * See general notes for xeqScheduleProlongation().
+ *
+ * @return refine schedule for filling ghost data from coarser level
+ * and physical bc.
+ */
+ void
+ xeqScheduleGhostFill(
+ int dst_id,
+ int dest_ln);
+
+ /*!
+ * @brief Execute schedule for filling ghosts on the specified
+ * level or reregister an existing one.
+ * This version does not get data from coarser levels.
+ *
+ * See general notes for xeqScheduleProlongation().
+ *
+ * This function is used for the bottom solve level, since it does
+ * not access data from any coarser level. (Ghost data obtained
+ * from coarser level must have been placed there before solve begins!)
+ *
+ * @return refine schedule for filling ghost data from same level
+ * and physical bc.
+ */
+ void
+ xeqScheduleGhostFillNoCoarse(
+ int dst_id,
+ int dest_ln);
+
+ //@}
+
+ //! @brief Return the patch data index for cell scratch data.
+ int
+ registerCellScratch() const;
+ //! @brief Return the patch data index for flux scratch data.
+ int
+ registerFluxScratch() const;
+ //! @brief Return the patch data index for outerflux scratch data.
+ int
+ registerOfluxScratch() const;
+
+ //! @brief Free static variables at shutdown time.
+ static void
+ finalizeCallback();
+
+ /*!
+ * @brief Object dimension.
+ */
+ const tbox::Dimension d_dim;
+
+ /*!
+ * @brief Object name.
+ */
+ std::string d_object_name;
+
+ //@{ @name Hierarchy-dependent objects.
+
+ /*!
+ * @brief Reference hierarchy
+ *
+ * This variable is non-null between the initializeOperatorState()
+ * and deallocateOperatorState() calls. It is not truly needed,
+ * because the hierarchy is obtainable through variables in most
+ * function argument lists. We use it to enforce working on one
+ * hierarchy at a time.
+ */
+ tbox::Pointer<hier::PatchHierarchy> d_hierarchy;
+
+ /*!
+ * @brief Coarsest level for solve.
+ */
+ int d_ln_min;
+
+ /*!
+ * @brief Finest level for solve.
+ */
+ int d_ln_max;
+
+ /*!
+ * @brief Description of coarse-fine boundaries.
+ *
+ * There is one coarse-fine boundary object for each level.
+ * d_coarse_fine_boundary[i] is the description of
+ * the coarse-fine boundary between level i and level i-1.
+ * The coarse-fine boundary does not exist at the coarsest level,
+ * although the hier::CoarseFineBoundary object still exists (it
+ * should not contain any boxes).
+ *
+ * This array is initialized in initializeOperatorState() and
+ * deallocated in deallocateOperatorState(). When allocated,
+ * it is allocated for the index range [0,d_ln_max], though
+ * the range [0,d_ln_min-1] is not used. This is okay because
+ * hier::CoarseFineBoundary is a light object before
+ * it is set for a level.
+ */
+ tbox::Array<tbox::Pointer<hier::CoarseFineBoundary> > d_cf_boundary;
+
+ //@}
+
+ //@{
+ /*!
+ * @name Private state variables for solution process.
+ */
+
+ /*!
+ * @brief Scalar Stokes equations specifications.
+ * @see setStokesSpecifications().
+ */
+ StokesSpecifications d_stokes_spec;
+
+ /*!
+ * @brief Smoothing choice.
+ * @see setSmoothingChoice.
+ */
+ std::string d_smoothing_choice;
+
+ /*!
+ * @brief Coarse level solver.
+ * @see setCoarsestLevelSolverChoice
+ */
+ std::string d_coarse_solver_choice;
+
+ /*!
+ * @brief Coarse-fine discretization method.
+ * @see setCoarseFineDiscretization().
+ */
+ std::string d_cf_discretization;
+
+ /*!
+ * @brief Coarse-fine discretization method.
+ *
+ * The name of the refinement operator used to prolong the
+ * coarse grid correction.
+ *
+ * @see setProlongationMethod()
+ */
+ std::string d_prolongation_method;
+
+ /*!
+ * @brief Tolerance specified to coarse solver
+ * @see setCoarsestLevelSolverTolerance()
+ */
+ double d_coarse_solver_tolerance;
+
+ /*!
+ * @brief Coarse level solver iteration limit.
+ * @see setCoarsestLevelSolverMaxIterations()
+ */
+ int d_coarse_solver_max_iterations;
+
+ /*!
+ * @brief Residual tolerance to govern smoothing.
+ *
+ * When we use one of the internal error smoothing functions
+ * and want to terminate the smoothing sweeps at a certain
+ * level of residual, this will be set to > 0. If it is
+ * < 0, the smoothing function effectively ignores it.
+ *
+ * This variable is needed because some coarse-level solver
+ * simply runs the smoothing function until convergence.
+ * It sets this variable to > 0, calls the smoothing function,
+ * then resets it to < 0.
+ */
+ double d_residual_tolerance_during_smoothing;
+
+ /*!
+ * @brief Id of the flux.
+ *
+ * If set to -1, create and delete storage space on the fly.
+ * Else, user has provided space for flux.
+ *
+ * @see setFluxId
+ */
+ int d_flux_id;
+
+#ifdef HAVE_HYPRE
+ /*!
+ * @brief HYPRE coarse-level solver object.
+ */
+ CellStokesHypreSolver d_hypre_solver;
+#endif
+
+ /*!
+ * @brief Externally provided physical boundary condition object.
+ *
+ * see setPhysicalBcCoefObject()
+ */
+ const RobinBcCoefStrategy* d_physical_bc_coef;
+
+ //@}
+
+ //@{ @name Internal context and scratch data
+
+ static tbox::Pointer<pdat::CellVariable<double> >
+ s_cell_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+
+ static tbox::Pointer<pdat::SideVariable<double> >
+ s_flux_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+
+ static tbox::Pointer<pdat::OutersideVariable<double> >
+ s_oflux_scratch_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+
+ /*!
+ * @brief Default context of internally maintained hierarchy data.
+ */
+ tbox::Pointer<hier::VariableContext> d_context;
+
+ /*!
+ * @brief ID of the solution-like scratch data.
+ *
+ * Set in constructor and never changed.
+ * Corresponds to a pdat::CellVariable<double> named
+ * @c d_object_name+"::cell_scratch".
+ * Scratch data is allocated and removed as needed
+ * to reduce memory usage.
+ */
+ int d_cell_scratch_id;
+
+ /*!
+ * @brief ID of the side-centered scratch data.
+ *
+ * Set in constructor and never changed.
+ * Corresponds to a pdat::SideVariable<double> named
+ * @c d_object_name+"::flux_scratch".
+ *
+ * This data is allocated only as needed and deallocated
+ * immediately after use.
+ */
+ int d_flux_scratch_id;
+
+ /*!
+ * @brief ID of the outerside-centered scratch data.
+ *
+ * Set in constructor and never changed.
+ * Corresponds to a pdat::OutersideVariable<double> named
+ * @c d_object_name+"::oflux_scratch".
+ */
+ int d_oflux_scratch_id;
+
+ //@}
+
+ //@{
+ /*!
+ * @name Various refine and coarsen objects used internally.
+ */
+
+ //! @brief Error prolongation (refinement) operator.
+ tbox::Pointer<xfer::RefineOperator> d_prolongation_refine_operator;
+ tbox::Pointer<xfer::RefineAlgorithm> d_prolongation_refine_algorithm;
+ tbox::Array<tbox::Pointer<xfer::RefineSchedule> >
+ d_prolongation_refine_schedules;
+
+ //! @brief Solution restriction (coarsening) operator.
+ tbox::Pointer<xfer::CoarsenOperator> d_urestriction_coarsen_operator;
+ tbox::Pointer<xfer::CoarsenAlgorithm> d_urestriction_coarsen_algorithm;
+ tbox::Array<tbox::Pointer<xfer::CoarsenSchedule> >
+ d_urestriction_coarsen_schedules;
+
+ //! @brief Residual restriction (coarsening) operator.
+ tbox::Pointer<xfer::CoarsenOperator> d_rrestriction_coarsen_operator;
+ tbox::Pointer<xfer::CoarsenAlgorithm> d_rrestriction_coarsen_algorithm;
+ tbox::Array<tbox::Pointer<xfer::CoarsenSchedule> >
+ d_rrestriction_coarsen_schedules;
+
+ //! @brief Coarsen operator for outerflux-to-flux
+ tbox::Pointer<xfer::CoarsenOperator> d_flux_coarsen_operator;
+ tbox::Pointer<xfer::CoarsenAlgorithm> d_flux_coarsen_algorithm;
+ tbox::Array<tbox::Pointer<xfer::CoarsenSchedule> >
+ d_flux_coarsen_schedules;
+
+ //! @brief Refine operator for cell-like data from coarser level.
+ tbox::Pointer<xfer::RefineOperator> d_ghostfill_refine_operator;
+ tbox::Pointer<xfer::RefineAlgorithm> d_ghostfill_refine_algorithm;
+ tbox::Array<tbox::Pointer<xfer::RefineSchedule> >
+ d_ghostfill_refine_schedules;
+
+ //! @brief Refine operator for cell-like data from same level.
+ tbox::Pointer<xfer::RefineOperator> d_ghostfill_nocoarse_refine_operator;
+ tbox::Pointer<xfer::RefineAlgorithm> d_ghostfill_nocoarse_refine_algorithm;
+ tbox::Array<tbox::Pointer<xfer::RefineSchedule> >
+ d_ghostfill_nocoarse_refine_schedules;
+
+ //@}
+
+ /*!
+ * @brief Utility object employed in setting ghost cells and providing
+ * xfer::RefinePatchStrategy implementation.
+ *
+ * Since this class deals only in scalar variables having
+ * Robin boundary conditions, we take advantage of the corresponding
+ * implementation in CartesianRobinBcHelper. Whenever
+ * we need an implementation of xfer::RefinePatchStrategy,
+ * this object is used. Note that in the code, before we
+ * use this object to set ghost cell values, directly or
+ * indirectly by calling xfer::RefineSchedule::fillData(),
+ * we must tell d_bc_helper the patch data index we want
+ * to set and whether we are setting data with homogeneous
+ * boundary condition.
+ */
+ CartesianRobinBcHelper d_bc_helper;
+
+ //@{
+ /*!
+ * @name Non-essential objects used in outputs and debugging.
+ */
+
+ /*!
+ * @brief Logging flag.
+ */
+ bool d_enable_logging;
+
+ /*!
+ * @brief Preconditioner using this object.
+ *
+ * See setPreconditioner().
+ */
+ const FACPreconditioner* d_preconditioner;
+
+ /*!
+ * @brief Hierarchy cell operator used in debugging.
+ */
+ tbox::Pointer<math::HierarchyCellDataOpsReal<double> > d_hopscell;
+
+ /*!
+ * @brief Hierarchy side operator used in debugging.
+ */
+ tbox::Pointer<math::HierarchySideDataOpsReal<double> > d_hopsside;
+
+ /*!
+ * @brief Timers for performance measurement.
+ */
+ tbox::Pointer<tbox::Timer> t_restrict_solution;
+ tbox::Pointer<tbox::Timer> t_restrict_residual;
+ tbox::Pointer<tbox::Timer> t_prolong;
+ tbox::Pointer<tbox::Timer> t_smooth_error;
+ tbox::Pointer<tbox::Timer> t_solve_coarsest;
+ tbox::Pointer<tbox::Timer> t_compute_composite_residual;
+ tbox::Pointer<tbox::Timer> t_compute_residual_norm;
+
+ static tbox::StartupShutdownManager::Handler
+ s_finalize_handler;
+};
+
+}
+}
+
+#ifdef SAMRAI_INLINE
+#include "CellStokesFACOps.I"
+#endif
+
+#endif // included_solv_CellStokesFACOps
diff -r cd1d77fa073e -r fe9d63509c19 CellStokesFACSolver.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/CellStokesFACSolver.C Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,585 @@
+/*************************************************************************
+ *
+ * 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: High-level solver (wrapper) for scalar stokes equation.
+ *
+ ************************************************************************/
+#ifndef included_solv_CellStokesFACSolver_C
+#define included_solv_CellStokesFACSolver_C
+
+#include "SAMRAI/pdat/CellVariable.h"
+#include "CellStokesFACSolver.h"
+#include "SAMRAI/tbox/PIO.h"
+#include "SAMRAI/tbox/Utilities.h"
+#include "SAMRAI/tbox/StartupShutdownManager.h"
+
+#include IOMANIP_HEADER_FILE
+
+#ifndef SAMRAI_INLINE
+#include "CellStokesFACSolver.I"
+#endif
+
+namespace SAMRAI {
+namespace solv {
+
+/*
+ *************************************************************************
+ * *
+ * Initialize the static data members. *
+ * *
+ *************************************************************************
+ */
+
+bool CellStokesFACSolver::s_initialized = 0;
+int CellStokesFACSolver::s_weight_id[SAMRAI::tbox::Dimension::
+ MAXIMUM_DIMENSION_VALUE];
+int CellStokesFACSolver::s_instance_counter[SAMRAI::tbox::Dimension::
+ MAXIMUM_DIMENSION_VALUE];
+
+/*
+ *************************************************************************
+ * *
+ * Constructor sets uninitialized solver state. *
+ * Set default iteration and convergence parameters. *
+ * *
+ * By default settings: *
+ * - Stokes equation specified has D=1, C=0. *
+ * - State is uninitialized *
+ * - Logging is disabled *
+ * - Context for internal data is set based on object name. *
+ * *
+ *************************************************************************
+ */
+
+CellStokesFACSolver::CellStokesFACSolver(
+ const tbox::Dimension& dim,
+ const std::string& object_name,
+ tbox::Pointer<tbox::Database> database):
+ d_dim(dim),
+ d_object_name(object_name),
+ d_stokes_spec(object_name + "::stokes_spec"),
+ d_fac_ops(d_dim, object_name + "::fac_ops"),
+ d_fac_precond(object_name + "::fac_precond", d_fac_ops),
+ d_bc_object(NULL),
+ d_simple_bc(d_dim, object_name + "::bc"),
+ d_hierarchy(NULL),
+ d_ln_min(-1),
+ d_ln_max(-1),
+ d_context(hier::VariableDatabase::getDatabase()
+ ->getContext(object_name + "::CONTEXT")),
+ d_uv(NULL),
+ d_fv(NULL),
+ d_solver_is_initialized(false),
+ d_enable_logging(false)
+{
+
+ if (!s_initialized) {
+ initializeStatics();
+ }
+
+ setMaxCycles(10);
+ setResidualTolerance(1e-6);
+ setPresmoothingSweeps(1);
+ setPostsmoothingSweeps(1);
+ setCoarseFineDiscretization("Ewing");
+#ifdef HAVE_HYPRE
+ setCoarsestLevelSolverChoice("hypre");
+ setCoarsestLevelSolverTolerance(1e-10);
+ setCoarsestLevelSolverMaxIterations(20);
+ setUseSMG(true);
+#else
+ setCoarsestLevelSolverChoice("redblack");
+ setCoarsestLevelSolverTolerance(1e-8);
+ setCoarsestLevelSolverMaxIterations(500);
+#endif
+
+ /*
+ * Construct integer tag variables and add to variable database. Note that
+ * variables and patch data indices are shared among all instances.
+ * The VariableDatabase holds the variables, once contructed and
+ * registered via the VariableDatabase::registerInternalSAMRAIVariable()
+ * function call. Note that variables are registered and patch data indices
+ * are made only for the first time through the constructor.
+ */
+ hier::VariableDatabase* var_db = hier::VariableDatabase::getDatabase();
+
+ static std::string weight_variable_name("CellStokesFACSolver_weight");
+
+ tbox::Pointer<pdat::CellVariable<double> > weight = var_db->getVariable(
+ weight_variable_name);
+ if (weight.isNull()) {
+ weight = new pdat::CellVariable<double>(d_dim, weight_variable_name, 1);
+ }
+
+ if (s_weight_id[d_dim.getValue() - 1] < 0) {
+ s_weight_id[d_dim.getValue() - 1] = var_db->registerInternalSAMRAIVariable(
+ weight,
+ hier::IntVector::getZero(d_dim));
+ }
+
+ /*
+ * The default RobinBcCoefStrategy used,
+ * SimpleCellRobinBcCoefs only works with constant refine
+ * for prolongation. So we use constant refinement
+ * for prolongation by default.
+ */
+ setProlongationMethod("CONSTANT_REFINE");
+
+ /*
+ * The FAC operator optionally uses the preconditioner
+ * to get data for logging.
+ */
+ d_fac_ops.setPreconditioner((const FACPreconditioner *)(&d_fac_precond));
+
+ if (database) {
+ getFromInput(database);
+ }
+
+ s_instance_counter[d_dim.getValue() - 1]++;
+}
+
+/*
+ *************************************************************************
+ * *
+ * Destructor for CellStokesFACSolver. *
+ * Deallocate internal data. *
+ * *
+ *************************************************************************
+ */
+
+CellStokesFACSolver::~CellStokesFACSolver()
+{
+ s_instance_counter[d_dim.getValue() - 1]--;
+
+ deallocateSolverState();
+
+ if (s_instance_counter[d_dim.getValue() - 1] == 0) {
+ hier::VariableDatabase::getDatabase()->
+ removeInternalSAMRAIVariablePatchDataIndex(s_weight_id[d_dim.getValue() - 1]);
+ s_weight_id[d_dim.getValue() - 1] = -1;
+ }
+}
+
+/*
+ ********************************************************************
+ * Set state from database *
+ * *
+ * Do not allow FAC preconditioner and Stokes FAC operators to be *
+ * set from database, as that may cause them to be inconsistent *
+ * with this object if user does not coordinate the inputs *
+ * correctly. This is also why we don't allow direct access to *
+ * those objects. The responsibility for maintaining consistency *
+ * lies in the public functions to set parameters, so use them *
+ * instead of setting the parameters directly in this function. *
+ ********************************************************************
+ */
+
+void CellStokesFACSolver::getFromInput(
+ tbox::Pointer<tbox::Database> database)
+{
+ if (database) {
+ if (database->isBool("enable_logging")) {
+ bool logging = database->getBool("enable_logging");
+ enableLogging(logging);
+ }
+ if (database->isInteger("max_cycles")) {
+ int max_cycles = database->getInteger("max_cycles");
+ setMaxCycles(max_cycles);
+ }
+ if (database->isDouble("residual_tol")) {
+ double residual_tol = database->getDouble("residual_tol");
+ setResidualTolerance(residual_tol);
+ }
+ if (database->isInteger("num_pre_sweeps")) {
+ int num_pre_sweeps = database->getInteger("num_pre_sweeps");
+ setPresmoothingSweeps(num_pre_sweeps);
+ }
+ if (database->isInteger("num_post_sweeps")) {
+ int num_post_sweeps = database->getInteger("num_post_sweeps");
+ setPostsmoothingSweeps(num_post_sweeps);
+ }
+ if (database->isString("coarse_fine_discretization")) {
+ std::string s = database->getString("coarse_fine_discretization");
+ setCoarseFineDiscretization(s);
+ }
+ if (database->isString("prolongation_method")) {
+ std::string s = database->getString("prolongation_method");
+ setProlongationMethod(s);
+ }
+ if (database->isString("coarse_solver_choice")) {
+ std::string s = database->getString("coarse_solver_choice");
+ setCoarsestLevelSolverChoice(s);
+ }
+ if (database->isDouble("coarse_solver_tolerance")) {
+ double tol = database->getDouble("coarse_solver_tolerance");
+ setCoarsestLevelSolverTolerance(tol);
+ }
+ if (database->isInteger("coarse_solver_max_iterations")) {
+ int itr = database->getInteger("coarse_solver_max_iterations");
+ setCoarsestLevelSolverMaxIterations(itr);
+ }
+#ifdef HAVE_HYPRE
+ if (database->isBool("use_smg")) {
+ bool smg = database->getBool("use_smg");
+ setUseSMG(smg);
+ }
+#endif
+ }
+}
+
+/*
+ *************************************************************************
+ * *
+ * Prepare internal data for solve. *
+ * Allocate scratch data. Create vectors for u and f *
+ * required by the FACPreconditioner interface. *
+ * Set up internal boundary condition object. *
+ * Share data to coordinate with FAC preconditioner and *
+ * Stokes FAC operator. *
+ * *
+ *************************************************************************
+ */
+
+void CellStokesFACSolver::initializeSolverState(
+ const int solution,
+ const int rhs,
+ tbox::Pointer<hier::PatchHierarchy> hierarchy,
+ const int coarse_level,
+ const int fine_level)
+{
+ TBOX_ASSERT(!hierarchy.isNull());
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS1(d_dim, *hierarchy);
+
+ if (d_bc_object == NULL) {
+ TBOX_ERROR(
+ d_object_name << ": No BC coefficient strategy object!\n"
+ <<
+ "Use either setBoundaries or setPhysicalBcCoefObject\n"
+ << "to specify the boundary conidition.\n");
+ }
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (solution < 0 || rhs < 0) {
+ TBOX_ERROR(d_object_name << ": Bad patch data id.\n");
+ }
+#endif
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (!hierarchy) {
+ TBOX_ERROR(d_object_name << ": NULL hierarchy pointer not allowed\n"
+ << "in inititialization.");
+ }
+#endif
+ d_hierarchy = hierarchy;
+
+ d_ln_min = coarse_level;
+ d_ln_max = fine_level;
+ if (d_ln_min == -1) {
+ d_ln_min = 0;
+ }
+ if (d_ln_max == -1) {
+ d_ln_max = d_hierarchy->getFinestLevelNumber();
+ }
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (d_ln_min < 0 || d_ln_max < 0 || d_ln_min > d_ln_max) {
+ TBOX_ERROR(d_object_name << ": Bad range of levels in\n"
+ << "inititialization.\n");
+ }
+#endif
+
+ int ln;
+ for (ln = d_ln_min; ln <= d_ln_max; ++ln) {
+ d_hierarchy->getPatchLevel(ln)->allocatePatchData(s_weight_id[d_dim.getValue() - 1]);
+ }
+
+ d_fac_ops.computeVectorWeights(d_hierarchy,
+ s_weight_id[d_dim.getValue() - 1],
+ d_ln_min,
+ d_ln_max);
+
+ if (d_bc_object == &d_simple_bc) {
+ d_simple_bc.setHierarchy(d_hierarchy,
+ d_ln_min,
+ d_ln_max);
+ if (d_stokes_spec.dIsConstant()) {
+ d_simple_bc.setDiffusionCoefConstant(d_stokes_spec.getDConstant());
+ } else {
+ d_simple_bc.setDiffusionCoefId(d_stokes_spec.getDPatchDataId());
+ }
+ }
+
+ d_fac_ops.setStokesSpecifications(d_stokes_spec);
+
+ createVectorWrappers(solution, rhs);
+
+ d_fac_precond.initializeSolverState(*d_uv, *d_fv);
+
+ d_solver_is_initialized = true;
+}
+
+void CellStokesFACSolver::deallocateSolverState()
+{
+ if (d_hierarchy) {
+
+ d_fac_precond.deallocateSolverState();
+
+ /*
+ * Delete internally managed data.
+ */
+ int ln;
+ for (ln = d_ln_min; ln <= d_ln_max; ++ln) {
+ d_hierarchy->getPatchLevel(ln)->deallocatePatchData(s_weight_id[d_dim.getValue()
+ - 1]);
+ }
+
+ d_hierarchy.setNull();
+ d_ln_min = -1;
+ d_ln_max = -1;
+ d_solver_is_initialized = false;
+
+ destroyVectorWrappers();
+
+ }
+}
+
+/*
+ *************************************************************************
+ * Enable logging and propagate logging flag to major components. *
+ *************************************************************************
+ */
+
+void CellStokesFACSolver::enableLogging(
+ bool logging)
+{
+ d_enable_logging = logging;
+ d_fac_precond.enableLogging(d_enable_logging);
+ d_fac_ops.enableLogging(d_enable_logging);
+}
+
+void CellStokesFACSolver::setBoundaries(
+ const std::string& boundary_type,
+ const int fluxes,
+ const int flags,
+ int* bdry_types)
+{
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (d_bc_object != NULL && d_bc_object != &d_simple_bc) {
+ TBOX_ERROR(
+ d_object_name << ": Bad attempt to set boundary condition\n"
+ <<
+ "by using default bc object after it has been overriden.\n");
+ }
+#endif
+ d_simple_bc.setBoundaries(boundary_type,
+ fluxes,
+ flags,
+ bdry_types);
+ d_bc_object = &d_simple_bc;
+ d_fac_ops.setPhysicalBcCoefObject(d_bc_object);
+}
+
+void CellStokesFACSolver::setBcObject(
+ const RobinBcCoefStrategy* bc_object)
+{
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (!bc_object) {
+ TBOX_ERROR(d_object_name << ": NULL pointer for boundary condition\n"
+ << "object.\n");
+ }
+#endif
+ d_bc_object = bc_object;
+ d_fac_ops.setPhysicalBcCoefObject(d_bc_object);
+}
+
+/*
+ *************************************************************************
+ * *
+ * Solve the linear system and report whether iteration converged. *
+ * *
+ * This version is for an initialized solver state. *
+ * Before solving, set the final piece of the boundary condition, *
+ * which is not known until now, and initialize some internal *
+ * solver quantities. *
+ * *
+ *************************************************************************
+ */
+
+bool CellStokesFACSolver::solveSystem(
+ const int u,
+ const int f)
+{
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (!d_solver_is_initialized) {
+ TBOX_ERROR(
+ d_object_name << ".solveSystem(int,int): uninitialized\n"
+ <<
+ "solver state. You must call initializeSolverState()\n"
+ <<
+ "before using this function. Or you can use\n"
+ <<
+ "solveSystem(int,int,...) to initialize the solver,\n"
+ << "solve and deallocate the solver.\n");
+ }
+ if (u < 0 || f < 0) {
+ TBOX_ERROR(d_object_name << ": Bad patch data id.\n");
+ }
+#endif
+ if (d_bc_object == &d_simple_bc) {
+ /*
+ * Knowing that we are using the SimpelCellRobinBcCoefsX
+ * implementation of RobinBcCoefStrategy, we must save
+ * the ghost data in u before solving.
+ * The solver overwrites it, but SimpleCellRobinBcCoefs
+ * needs to get to access it repeatedly.
+ */
+ d_simple_bc.cacheDirichletData(u);
+ }
+
+ createVectorWrappers(u, f);
+ bool solver_rval;
+ solver_rval = d_fac_precond.solveSystem(*d_uv, *d_fv);
+
+ if (d_bc_object == &d_simple_bc) {
+ /*
+ * Restore the Dirichlet cell data that were overwritten by the
+ * solve process. We do this to be backward compatible with the
+ * user code.
+ */
+ d_simple_bc.restoreDirichletData(u);
+ }
+
+ return solver_rval;
+}
+
+/*
+ *************************************************************************
+ * *
+ * Solve the linear system and report whether iteration converged. *
+ * *
+ * This version is for an uninitialized solver state. *
+ * 1. Initialize the (currently uninitialized) solver state. *
+ * 2. Solve. *
+ * 3. Deallocate the solver state. *
+ * *
+ *************************************************************************
+ */
+
+bool CellStokesFACSolver::solveSystem(
+ const int u,
+ const int f,
+ tbox::Pointer<hier::PatchHierarchy> hierarchy,
+ int coarse_ln,
+ int fine_ln)
+{
+ TBOX_ASSERT(!hierarchy.isNull());
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS1(d_dim, *hierarchy);
+
+ if (d_enable_logging) {
+ tbox::plog << "CellStokesFACSolver::solveSystem (" << d_object_name
+ << ")\n";
+ d_stokes_spec.printClassData(tbox::plog);
+ }
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (d_solver_is_initialized) {
+ TBOX_ERROR(
+ d_object_name << ".solveSystem(int,int,...): initialized\n"
+ <<
+ "solver state. This function can only used when the\n"
+ <<
+ "solver state is uninitialized. You should deallocate\n"
+ <<
+ "the solver state or use solveSystem(int,int).\n");
+ }
+ if (!hierarchy) {
+ TBOX_ERROR(d_object_name << ".solveSystem(): Null hierarchy\n"
+ << "specified.\n");
+ }
+#endif
+ initializeSolverState(u, f, hierarchy, coarse_ln, fine_ln);
+
+ bool solver_rval;
+ solver_rval = solveSystem(u, f);
+
+ deallocateSolverState();
+
+ return solver_rval;
+}
+
+void CellStokesFACSolver::createVectorWrappers(
+ int u,
+ int f) {
+
+ hier::VariableDatabase& vdb(*hier::VariableDatabase::getDatabase());
+ tbox::Pointer<hier::Variable> variable;
+
+ if (!d_uv || d_uv->getComponentDescriptorIndex(0) != u) {
+ d_uv.setNull();
+ d_uv = new SAMRAIVectorReal<double>(d_object_name + "::uv",
+ d_hierarchy,
+ d_ln_min,
+ d_ln_max);
+ vdb.mapIndexToVariable(u, variable);
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (!variable) {
+ TBOX_ERROR(d_object_name << ": No variable for patch data index "
+ << u << "\n");
+ }
+ tbox::Pointer<pdat::CellVariable<double> > cell_variable = variable;
+ if (!cell_variable) {
+ TBOX_ERROR(d_object_name << ": hier::Patch data index " << u
+ << " is not a cell-double variable.\n");
+ }
+#endif
+ d_uv->addComponent(variable, u, s_weight_id[d_dim.getValue() - 1]);
+ }
+
+ if (!d_fv || d_fv->getComponentDescriptorIndex(0) != f) {
+ d_fv.setNull();
+ d_fv = new SAMRAIVectorReal<double>(d_object_name + "::fv",
+ d_hierarchy,
+ d_ln_min,
+ d_ln_max);
+ vdb.mapIndexToVariable(f, variable);
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (!variable) {
+ TBOX_ERROR(d_object_name << ": No variable for patch data index "
+ << f << "\n");
+ }
+ tbox::Pointer<pdat::CellVariable<double> > cell_variable = variable;
+ if (!cell_variable) {
+ TBOX_ERROR(d_object_name << ": hier::Patch data index " << f
+ << " is not a cell-double variable.\n");
+ }
+#endif
+ d_fv->addComponent(variable, f, s_weight_id[d_dim.getValue() - 1]);
+ }
+}
+
+/*
+ ***********************************************************************
+ * Delete the vector wrappers. Do not freeVectorComponents because *
+ * we do not control their data allocation. The user does that. *
+ ***********************************************************************
+ */
+void CellStokesFACSolver::destroyVectorWrappers() {
+ d_uv.setNull();
+ d_fv.setNull();
+}
+
+void CellStokesFACSolver::initializeStatics() {
+
+ for (int d = 0; d < SAMRAI::tbox::Dimension::MAXIMUM_DIMENSION_VALUE; ++d) {
+ s_weight_id[d] = -1;
+ s_instance_counter[d] = -1;
+ }
+
+ s_initialized = 1;
+}
+
+}
+}
+#endif
diff -r cd1d77fa073e -r fe9d63509c19 CellStokesFACSolver.I
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/CellStokesFACSolver.I Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,142 @@
+/*************************************************************************
+ *
+ * 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: High-level solver (wrapper) for scalar stokes equation.
+ *
+ ************************************************************************/
+namespace SAMRAI {
+namespace solv {
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setDPatchDataId(
+ int id) {
+ d_stokes_spec.setDPatchDataId(id);
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setDConstant(
+ double scalar) {
+ d_stokes_spec.setDConstant(scalar);
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setCPatchDataId(
+ int id) {
+ d_stokes_spec.setCPatchDataId(id);
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setCConstant(
+ double scalar) {
+// Disable Intel warning on real comparison
+#ifdef __INTEL_COMPILER
+#pragma warning (disable:1572)
+#endif
+ if (scalar == 0.0) {
+ d_stokes_spec.setCZero();
+ } else {
+ d_stokes_spec.setCConstant(scalar);
+ }
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setProlongationMethod(
+ const std::string& prolongation_method)
+{
+ d_fac_ops.setProlongationMethod(prolongation_method);
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setCoarsestLevelSolverChoice(
+ const std::string& choice)
+{
+ d_fac_ops.setCoarsestLevelSolverChoice(choice);
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setCoarsestLevelSolverTolerance(
+ double tol)
+{
+ d_fac_ops.setCoarsestLevelSolverTolerance(tol);
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setCoarsestLevelSolverMaxIterations(
+ int max_iterations)
+{
+ d_fac_ops.setCoarsestLevelSolverMaxIterations(max_iterations);
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setCoarseFineDiscretization(
+ const std::string& coarsefine_method)
+{
+ d_fac_ops.setCoarseFineDiscretization(coarsefine_method);
+}
+
+#ifdef HAVE_HYPRE
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setUseSMG(
+ bool use_smg)
+{
+ if (d_solver_is_initialized) {
+ TBOX_ERROR(
+ d_object_name << ": setUseSMG(bool) may NOT be called\n"
+ <<
+ "while the solver state is initialized, as that\n"
+ << "would lead to a corrupted solver state.\n");
+ }
+ d_fac_ops.setUseSMG(use_smg);
+}
+#endif
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setPresmoothingSweeps(
+ int num_pre_sweeps) {
+ d_fac_precond.setPresmoothingSweeps(num_pre_sweeps);
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setPostsmoothingSweeps(
+ int num_post_sweeps) {
+ d_fac_precond.setPostsmoothingSweeps(num_post_sweeps);
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setMaxCycles(
+ int max_cycles) {
+ d_fac_precond.setMaxCycles(max_cycles);
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::setResidualTolerance(
+ double residual_tol) {
+ d_fac_precond.setResidualTolerance(residual_tol);
+}
+
+SAMRAI_INLINE_KEYWORD
+int CellStokesFACSolver::getNumberOfIterations() const
+{
+ return d_fac_precond.getNumberOfIterations();
+}
+
+SAMRAI_INLINE_KEYWORD
+double CellStokesFACSolver::getResidualNorm() const
+{
+ return d_fac_precond.getResidualNorm();
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesFACSolver::getConvergenceFactors(
+ double& avg_factor,
+ double& final_factor)
+const
+{
+ d_fac_precond.getConvergenceFactors(avg_factor, final_factor);
+}
+
+}
+}
diff -r cd1d77fa073e -r fe9d63509c19 CellStokesFACSolver.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/CellStokesFACSolver.h Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,660 @@
+/*************************************************************************
+ *
+ * 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: High-level solver (wrapper) for scalar stokes equation.
+ *
+ ************************************************************************/
+#ifndef included_solv_CellStokesFACSolver
+#define included_solv_CellStokesFACSolver
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#include "SAMRAI/solv/FACPreconditioner.h"
+#include "CellStokesFACOps.h"
+#include "StokesSpecifications.h"
+#include "SAMRAI/solv/SimpleCellRobinBcCoefs.h"
+#include "SAMRAI/tbox/Database.h"
+#include "SAMRAI/tbox/Pointer.h"
+
+namespace SAMRAI {
+namespace solv {
+
+/*!
+ * @brief Class for solving scalar Stokes's equation on SAMR grid,
+ * wrapping up lower-level components (FAC cycling, Stokes equation
+ * operations and boundary conditions) in a single high-level interface.
+ *
+ * Note: this class provides a backward-compatible interface to
+ * the soon-to-be obsolete StokesHierarchySolver<DIM> class.
+ * Although this class hides the lower-level components (FAC cycling,
+ * Stokes equation operations and boundary conditions), it is
+ * perfectly acceptable to use those lower-level components directly.
+ *
+ * We solve the equation
+ * div(D grad(u)) + Cu = f
+ * where D is a side-centered array and C is a cell-centered array.
+ * u and f are also cell-centered.
+ * Boundary conditions supported are Dirichlet, Neumann and mixed
+ * (Dirichlet on some faces and Neumann on others).
+ *
+ * This class is a wrapper, providing a single class that coordinates
+ * three major components: the FAC solver, the cell-centered Stokes
+ * FAC operator and a default Robin bc coefficient implelemtation.
+ * It is perfectly acceptable to use those classes outside of this
+ * class.
+ *
+ * The underlying solver is an FAC solver using cell-centered
+ * discretization. The difference scheme is second-order
+ * central-difference. On coarse-fine boundaries within the
+ * solution levels, the composite grid operator uses, by default,
+ * the discretization method of Ewing, Lazarov and Vassilevski
+ * ("Local Refinement Techniques for Elliptic Problems on
+ * Cell-Centered Grids, I. Error Analysis", Mathematics of
+ * Computation, Vol. 56, No. 194, April 1991, pp. 437-461).
+ *
+ * Typical use of this class is:
+ * -# Construct a CellStokesFACSolver object, providing it
+ * the hierarchy and range of levels participating in the solve.
+ * -# Set the parameters C and D using the functions named @c setC...
+ * and @c setD... By default, D=1 and C=0 everywhere.
+ * -# Call setBoundaries() to state the types boundary conditions,
+ * along with supplemental data for setting those boundary
+ * conditions.
+ * -# Call initializeSolverState() to set up information
+ * internal to the solver. This is step is not required
+ * but will save setup costs if you are making multiple
+ * solves. This commits the object to the current hierarchy state
+ * and the specific @em types of boundary conditions you selected,
+ * It does NOT commit to the specific @em values of the boundary
+ * condition. A hierarchy change (through adaption or other means)
+ * invalidates the state, thus you must reinitialize or
+ * deallocateSolverState() the state before another solve.
+ * -# Solve the equation with solveSystem(). You provide the
+ * patch data indices for the solution u and the right hand
+ * side f. u must have at least one ghost cell and where
+ * a Dirichlet boundary condition applies, those cells
+ * must be set to the value on the boundary. If only Neumann
+ * boundary conditions are used, the ghost cell values
+ * do not matter.
+ * -# Call deallocateSolverState() to free up internal resources,
+ * if initializeSolverState() was called before the solve.
+ *
+ * After the solve, information on the solve can be obtained
+ * by calling one of these functions:
+ * - getNumberOfIterations() gives the number of FAC cycles used.
+ * - getConvergenceFactors() gives the average and final convergence
+ * factors for the solve.
+ * - getResidualNorm() gives the final residual
+ *
+ * Finer solver controls can be set using the functions in this class.
+ *
+ * Object of this class can be set using input databases.
+ * The following parameters can be set. Each is shown with its
+ * default value in the case where hypre is used.
+ * @verbatim
+ * enable_logging = TRUE // Bool flag to switch logging on/off
+ * max_cycles = 10 // Integer number of max FAC cycles to use
+ * residual_tol = 1.e-6 // Residual tolerance to solve for
+ * num_pre_sweeps = 1 // Number of presmoothing sweeps to use
+ * num_post_sweeps = 1 // Number of postsmoothing sweeps to use
+ * coarse_fine_discretization = "Ewing" // Name of coarse-fine discretization
+ * prolongation_method = "CONSTANT_REFINE" // Name of prolongation method
+ * coarse_solver_choice = "hypre" // Name of coarse level solver
+ * coarse_solver_tolerance = 1e-10 // Coarse level tolerance
+ * coarse_solver_max_iterations = 20 // Coarse level max iterations
+ * use_smg = "FALSE" // Whether to use hypre's smg solver
+ * // (alternative is the pfmg solver)
+ * @endverbatim
+ *
+ */
+class CellStokesFACSolver
+{
+
+public:
+ /*!
+ * @brief Construct a solver.
+ *
+ * If the database is not NULL, initial settings will be set
+ * using the database.
+ * The solver is uninitialized until initializeSolverState()
+ * is called.
+ *
+ * @param object_name Name of object used in outputs
+ * @param database tbox::Database for initialization (may be NULL)
+ */
+ CellStokesFACSolver(
+ const tbox::Dimension& dim,
+ const std::string& object_name,
+ tbox::Pointer<tbox::Database> database =
+ tbox::Pointer<tbox::Database>());
+
+ /*!
+ * @brief Destructor.
+ */
+ ~CellStokesFACSolver(
+ void);
+
+ /*!
+ * @brief Enable logging.
+ *
+ * To disable, pass in @c false.
+ */
+ void
+ enableLogging(
+ bool logging);
+
+ /*!
+ * @brief Solve Stokes's equation, assuming an uninitialized
+ * solver state.
+ *
+ * Here, u is the "solution" patch data index and f is the
+ * right hand side patch data index.
+ * The return value is true if the solver converged and false otherwise.
+ *
+ * This function is a wrapper.
+ * It simply initializes the solver state, call the
+ * solveSystem(const int,const int) for the initialized solver then
+ * deallocates the solver state.
+ *
+ * Upon return from this function,
+ * solution will contain the result of the solve.
+ *
+ * See initializeSolverState() for opportunities to save overhead
+ * when using multiple consecutive solves.
+ *
+ * @see solveSystem(const int,const int)
+ *
+ * @param solution hier::Patch data index for solution u
+ * @param rhs hier::Patch data index for right hand side f
+ * @param hierarchy The patch hierarchy to solve on
+ * @param coarse_ln The coarsest level in the solve.
+ * @param fine_ln The finest level in the solve.
+ *
+ * @return whether solver converged to specified level
+ *
+ * @see initializeSolverState
+ */
+ bool
+ solveSystem(
+ const int solution,
+ const int rhs,
+ tbox::Pointer<hier::PatchHierarchy> hierarchy,
+ int coarse_ln = -1,
+ int fine_ln = -1);
+
+ /*!
+ * @brief Solve Stokes's equation using the current solver state
+ * set by initializeSolverState().
+ *
+ * When the solver state has been initialized, this function may
+ * be called repeadedly with different values on the rhs.
+ * There is some cost savings for multiple solves when this
+ * is done.
+ *
+ * Before calling this function, the solution and
+ * right-hand-side quantities should be set properly by the user
+ * on all patch interiors on the range of levels covered by the
+ * FAC iteration. All data for these patch data index should be allocated.
+ * Thus, the user is responsible for managing the
+ * storage for the solution and right-hand-side.
+ *
+ * @return whether solver converged to specified level
+ *
+ * @see solveSystem( const int, const int, tbox::Pointer< hier::PatchHierarchy >, int, int);
+ */
+ bool
+ solveSystem(
+ const int solution,
+ const int rhs);
+
+ /*!
+ * @brief Specify the boundary conditions that are to be used at the
+ * physical domain boundary.
+ *
+ * This method is used to set up the default SimpleCellRobinBcCoefs
+ * object for specifying boundary conditions. Note that you may
+ * alternatively provide your own implementation of the Robin
+ * boundary condition coefficients using the setBcObject() method.
+ *
+ * The boundary conditions specified as the
+ * std::string argument "boundary_type." The boundary type argument can be
+ * "Dirichlet", "Neumann", or "Mixed".
+ *
+ * If using Dirichlet boundary conditions, then before the solver is
+ * called, the storage for the unknown u
+ * must have a mapped_box_level of ghost cells at least one cell wide that includes
+ * the Dirichlet boundary values.
+ *
+ * If using Neumann boundary conditions, then before the solver is called,
+ * the outerface boundary flux data must be set for the Neumann conditions.
+ * The fluxes argument gives the patch data index of this flux
+ * data.
+ *
+ * The mixed boundary type is for a mixture of Dirichlet and Neumann
+ * boundary conditions are used at the physical domain boundary.
+ * The fluxes argument gives the patch data index of the outerface data
+ * that specifies the flux data for the Neumann conditions. The flags
+ * array is an outerface data array of integer flags that specifies whether
+ * Dirichlet (flag == zero) or Neumann (flag == one) conditions are to be
+ * used at a particular cell boundary face. Note that the flag data must
+ * be set before the matrix entries can be computed and the flux data
+ * must be set before the solver is called. The bdry_types argument can
+ * be used if the boundary conditions are mixed but one or more of the
+ * faces of the physical boundary are entirely either Dirichlet or
+ * Neumann boundaries. The bdry_types argument should be an array of
+ * 2*DIM integers, specifying the boundary conditions on each side of
+ * the physical domain. It should be ordered {x_lo, x_hi, y_lo, y_hi,
+ * z_lo, z_hi}, with the values for each face being 0 for Dirichlet
+ * conditions, 1 for Neumann conditions, and 2 for mixed boundary
+ * conditions. The bdry_type argument is never required, but if used
+ * it can sometimes make the StokesHYPRESolver class more efficient.
+ */
+
+ void
+ setBoundaries(
+ const std::string& boundary_type,
+ const int fluxes = -1,
+ const int flags = -1,
+ int* bdry_types = NULL);
+
+ /*!
+ * @brief Override internal implementation to set boundary condition
+ * coefficients with user-provided implementation.
+ *
+ * This function is used to override the default internal
+ * object for setting Robin boundary condition coefficients.
+ * You should override when you need to avoid the limitations
+ * of the SimpleCellRobinBcCoefs class or you prefer to
+ * use your own implementation.
+ *
+ * Note that an important limitation of the SimpleCellRobinBcCoefs
+ * class is the inability to support linear interpolation in
+ * the prolongation step.
+ *
+ * Once the boundary condition object is overwritten by this
+ * method, you must no longer call the setBoundaries() method.
+ */
+ void
+ setBcObject(
+ const RobinBcCoefStrategy* bc_object);
+
+ //!@{ @name Specifying PDE parameters
+
+ /*!
+ * @brief Set the patch data index for variable D.
+ *
+ * In addition, disregard any previous D
+ * specified by setDConstant().
+ */
+ void
+ setDPatchDataId(
+ int id);
+
+ /*!
+ * @brief Set the scalar value variable D.
+ *
+ * In addition, disregard any previous D
+ * specified by setDPatchDataId().
+ */
+ void
+ setDConstant(
+ double scalar);
+
+ /*!
+ * @brief Set the scalar value variable C.
+ *
+ * In addition, disregard any previous C
+ * specified by setCConstant().
+ */
+ void
+ setCPatchDataId(
+ int id);
+
+ /*!
+ * @brief Set the patch data index for variable C.
+ *
+ * In addition, disregard any previous C
+ * specified by setCConstant().
+ */
+ void
+ setCConstant(
+ double scalar);
+
+ //@}
+
+ //@{ @name Functions for setting solver mathematic algorithm controls
+
+ /*!
+ * @brief Set coarse level solver.
+ *
+ * Select from these:
+ * - @c "redblack"
+ * - @c "hypre" (only if the HYPRE library is available).
+ */
+ void
+ setCoarsestLevelSolverChoice(
+ const std::string& choice);
+
+ /*!
+ * @brief Set tolerance for coarse level solve.
+ *
+ * If the coarse level solver requires a tolerance
+ * (currently, they all do), the specified value is used.
+ */
+ void
+ setCoarsestLevelSolverTolerance(
+ double tol);
+
+ /*!
+ * @brief Set max iterations for coarse level solve.
+ *
+ * If the coarse level solver requires a max iteration limit
+ * (currently, they all do), the specified value is used.
+ */
+ void
+ setCoarsestLevelSolverMaxIterations(
+ int max_iterations);
+
+#ifdef HAVE_HYPRE
+ /*!
+ * @brief Set whether to use HYPRe's PFMG algorithm instead of the
+ * SMG algorithm.
+ *
+ * The flag is used to select which of HYPRE's linear solver algorithms
+ * to use if true, the semicoarsening multigrid algorithm is used, and if
+ * false, the ``PF'' multigrid algorithm is used.
+ * By default, the SMG algorithm is used.
+ *
+ * This setting has effect only when HYPRe is chosen for the coarsest
+ * level solver. See setCoarsestLevelSolverChoice().
+ *
+ * Changing the algorithm must be done before setting up the matrix
+ * coefficients.
+ */
+ void
+ setUseSMG(
+ bool use_smg);
+#endif
+
+ /*!
+ * @brief Set the coarse-fine boundary discretization method.
+ *
+ * Specify the @c op_name std::string which will be passed to
+ * xfer::Geometry::lookupRefineOperator() to get the operator
+ * for setting fine grid ghost cells from the coarse grid.
+ * Note that chosing this operator implicitly choses the
+ * discretization method at the coarse-fine boundary.
+ *
+ * There is one important instance where this std::string is
+ * @em not passed to xfer::Geometry::lookupRefineOperator().
+ * If this variable is set to "Ewing", a constant refinement
+ * method is used along with Ewing's correction.
+ * For a reference to the correction method, see
+ * "Local Refinement Techniques for Elliptic Problems on Cell-Centered
+ * Grids, I. Error Analysis", Mathematics of Computation, Vol. 56, No. 194,
+ * April 1991, pp. 437-461.
+ *
+ * @param coarsefine_method String selecting the coarse-fine discretization method.
+ */
+ void
+ setCoarseFineDiscretization(
+ const std::string& coarsefine_method);
+
+ /*!
+ * @brief Set the name of the prolongation method.
+ *
+ * Specify the @c op_name std::string which will be passed to
+ * xfer::Geometry::lookupRefineOperator() to get the operator
+ * for prolonging the coarse-grid correction.
+ *
+ * By default, "CONSTANT_REFINE" is used. "LINEAR_REFINE" seems to
+ * to lead to faster convergence, but it does NOT satisfy the Galerkin
+ * condition.
+ *
+ * Prolonging using linear refinement requires a Robin bc
+ * coefficient implementation that is capable of delivering
+ * coefficients for non-hierarchy data, because linear refinement
+ * requires boundary conditions to be set on temporary levels.
+ *
+ * @param prolongation_method String selecting the coarse-fine discretization method.
+ */
+ void
+ setProlongationMethod(
+ const std::string& prolongation_method);
+
+ /*!
+ * @brief Set the number of pre-smoothing sweeps during
+ * FAC iteration process.
+ *
+ * Presmoothing is applied during the fine-to-coarse phase of the
+ * iteration. The default is to use one sweep.
+ *
+ * @param num_pre_sweeps Number of presmoothing sweeps
+ */
+ void
+ setPresmoothingSweeps(
+ int num_pre_sweeps);
+
+ /*!
+ * @brief Set the number of post-smoothing sweeps during
+ * FAC iteration process.
+ *
+ * Postsmoothing is applied during the coarse-to-fine phase of the
+ * iteration. The default is to use one sweep.
+ *
+ * @param num_post_sweeps Number of postsmoothing sweeps
+ */
+ void
+ setPostsmoothingSweeps(
+ int num_post_sweeps);
+
+ /*!
+ * @brief Set the max number of iterations (cycles) to use per solve.
+ */
+ void
+ setMaxCycles(
+ int max_cycles);
+
+ /*!
+ * @brief Set the residual tolerance for stopping.
+ *
+ * If you want the prescribed maximum number of cycles to always be taken,
+ * set the residual tolerance to a negative number.
+ */
+ void
+ setResidualTolerance(
+ double residual_tol);
+
+ //@}
+
+ /*!
+ * @brief Prepare the solver's internal state for solving
+ *
+ * In the interest of efficiency, this class may prepare and
+ * cache some hierarchy-dependent objects. Though it is not required,
+ * initializing the solver state makes for greater efficiency
+ * when you are doing multiple solves on the same system of
+ * equation. If you do not initialize the state, it is initialized
+ * and deallocated each time you call solveSystem(const int, const int).
+ * The state must be reinitialized if the hierarchy or a boundary
+ * condition type changes.
+ *
+ * To unset the data set in this function,
+ * see deallocateSolverState().
+ *
+ * The @c solution and @c rhs patch data indices in the argument
+ * list are used to determine the @em form of the data you
+ * plan to use in the solve. They need not be the same data
+ * you solve on, but they should be similar. Both must represent
+ * cell-centered double data. The solution must have at least one
+ * ghost cell width, though this is not checked in the initialize
+ * phase, because data is not required yet.
+ *
+ * @param solution solution patch data index for u
+ * @param rhs right hand side patch data index for f
+ * @param hierarchy The patch hierarchy to solve on
+ * @param coarse_level The coarsest level in the solve
+ * @param fine_level The finest level in the solve
+ */
+ void
+ initializeSolverState(
+ const int solution,
+ const int rhs,
+ tbox::Pointer<hier::PatchHierarchy> hierarchy,
+ const int coarse_level = -1,
+ const int fine_level = -1);
+
+ /*!
+ * @brief Remove the solver's internal state data
+ *
+ * Remove all hierarchy-dependent data set by initializeSolverState.
+ * It is safe to call deallocateSolverState() even state is already
+ * deallocated, but nothing is done in that case.
+ *
+ * @see initializeSolverState()
+ */
+ void
+ deallocateSolverState();
+
+ //@{
+ //! @name Functions to get data on last solve.
+
+ /*!
+ * @brief Return FAC iteration count from last (or current
+ * if there is one) FAC iteration process.
+ */
+ int
+ getNumberOfIterations() const;
+
+ /*!
+ * @brief Get average convergance rate and convergence rate of
+ * the last (or current if there is one) FAC solve.
+ *
+ * @param avg_factor average convergence factor over current FAC cycles
+ * @param final_factor convergence factor of the last FAC cycle
+ */
+ void
+ getConvergenceFactors(
+ double& avg_factor,
+ double& final_factor) const;
+
+ /*!
+ * @brief Return residual norm from the just-completed FAC iteration.
+ *
+ * The norm return value is computed as the maximum norm over all
+ * patch levels involved in the solve. The value corresponds to the
+ * norm applied in the user-defined residual computation.
+ *
+ * The latest computed norm is the one returned.
+ */
+ double
+ getResidualNorm() const;
+
+ //@}
+
+private:
+ /*!
+ * @brief Set state using database
+ *
+ * See the class description for the parameters that can be set
+ * from a database.
+ *
+ * @param database Input database. If a NULL pointer is given,
+ * nothing is done.
+ */
+ void
+ getFromInput(
+ tbox::Pointer<tbox::Database> database);
+
+ /*
+ * @brief Set @c d_uv and @c d_fv to vectors wrapping the data
+ * specified by patch data indices u and f.
+ */
+ void
+ createVectorWrappers(
+ int u,
+ int f);
+
+ /*
+ * @brief Destroy vector wrappers referenced to by @c d_uv and @c d_fv.
+ */
+ void
+ destroyVectorWrappers();
+
+ /*
+ * @brief Initialize static members
+ */
+ static void
+ initializeStatics();
+
+ const tbox::Dimension d_dim;
+
+ /*!
+ * @brief Object name.
+ */
+ std::string d_object_name;
+
+ /*!
+ * @brief Object holding the specifications of the Stokes equation.
+ */
+ StokesSpecifications d_stokes_spec;
+
+ /*!
+ * @brief FAC operator implementation corresponding to cell-centered
+ * Stokes discretization.
+ */
+ CellStokesFACOps d_fac_ops;
+
+ /*!
+ * @brief FAC preconditioner algorithm.
+ */
+ FACPreconditioner d_fac_precond;
+
+ /*!
+ * @brief Robin bc object in use.
+ */
+ const RobinBcCoefStrategy* d_bc_object;
+
+ /*
+ * @brief Default implementation of RobinBcCoefStrategy
+ */
+ SimpleCellRobinBcCoefs d_simple_bc;
+
+ tbox::Pointer<hier::PatchHierarchy> d_hierarchy;
+ int d_ln_min;
+ int d_ln_max;
+
+ /*!
+ * @brief Context for all internally maintained data.
+ */
+ tbox::Pointer<hier::VariableContext> d_context;
+ /*
+ * @brief Vector wrapper for solution.
+ * @see createVectorWrappers(), destroyVectorWrappers()
+ */
+ tbox::Pointer<SAMRAIVectorReal<double> > d_uv;
+ /*
+ * @brief Vector wrapper for source.
+ * @see createVectorWrappers(), destroyVectorWrappers()
+ */
+ tbox::Pointer<SAMRAIVectorReal<double> > d_fv;
+
+ bool d_solver_is_initialized;
+ bool d_enable_logging;
+
+ static bool s_initialized;
+ static int s_weight_id[SAMRAI::tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+ static int s_instance_counter[SAMRAI::tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+};
+
+}
+}
+
+#ifdef SAMRAI_INLINE
+#include "CellStokesFACSolver.I"
+#endif
+
+#endif // included_solv_CellStokesFACSolver
diff -r cd1d77fa073e -r fe9d63509c19 CellStokesHypreSolver.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/CellStokesHypreSolver.C Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,1590 @@
+/*************************************************************************
+ *
+ * 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: Hypre solver interface for diffusion-like elliptic problems.
+ *
+ ************************************************************************/
+#ifndef included_solv_CellStokesHypreSolver_C
+#define included_solv_CellStokesHypreSolver_C
+
+#include "CellStokesHypreSolver.h"
+
+#ifdef HAVE_HYPRE
+
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/geom/CartesianGridGeometry.h"
+#include "SAMRAI/math/ArrayDataBasicOps.h"
+#include "SAMRAI/math/PatchSideDataBasicOps.h"
+#include "SAMRAI/pdat/ArrayData.h"
+#include "SAMRAI/pdat/CellIndex.h"
+#include "SAMRAI/pdat/CellIterator.h"
+#include "SAMRAI/pdat/FaceIndex.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/pdat/SideIndex.h"
+#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/pdat/OuterfaceData.h"
+#include "SAMRAI/pdat/OutersideData.h"
+#include "SAMRAI/hier/BoundaryBoxUtils.h"
+#include "SAMRAI/hier/VariableDatabase.h"
+#include "SAMRAI/tbox/MathUtilities.h"
+#include "SAMRAI/tbox/SAMRAI_MPI.h"
+#include "SAMRAI/tbox/SAMRAIManager.h"
+#include "SAMRAI/tbox/PIO.h"
+#include "SAMRAI/tbox/Timer.h"
+#include "SAMRAI/tbox/TimerManager.h"
+#include "SAMRAI/tbox/StartupShutdownManager.h"
+#include "SAMRAI/tbox/Utilities.h"
+
+#include <cstdlib>
+
+#ifndef SAMRAI_INLINE
+#include "CellStokesHypreSolver.I"
+#endif
+
+extern "C" {
+
+#ifdef __INTEL_COMPILER
+#pragma warning (disable:1419)
+#endif
+
+void F77_FUNC(compdiagvariablec2d, COMPDIAGVARIABLEC2D) (
+ double* diag,
+ const double* c,
+ const double* offdiagi,
+ const double* offdiagj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* cscale,
+ const double* dscale);
+void F77_FUNC(compdiagscalarc2d, COMPDIAGSCALARC2D) (
+ double* diag,
+ const double* c,
+ const double* offdiagi,
+ const double* offdiagj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* cscale,
+ const double* dscale);
+void F77_FUNC(compdiagzeroc2d, COMPDIAGZEROC2D) (
+ double* diag,
+ const double* offdiagi,
+ const double* offdiagj,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const double* cscale,
+ const double* dscale);
+void F77_FUNC(adjbdry2d, ADJBDRY2D) (
+ double* diag,
+ const double* offdiagi,
+ const double* offdiagj,
+ const int* pifirst, const int* pilast,
+ const int* pjfirst, const int* pjlast,
+ const double* acoef,
+ const double* bcoef,
+ const int* aifirst, const int* ailast,
+ const int* ajfirst, const int* ajlast,
+ const double* Ak0,
+ const int* kifirst, const int* kilast,
+ const int* kjfirst, const int* kjlast,
+ const int* lower, const int* upper,
+ const int* location,
+ const double* h);
+void F77_FUNC(adjbdryconstoffdiags2d, ADJBDRYCONSTOFFDIAGS2D) (
+ double* diag,
+ const double* offdiag,
+ const int* pifirst,
+ const int* pilast,
+ const int* pjfirst,
+ const int* pjlast,
+ const double* acoef,
+ const int* aifirst,
+ const int* ailast,
+ const int* ajfirst,
+ const int* ajlast,
+ const double* Ak0,
+ const int* kifirst,
+ const int* kilast,
+ const int* kjfirst,
+ const int* kjlast,
+ const int* lower, const int* upper,
+ const int* location,
+ const double* h);
+void F77_FUNC(adjustrhs2d, ADJUSTRHS2D) (double* rhs,
+ const int* rifirst,
+ const int* rilast,
+ const int* rjfirst,
+ const int* rjlast,
+ const double* Ak0,
+ const int* kifirst,
+ const int* kilast,
+ const int* kjfirst,
+ const int* kjlast,
+ const double* gcoef,
+ const int* aifirst,
+ const int* ailast,
+ const int* ajfirst,
+ const int* ajlast,
+ const int* lower, const int* upper,
+ const int* location);
+
+void F77_FUNC(compdiagvariablec3d, COMPDIAGVARIABLEC3D) (
+ double* diag,
+ const double* c,
+ const double* offdiagi,
+ const double* offdiagj,
+ const double* offdiagk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* cscale,
+ const double* dscale);
+void F77_FUNC(compdiagscalarc3d, COMPDIAGSCALARC3D) (
+ double* diag,
+ const double* c,
+ const double* offdiagi,
+ const double* offdiagj,
+ const double* offdiagk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* cscale,
+ const double* dscale);
+void F77_FUNC(compdiagzeroc3d, COMPDIAGZEROC3D) (
+ double* diag,
+ const double* offdiagi,
+ const double* offdiagj,
+ const double* offdiagk,
+ const int* ifirst,
+ const int* ilast,
+ const int* jfirst,
+ const int* jlast,
+ const int* kfirst,
+ const int* klast,
+ const double* cscale,
+ const double* dscale);
+void F77_FUNC(adjbdry3d, ADJBDRY3D) (
+ double* diag,
+ const double* offdiagi,
+ const double* offdiagj,
+ const double* offdiagk,
+ const int* pifirst,
+ const int* pilast,
+ const int* pjfirst,
+ const int* pjlast,
+ const int* pkfirst,
+ const int* pklast,
+ const double* acoef,
+ const double* bcoef,
+ const int* aifirst,
+ const int* ailast,
+ const int* ajfirst,
+ const int* ajlast,
+ const int* akfirst,
+ const int* aklast,
+ const double* Ak0,
+ const int* kifirst,
+ const int* kilast,
+ const int* kjfirst,
+ const int* kjlast,
+ const int* kkfirst,
+ const int* kklast,
+ const int* lower, const int* upper,
+ const int* location,
+ const double* h);
+void F77_FUNC(adjbdryconstoffdiags3d, ADJBDRYCONSTOFFDIAGS3D) (
+ double* diag,
+ const double* offdiag,
+ const int* pifirst,
+ const int* pilast,
+ const int* pjfirst,
+ const int* pjlast,
+ const int* pkfirst,
+ const int* pklast,
+ const double* acoef,
+ const int* aifirst,
+ const int* ailast,
+ const int* ajfirst,
+ const int* ajlast,
+ const int* akfirst,
+ const int* aklast,
+ const double* Ak0,
+ const int* kifirst,
+ const int* kilast,
+ const int* kjfirst,
+ const int* kjlast,
+ const int* kkfirst,
+ const int* kklast,
+ const int* lower, const int* upper,
+ const int* location,
+ const double* h);
+void F77_FUNC(adjustrhs3d, ADJUSTRHS3D) (double* rhs,
+ const int* rifirst,
+ const int* rilast,
+ const int* rjfirst,
+ const int* rjlast,
+ const int* rkfirst,
+ const int* rklast,
+ const double* Ak0,
+ const int* kifirst,
+ const int* kilast,
+ const int* kjfirst,
+ const int* kjlast,
+ const int* kkfirst,
+ const int* kklast,
+ const double* gcoef,
+ const int* aifirst,
+ const int* ailast,
+ const int* ajfirst,
+ const int* ajlast,
+ const int* akfirst,
+ const int* aklast,
+ const int* lower, const int* upper,
+ const int* location);
+
+}
+
+namespace SAMRAI {
+namespace solv {
+
+tbox::Pointer<pdat::OutersideVariable<double> >
+CellStokesHypreSolver::s_Ak0_var[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+
+tbox::StartupShutdownManager::Handler CellStokesHypreSolver::s_finalize_handler(
+ 0,
+ 0,
+ 0,
+ CellStokesHypreSolver::finalizeCallback,
+ tbox::StartupShutdownManager::priorityVariables);
+
+/*
+ *************************************************************************
+ * Constructor *
+ *************************************************************************
+ */
+
+CellStokesHypreSolver::CellStokesHypreSolver(
+ const tbox::Dimension& dim,
+ const std::string& object_name,
+ tbox::Pointer<tbox::Database> database):
+ d_dim(dim),
+ d_object_name(object_name),
+ d_hierarchy(NULL),
+ d_ln(-1),
+ d_context(hier::VariableDatabase::getDatabase()->
+ getContext(object_name + "::context")),
+ d_cf_boundary(),
+ d_physical_bc_coef_strategy(&d_physical_bc_simple_case),
+ d_physical_bc_variable(NULL),
+ d_physical_bc_simple_case(dim, d_object_name + "::simple bc"),
+ d_cf_bc_coef(dim, object_name + "::coarse-fine bc coefs"),
+ d_coarsefine_bc_variable(NULL),
+ d_Ak0_id(-1),
+ d_soln_depth(0),
+ d_rhs_depth(0),
+ d_max_iterations(10),
+ d_relative_residual_tol(1e-10),
+ d_number_iterations(-1),
+ d_num_pre_relax_steps(1),
+ d_num_post_relax_steps(1),
+ d_relative_residual_norm(-1.0),
+ d_use_smg(false),
+ d_grid(NULL),
+ d_stencil(NULL),
+ d_matrix(NULL),
+ d_linear_rhs(NULL),
+ d_linear_sol(NULL),
+ d_mg_data(NULL),
+ d_print_solver_info(false)
+{
+ if (d_dim == tbox::Dimension(1) || d_dim > tbox::Dimension(3)) {
+ TBOX_ERROR(" CellStokesHypreSolver : DIM == 1 or > 3 not implemented");
+ }
+
+ t_solve_system = tbox::TimerManager::getManager()->
+ getTimer("solv::CellStokesHypreSolver::solveSystem()");
+ t_set_matrix_coefficients = tbox::TimerManager::getManager()->
+ getTimer("solv::CellStokesHypreSolver::setMatrixCoefficients()");
+
+ hier::VariableDatabase* vdb = hier::VariableDatabase::getDatabase();
+ if (s_Ak0_var[d_dim.getValue() - 1].isNull()) {
+ s_Ak0_var[d_dim.getValue() - 1] = new
+ pdat::OutersideVariable<double>(d_dim, d_object_name + "::Ak0", 1);
+ }
+ d_Ak0_id =
+ vdb->registerVariableAndContext(s_Ak0_var[d_dim.getValue() - 1],
+ d_context,
+ hier::IntVector::getZero(d_dim));
+ if (database) {
+ getFromInput(database);
+ }
+}
+
+/*
+ ********************************************************************
+ * Set state from database *
+ ********************************************************************
+ */
+
+void CellStokesHypreSolver::getFromInput(
+ tbox::Pointer<tbox::Database> database)
+{
+ if (database) {
+ d_print_solver_info = database->getBoolWithDefault("print_solver_info",
+ d_print_solver_info);
+ d_max_iterations = database->getIntegerWithDefault("max_iterations",
+ d_max_iterations);
+ d_relative_residual_tol = database->getDoubleWithDefault(
+ "relative_residual_tol",
+ d_relative_residual_tol);
+ if (database->isDouble("residual_tol")) {
+ TBOX_ERROR("CellStokesHypreSolver input error.\n"
+ << "The parameter 'residual_tol' has been replaced\n"
+ << "by 'relative_residual_tol' to be more descriptive.\n"
+ << "Please change the parameter name in the input database.");
+ }
+ d_num_pre_relax_steps =
+ database->getIntegerWithDefault("num_pre_relax_steps",
+ d_num_pre_relax_steps);
+ if (d_num_pre_relax_steps < 0) {
+ TBOX_ERROR(d_object_name << ": Number of relaxation steps must be\n"
+ << "non-negative.\n");
+ }
+ d_num_post_relax_steps =
+ database->getIntegerWithDefault("num_post_relax_steps",
+ d_num_post_relax_steps);
+ if (d_num_post_relax_steps < 0) {
+ TBOX_ERROR(d_object_name << ": Number of relaxation steps must be\n"
+ << "non-negative.\n");
+ }
+ if (database->isBool("use_smg")) {
+ bool use_smg = database->getBool("use_smg");
+ if (use_smg != d_use_smg) {
+ setUseSMG(use_smg);
+ }
+ }
+ }
+}
+
+/*
+ ********************************************************************
+ * Initialize internal data for a given hierarchy level *
+ * After setting internal data, propagate the information *
+ * to the major algorithm objects. Allocate data for *
+ * storing boundary condition-dependent quantities for *
+ * adding to souce term before solving. *
+ ********************************************************************
+ */
+
+void CellStokesHypreSolver::initializeSolverState(
+ tbox::Pointer<hier::PatchHierarchy> hierarchy,
+ int ln)
+{
+ TBOX_ASSERT(!hierarchy.isNull());
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS1(d_dim, *hierarchy);
+
+ deallocateSolverState();
+
+ d_hierarchy = hierarchy;
+ d_ln = ln;
+
+ hier::IntVector max_gcw(d_dim, 1);
+ d_cf_boundary = new hier::CoarseFineBoundary(*d_hierarchy, d_ln, max_gcw);
+
+ d_physical_bc_simple_case.setHierarchy(d_hierarchy, d_ln, d_ln);
+
+ d_number_iterations = -1;
+ d_relative_residual_norm = -1.0;
+
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(d_ln);
+ level->allocatePatchData(d_Ak0_id);
+ allocateHypreData();
+}
+
+/*
+ ********************************************************************
+ * Deallocate data initialized by initializeSolverState *
+ ********************************************************************
+ */
+
+void CellStokesHypreSolver::deallocateSolverState()
+{
+ if (d_hierarchy.isNull()) return;
+
+ d_cf_boundary->clear();
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(d_ln);
+ level->deallocatePatchData(d_Ak0_id);
+ deallocateHypreData();
+ d_hierarchy.setNull();
+ d_ln = -1;
+}
+
+/*
+ *************************************************************************
+ * *
+ * Allocate the HYPRE data structures that depend only on the level *
+ * and will not change (grid, stencil, matrix, and vectors). *
+ * *
+ *************************************************************************
+ */
+void CellStokesHypreSolver::allocateHypreData()
+{
+ tbox::SAMRAI_MPI::Comm communicator = d_hierarchy->getDomainMappedBoxLevel().getMPI().getCommunicator();
+
+ /*
+ * Set up the grid data - only set grid data for local boxes
+ */
+
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(d_ln);
+ tbox::Pointer<geom::CartesianGridGeometry> grid_geometry =
+ d_hierarchy->getGridGeometry();
+ const hier::IntVector ratio = level->getRatioToLevelZero();
+ hier::IntVector periodic_shift =
+ grid_geometry->getPeriodicShift(ratio);
+
+ int periodic_flag[tbox::Dimension::MAXIMUM_DIMENSION_VALUE];
+ int d;
+ bool is_periodic = false;
+ for (d = 0; d < d_dim.getValue(); ++d) {
+ periodic_flag[d] = periodic_shift[d] != 0;
+ is_periodic = is_periodic || periodic_flag[d];
+ }
+
+ HYPRE_StructGridCreate(communicator, d_dim.getValue(), &d_grid);
+ for (hier::PatchLevel::Iterator p(level); p; p++) {
+ const hier::Box& box = (*p)->getBox();
+ hier::Index lower = box.lower();
+ hier::Index upper = box.upper();
+ HYPRE_StructGridSetExtents(d_grid, &lower[0], &upper[0]);
+ }
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (is_periodic) {
+ const hier::BoxArray& level_domain = level->getPhysicalDomain();
+ hier::Box domain_bound(level_domain[0]);
+ for (int i = 1; i < level_domain.size(); ++i) {
+ domain_bound.lower().min(level_domain[i].lower());
+ domain_bound.upper().max(level_domain[i].upper());
+ }
+ for (d = 0; d < d_dim.getValue(); ++d) {
+ if (periodic_flag[d] == true) {
+ int tmpi = 1;
+ unsigned int p_of_two;
+ for (p_of_two = 0; p_of_two < 8 * sizeof(p_of_two) - 1;
+ ++p_of_two) {
+ if (tmpi == domain_bound.numberCells(d)) {
+ break;
+ }
+ if (tmpi > domain_bound.numberCells(d)) {
+ TBOX_ERROR(
+ d_object_name << ": Hypre currently requires\n"
+ <<
+ "that grid size in periodic directions be\n"
+ <<
+ "powers of two. (This requirement may go\n"
+ <<
+ "away in future versions of hypre.)\n"
+ << "Size problem in direction "
+ << d << "\n"
+ << "Domain bound is "
+ << domain_bound << ",\n"
+ << "Size of "
+ << domain_bound.numberCells() << "\n");
+ }
+ tmpi = tmpi ? tmpi << 1 : 1;
+ }
+ }
+ }
+ }
+#endif
+
+ HYPRE_StructGridSetPeriodic(d_grid, &periodic_shift[0]);
+ HYPRE_StructGridAssemble(d_grid);
+
+ {
+ /*
+ * Allocate stencil data and set stencil offsets
+ */
+
+ if (d_dim == tbox::Dimension(1)) {
+ const int stencil_size = 2;
+ int stencil_offsets[2][1] = {
+ { -1 }, { 0 }
+ };
+ HYPRE_StructStencilCreate(d_dim.getValue(), stencil_size, &d_stencil);
+ for (int s = 0; s < stencil_size; s++) {
+ HYPRE_StructStencilSetElement(d_stencil, s,
+ stencil_offsets[s]);
+ }
+ } else if (d_dim == tbox::Dimension(2)) {
+ const int stencil_size = 3;
+ int stencil_offsets[3][2] = {
+ { -1, 0 }, { 0, -1 }, { 0, 0 }
+ };
+ HYPRE_StructStencilCreate(d_dim.getValue(), stencil_size, &d_stencil);
+ for (int s = 0; s < stencil_size; s++) {
+ HYPRE_StructStencilSetElement(d_stencil, s,
+ stencil_offsets[s]);
+ }
+ } else if (d_dim == tbox::Dimension(3)) {
+ const int stencil_size = 4;
+ int stencil_offsets[4][3] = {
+ { -1, 0, 0 }, { 0, -1, 0 }, { 0, 0, -1 }, { 0, 0, 0 }
+ };
+ HYPRE_StructStencilCreate(d_dim.getValue(), stencil_size, &d_stencil);
+ for (int s = 0; s < stencil_size; s++) {
+ HYPRE_StructStencilSetElement(d_stencil, s,
+ stencil_offsets[s]);
+ }
+ }
+ }
+
+ {
+ int full_ghosts1[2 * 3] = { 1, 1, 0, 0, 0, 0 };
+ int no_ghosts1[2 * 3] = { 0, 0, 0, 0, 0, 0 };
+
+ int full_ghosts2[2 * 3] = { 1, 1, 1, 1, 0, 0 };
+ int no_ghosts2[2 * 3] = { 0, 0, 0, 0, 0, 0 };
+
+ int full_ghosts3[2 * 3] = { 1, 1, 1, 1, 1, 1 };
+ int no_ghosts3[2 * 3] = { 0, 0, 0, 0, 0, 0 };
+
+ /*
+ * Allocate the structured matrix
+ */
+
+ int* full_ghosts = NULL;
+ int* no_ghosts = NULL;
+
+ if (d_dim == tbox::Dimension(1)) {
+ full_ghosts = full_ghosts1;
+ no_ghosts = no_ghosts1;
+ } else if (d_dim == tbox::Dimension(2)) {
+ full_ghosts = full_ghosts2;
+ no_ghosts = no_ghosts2;
+ } else if (d_dim == tbox::Dimension(3)) {
+ full_ghosts = full_ghosts3;
+ no_ghosts = no_ghosts3;
+ } else {
+ TBOX_ERROR(
+ "CellStokesHypreSolver does not yet support dimension " << d_dim);
+ }
+
+ HYPRE_StructMatrixCreate(communicator,
+ d_grid,
+ d_stencil,
+ &d_matrix);
+ HYPRE_StructMatrixSetNumGhost(d_matrix, full_ghosts);
+ HYPRE_StructMatrixSetSymmetric(d_matrix, 1);
+ HYPRE_StructMatrixInitialize(d_matrix);
+
+ HYPRE_StructVectorCreate(communicator,
+ d_grid,
+ &d_linear_rhs);
+ HYPRE_StructVectorSetNumGhost(d_linear_rhs, no_ghosts);
+ HYPRE_StructVectorInitialize(d_linear_rhs);
+
+ HYPRE_StructVectorCreate(communicator,
+ d_grid,
+ &d_linear_sol);
+ HYPRE_StructVectorSetNumGhost(d_linear_sol, full_ghosts);
+ HYPRE_StructVectorInitialize(d_linear_sol);
+ }
+}
+
+/*
+ *************************************************************************
+ * *
+ * The destructor deallocates solver data. *
+ * *
+ *************************************************************************
+ */
+
+CellStokesHypreSolver::~CellStokesHypreSolver()
+{
+ deallocateHypreData();
+
+ if (!d_hierarchy.isNull()) {
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(0);
+ level->deallocatePatchData(d_Ak0_id);
+ }
+ hier::VariableDatabase* vdb =
+ hier::VariableDatabase::getDatabase();
+ vdb->removePatchDataIndex(d_Ak0_id);
+}
+
+/*
+ *************************************************************************
+ * *
+ * Deallocate HYPRE data and solver. HYPRE requires that we *
+ * check whether HYPRE has already deallocated this data. *
+ * Note that the HYPRE solver, d_mg_data, was created at *
+ * the end of setMatrixCoefficients. *
+ * *
+ *************************************************************************
+ */
+
+void CellStokesHypreSolver::deallocateHypreData()
+{
+ if (d_stencil) {
+ HYPRE_StructStencilDestroy(d_stencil);
+ d_stencil = NULL;
+ }
+ if (d_grid) {
+ HYPRE_StructGridDestroy(d_grid);
+ d_grid = NULL;
+ }
+ if (d_matrix) {
+ HYPRE_StructMatrixDestroy(d_matrix);
+ d_matrix = NULL;
+ }
+ if (d_linear_rhs) {
+ HYPRE_StructVectorDestroy(d_linear_rhs);
+ d_linear_rhs = NULL;
+ }
+ if (d_linear_sol) {
+ HYPRE_StructVectorDestroy(d_linear_sol);
+ d_linear_sol = NULL;
+ }
+ destroyHypreSolver();
+}
+
+/*
+ *************************************************************************
+ * *
+ * Copy data into the HYPRE vector structures. *
+ * *
+ *************************************************************************
+ */
+
+void CellStokesHypreSolver::copyToHypre(
+ HYPRE_StructVector vector,
+ pdat::CellData<double>& src,
+ int depth,
+ const hier::Box& box)
+{
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS2(d_dim, src, box);
+
+ for (pdat::CellIterator c(box); c; c++) {
+ hier::IntVector ic = c();
+ HYPRE_StructVectorSetValues(vector, &ic[0], src(c(), depth));
+ }
+}
+
+/*
+ *************************************************************************
+ * *
+ * Copy data out of the HYPRE vector structures. *
+ * *
+ *************************************************************************
+ */
+
+void CellStokesHypreSolver::copyFromHypre(
+ pdat::CellData<double>& dst,
+ int depth,
+ HYPRE_StructVector vector,
+ const hier::Box box)
+{
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS2(d_dim, dst, box);
+
+ for (pdat::CellIterator c(box); c; c++) {
+ double value;
+ hier::IntVector ic = c();
+ HYPRE_StructVectorGetValues(vector, &ic[0], &value);
+ dst(c(), depth) = value;
+ }
+}
+
+/*
+ *************************************************************************
+ * *
+ * Set the matrix coefficients for the linear system. *
+ * The matrix coefficients are dependent on the problem *
+ * specification described by the StokesSpecificiations *
+ * object and by the boundary condition. *
+ * *
+ *************************************************************************
+ */
+
+void CellStokesHypreSolver::setMatrixCoefficients(
+ const StokesSpecifications& spec)
+{
+ if (d_physical_bc_coef_strategy == NULL) {
+ TBOX_ERROR(
+ d_object_name << ": No BC coefficient strategy object!\n"
+ <<
+ "Use either setBoundaries or setPhysicalBcCoefObject\n"
+ <<
+ "to specify the boundary conidition. Do it before\n"
+ << "calling setMatrixCoefficients.");
+ }
+
+ t_set_matrix_coefficients->start();
+
+ int i = 0;
+
+ tbox::Pointer<pdat::CellData<double> > C_data;
+ tbox::Pointer<pdat::SideData<double> > D_data;
+
+ /*
+ * Some computations can be done using high-level math objects.
+ * Define the math objects.
+ */
+ math::ArrayDataBasicOps<double> array_math;
+ math::PatchSideDataBasicOps<double> patch_side_math;
+
+ /*
+ * The value of the ghost cell based on the Robin boundary condition
+ * can be written as the sum of a constant, k0, plus a multiple of the
+ * internal cell value, k1*ui. k1*ui depends on the value of u so it
+ * contributes to the product Au,
+ * while the constant k0 contributes the right hand side f.
+ * We save Ak0 = A*k0(a) to add to f when solving.
+ * We assume unit g here because we will multiply it in just before
+ * solving, thus allowing everything that does not affect A to change
+ * from solve to solve.
+ */
+ tbox::Pointer<pdat::OutersideData<double> > Ak0;
+
+ /*
+ * Loop over patches and set matrix entries for each patch.
+ */
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(d_ln);
+ const hier::IntVector no_ghosts(d_dim, 0);
+ for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
+
+ hier::Patch& patch = **pi;
+
+ tbox::Pointer<geom::CartesianPatchGeometry> pg =
+ patch.getPatchGeometry();
+
+ const double* h = pg->getDx();
+
+ const hier::Box patch_box = patch.getBox();
+ const hier::Index patch_lo = patch_box.lower();
+ const hier::Index patch_up = patch_box.upper();
+
+ if (!spec.cIsZero() && !spec.cIsConstant()) {
+ C_data = patch.getPatchData(spec.getCPatchDataId());
+ if (C_data.isNull()) {
+ TBOX_ERROR(d_object_name << ": Invalid cell variable index "
+ << spec.getCPatchDataId()
+ << " for the C parameter. It is not\n"
+ << "cell-centered double data.");
+ }
+ }
+
+ if (!spec.dIsConstant()) {
+ D_data = patch.getPatchData(spec.getDPatchDataId());
+ if (D_data.isNull()) {
+ TBOX_ERROR(d_object_name << ": Invalid cell variable index "
+ << spec.getDPatchDataId()
+ <<
+ " for diffusion coefficient. It is not\n"
+ << "side-centered double data.");
+ }
+ }
+
+ Ak0 = patch.getPatchData(d_Ak0_id);
+
+ Ak0->fillAll(0.0);
+
+ pdat::CellData<double> diagonal(patch_box, 1, no_ghosts);
+
+ /*
+ * Set diagonals to zero so we can accumulate to it.
+ * Accumulation is used at boundaries to shift weights
+ * for ghost cells onto the diagonal.
+ */
+ diagonal.fillAll(0.0);
+
+ const tbox::Pointer<geom::CartesianPatchGeometry>
+ geometry = patch.getPatchGeometry();
+
+ const hier::Index ifirst = patch_box.lower();
+ const hier::Index ilast = patch_box.upper();
+
+ /*
+ * Storage for off-diagonal entries,
+ * which can be variable or constant.
+ */
+ pdat::SideData<double> off_diagonal(patch_box, 1, no_ghosts);
+
+ /*
+ * Compute all off-diagonal entries with no regard to BCs.
+ * These off-diagonal entries are simply D/(h*h), according
+ * to our central difference formula.
+ */
+ if (spec.dIsConstant()) {
+ for (i = 0; i < d_dim.getValue(); ++i) {
+ double dhh = spec.getDConstant() / (h[i] * h[i]);
+ pdat::ArrayData<double>& off_diag_array(off_diagonal.getArrayData(i));
+ off_diag_array.fill(dhh);
+ }
+ } else {
+ for (i = 0; i < d_dim.getValue(); ++i) {
+ hier::Box sbox(patch_box);
+ sbox.growUpper(i, 1);
+ array_math.scale(off_diagonal.getArrayData(i),
+ 1.0 / (h[i] * h[i]),
+ D_data->getArrayData(i),
+ sbox);
+ }
+ }
+
+ /*
+ * Compute diagonal entries using off-diagonal contributions.
+ */
+ if (spec.cIsZero()) {
+ computeDiagonalEntries(diagonal,
+ off_diagonal,
+ patch_box);
+ } else if (spec.cIsConstant()) {
+ computeDiagonalEntries(diagonal,
+ spec.getCConstant(),
+ off_diagonal,
+ patch_box);
+ } else {
+ computeDiagonalEntries(diagonal,
+ *C_data,
+ off_diagonal,
+ patch_box);
+ }
+
+ /*
+ * Walk physical domain boundaries and adjust off-diagonals
+ * before computation of diagonal entries.
+ * The exterior cell's value is
+ * uo = ( h*gamma + ui*(beta-h*alpha/2) )/( beta+h*alpha/2 )
+ * = k0 + k1*ui
+ * where k0 = h*gamma/( beta+h*alpha/2 )
+ * k1 = ( beta-h*alpha/2 )/( beta+h*alpha/2 )
+ * Split coupling between interior-exterior cells
+ * into two parts: interior-interior coupling (k1)
+ * and rhs contribution (k0).
+ */
+ {
+ const tbox::Array<hier::BoundaryBox>& surface_boxes =
+ pg->getCodimensionBoundaries(1);
+ const int n_bdry_boxes = surface_boxes.getSize();
+ for (int n = 0; n < n_bdry_boxes; ++n) {
+
+ const hier::BoundaryBox& boundary_box = surface_boxes[n];
+ if (boundary_box.getBoundaryType() != 1) {
+ TBOX_ERROR(
+ d_object_name << ": Illegal boundary type in "
+ <<
+ "CellStokesHypreSolver::setMatrixCoefficients\n");
+ }
+ const hier::BoundaryBoxUtils bbu(boundary_box);
+ const int location_index = boundary_box.getLocationIndex();
+ const hier::BoundaryBox trimmed_boundary_box =
+ bbu.trimBoundaryBox(patch.getBox());
+ const hier::Box bccoef_box =
+ bbu.getSurfaceBoxFromBoundaryBox();
+ tbox::Pointer<pdat::ArrayData<double> >
+ acoef_data(new pdat::ArrayData<double>(bccoef_box, 1));
+ tbox::Pointer<pdat::ArrayData<double> >
+ bcoef_data(new pdat::ArrayData<double>(bccoef_box, 1));
+ tbox::Pointer<pdat::ArrayData<double> >
+ gcoef_data(NULL);
+ static const double fill_time = 0.0;
+ d_physical_bc_coef_strategy->setBcCoefs(acoef_data,
+ bcoef_data,
+ gcoef_data,
+ d_physical_bc_variable,
+ patch,
+ boundary_box,
+ fill_time);
+ pdat::ArrayData<double>& Ak0_data =
+ Ak0->getArrayData(location_index / 2,
+ location_index % 2);
+ adjustBoundaryEntries(diagonal,
+ off_diagonal,
+ patch_box,
+ *acoef_data,
+ *bcoef_data,
+ bccoef_box,
+ Ak0_data,
+ trimmed_boundary_box,
+ h);
+ }
+ }
+
+ /*
+ * Walk coarse-fine boundaries and adjust off-diagonals
+ * according data in ghost cells.
+ */
+ if (d_ln > 0) {
+ /*
+ * There are potentially coarse-fine boundaries to deal with.
+ */
+
+ tbox::Array<hier::BoundaryBox> surface_boxes;
+
+ if (d_dim == tbox::Dimension(2)) {
+ surface_boxes = d_cf_boundary->getEdgeBoundaries(pi->getGlobalId());
+ } else if (d_dim == tbox::Dimension(3)) {
+ surface_boxes = d_cf_boundary->getFaceBoundaries(pi->getGlobalId());
+ }
+
+ const int n_bdry_boxes = surface_boxes.getSize();
+ for (int n = 0; n < n_bdry_boxes; ++n) {
+
+ const hier::BoundaryBox& boundary_box = surface_boxes[n];
+ if (boundary_box.getBoundaryType() != 1) {
+ TBOX_ERROR(
+ d_object_name << ": Illegal boundary type in "
+ <<
+ "CellStokesHypreSolver::setMatrixCoefficients\n");
+ }
+ const int location_index = boundary_box.getLocationIndex();
+ const hier::BoundaryBoxUtils bbu(boundary_box);
+ const hier::BoundaryBox trimmed_boundary_box =
+ bbu.trimBoundaryBox(patch.getBox());
+ const hier::Box bccoef_box =
+ bbu.getSurfaceBoxFromBoundaryBox();
+ tbox::Pointer<pdat::ArrayData<double> >
+ acoef_data(new pdat::ArrayData<double>(bccoef_box, 1));
+ tbox::Pointer<pdat::ArrayData<double> >
+ bcoef_data(new pdat::ArrayData<double>(bccoef_box, 1));
+ tbox::Pointer<pdat::ArrayData<double> >
+ gcoef_data(NULL);
+ static const double fill_time = 0.0;
+ /*
+ * Reset invalid ghost data id to help detect use in setBcCoefs.
+ */
+ d_cf_bc_coef.setGhostDataId(-1, hier::IntVector::getZero(d_dim));
+ d_cf_bc_coef.setBcCoefs(acoef_data,
+ bcoef_data,
+ gcoef_data,
+ d_coarsefine_bc_variable,
+ patch,
+ boundary_box,
+ fill_time);
+ pdat::ArrayData<double>& Ak0_data =
+ Ak0->getArrayData(location_index / 2,
+ location_index % 2);
+ adjustBoundaryEntries(diagonal,
+ off_diagonal,
+ patch_box,
+ *acoef_data,
+ *bcoef_data,
+ bccoef_box,
+ Ak0_data,
+ trimmed_boundary_box,
+ h);
+ }
+ }
+
+ /*
+ * Copy matrix entries to HYPRE matrix structure. Note that
+ * we translate our temporary diagonal/off-diagonal storage into the
+ * HYPRE symmetric storage scheme for the stencil specified earlier.
+ */
+ const int stencil_size = d_dim.getValue() + 1;
+ int stencil_indices[stencil_size];
+ double mat_entries[stencil_size];
+
+ for (i = 0; i < stencil_size; i++) stencil_indices[i] = i;
+
+ pdat::CellIterator ic(patch_box);
+
+ /*
+ * To do: This loop uses inefficient high-level syntax.
+ * See if it can be replaced by a Fortran loop or if we
+ * can set matrix entries for an entire box at once.
+ */
+ for ( ; ic; ic++) {
+
+ hier::IntVector icell = ic();
+ pdat::SideIndex ixlower(ic(),
+ pdat::SideIndex::X,
+ pdat::SideIndex::Lower);
+ mat_entries[0] = (off_diagonal)(ixlower);
+
+ if (d_dim > tbox::Dimension(1)) {
+ pdat::SideIndex iylower(ic(),
+ pdat::SideIndex::Y,
+ pdat::SideIndex::Lower);
+ mat_entries[1] = (off_diagonal)(iylower);
+ }
+
+ if (d_dim > tbox::Dimension(2)) {
+ pdat::SideIndex izlower(ic(),
+ pdat::SideIndex::Z,
+ pdat::SideIndex::Lower);
+ // The "funny" indexing prevents a warning when compiling for
+ // DIM < 2. This code is only reached if DIM > 2 when
+ // executing.
+ mat_entries[d_dim.getValue() > 2 ? 2 : 0] = (off_diagonal)(izlower);
+ }
+
+ mat_entries[d_dim.getValue()] = (diagonal)(ic());
+ HYPRE_StructMatrixSetValues(d_matrix, &icell[0],
+ stencil_size, stencil_indices,
+ mat_entries);
+ } // end cell loop
+
+ } // end patch loop
+
+ if (d_print_solver_info) {
+ HYPRE_StructMatrixPrint("mat_bA.out", d_matrix, 1);
+ }
+
+ HYPRE_StructMatrixAssemble(d_matrix);
+
+ if (d_print_solver_info) {
+ HYPRE_StructMatrixPrint("mat_aA.out", d_matrix, 1);
+ }
+
+ t_set_matrix_coefficients->stop();
+
+ setupHypreSolver();
+}
+
+/*
+ **********************************************************************
+ * Add g*A*k0(a) from physical boundaries to rhs. *
+ * This operation is done for physical as well as cf boundaries, *
+ * so it is placed in a function. *
+ **********************************************************************
+ */
+
+void CellStokesHypreSolver::add_gAk0_toRhs(
+ const hier::Patch& patch,
+ const tbox::Array<hier::BoundaryBox>& bdry_boxes,
+ const RobinBcCoefStrategy* robin_bc_coef,
+ pdat::CellData<double>& rhs)
+{
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS2(d_dim, patch, rhs);
+
+ /*
+ * g*A*k0(a) is the storage for adjustments to be made to the rhs
+ * when we solve. This is the value of the weight of the ghost cell
+ * value for the interior cell, times k0. It is independent of u,
+ * and so is moved to the rhs. Before solving, g*A*k0(a) is added
+ * to rhs.
+ */
+ tbox::Pointer<pdat::OutersideData<double> > Ak0;
+
+ tbox::Pointer<geom::CartesianPatchGeometry> pg =
+ patch.getPatchGeometry();
+
+ Ak0 = patch.getPatchData(d_Ak0_id);
+
+ const int n_bdry_boxes = bdry_boxes.getSize();
+ for (int n = 0; n < n_bdry_boxes; ++n) {
+
+ const hier::BoundaryBox& boundary_box = bdry_boxes[n];
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (boundary_box.getBoundaryType() != 1) {
+ TBOX_ERROR(d_object_name << ": Illegal boundary type in "
+ << "CellStokesHypreSolver::add_gAk0_toRhs\n");
+ }
+#endif
+ const int location_index = boundary_box.getLocationIndex();
+ const hier::BoundaryBoxUtils bbu(boundary_box);
+ const hier::BoundaryBox trimmed_boundary_box =
+ bbu.trimBoundaryBox(patch.getBox());
+ const hier::Index& lower = trimmed_boundary_box.getBox().lower();
+ const hier::Index& upper = trimmed_boundary_box.getBox().upper();
+ const hier::Box& rhsbox = rhs.getArrayData().getBox();
+ const hier::Box& Ak0box = Ak0->getArrayData(location_index / 2,
+ location_index % 2).getBox();
+ const hier::Box bccoef_box = bbu.getSurfaceBoxFromBoundaryBox();
+ tbox::Pointer<pdat::ArrayData<double> >
+ acoef_data(NULL);
+ tbox::Pointer<pdat::ArrayData<double> >
+ bcoef_data(NULL);
+ tbox::Pointer<pdat::ArrayData<double> >
+ gcoef_data(new pdat::ArrayData<double>(bccoef_box, 1));
+ static const double fill_time = 0.0;
+ robin_bc_coef->setBcCoefs(acoef_data,
+ bcoef_data,
+ gcoef_data,
+ d_physical_bc_variable,
+ patch,
+ boundary_box,
+ fill_time);
+ /*
+ * Nomenclature for indices: cel=first-cell, gho=ghost,
+ * beg=beginning, end=ending.
+ */
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(adjustrhs2d, ADJUSTRHS2D) (rhs.getPointer(d_rhs_depth),
+ &rhsbox.lower()[0],
+ &rhsbox.upper()[0],
+ &rhsbox.lower()[1],
+ &rhsbox.upper()[1],
+ Ak0->getPointer(location_index / 2, location_index % 2),
+ &Ak0box.lower()[0],
+ &Ak0box.upper()[0],
+ &Ak0box.lower()[1],
+ &Ak0box.upper()[1],
+ gcoef_data->getPointer(),
+ &bccoef_box.lower()[0],
+ &bccoef_box.upper()[0],
+ &bccoef_box.lower()[1],
+ &bccoef_box.upper()[1],
+ &lower[0], &upper[0],
+ &location_index);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(adjustrhs3d, ADJUSTRHS3D) (rhs.getPointer(d_rhs_depth),
+ &rhsbox.lower()[0],
+ &rhsbox.upper()[0],
+ &rhsbox.lower()[1],
+ &rhsbox.upper()[1],
+ &rhsbox.lower()[2],
+ &rhsbox.upper()[2],
+ Ak0->getPointer(location_index / 2, location_index % 2),
+ &Ak0box.lower()[0],
+ &Ak0box.upper()[0],
+ &Ak0box.lower()[1],
+ &Ak0box.upper()[1],
+ &Ak0box.lower()[2],
+ &Ak0box.upper()[2],
+ gcoef_data->getPointer(),
+ &bccoef_box.lower()[0],
+ &bccoef_box.upper()[0],
+ &bccoef_box.lower()[1],
+ &bccoef_box.upper()[1],
+ &bccoef_box.lower()[2],
+ &bccoef_box.upper()[2],
+ &lower[0], &upper[0],
+ &location_index);
+ }
+ }
+}
+
+/*
+ *************************************************************************
+ * Create the hypre solver and set it according to the current state. *
+ *************************************************************************
+ */
+void CellStokesHypreSolver::setupHypreSolver()
+{
+ TBOX_ASSERT(d_mg_data == NULL);
+
+ tbox::SAMRAI_MPI::Comm communicator = d_hierarchy->getDomainMappedBoxLevel().getMPI().getCommunicator();
+
+ if (d_use_smg) {
+ HYPRE_StructSMGCreate(communicator, &d_mg_data);
+ HYPRE_StructSMGSetMemoryUse(d_mg_data, 0);
+ HYPRE_StructSMGSetMaxIter(d_mg_data, d_max_iterations);
+ HYPRE_StructSMGSetTol(d_mg_data, d_relative_residual_tol);
+ HYPRE_StructSMGSetLogging(d_mg_data, 1);
+ HYPRE_StructSMGSetNumPreRelax(d_mg_data,
+ d_num_pre_relax_steps);
+ HYPRE_StructSMGSetNumPostRelax(d_mg_data,
+ d_num_post_relax_steps);
+ HYPRE_StructSMGSetup(d_mg_data,
+ d_matrix,
+ d_linear_rhs,
+ d_linear_sol);
+ } else {
+ HYPRE_StructPFMGCreate(communicator, &d_mg_data);
+ HYPRE_StructPFMGSetMaxIter(d_mg_data, d_max_iterations);
+ HYPRE_StructPFMGSetTol(d_mg_data, d_relative_residual_tol);
+ HYPRE_StructPFMGSetLogging(d_mg_data, 1);
+ HYPRE_StructPFMGSetNumPreRelax(d_mg_data,
+ d_num_pre_relax_steps);
+ HYPRE_StructPFMGSetNumPostRelax(d_mg_data,
+ d_num_post_relax_steps);
+ HYPRE_StructPFMGSetup(d_mg_data,
+ d_matrix,
+ d_linear_rhs,
+ d_linear_sol);
+ }
+}
+
+void CellStokesHypreSolver::destroyHypreSolver()
+{
+ if (d_mg_data != NULL) {
+ if (d_use_smg) {
+ HYPRE_StructSMGDestroy(d_mg_data);
+ } else {
+ HYPRE_StructPFMGDestroy(d_mg_data);
+ }
+ d_mg_data = NULL;
+ }
+}
+
+/*
+ *************************************************************************
+ * *
+ * Solve the linear system. This routine assumes that the boundary *
+ * conditions and the matrix coefficients have been specified. *
+ * *
+ *************************************************************************
+ */
+
+int CellStokesHypreSolver::solveSystem(
+ const int u,
+ const int f,
+ bool homogeneous_bc)
+{
+ if (d_physical_bc_coef_strategy == NULL) {
+ TBOX_ERROR(
+ d_object_name << ": No BC coefficient strategy object!\n"
+ <<
+ "Use either setBoundaries or setPhysicalBcCoefObject\n"
+ <<
+ "to specify the boundary conidition. Do it before\n"
+ << "calling solveSystem.");
+ }
+ // Tracer t("CellStokesHypreSolver::solveSystem");
+
+ t_solve_system->start();
+
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(d_ln);
+#ifdef DEBUG_CHECK_ASSERTIONS
+ TBOX_ASSERT(u >= 0);
+ TBOX_ASSERT(
+ u < level->getPatchDescriptor()->getMaxNumberRegisteredComponents());
+ TBOX_ASSERT(f >= 0);
+ TBOX_ASSERT(
+ f < level->getPatchDescriptor()->getMaxNumberRegisteredComponents());
+#endif
+
+ if (d_physical_bc_coef_strategy == &d_physical_bc_simple_case) {
+ /*
+ * If we are using the simple bc implementation, the final piece
+ * of information it requires is the Dirichlet boundary value
+ * set in the ghost cells. Now that we have the ghost cell data,
+ * we can complete the boundary condition setup.
+ */
+ d_physical_bc_simple_case.cacheDirichletData(u);
+ }
+
+ /*
+ * Modify right-hand-side to account for boundary conditions and
+ * copy solution and right-hand-side to HYPRE structures.
+ */
+
+ const hier::IntVector no_ghosts(d_dim, 0);
+ const hier::IntVector ghosts(d_dim, 1);
+
+ /*
+ * At coarse-fine boundaries, we expect ghost cells to have correct
+ * values to be used in our bc, so u provides the ghost cell data.
+ * Assume that the user only provided data for the immediate first
+ * ghost cell, so pass zero for the number of extensions fillable.
+ */
+ d_cf_bc_coef.setGhostDataId(u, hier::IntVector::getZero(d_dim));
+
+ for (hier::PatchLevel::Iterator p(level); p; p++) {
+ tbox::Pointer<hier::Patch> patch = *p;
+
+ const hier::Box box = patch->getBox();
+
+ /*
+ * Set up variable data needed to prepare linear system solver.
+ */
+ tbox::Pointer<pdat::CellData<double> > u_data_ = patch->getPatchData(u);
+#ifdef DEBUG_CHECK_ASSERTIONS
+ TBOX_ASSERT(!u_data_.isNull());
+#endif
+ pdat::CellData<double>& u_data = *u_data_;
+ pdat::CellData<double> rhs_data(box, 1, no_ghosts);
+
+ /*
+ * Copy rhs and solution from the hierarchy into HYPRE structures.
+ * For rhs, add in the contribution from boundary conditions, if
+ * needed. If boundary condition is homogenous, this only adds
+ * zero, so we skip it.
+ */
+ copyToHypre(d_linear_sol, u_data, d_soln_depth, box);
+ rhs_data.copy(*(patch->getPatchData(f)));
+ if (!homogeneous_bc) {
+ /*
+ * Add g*A*k0(a) from physical and coarse-fine boundaries to rhs.
+ */
+ add_gAk0_toRhs(*patch,
+ patch->getPatchGeometry()->getCodimensionBoundaries(1),
+ d_physical_bc_coef_strategy,
+ rhs_data);
+ add_gAk0_toRhs(*patch,
+ d_cf_boundary->getBoundaries(patch->getGlobalId(), 1),
+ &d_cf_bc_coef,
+ rhs_data);
+ }
+ copyToHypre(d_linear_rhs, rhs_data, d_rhs_depth, box);
+
+ } // end patch loop
+
+ /*
+ * Reset invalid ghost data id to help detect erroneous further use.
+ */
+ d_cf_bc_coef.setGhostDataId(-1, hier::IntVector::getZero(d_dim));
+
+ /*
+ * Finish assembly of the vectors
+ */
+ HYPRE_StructVectorAssemble(d_linear_sol);
+
+ HYPRE_StructVectorAssemble(d_linear_rhs);
+
+ /*
+ * Solve the system - zero means convergence
+ * Solve takes the same arguments as Setup
+ */
+
+ if (d_print_solver_info) {
+ HYPRE_StructVectorPrint("sol0.out", d_linear_sol, 1);
+ HYPRE_StructMatrixPrint("mat0.out", d_matrix, 1);
+ HYPRE_StructVectorPrint("rhs.out", d_linear_rhs, 1);
+ }
+
+ if (d_use_smg) {
+ // HYPRE_StructSMGSetMaxIter(d_mg_data, d_max_iterations);
+ HYPRE_StructSMGSetTol(d_mg_data, d_relative_residual_tol);
+ /* converge = */ HYPRE_StructSMGSolve(d_mg_data,
+ d_matrix,
+ d_linear_rhs,
+ d_linear_sol);
+ } else {
+ // HYPRE_StructPFMGSetMaxIter(d_mg_data, d_max_iterations);
+ HYPRE_StructPFMGSetTol(d_mg_data, d_relative_residual_tol);
+ /* converge = */ HYPRE_StructPFMGSolve(d_mg_data,
+ d_matrix,
+ d_linear_rhs,
+ d_linear_sol);
+ }
+
+ if (d_print_solver_info) {
+ HYPRE_StructMatrixPrint("mat.out", d_matrix, 1);
+ HYPRE_StructVectorPrint("sol.out", d_linear_sol, 1);
+ }
+
+ if (d_use_smg) {
+ HYPRE_StructSMGGetNumIterations(d_mg_data,
+ &d_number_iterations);
+ HYPRE_StructSMGGetFinalRelativeResidualNorm(d_mg_data,
+ &d_relative_residual_norm);
+ } else {
+ HYPRE_StructPFMGGetNumIterations(d_mg_data,
+ &d_number_iterations);
+ HYPRE_StructPFMGGetFinalRelativeResidualNorm(d_mg_data,
+ &d_relative_residual_norm);
+ }
+
+ /*
+ * Pull the solution vector out of the HYPRE structures
+ */
+ for (hier::PatchLevel::Iterator ip(level); ip; ip++) {
+ tbox::Pointer<hier::Patch> patch = *ip;
+ tbox::Pointer<pdat::CellData<double> > u_data_ = patch->getPatchData(u);
+ pdat::CellData<double>& u_data = *u_data_;
+ copyFromHypre(u_data,
+ d_soln_depth,
+ d_linear_sol,
+ patch->getBox());
+ }
+
+ t_solve_system->stop();
+
+ return d_relative_residual_norm <= d_relative_residual_tol;
+}
+
+void CellStokesHypreSolver::computeDiagonalEntries(
+ pdat::CellData<double>& diagonal,
+ const pdat::CellData<double>& C_data,
+ const pdat::SideData<double>& off_diagonal,
+ const hier::Box& patch_box)
+{
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(d_dim,
+ diagonal,
+ C_data,
+ off_diagonal,
+ patch_box);
+
+ const hier::Index patch_lo = patch_box.lower();
+ const hier::Index patch_up = patch_box.upper();
+ const double c = 1.0, d = 1.0;
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(compdiagvariablec2d, COMPDIAGVARIABLEC2D) (diagonal.getPointer(),
+ C_data.getPointer(),
+ off_diagonal.getPointer(0),
+ off_diagonal.getPointer(1),
+ &patch_lo[0], &patch_up[0],
+ &patch_lo[1], &patch_up[1],
+ &c, &d);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(compdiagvariablec3d, COMPDIAGVARIABLEC3D) (diagonal.getPointer(),
+ C_data.getPointer(),
+ off_diagonal.getPointer(0),
+ off_diagonal.getPointer(1),
+ off_diagonal.getPointer(2),
+ &patch_lo[0], &patch_up[0],
+ &patch_lo[1], &patch_up[1],
+ &patch_lo[2], &patch_up[2],
+ &c, &d);
+ }
+}
+
+void CellStokesHypreSolver::computeDiagonalEntries(
+ pdat::CellData<double>& diagonal,
+ const double C,
+ const pdat::SideData<double>& off_diagonal,
+ const hier::Box& patch_box)
+{
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS3(d_dim, diagonal, off_diagonal, patch_box);
+
+ const hier::Index patch_lo = patch_box.lower();
+ const hier::Index patch_up = patch_box.upper();
+ const double c = 1.0, d = 1.0;
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(compdiagscalarc2d, COMPDIAGSCALARC2D) (diagonal.getPointer(),
+ &C,
+ off_diagonal.getPointer(0),
+ off_diagonal.getPointer(1),
+ &patch_lo[0], &patch_up[0],
+ &patch_lo[1], &patch_up[1],
+ &c, &d);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(compdiagscalarc3d, COMPDIAGSCALARC3D) (diagonal.getPointer(),
+ &C,
+ off_diagonal.getPointer(0),
+ off_diagonal.getPointer(1),
+ off_diagonal.getPointer(2),
+ &patch_lo[0], &patch_up[0],
+ &patch_lo[1], &patch_up[1],
+ &patch_lo[2], &patch_up[2],
+ &c, &d);
+ } else {
+ TBOX_ERROR("CellStokesHypreSolver error...\n"
+ << "DIM > 3 not supported." << std::endl);
+ }
+}
+
+void CellStokesHypreSolver::computeDiagonalEntries(
+ pdat::CellData<double>& diagonal,
+ const pdat::SideData<double>& off_diagonal,
+ const hier::Box& patch_box)
+{
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS3(d_dim, diagonal, off_diagonal, patch_box);
+
+ const hier::Index patch_lo = patch_box.lower();
+ const hier::Index patch_up = patch_box.upper();
+ const double c = 1.0, d = 1.0;
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(compdiagzeroc2d, COMPDIAGZEROC2D) (diagonal.getPointer(),
+ off_diagonal.getPointer(0),
+ off_diagonal.getPointer(1),
+ &patch_lo[0], &patch_up[0],
+ &patch_lo[1], &patch_up[1],
+ &c, &d);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(compdiagzeroc3d, COMPDIAGZEROC3D) (diagonal.getPointer(),
+ off_diagonal.getPointer(0),
+ off_diagonal.getPointer(1),
+ off_diagonal.getPointer(2),
+ &patch_lo[0], &patch_up[0],
+ &patch_lo[1], &patch_up[1],
+ &patch_lo[2], &patch_up[2],
+ &c, &d);
+ } else {
+ TBOX_ERROR("CellStokesHypreSolver error...\n"
+ << "DIM > 3 not supported." << std::endl);
+ }
+}
+
+void CellStokesHypreSolver::adjustBoundaryEntries(
+ pdat::CellData<double>& diagonal,
+ const pdat::SideData<double>& off_diagonal,
+ const hier::Box& patch_box,
+ const pdat::ArrayData<double>& acoef_data,
+ const pdat::ArrayData<double>& bcoef_data,
+ const hier::Box bccoef_box,
+ pdat::ArrayData<double>& Ak0_data,
+ const hier::BoundaryBox& trimmed_boundary_box,
+ const double h[tbox::Dimension::MAXIMUM_DIMENSION_VALUE])
+{
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS8(d_dim, diagonal, off_diagonal, patch_box,
+ acoef_data, bcoef_data,
+ bccoef_box, Ak0_data, trimmed_boundary_box);
+
+ const hier::Index patch_lo = patch_box.lower();
+ const hier::Index patch_up = patch_box.upper();
+ const int location_index = trimmed_boundary_box.getLocationIndex();
+ const hier::Index& lower = trimmed_boundary_box.getBox().lower();
+ const hier::Index& upper = trimmed_boundary_box.getBox().upper();
+ const hier::Box& Ak0_box = Ak0_data.getBox();
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(adjbdry2d, ADJBDRY2D) (diagonal.getPointer(),
+ off_diagonal.getPointer(0),
+ off_diagonal.getPointer(1),
+ &patch_lo[0], &patch_up[0],
+ &patch_lo[1], &patch_up[1],
+ acoef_data.getPointer(),
+ bcoef_data.getPointer(),
+ &bccoef_box.lower()[0],
+ &bccoef_box.upper()[0],
+ &bccoef_box.lower()[1],
+ &bccoef_box.upper()[1],
+ Ak0_data.getPointer(),
+ &Ak0_box.lower()[0],
+ &Ak0_box.upper()[0],
+ &Ak0_box.lower()[1],
+ &Ak0_box.upper()[1],
+ &lower[0], &upper[0],
+ &location_index, h);
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(adjbdry3d, ADJBDRY3D) (diagonal.getPointer(),
+ off_diagonal.getPointer(0),
+ off_diagonal.getPointer(1),
+ off_diagonal.getPointer(2),
+ &patch_lo[0], &patch_up[0],
+ &patch_lo[1], &patch_up[1],
+ &patch_lo[2], &patch_up[2],
+ acoef_data.getPointer(),
+ bcoef_data.getPointer(),
+ &bccoef_box.lower()[0],
+ &bccoef_box.upper()[0],
+ &bccoef_box.lower()[1],
+ &bccoef_box.upper()[1],
+ &bccoef_box.lower()[2],
+ &bccoef_box.upper()[2],
+ Ak0_data.getPointer(),
+ &Ak0_box.lower()[0],
+ &Ak0_box.upper()[0],
+ &Ak0_box.lower()[1],
+ &Ak0_box.upper()[1],
+ &Ak0_box.lower()[2],
+ &Ak0_box.upper()[2],
+ &lower[0], &upper[0],
+ &location_index, h);
+ } else {
+ TBOX_ERROR("CellStokesHypreSolver error...\n"
+ << "DIM > 3 not supported." << std::endl);
+ }
+}
+
+void
+CellStokesHypreSolver::finalizeCallback()
+{
+ for (int d = 0; d < tbox::Dimension::MAXIMUM_DIMENSION_VALUE; ++d) {
+ s_Ak0_var[d].setNull();
+ }
+}
+
+}
+}
+
+#endif
+#endif
diff -r cd1d77fa073e -r fe9d63509c19 CellStokesHypreSolver.I
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/CellStokesHypreSolver.I Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,121 @@
+/*************************************************************************
+ *
+ * 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: Level solver for diffusion-like elliptic problems.
+ *
+ ************************************************************************/
+namespace SAMRAI {
+namespace solv {
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesHypreSolver::setSolnIdDepth(
+ const int depth) {
+ d_soln_depth = depth;
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesHypreSolver::setRhsIdDepth(
+ const int depth) {
+ d_rhs_depth = depth;
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesHypreSolver::setUseSMG(
+ bool use_smg) {
+ d_use_smg = use_smg;
+}
+
+/*
+ ********************************************************************
+ * Specify bc using the default internal bc coefficient object. *
+ * Clear up data supporting external bc coefficient setter. *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesHypreSolver::setBoundaries(
+ const std::string& boundary_type,
+ const int fluxes,
+ const int flags,
+ int* bdry_types)
+{
+ d_physical_bc_simple_case.setBoundaries(boundary_type,
+ fluxes,
+ flags,
+ bdry_types);
+ d_physical_bc_coef_strategy = &d_physical_bc_simple_case;
+ d_physical_bc_variable.setNull();
+}
+
+/*
+ ********************************************************************
+ * Set the physical boundary condition object. *
+ ********************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesHypreSolver::setPhysicalBcCoefObject(
+ const RobinBcCoefStrategy* physical_bc_coef_strategy,
+ const tbox::Pointer<hier::Variable> variable)
+{
+ d_physical_bc_coef_strategy = physical_bc_coef_strategy;
+ d_physical_bc_variable = variable;
+}
+
+SAMRAI_INLINE_KEYWORD
+int CellStokesHypreSolver::getNumberOfIterations() const
+{
+ return d_number_iterations;
+}
+
+SAMRAI_INLINE_KEYWORD
+double CellStokesHypreSolver::getRelativeResidualNorm() const
+{
+ return d_relative_residual_norm;
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesHypreSolver::setNumPreRelaxSteps(
+ const int steps)
+{
+#ifdef DEBUG_CHECK_ASSERTIONS
+ TBOX_ASSERT(!d_hierarchy.isNull());
+#endif
+ d_num_pre_relax_steps = steps;
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesHypreSolver::setNumPostRelaxSteps(
+ const int steps)
+{
+#ifdef DEBUG_CHECK_ASSERTIONS
+ TBOX_ASSERT(!d_hierarchy.isNull());
+#endif
+ d_num_post_relax_steps = steps;
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesHypreSolver::setPrintSolverInfo(
+ const bool print)
+{
+ d_print_solver_info = print;
+}
+
+SAMRAI_INLINE_KEYWORD
+void CellStokesHypreSolver::setStoppingCriteria(
+ const int max_iterations,
+ const double residual_tol)
+{
+#ifdef DEBUG_CHECK_ASSERTIONS
+ TBOX_ASSERT(max_iterations >= 0);
+ TBOX_ASSERT(residual_tol >= 0.0);
+#endif
+ d_max_iterations = max_iterations;
+ d_relative_residual_tol = residual_tol;
+}
+
+}
+}
diff -r cd1d77fa073e -r fe9d63509c19 CellStokesHypreSolver.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/CellStokesHypreSolver.h Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,617 @@
+/*************************************************************************
+ *
+ * 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: Hypre solver interface for diffusion-like elliptic problems.
+ *
+ ************************************************************************/
+#ifndef included_solv_CellStokesHypreSolver
+#define included_solv_CellStokesHypreSolver
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#ifdef HAVE_HYPRE
+
+#ifndef included_HYPRE_struct_ls
+/*
+ * This might break things if F77_FUNC_ is different for hypre vs
+ * SAMRAI autoconf detection. But then C/C++ macros are totally
+ * broken due to namespace collision as this example highlights so
+ * resorting to hacks are necessary.
+ */
+#ifdef F77_FUNC_
+#undef F77_FUNC_
+#endif
+#include "HYPRE_struct_ls.h"
+#define included_HYPRE_struct_ls
+#endif
+
+#include "SAMRAI/solv/GhostCellRobinBcCoefs.h"
+#include "SAMRAI/solv/RobinBcCoefStrategy.h"
+#include "StokesSpecifications.h"
+#include "SAMRAI/solv/SimpleCellRobinBcCoefs.h"
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/pdat/OutersideVariable.h"
+#include "SAMRAI/hier/BoxList.h"
+#include "SAMRAI/hier/CoarseFineBoundary.h"
+#include "SAMRAI/hier/PatchHierarchy.h"
+#include "SAMRAI/hier/PatchLevel.h"
+#include "SAMRAI/hier/VariableContext.h"
+#include "SAMRAI/tbox/Database.h"
+#include "SAMRAI/tbox/Pointer.h"
+
+#include <string>
+
+namespace SAMRAI {
+namespace solv {
+
+/*!
+ * @brief Use the HYPRE preconditioner library to solve (the cell-centered)
+ * Stokes's equation on a single level in a hierarchy.
+ *
+ * Class CellStokesHypreSolver uses the HYPRE preconditioner library
+ * to solve linear equations of the form
+ * @f$ \nabla ( D \nabla u ) + C u = f @f$, where
+ * C is a cell-centered array, D is a face-centered array,
+ * and u and f are cell-centered arrays
+ * (see StokesSpecifications).
+ * The discretization is the standard second order
+ * finite difference stencil.
+ *
+ * Robin boundary conditions are used through the
+ * interface class RobinBcCoefStrategy.
+ * Periodic boundary conditions are not supported yet.
+ *
+ * The user must perform the following steps to use
+ * CellStokesHypreSolver:
+ * - Create a CellStokesHypreSolver object.
+ * - Initialize CellStokesHypreSolver object with a patch hierarchy,
+ * using the function initializeSolverState().
+ * - Use the functions setPhysicalBcCoefObject()
+ * to provide implementations of RobinBcCoefStrategy.
+ * (For most problems you can probably find a suitable
+ * implementation to use without implementing the
+ * strategy yourself. See for example
+ * SimpleCellRobinBcCoefs and GhostCellRobinBcCoefs.)
+ * - Set the matrix coefficients in the linear system,
+ * using the function setMatrixCoefficients().
+ * - Specify the stopping criteria using setStoppingCriteria().
+ * - Solve the linear system, passing in u and f as the patch
+ * indices of the solution and the right hand side, respectively.
+ *
+ * Sample parameters for initialization from database (and their
+ * default values):
+ * @verbatim
+ * print_solver_info = FALSE // Whether to print some data for debugging
+ * max_iterations = 10 // Max iterations used by Hypre
+ * relative_residual_tol = 1.0e-8 // Residual tolerance used by Hypre
+ * num_pre_relax_steps = 1 // # of presmoothing steps used by Hypre
+ * num_post_relax_steps = 1 // # of postsmoothing steps used by Hypre
+ * use_smg = FALSE // Whether to use hypre's smg solver
+ * // (alternative is the pfmg solver)
+ * @endverbatim
+ */
+
+class CellStokesHypreSolver
+{
+public:
+ /*!
+ * @brief Constructor.
+ *
+ * @param object_name Name of object.
+ * @param database tbox::Database for input.
+ */
+ CellStokesHypreSolver(
+ const tbox::Dimension& dim,
+ const std::string& object_name,
+ tbox::Pointer<tbox::Database> database =
+ tbox::Pointer<tbox::Database>(NULL));
+
+ /*!
+ * The Stokes destructor releases all internally managed data.
+ */
+ ~CellStokesHypreSolver();
+
+ /*!
+ * @brief Initialize to a given hierarchy.
+ *
+ * Initializer Stokes solver for a patch level in a hierarchy.
+ *
+ * @param hierarchy Hierarchy
+ * @param ln Level number
+ */
+ void
+ initializeSolverState(
+ tbox::Pointer<hier::PatchHierarchy> hierarchy,
+ int ln = 0);
+
+ /*!
+ * @brief Reset to an uninitialized state.
+ */
+ void
+ deallocateSolverState();
+
+ /*!
+ * @brief Set the matrix coefficients
+ *
+ * For information describing the Stokes equation parameters,
+ * see the light-weight StokesSpecifications class where
+ * you set the values of C and D.
+ *
+ * This method must be called before solveSystem().
+ */
+ void
+ setMatrixCoefficients(
+ const StokesSpecifications& spec);
+
+ /*!
+ * @brief Set default depth of the solution data involved in the solve.
+ *
+ * If the solution data has multiple depths,
+ * the solver uses just one depth at a time.
+ * The default depth is the first depth.
+ * Use this function to change it.
+ * This is not used to set the depth of the data (which is not
+ * controled by this class) but the depth used in the solve.
+ *
+ * Changing the depth after setting up the matrix is permissible,
+ * as the solution data does not affect the matrix.
+ */
+ void
+ setSolnIdDepth(
+ const int depth);
+
+ /*!
+ * @brief Set default depth of the rhs data involved in the solve.
+ *
+ * If the rhs data has multiple depths,
+ * the solver uses just one depth at a time.
+ * The default depth is the first depth.
+ * Use this function to change it.
+ * This is not used to set the depth of the data (which is not
+ * controled by this class) but the depth used in the solve.
+ *
+ * Changing the depth after setting up the matrix is permissible,
+ * as the rhs data does not affect the matrix.
+ */
+ void
+ setRhsIdDepth(
+ const int depth);
+
+ /*!
+ * @brief Set the stopping criteria (max iterations and residual
+ * tolerance) for the linear solver.
+ *
+ * @param max_iterations gives the maximum number of iterations
+ * @param relative_residual_tol the maximum error tolerance
+ */
+ void
+ setStoppingCriteria(
+ const int max_iterations = 10,
+ const double relative_residual_tol = 1.0e-6);
+
+ /*!
+ * @brief Solve the linear system Au=f.
+ *
+ * The solution u and the right hand side f are
+ * specified via patch indices on the patch hierarchy.
+ *
+ * Member functions getNumberOfIterations() return the iterations
+ * from the solver.
+ * Note that the matrix coefficients and boundary condition object
+ * must have been set up before this function is called.
+ * As long as the matrix coefficients do not change,
+ * this routine may be called repeatedly to solve any number of linear
+ * systems (with the right-hand side varying).
+ * If the boundary conditions or matrix coefficients are changed
+ * then function setMatrixCoefficients() must be called again.
+ *
+ * When computing the matrix coefficients in setMatrixCoefficients(),
+ * the inhomogeneous portion of the boundary condition (constant
+ * terms, independent of u and thus having no effect on the matrix)
+ * are saved and added to the source term, f,
+ * before performing the matrix solve. In some situations, it may be
+ * useful to not add the inhomogeneous portion to f. The flag argument
+ * @c homoegneous_bc is used for this. (This is a sort of optimization,
+ * to avoid having to re-call setMatrixCoefficients() to change the
+ * inhomogeneous portion.)
+ *
+ * @param u Descriptor of cell-centered unknown variable.
+ * @param f Descriptor of cell-centered source variable.
+ * @param homogeneous_bc Whether homogeneous boundary conditions
+ * are assumed.
+ *
+ * @return whether solver converged to specified level
+ */
+ int
+ solveSystem(
+ const int u,
+ const int f,
+ bool homogeneous_bc = false);
+
+ /*!
+ * @brief Return the number of iterations taken by the solver to converge.
+ *
+ * @return number of iterations taken by the solver to converge
+ */
+ int
+ getNumberOfIterations() const;
+
+ /*!
+ * @brief Set the number of pre-relax steps used by the Hypre solve.
+ */
+ void
+ setNumPreRelaxSteps(
+ const int steps);
+
+ /*!
+ * @brief Set the number of post-relax steps used by the Hypre solve.
+ */
+ void
+ setNumPostRelaxSteps(
+ const int steps);
+
+ /*!
+ * @brief Return the final residual norm returned by the Hypre solve.
+ * @return final residual norm returned by the Hypre solve.
+ */
+ double
+ getRelativeResidualNorm() const;
+
+ /*!
+ * @brief Set whether to use Hypre's PFMG algorithm instead of the
+ * SMG algorithm.
+ *
+ * The flag is used to select which of HYPRE's linear solver algorithms
+ * to use if true, the semicoarsening multigrid algorithm is used, and if
+ * false, the "PF" multigrid algorithm is used.
+ * By default, the SMG algorithm is used.
+ *
+ * Changing the algorithm must be done before setting up the matrix
+ * coefficients.
+ */
+ void
+ setUseSMG(
+ bool use_smg);
+
+ /*!
+ * @brief Specify boundary condition directly, without using
+ * a RobinBcCoefStrategy object.
+ *
+ * Use @em either setBoundaries() @em or setPhysicalBcCoefObject(),
+ * but not both.
+ *
+ * A SimpleCelBcCoef object is used to interpret and implement
+ * the specified boundary conditions.
+ * See SimpleCellRobinBcCoefs::setBoundaries()
+ * for an explanation of the arguments.
+ */
+ void
+ setBoundaries(
+ const std::string& boundary_type,
+ const int fluxes = -1,
+ const int flags = -1,
+ int* bdry_types = NULL);
+
+ /*!
+ * @brief Specify boundary condition through the use of a
+ * Robin boundary condition object.
+ *
+ * Use @em either setBoundaries() @em or setPhysicalBcCoefObject(),
+ * but not both.
+ *
+ * The Robin boundary condition object is used when setting
+ * the matrix coefficient and when solving the system.
+ * If your boundary conditions are fixed values at ghost
+ * cell centers, use the GhostCellRobinBcCoefs
+ * implementation of the RobinBcCoefStrategy strategy.
+ *
+ * @param physical_bc_coef_strategy tbox::Pointer a concrete
+ * implementation of the Robin bc strategy.
+ * @param variable hier::Variable pointer to be passed
+ * to RobinBcCoefStrategy::setBcCoefs(),
+ * but otherwise unused by this class.
+ */
+ void
+ setPhysicalBcCoefObject(
+ const RobinBcCoefStrategy* physical_bc_coef_strategy,
+ const tbox::Pointer<hier::Variable> variable =
+ tbox::Pointer<hier::Variable>(NULL));
+
+ /*!
+ * @brief Set the flag for printing solver information.
+ *
+ * This optional function is used primarily for debugging.
+ *
+ * If set true, it will print the HYPRE matrix information
+ * to the following files:
+ *
+ * - mat_bA.out - before setting matrix coefficients in matrix assemble
+ * - mat_aA.out - after setting matrix coefficients in matrix assemble
+ * - sol0.out - u before solve (i.e. for system Au = b)
+ * - sol.out - u after solve
+ * - mat0.out - A before solve
+ * - mat.out - A after solve
+ * - rhs.out - b before and after solve
+ *
+ * If this method is not called, or the flag is set false, no printing
+ * will occur.
+ */
+ void
+ setPrintSolverInfo(
+ const bool print);
+
+private:
+ /*!
+ * @brief Set state using database
+ *
+ * See the class description for the parameters that can be set
+ * from a database.
+ *
+ * @param database Input database. If a NULL pointer is given,
+ * nothing is done.
+ */
+ void
+ getFromInput(
+ tbox::Pointer<tbox::Database> database);
+
+ void
+ setupHypreSolver();
+ void
+ destroyHypreSolver();
+ void
+ allocateHypreData();
+ void
+ deallocateHypreData();
+
+ void
+ copyToHypre(
+ HYPRE_StructVector vector,
+ pdat::CellData<double>& src,
+ int depth,
+ const hier::Box& box);
+ void
+ copyFromHypre(
+ pdat::CellData<double>& dst,
+ int depth,
+ HYPRE_StructVector vector,
+ const hier::Box box);
+
+ /*!
+ * @brief Add g*A*k0(a) from boundaries to rhs.
+ *
+ * Move the constant portion of the boundary condition
+ * contribution to the right hand side and add it to the existing rhs.
+ * This operation is done for physical as well as cf boundaries,
+ * so it is placed in a function.
+ *
+ * The boundary boxes given must be to either the physical
+ * boundary or coarse-fine boundary for the patch. The
+ * bc coefficient implementation should correspond to the
+ * boundary being worked on.
+ */
+ void
+ add_gAk0_toRhs(
+ const hier::Patch& patch,
+ const tbox::Array<hier::BoundaryBox>& bdry_boxes,
+ const RobinBcCoefStrategy* robin_bc_coef,
+ pdat::CellData<double>& rhs);
+
+ //@{
+
+ /*!
+ * @name Dimension-independent functions to organize Fortran interface.
+ */
+
+ //! @brief Compute diagonal entries of the matrix when C is variable.
+ void
+ computeDiagonalEntries(
+ pdat::CellData<double>& diagonal,
+ const pdat::CellData<double>& C_data,
+ const pdat::SideData<double>& variable_off_diagonal,
+ const hier::Box& patch_box);
+ //! @brief Compute diagonal entries of the matrix when C is constant.
+ void
+ computeDiagonalEntries(
+ pdat::CellData<double>& diagonal,
+ const double C,
+ const pdat::SideData<double>& variable_off_diagonal,
+ const hier::Box& patch_box);
+ //! @brief Compute diagonal entries of the matrix when C is zero.
+ void
+ computeDiagonalEntries(
+ pdat::CellData<double>& diagonal,
+ const pdat::SideData<double>& variable_off_diagonal,
+ const hier::Box& patch_box);
+ /*!
+ * @brief Adjust boundary entries for variable off-diagonals.
+ *
+ * At the same time, save information that are needed to adjust
+ * the rhs.
+ */
+ void
+ adjustBoundaryEntries(
+ pdat::CellData<double>& diagonal,
+ const pdat::SideData<double>& variable_off_diagonal,
+ const hier::Box& patch_box,
+ const pdat::ArrayData<double>& acoef_data,
+ const pdat::ArrayData<double>& bcoef_data,
+ const hier::Box bccoef_box,
+ pdat::ArrayData<double>& Ak0_data,
+ const hier::BoundaryBox& trimmed_boundary_box,
+ const double h[tbox::Dimension::MAXIMUM_DIMENSION_VALUE]);
+
+ //@}
+
+ //! @brief Free static variables at shutdown time.
+ static void
+ finalizeCallback();
+
+ /*!
+ * @brief Object dimension.
+ */
+ const tbox::Dimension d_dim;
+
+ /*!
+ * @brief Object name.
+ */
+ std::string d_object_name;
+
+ /*!
+ * @brief Associated hierarchy.
+ */
+ tbox::Pointer<hier::PatchHierarchy> d_hierarchy;
+
+ /*!
+ * @brief Associated level number.
+ *
+ * Currently, this must be level number 0.
+ */
+ int d_ln;
+
+ /*!
+ * @brief Scratch context for this object.
+ */
+ tbox::Pointer<hier::VariableContext> d_context;
+
+ //@{ @name Boundary condition handling
+
+ /*!
+ * @brief The coarse-fine boundary description for level d_ln.
+ *
+ * The coarse-fine boundary is computed when the operator
+ * state is initialized. It is used to allow solves on
+ * levels that are not the coarsest in the hierarchy.
+ */
+ tbox::Pointer<hier::CoarseFineBoundary> d_cf_boundary;
+
+ /*!
+ * @brief Robin boundary coefficient object for physical
+ * boundaries.
+ *
+ * If d_physical_bc_coef_strategy is set, use it, otherwise,
+ * use d_physical_bc_simple_case.
+ */
+ const RobinBcCoefStrategy* d_physical_bc_coef_strategy;
+ tbox::Pointer<hier::Variable> d_physical_bc_variable;
+
+ /*!
+ * @brief Implementation of Robin boundary conefficients
+ * for the case of simple boundary conditions.
+ */
+ SimpleCellRobinBcCoefs d_physical_bc_simple_case;
+
+ /*!
+ * @brief Robin boundary coefficient object for coarse-fine
+ * boundaries.
+ *
+ * This is a GhostCellRobinBcCoefs object because we
+ * expect the users to have the correct ghost cell values
+ * in the coarse-fine boundaries before solving.
+ */
+ GhostCellRobinBcCoefs d_cf_bc_coef;
+ tbox::Pointer<hier::Variable> d_coarsefine_bc_variable;
+
+ //@}
+
+ /*!
+ * @brief hier::Patch index of A*k0(a) quantity
+ *
+ * A*k0(a) is the quantity that is saved for
+ * later adding to the rhs.
+ *
+ * The Robin bc is expressed by the coefficients a and g
+ * on the boundary (see RobinBcCoefStrategy).
+ * This class uses a central difference approximation of
+ * the Robin bc, which results in the value at a ghost cell,
+ * uo, being writen as uo = g*k0(a) + k1(a)*ui, where ui is
+ * the first interior cell value, k0 and k1 depend on a as
+ * indicated.
+ *
+ * In setting up the Au=f system, the contribution of k1(a)*ui
+ * is incorporated into the product Au. The contribution of
+ * A*g*k0(a) should be moved to the right hand side and saved for
+ * later adding to f. However, the value of g is not provided
+ * until solve time. Therefore, we save just A*k0(a) at the
+ * patch data index d_Ak0_id.
+ */
+ int d_Ak0_id;
+
+ static tbox::Pointer<pdat::OutersideVariable<double> > s_Ak0_var[tbox::
+ Dimension::
+ MAXIMUM_DIMENSION_VALUE];
+
+ /*!
+ * @brief Depth of the solution variable.
+ */
+ int d_soln_depth;
+
+ /*!
+ * @brief Depth of the rhs variable.
+ */
+ int d_rhs_depth;
+
+ int d_max_iterations;
+ double d_relative_residual_tol;
+
+ int d_number_iterations; // iterations in solver
+ int d_num_pre_relax_steps; // pre-relax steps in solver
+ int d_num_post_relax_steps; // post-relax steps in solver
+ double d_relative_residual_norm; // norm from solver
+
+ /*@
+ * @brief Flag to use SMG or PFMG (default)
+ */
+ bool d_use_smg;
+
+ //@{
+ //! @name Hypre object
+ //! @brief HYPRE grid
+ HYPRE_StructGrid d_grid;
+ //! @brief HYPRE stencil
+ HYPRE_StructStencil d_stencil;
+ //! @brief HYPRE structured matrix
+ HYPRE_StructMatrix d_matrix;
+ //! @brief Hypre RHS vector for linear solves
+ HYPRE_StructVector d_linear_rhs;
+ //! @brief Hypre solution vector
+ HYPRE_StructVector d_linear_sol;
+ //! @brief Hypre SMG solver data
+ HYPRE_StructSolver d_mg_data;
+ //@}
+
+ //@{
+
+ //! @name Variables for debugging and analysis.
+
+ /*!
+ * @brief Flag to print solver info
+ *
+ * See setPrintSolverInfo().
+ */
+ bool d_print_solver_info;
+
+ //@}
+
+ /*!
+ * @brief Timers for performance measurement.
+ */
+ tbox::Pointer<tbox::Timer> t_solve_system;
+ tbox::Pointer<tbox::Timer> t_set_matrix_coefficients;
+
+ static tbox::StartupShutdownManager::Handler s_finalize_handler;
+};
+
+}
+} // namespace SAMRAI
+
+#ifdef SAMRAI_INLINE
+#include "CellStokesHypreSolver.I"
+#endif
+
+#endif
+
+#endif
diff -r cd1d77fa073e -r fe9d63509c19 FACStokes.C
--- a/FACStokes.C Thu Dec 30 20:01:24 2010 -0800
+++ b/FACStokes.C Fri Dec 31 07:12:28 2010 -0800
@@ -16,7 +16,7 @@
#include "SAMRAI/pdat/CellData.h"
#include "SAMRAI/math/HierarchyCellDataOpsReal.h"
#include "SAMRAI/pdat/SideData.h"
-#include "SAMRAI/solv/PoissonSpecifications.h"
+#include "StokesSpecifications.h"
#include "SAMRAI/tbox/Utilities.h"
#include "SAMRAI/hier/Variable.h"
#include "SAMRAI/hier/VariableDatabase.h"
@@ -59,7 +59,7 @@ FACStokes::FACStokes(
d_dim(dim),
d_hierarchy(NULL),
d_stokes_fac_solver((d_dim),
- object_name + "::poisson_hypre",
+ object_name + "::stokes_hypre",
(!database.isNull() &&
database->isDatabase("fac_solver")) ?
database->getDatabase("fac_solver"):
diff -r cd1d77fa073e -r fe9d63509c19 FACStokes.h
--- a/FACStokes.h Thu Dec 30 20:01:24 2010 -0800
+++ b/FACStokes.h Fri Dec 31 07:12:28 2010 -0800
@@ -10,7 +10,7 @@
#ifndef included_FACStokes
#define included_FACStokes
-#include "SAMRAI/solv/CellPoissonFACSolver.h"
+#include "CellStokesFACSolver.h"
#include "SAMRAI/pdat/CellVariable.h"
#include "SAMRAI/tbox/Database.h"
#include "SAMRAI/hier/Box.h"
@@ -175,7 +175,7 @@ private:
/*!
* @brief FAC stokes solver.
*/
- solv::CellPoissonFACSolver d_stokes_fac_solver;
+ solv::CellStokesFACSolver d_stokes_fac_solver;
/*!
* @brief Boundary condition coefficient implementation.
diff -r cd1d77fa073e -r fe9d63509c19 Makefile
--- a/Makefile Thu Dec 30 20:01:24 2010 -0800
+++ b/Makefile Fri Dec 31 07:12:28 2010 -0800
@@ -24,7 +24,8 @@ NUM_TESTS = 2
TEST_NPROCS = 0,2
-CXX_OBJS = main.o FACStokes.o
+CXX_OBJS = main.o FACStokes.o CellStokesFACOps.o \
+ CellStokesHypreSolver.o StokesSpecifications.o CellStokesFACSolver.o
F_OBJS = facpoisson2d.o facpoisson3d.o
main: $(CXX_OBJS) $(F_OBJS) $(LIBSAMRAIDEPEND)
diff -r cd1d77fa073e -r fe9d63509c19 StokesSpecifications.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesSpecifications.C Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,40 @@
+/*************************************************************************
+ *
+ * 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: Specifications for the scalar Stokes equation
+ *
+ ************************************************************************/
+#include "StokesSpecifications.h"
+
+#ifndef SAMRAI_INLINE
+#include "StokesSpecifications.I"
+#endif
+
+namespace SAMRAI {
+namespace solv {
+
+void StokesSpecifications::printClassData(
+ std::ostream& stream) const
+{
+ stream << "StokesSpecifications " << d_object_name << "\n"
+ << " D is ";
+ if (d_D_id != -1) {
+ stream << "variable with patch id " << d_D_id << "\n";
+ } else {
+ stream << "constant with value " << d_D_constant << "\n";
+ }
+ stream << " C is ";
+ if (d_C_zero) {
+ stream << "zero\n";
+ } else if (d_C_id != -1) {
+ stream << "variable with patch id " << d_C_id << "\n";
+ } else {
+ stream << "constant with value " << d_C_constant << "\n";
+ }
+}
+
+} // namespace solv
+} // namespace SAMRAI
diff -r cd1d77fa073e -r fe9d63509c19 StokesSpecifications.I
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesSpecifications.I Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,259 @@
+/*************************************************************************
+ *
+ * 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: Specifications for the scalar Stokes equation
+ *
+ ************************************************************************/
+#include "SAMRAI/tbox/Utilities.h"
+#include "SAMRAI/tbox/Pointer.h"
+
+namespace SAMRAI {
+namespace solv {
+
+/*
+ *******************************************************************
+ * Default constructor *
+ *******************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+StokesSpecifications::StokesSpecifications(
+ const std::string& object_name):d_object_name(object_name),
+ d_D_id(-1),
+ d_D_constant(1.0),
+ d_C_zero(true),
+ d_C_id(-1),
+ d_C_constant(0.0) {
+}
+
+/*
+ *******************************************************************
+ * Copy constructor *
+ *******************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+StokesSpecifications::StokesSpecifications(
+ const std::string& object_name,
+ const StokesSpecifications& r):d_object_name(object_name),
+ d_D_id(r.d_D_id),
+ d_D_constant(r.d_D_constant),
+ d_C_zero(r.d_C_zero),
+ d_C_id(r.d_C_id),
+ d_C_constant(r.d_C_constant) {
+}
+
+/*
+ *******************************************************************
+ * Destructor (does nothing). *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+StokesSpecifications::~StokesSpecifications() {
+}
+
+/*
+ *******************************************************************
+ * Assignment operator *
+ *******************************************************************
+ */
+
+SAMRAI_INLINE_KEYWORD
+const StokesSpecifications
+& StokesSpecifications::operator = (
+ const StokesSpecifications& r) {
+ d_D_id = r.d_D_id;
+ d_D_constant = r.d_D_constant;
+ d_C_zero = r.d_C_zero;
+ d_C_id = r.d_C_id;
+ d_C_constant = r.d_C_constant;
+ return *this;
+}
+
+/*
+ *******************************************************************
+ * Set the patch data index for variable D. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+void StokesSpecifications::setDPatchDataId(
+ int id) {
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (id < 0) {
+ TBOX_ERROR(d_object_name << ": Invalid patch data id.\n");
+ }
+#endif
+ d_D_id = id;
+ d_D_constant = 0.0;
+}
+
+/*
+ *******************************************************************
+ * Set the constant value variable D. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+void StokesSpecifications::setDConstant(
+ double constant) {
+ d_D_id = -1;
+ d_D_constant = constant;
+}
+
+/*
+ *******************************************************************
+ * Whether D is variable. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+bool StokesSpecifications::dIsVariable() const {
+ return d_D_id != -1;
+}
+
+/*
+ *******************************************************************
+ * Whether D is constant. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+bool StokesSpecifications::dIsConstant() const {
+ return d_D_id == -1;
+}
+
+/*
+ *******************************************************************
+ * Get the patch data index for variable D. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+int StokesSpecifications::getDPatchDataId() const {
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (d_D_id == -1) {
+ TBOX_ERROR(d_object_name << ": D not prepresented by a patch data.\n");
+ }
+#endif
+ return d_D_id;
+}
+
+/*
+ *******************************************************************
+ * Get the constant D value. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+double StokesSpecifications::getDConstant() const {
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (d_D_id != -1) {
+ TBOX_ERROR(d_object_name << ": D not prepresented by a constant.\n");
+ }
+#endif
+ return d_D_constant;
+}
+
+/*
+ *******************************************************************
+ * Set the constant value variable C. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+void StokesSpecifications::setCPatchDataId(
+ int id) {
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (id < 0) {
+ TBOX_ERROR(d_object_name << ": Invalid patch data id.\n");
+ }
+#endif
+ d_C_zero = false;
+ d_C_id = id;
+ d_C_constant = 0.0;
+}
+
+/*
+ *******************************************************************
+ * Set the patch data index for variable C. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+void StokesSpecifications::setCConstant(
+ double constant) {
+ d_C_zero = false;
+ d_C_id = -1;
+ d_C_constant = constant;
+}
+
+/*
+ *******************************************************************
+ * Set the value of C to zero. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+void StokesSpecifications::setCZero() {
+ d_C_zero = true;
+ d_C_id = -1;
+ d_C_constant = 0.0;
+}
+
+/*
+ *******************************************************************
+ * Whether C is variable. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+bool StokesSpecifications::cIsVariable() const {
+ return d_C_id != -1;
+}
+
+/*
+ *******************************************************************
+ * Whether C is zero. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+bool StokesSpecifications::cIsZero() const {
+ return d_C_zero;
+}
+
+/*
+ *******************************************************************
+ * Whether C is constant. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+bool StokesSpecifications::cIsConstant() const {
+ return !d_C_zero && (d_C_id == -1);
+}
+
+/*
+ *******************************************************************
+ * Get the patch data index for variable C. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+int StokesSpecifications::getCPatchDataId() const {
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (d_C_id == -1) {
+ TBOX_ERROR(d_object_name << ": C not prepresented by a an index.\n");
+ }
+#endif
+ return d_C_id;
+}
+
+/*
+ *******************************************************************
+ * Get the constant C value. *
+ *******************************************************************
+ */
+SAMRAI_INLINE_KEYWORD
+double StokesSpecifications::getCConstant() const {
+#ifdef DEBUG_CHECK_ASSERTIONS
+ if (d_C_id != -1 || d_C_zero) {
+ TBOX_ERROR(d_object_name << ": C is not prepresented by a constant.\n");
+ }
+#endif
+ return d_C_constant;
+}
+
+} // namespace SAMRAI
+}
diff -r cd1d77fa073e -r fe9d63509c19 StokesSpecifications.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesSpecifications.h Fri Dec 31 07:12:28 2010 -0800
@@ -0,0 +1,256 @@
+/*************************************************************************
+ *
+ * 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: Specifications for the scalar Stokes equation
+ *
+ ************************************************************************/
+#ifndef included_solv_StokesSpecifications
+#define included_solv_StokesSpecifications
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#include <string>
+
+namespace SAMRAI {
+namespace solv {
+
+/*!
+ * @brief Light class holding specifications for cell-centered
+ * implementation of the scalar Stokes equation.
+ *
+ * The scalar Stokes equation is
+ * @f$ \nabla ( D \nabla u ) + C u = f @f$,
+ * where C is a scalar field, D is the diffusion coefficient.
+ * and u and f are scalar quantities.
+ *
+ * This class describes the things you can set: C, D.
+ *
+ * Note that the storage and alignment of u, f, C and D depend
+ * on the implementation of the solver. For example, if the
+ * solver is cell centered, u, f and C are cell-centered while
+ * D is side-centered.
+ */
+
+class StokesSpecifications
+{
+public:
+ /*!
+ * @brief Constructor.
+ *
+ * Sets the specifications to their default state:
+ * - C is zero
+ * - D is uniformly 1
+ *
+ * @param object_name Name of object.
+ */
+ StokesSpecifications(
+ const std::string& object_name);
+
+ /*!
+ * @brief Copy constructor.
+ */
+ StokesSpecifications(
+ const std::string& object_name,
+ const StokesSpecifications& r);
+
+ /*!
+ * @brief Destructor (does nothing).
+ */
+ virtual ~StokesSpecifications();
+
+ /*!
+ * @brief Assignment operator
+ *
+ * Assign everything except name.
+ */
+ const StokesSpecifications&
+ operator = (
+ const StokesSpecifications& r);
+
+ /*!
+ * @brief Print out class data.
+ */
+ virtual void
+ printClassData(
+ std::ostream& stream) const;
+
+ //@{
+ //! @name Functions for setting and getting D
+
+ /*!
+ * @brief Set the patch data index for variable D.
+ *
+ * In addition, disregard any previous value
+ * specified by setDConstant().
+ */
+ void
+ setDPatchDataId(
+ int id);
+
+ /*!
+ * @brief Set the constant value variable D.
+ *
+ * In addition, disregard any previous patch data index
+ * specified by setDPatchDataId().
+ */
+ void
+ setDConstant(
+ double constant);
+
+ /*!
+ * @brief Whether D is variable (described by a patch data id).
+ *
+ * @return True if D is variable, described by the patch data
+ * id given in setCPatchDataId().
+ */
+ bool
+ dIsVariable() const;
+
+ /*!
+ * @brief Whether D is constant.
+ *
+ * @return True if D is constant, as specified by setCConstant().
+ */
+ bool
+ dIsConstant() const;
+
+ /*!
+ * @brief Get D's patch data id
+ *
+ * Error if D is not represented by a patch data id.
+ *
+ * @return D's id
+ */
+ int
+ getDPatchDataId() const;
+
+ /*!
+ * @brief Get D constant value
+ *
+ * Error if D is not represented by a constant.
+ *
+ * @return D's constant value
+ */
+ double
+ getDConstant() const;
+
+ //@}
+
+ //@{
+ //! @name Functions for setting and getting C
+
+ /*!
+ * @brief Set the patch data index for C.
+ *
+ * In addition, disregard any previous values
+ * specified by setCConstant() or setCZero().
+ */
+ void
+ setCPatchDataId(
+ int id);
+
+ /*!
+ * @brief Set C to a constant.
+ *
+ * In addition, disregard any previous value
+ * specified by setCPatchDataId() or setCZero().
+ *
+ * If you want to set C to zero, use setCZero() instead.
+ * This allows solvers to take advantage of fact C is absent.
+ */
+ void
+ setCConstant(
+ double constant);
+
+ /*!
+ * @brief Set the value of C to zero.
+ *
+ * In addition, disregard any previous patch data index
+ * specified by setCPatchDataId() and any previous constant
+ * specified by setCConstant().
+ */
+ void
+ setCZero();
+
+ /*!
+ * @brief Whether C is variable (described by a patch data id).
+ *
+ * @return True if C is variable, described by the patch data
+ * id given in setCPatchDataId().
+ */
+ bool
+ cIsVariable() const;
+
+ /*!
+ * @brief Whether C is zero.
+ *
+ * As it pertains to what this function returns,
+ * C is zero @em only by calling setCZero().
+ * Calling setCConstant() does @em not make C zero,
+ * even if you pass in the value of zero.
+ *
+ * @return True if C is exactly zero, as set by setCZero().
+ */
+ bool
+ cIsZero() const;
+
+ /*!
+ * @brief Whether C is constant.
+ *
+ * As it pertains to what this function returns,
+ * C is constant @em only by calling setCConstant().
+ * Calling setCZero() does @em not make C a constant.
+ *
+ * @return True if C is constant, as specified by setCConstant().
+ */
+ bool
+ cIsConstant() const;
+
+ /*!
+ * @brief Get C's patch data id
+ *
+ * Error if C is not represented by a patch data id.
+ *
+ * @return C's patch data id
+ */
+ int
+ getCPatchDataId() const;
+
+ /*!
+ * @brief Get C as a constant value.
+ *
+ * Error if C is not represented by a constant.
+ *
+ * @return C's constant value
+ */
+ double
+ getCConstant() const;
+
+ //@}
+
+private:
+ /*!
+ * @brief Object name.
+ */
+ std::string d_object_name;
+
+ int d_D_id;
+ double d_D_constant;
+
+ bool d_C_zero;
+ int d_C_id;
+ double d_C_constant;
+
+};
+
+} // namespace SAMRAI
+}
+
+#ifdef SAMRAI_INLINE
+#include "StokesSpecifications.I"
+#endif
+
+#endif
More information about the CIG-COMMITS
mailing list