[cig-commits] commit: Rename Poisson to Stokes
Mercurial
hg at geodynamics.org
Fri Feb 25 14:12:18 PST 2011
changeset: 4:cd1d77fa073e
user: Walter Landry <wlandry at caltech.edu>
date: Thu Dec 30 20:01:24 2010 -0800
files: FACPoisson.C FACPoisson.h FACStokes.C FACStokes.h Makefile Makefile.depend example_inputs/const_refine.2d.input main.C
description:
Rename Poisson to Stokes
diff -r 390def91ff57 -r cd1d77fa073e FACPoisson.C
--- a/FACPoisson.C Wed Dec 29 16:12:07 2010 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,378 +0,0 @@
-/*************************************************************************
- *
- * This file is part of the SAMRAI distribution. For full copyright
- * information, see COPYRIGHT and COPYING.LESSER.
- *
- * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
- * Description: Numerical routines for example FAC Poisson solver
- *
- ************************************************************************/
-#include "FACPoisson.h"
-
-#include "SAMRAI/hier/IntVector.h"
-#include "SAMRAI/geom/CartesianGridGeometry.h"
-#include "SAMRAI/geom/CartesianPatchGeometry.h"
-#include "SAMRAI/solv/SimpleCellRobinBcCoefs.h"
-#include "SAMRAI/pdat/CellData.h"
-#include "SAMRAI/math/HierarchyCellDataOpsReal.h"
-#include "SAMRAI/pdat/SideData.h"
-#include "SAMRAI/solv/PoissonSpecifications.h"
-#include "SAMRAI/tbox/Utilities.h"
-#include "SAMRAI/hier/Variable.h"
-#include "SAMRAI/hier/VariableDatabase.h"
-
-extern "C" {
-void F77_FUNC(setexactandrhs2d, SETEXACTANDRHS2D) (const int & ifirst0,
- const int & ilast0,
- const int & ifirst1,
- const int & ilast1,
- double* exact,
- double* rhs,
- const double* dx,
- const double* xlower);
-
-void F77_FUNC(setexactandrhs3d, SETEXACTANDRHS3D) (const int & ifirst0,
- const int & ilast0,
- const int & ifirst1,
- const int & ilast1,
- const int & ifirst2,
- const int & ilast2,
- double* exact,
- double* rhs,
- const double* dx,
- const double* xlower);
-}
-
-namespace SAMRAI {
-
-/*
- *************************************************************************
- * Constructor creates a unique context for the object and register *
- * all its internal variables with the variable database. *
- *************************************************************************
- */
-FACPoisson::FACPoisson(
- const std::string& object_name,
- const tbox::Dimension& dim,
- tbox::Pointer<tbox::Database> database):
- d_object_name(object_name),
- d_dim(dim),
- d_hierarchy(NULL),
- d_poisson_fac_solver((d_dim),
- object_name + "::poisson_hypre",
- (!database.isNull() &&
- database->isDatabase("fac_solver")) ?
- database->getDatabase("fac_solver"):
- tbox::Pointer<tbox::Database>(NULL)),
- d_bc_coefs(d_dim,
- object_name + "::bc_coefs",
- (!database.isNull() &&
- database->isDatabase("bc_coefs")) ?
- database->getDatabase("bc_coefs"):
- tbox::Pointer<tbox::Database>(NULL)),
- d_context()
-{
-
- hier::VariableDatabase* vdb =
- hier::VariableDatabase::getDatabase();
-
- /*
- * Get a unique context for variables owned by this object.
- */
- d_context = vdb->getContext(d_object_name + ":Context");
-
- /*
- * Register variables with hier::VariableDatabase
- * and get the descriptor indices for those variables.
- */
-
- tbox::Pointer<pdat::CellVariable<double> > comp_soln(
- new pdat::CellVariable<double>(dim, object_name + ":computed solution", 1));
- d_comp_soln_id =
- vdb->registerVariableAndContext(
- comp_soln,
- d_context,
- hier::IntVector(dim, 1) /* ghost cell width is 1 for stencil widths */);
-
- tbox::Pointer<pdat::CellVariable<double> > exact_solution(
- new pdat::CellVariable<double>(dim, object_name + ":exact solution"));
- d_exact_id =
- vdb->registerVariableAndContext(
- exact_solution,
- d_context,
- hier::IntVector(dim, 1) /* ghost cell width is 1 in case needed */);
-
- tbox::Pointer<pdat::CellVariable<double> > rhs_variable(
- new pdat::CellVariable<double>(
- dim,
- object_name
- + ":linear system right hand side"));
- d_rhs_id =
- vdb->registerVariableAndContext(
- rhs_variable,
- d_context,
- hier::IntVector(dim, 0) /* ghost cell width is 0 */);
-
- /*
- * Specify an implementation of solv::RobinBcCoefStrategy for the solver to use.
- * We use the implementation solv::LocationIndexRobinBcCoefs, but other
- * implementations are possible, including user-implemented.
- */
- d_poisson_fac_solver.setBcObject(&d_bc_coefs);
-}
-
-/*
- *************************************************************************
- * Destructor does nothing interesting *
- *************************************************************************
- */
-FACPoisson::~FACPoisson()
-{
-}
-
-/*
- *************************************************************************
- * Initialize data on a level. *
- * *
- * Allocate the solution, exact solution and rhs memory. *
- * Fill the rhs and exact solution. *
- *************************************************************************
- */
-void FACPoisson::initializeLevelData(
- const tbox::Pointer<hier::BasePatchHierarchy> patch_hierarchy,
- const int level_number,
- const double init_data_time,
- const bool can_be_refined,
- const bool initial_time,
- const tbox::Pointer<hier::BasePatchLevel> old_level,
- const bool allocate_data)
-{
-
- (void)init_data_time;
- (void)can_be_refined;
- (void)initial_time;
- (void)old_level;
-
- tbox::Pointer<hier::PatchHierarchy> hierarchy = patch_hierarchy;
- tbox::Pointer<geom::CartesianGridGeometry> grid_geom =
- hierarchy->getGridGeometry();
-
- tbox::Pointer<hier::PatchLevel> level =
- hierarchy->getPatchLevel(level_number);
-
- if (allocate_data) {
- level->allocatePatchData(d_comp_soln_id);
- level->allocatePatchData(d_rhs_id);
- level->allocatePatchData(d_exact_id);
- }
-
- /*
- * Initialize data in all patches in the level.
- */
- hier::PatchLevel::Iterator pi(*level);
- for (pi.initialize(*level); pi; pi++) {
-
- tbox::Pointer<hier::Patch> patch = *pi;
- if (patch.isNull()) {
- TBOX_ERROR(d_object_name
- << ": Cannot find patch. Null patch pointer.");
- }
- hier::Box pbox = patch->getBox();
- tbox::Pointer<geom::CartesianPatchGeometry> patch_geom =
- patch->getPatchGeometry();
-
- tbox::Pointer<pdat::CellData<double> > exact_data =
- patch->getPatchData(d_exact_id);
- tbox::Pointer<pdat::CellData<double> > rhs_data =
- patch->getPatchData(d_rhs_id);
-
- /*
- * Set source function and exact solution.
- */
- if (d_dim == tbox::Dimension(2)) {
- F77_FUNC(setexactandrhs2d, SETEXACTANDRHS2D) (
- pbox.lower()[0],
- pbox.upper()[0],
- pbox.lower()[1],
- pbox.upper()[1],
- exact_data->getPointer(),
- rhs_data->getPointer(),
- patch_geom->getDx(),
- patch_geom->getXLower());
- } else if (d_dim == tbox::Dimension(3)) {
- F77_FUNC(setexactandrhs3d, SETEXACTANDRHS3D) (
- pbox.lower()[0],
- pbox.upper()[0],
- pbox.lower()[1],
- pbox.upper()[1],
- pbox.lower()[2],
- pbox.upper()[2],
- exact_data->getPointer(),
- rhs_data->getPointer(),
- patch_geom->getDx(),
- patch_geom->getXLower());
- }
-
- } // End patch loop.
-}
-
-/*
- *************************************************************************
- * Reset the hierarchy-dependent internal information. *
- *************************************************************************
- */
-void FACPoisson::resetHierarchyConfiguration(
- tbox::Pointer<hier::BasePatchHierarchy> new_hierarchy,
- int coarsest_level,
- int finest_level)
-{
- (void)coarsest_level;
- (void)finest_level;
-
- d_hierarchy = new_hierarchy;
-}
-
-/*
- *************************************************************************
- * Set up the initial guess and problem parameters *
- * and solve the Poisson problem. We explicitly initialize and *
- * deallocate the solver state in this example. *
- *************************************************************************
- */
-int FACPoisson::solvePoisson()
-{
-
- if (d_hierarchy.isNull()) {
- TBOX_ERROR(d_object_name
- << "Cannot solve using an uninitialized object.\n");
- }
-
- int ln;
- /*
- * Fill in the initial guess.
- */
- for (ln = 0; ln <= d_hierarchy->getFinestLevelNumber(); ++ln) {
- tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
- hier::PatchLevel::Iterator ip(*level);
- for ( ; ip; ip++) {
- tbox::Pointer<hier::Patch> patch = *ip;
- tbox::Pointer<pdat::CellData<double> > data = patch->getPatchData(
- d_comp_soln_id);
- data->fill(0.0);
- }
- }
-
- /*
- * Set the parameters for the Poisson equation.
- * See classes solv::CellPoissonFACSolver or
- * solv::PoissonSpecifications.
- * (D is the diffusion coefficient.
- * C is the source term which is not used in this example.)
- */
- d_poisson_fac_solver.setDConstant(1.0);
- d_poisson_fac_solver.setCConstant(0.0);
-
- d_poisson_fac_solver.initializeSolverState(
- d_comp_soln_id,
- d_rhs_id,
- d_hierarchy,
- 0,
- d_hierarchy->getFinestLevelNumber());
-
- tbox::plog << "solving..." << std::endl;
- int solver_ret;
- solver_ret = d_poisson_fac_solver.solveSystem(d_comp_soln_id,
- d_rhs_id);
- /*
- * Present data on the solve.
- */
- double avg_factor, final_factor;
- d_poisson_fac_solver.getConvergenceFactors(avg_factor, final_factor);
- tbox::plog << "\t" << (solver_ret ? "" : "NOT ") << "converged " << "\n"
- << " iterations: "
- << d_poisson_fac_solver.getNumberOfIterations() << "\n"
- << " residual: "<< d_poisson_fac_solver.getResidualNorm()
- << "\n"
- << " average convergence: "<< avg_factor << "\n"
- << " final convergence: "<< final_factor << "\n"
- << std::flush;
-
- d_poisson_fac_solver.deallocateSolverState();
-
- return 0;
-}
-
-#ifdef HAVE_HDF5
-/*
- *************************************************************************
- * Set up external plotter to plot internal data from this class. *
- * Register variables appropriate for plotting. *
- *************************************************************************
- */
-int FACPoisson::setupPlotter(
- appu::VisItDataWriter& plotter) const {
- if (d_hierarchy.isNull()) {
- TBOX_ERROR(d_object_name << ": No hierarchy in\n"
- << " FACPoisson::setupPlotter\n"
- << "The hierarchy must be set before calling\n"
- << "this function.\n");
- }
- plotter.registerPlotQuantity("Computed solution",
- "SCALAR",
- d_comp_soln_id);
- plotter.registerDerivedPlotQuantity("Error",
- "SCALAR",
- (appu::VisDerivedDataStrategy *)this);
- plotter.registerPlotQuantity("Exact solution",
- "SCALAR",
- d_exact_id);
- plotter.registerPlotQuantity("Poisson source",
- "SCALAR",
- d_rhs_id);
-
- return 0;
-}
-#endif
-
-/*
- *************************************************************************
- * Write derived data to the given stream. *
- *************************************************************************
- */
-bool FACPoisson::packDerivedDataIntoDoubleBuffer(
- double* buffer,
- const hier::Patch& patch,
- const hier::Box& region,
- const std::string& variable_name,
- int depth_id) const
-{
- (void)depth_id;
-
- pdat::CellData<double>::Iterator icell(region);
-
- if (variable_name == "Error") {
- tbox::Pointer<pdat::CellData<double> > current_solution_ =
- patch.getPatchData(d_comp_soln_id);
- tbox::Pointer<pdat::CellData<double> > exact_solution_ =
- patch.getPatchData(d_exact_id);
- pdat::CellData<double>& current_solution = *current_solution_;
- pdat::CellData<double>& exact_solution = *exact_solution_;
- for ( ; icell; icell++) {
- double diff = (current_solution(*icell) - exact_solution(*icell));
- *buffer = diff;
- buffer = buffer + 1;
- }
- } else {
- // Did not register this name.
- TBOX_ERROR(
- "Unregistered variable name '" << variable_name << "' in\n"
- <<
- "FACPoissonX::packDerivedDataIntoDoubleBuffer");
-
- }
- // Return true if this patch has derived data on it.
- // False otherwise.
- return true;
-}
-
-}
diff -r 390def91ff57 -r cd1d77fa073e FACPoisson.h
--- a/FACPoisson.h Wed Dec 29 16:12:07 2010 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,211 +0,0 @@
-/*************************************************************************
- *
- * This file is part of the SAMRAI distribution. For full copyright
- * information, see COPYRIGHT and COPYING.LESSER.
- *
- * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
- * Description: Numerical routines for example FAC Poisson solver
- *
- ************************************************************************/
-#ifndef included_FACPoisson
-#define included_FACPoisson
-
-#include "SAMRAI/solv/CellPoissonFACSolver.h"
-#include "SAMRAI/pdat/CellVariable.h"
-#include "SAMRAI/tbox/Database.h"
-#include "SAMRAI/hier/Box.h"
-#include "SAMRAI/solv/LocationIndexRobinBcCoefs.h"
-#include "SAMRAI/hier/Patch.h"
-#include "SAMRAI/hier/PatchHierarchy.h"
-#include "SAMRAI/hier/PatchLevel.h"
-#include "SAMRAI/tbox/Pointer.h"
-#include "SAMRAI/pdat/SideVariable.h"
-#include "SAMRAI/mesh/StandardTagAndInitStrategy.h"
-#include "SAMRAI/hier/VariableContext.h"
-#include "SAMRAI/appu/VisDerivedDataStrategy.h"
-#include "SAMRAI/appu/VisItDataWriter.h"
-
-namespace SAMRAI {
-
-/*!
- * @brief Class to solve a sample Poisson equation on a SAMR grid.
- *
- * This class demonstrates how use the FAC Poisson solver
- * class to solve Poisson's equation on a single level
- * within a hierarchy.
- *
- * We set up and solve the following problem:
- *
- * 2d: div(grad(u)) = -2 (pi^2) sin(pi x) sin(pi y)
- *
- * 3d: div(grad(u)) = -3 (pi^2) sin(pi x) sin(pi y) sin(pi z)
- *
- * which has the exact solution
- *
- * 2d: u = sin(pi x) sin(pi y)
- *
- * 3d: u = sin(pi x) sin(pi y) sin(pi z)
- *
- * This class inherits and implements virtual functions from
- * - mesh::StandardTagAndInitStrategy to initialize data
- * on the SAMR grid.
- * - appu::VisDerivedDataStrategy to write out certain data
- * in a vis file, such as the error of the solution.
- *
- * Inputs: The only input parameter for this class is
- * "fac_poisson", the input database for the solv::CellPoissonFACSolver
- * object. See the documentation for solv::CellPoissonFACSolver
- * for its input parameters.
- */
-class FACPoisson:
- public mesh::StandardTagAndInitStrategy,
- public appu::VisDerivedDataStrategy
-{
-
-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 (may be NULL)
- */
- FACPoisson(
- const std::string& object_name,
- const tbox::Dimension& dim,
- tbox::Pointer<tbox::Database> database =
- tbox::Pointer<tbox::Database>(NULL));
-
- virtual ~FACPoisson();
-
- //@{ @name mesh::StandardTagAndInitStrategy virtuals
-
- /*!
- * @brief Allocate and initialize data for a new level
- * in the patch hierarchy.
- *
- * This is where you implement the code for initialize data on
- * the grid. All the information needed to initialize the grid
- * are in the arguments.
- *
- * @see mesh::StandardTagAndInitStrategy::initializeLevelData()
- */
- virtual void
- initializeLevelData(
- const tbox::Pointer<hier::BasePatchHierarchy> hierarchy,
- const int level_number,
- const double init_data_time,
- const bool can_be_refined,
- const bool initial_time,
- const tbox::Pointer<hier::BasePatchLevel> old_level,
- const bool allocate_data);
-
- /*!
- * @brief Reset any internal hierarchy-dependent information.
- */
- virtual void
- resetHierarchyConfiguration(
- tbox::Pointer<hier::BasePatchHierarchy> new_hierarchy,
- int coarsest_level,
- int finest_level);
-
- //@}
-
- //@{ @name appu::VisDerivedDataStrategy virtuals
-
- virtual bool
- packDerivedDataIntoDoubleBuffer(
- double* buffer,
- const hier::Patch& patch,
- const hier::Box& region,
- const std::string& variable_name,
- int depth_id) const;
-
- //@}
-
- /*!
- * @brief Solve using HYPRE Poisson solver
- *
- * Set up the linear algebra problem and use a
- * solv::CellPoissonFACSolver object to solve it.
- * -# Set initial guess
- * -# Set boundary conditions
- * -# Specify Poisson equation parameters
- * -# Call solver
- */
- int
- solvePoisson();
-
-#ifdef HAVE_HDF5
- /*!
- * @brief Set up external plotter to plot internal
- * data from this class.
- *
- * After calling this function, the external
- * data writer may be used to write the
- * viz file for this object.
- *
- * The internal hierarchy is used and must be
- * established before calling this function.
- * (This is commonly done by building a hierarchy
- * with the mesh::StandardTagAndInitStrategy virtual
- * functions implemented by this class.)
- *
- * @param viz_writer VisIt writer
- */
- int
- setupPlotter(
- appu::VisItDataWriter& plotter) const;
-#endif
-
-private:
- std::string d_object_name;
-
- const tbox::Dimension d_dim;
-
- tbox::Pointer<hier::PatchHierarchy> d_hierarchy;
-
- //@{
- /*!
- * @name Major algorithm objects.
- */
-
- /*!
- * @brief FAC poisson solver.
- */
- solv::CellPoissonFACSolver d_poisson_fac_solver;
-
- /*!
- * @brief Boundary condition coefficient implementation.
- */
- solv::LocationIndexRobinBcCoefs d_bc_coefs;
-
- //@}
-
- //@{
-
- /*!
- * @name Private state variables for solution.
- */
-
- /*!
- * @brief Context owned by this object.
- */
- tbox::Pointer<hier::VariableContext> d_context;
-
- /*!
- * @brief Descriptor indices of internal data.
- *
- * These are initialized in the constructor and never change.
- */
- int d_comp_soln_id, d_exact_id, d_rhs_id;
-
- //@}
-
-};
-
-}
-
-#endif // included_FACPoisson
diff -r 390def91ff57 -r cd1d77fa073e FACStokes.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/FACStokes.C Thu Dec 30 20:01:24 2010 -0800
@@ -0,0 +1,378 @@
+/*************************************************************************
+ *
+ * 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: Numerical routines for example FAC Stokes solver
+ *
+ ************************************************************************/
+#include "FACStokes.h"
+
+#include "SAMRAI/hier/IntVector.h"
+#include "SAMRAI/geom/CartesianGridGeometry.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/solv/SimpleCellRobinBcCoefs.h"
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/math/HierarchyCellDataOpsReal.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/solv/PoissonSpecifications.h"
+#include "SAMRAI/tbox/Utilities.h"
+#include "SAMRAI/hier/Variable.h"
+#include "SAMRAI/hier/VariableDatabase.h"
+
+extern "C" {
+void F77_FUNC(setexactandrhs2d, SETEXACTANDRHS2D) (const int & ifirst0,
+ const int & ilast0,
+ const int & ifirst1,
+ const int & ilast1,
+ double* exact,
+ double* rhs,
+ const double* dx,
+ const double* xlower);
+
+void F77_FUNC(setexactandrhs3d, SETEXACTANDRHS3D) (const int & ifirst0,
+ const int & ilast0,
+ const int & ifirst1,
+ const int & ilast1,
+ const int & ifirst2,
+ const int & ilast2,
+ double* exact,
+ double* rhs,
+ const double* dx,
+ const double* xlower);
+}
+
+namespace SAMRAI {
+
+/*
+ *************************************************************************
+ * Constructor creates a unique context for the object and register *
+ * all its internal variables with the variable database. *
+ *************************************************************************
+ */
+FACStokes::FACStokes(
+ const std::string& object_name,
+ const tbox::Dimension& dim,
+ tbox::Pointer<tbox::Database> database):
+ d_object_name(object_name),
+ d_dim(dim),
+ d_hierarchy(NULL),
+ d_stokes_fac_solver((d_dim),
+ object_name + "::poisson_hypre",
+ (!database.isNull() &&
+ database->isDatabase("fac_solver")) ?
+ database->getDatabase("fac_solver"):
+ tbox::Pointer<tbox::Database>(NULL)),
+ d_bc_coefs(d_dim,
+ object_name + "::bc_coefs",
+ (!database.isNull() &&
+ database->isDatabase("bc_coefs")) ?
+ database->getDatabase("bc_coefs"):
+ tbox::Pointer<tbox::Database>(NULL)),
+ d_context()
+{
+
+ hier::VariableDatabase* vdb =
+ hier::VariableDatabase::getDatabase();
+
+ /*
+ * Get a unique context for variables owned by this object.
+ */
+ d_context = vdb->getContext(d_object_name + ":Context");
+
+ /*
+ * Register variables with hier::VariableDatabase
+ * and get the descriptor indices for those variables.
+ */
+
+ tbox::Pointer<pdat::CellVariable<double> > comp_soln(
+ new pdat::CellVariable<double>(dim, object_name + ":computed solution", 1));
+ d_comp_soln_id =
+ vdb->registerVariableAndContext(
+ comp_soln,
+ d_context,
+ hier::IntVector(dim, 1) /* ghost cell width is 1 for stencil widths */);
+
+ tbox::Pointer<pdat::CellVariable<double> > exact_solution(
+ new pdat::CellVariable<double>(dim, object_name + ":exact solution"));
+ d_exact_id =
+ vdb->registerVariableAndContext(
+ exact_solution,
+ d_context,
+ hier::IntVector(dim, 1) /* ghost cell width is 1 in case needed */);
+
+ tbox::Pointer<pdat::CellVariable<double> > rhs_variable(
+ new pdat::CellVariable<double>(
+ dim,
+ object_name
+ + ":linear system right hand side"));
+ d_rhs_id =
+ vdb->registerVariableAndContext(
+ rhs_variable,
+ d_context,
+ hier::IntVector(dim, 0) /* ghost cell width is 0 */);
+
+ /*
+ * Specify an implementation of solv::RobinBcCoefStrategy for the solver to use.
+ * We use the implementation solv::LocationIndexRobinBcCoefs, but other
+ * implementations are possible, including user-implemented.
+ */
+ d_stokes_fac_solver.setBcObject(&d_bc_coefs);
+}
+
+/*
+ *************************************************************************
+ * Destructor does nothing interesting *
+ *************************************************************************
+ */
+FACStokes::~FACStokes()
+{
+}
+
+/*
+ *************************************************************************
+ * Initialize data on a level. *
+ * *
+ * Allocate the solution, exact solution and rhs memory. *
+ * Fill the rhs and exact solution. *
+ *************************************************************************
+ */
+void FACStokes::initializeLevelData(
+ const tbox::Pointer<hier::BasePatchHierarchy> patch_hierarchy,
+ const int level_number,
+ const double init_data_time,
+ const bool can_be_refined,
+ const bool initial_time,
+ const tbox::Pointer<hier::BasePatchLevel> old_level,
+ const bool allocate_data)
+{
+
+ (void)init_data_time;
+ (void)can_be_refined;
+ (void)initial_time;
+ (void)old_level;
+
+ tbox::Pointer<hier::PatchHierarchy> hierarchy = patch_hierarchy;
+ tbox::Pointer<geom::CartesianGridGeometry> grid_geom =
+ hierarchy->getGridGeometry();
+
+ tbox::Pointer<hier::PatchLevel> level =
+ hierarchy->getPatchLevel(level_number);
+
+ if (allocate_data) {
+ level->allocatePatchData(d_comp_soln_id);
+ level->allocatePatchData(d_rhs_id);
+ level->allocatePatchData(d_exact_id);
+ }
+
+ /*
+ * Initialize data in all patches in the level.
+ */
+ hier::PatchLevel::Iterator pi(*level);
+ for (pi.initialize(*level); pi; pi++) {
+
+ tbox::Pointer<hier::Patch> patch = *pi;
+ if (patch.isNull()) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot find patch. Null patch pointer.");
+ }
+ hier::Box pbox = patch->getBox();
+ tbox::Pointer<geom::CartesianPatchGeometry> patch_geom =
+ patch->getPatchGeometry();
+
+ tbox::Pointer<pdat::CellData<double> > exact_data =
+ patch->getPatchData(d_exact_id);
+ tbox::Pointer<pdat::CellData<double> > rhs_data =
+ patch->getPatchData(d_rhs_id);
+
+ /*
+ * Set source function and exact solution.
+ */
+ if (d_dim == tbox::Dimension(2)) {
+ F77_FUNC(setexactandrhs2d, SETEXACTANDRHS2D) (
+ pbox.lower()[0],
+ pbox.upper()[0],
+ pbox.lower()[1],
+ pbox.upper()[1],
+ exact_data->getPointer(),
+ rhs_data->getPointer(),
+ patch_geom->getDx(),
+ patch_geom->getXLower());
+ } else if (d_dim == tbox::Dimension(3)) {
+ F77_FUNC(setexactandrhs3d, SETEXACTANDRHS3D) (
+ pbox.lower()[0],
+ pbox.upper()[0],
+ pbox.lower()[1],
+ pbox.upper()[1],
+ pbox.lower()[2],
+ pbox.upper()[2],
+ exact_data->getPointer(),
+ rhs_data->getPointer(),
+ patch_geom->getDx(),
+ patch_geom->getXLower());
+ }
+
+ } // End patch loop.
+}
+
+/*
+ *************************************************************************
+ * Reset the hierarchy-dependent internal information. *
+ *************************************************************************
+ */
+void FACStokes::resetHierarchyConfiguration(
+ tbox::Pointer<hier::BasePatchHierarchy> new_hierarchy,
+ int coarsest_level,
+ int finest_level)
+{
+ (void)coarsest_level;
+ (void)finest_level;
+
+ d_hierarchy = new_hierarchy;
+}
+
+/*
+ *************************************************************************
+ * Set up the initial guess and problem parameters *
+ * and solve the Stokes problem. We explicitly initialize and *
+ * deallocate the solver state in this example. *
+ *************************************************************************
+ */
+int FACStokes::solveStokes()
+{
+
+ if (d_hierarchy.isNull()) {
+ TBOX_ERROR(d_object_name
+ << "Cannot solve using an uninitialized object.\n");
+ }
+
+ int ln;
+ /*
+ * Fill in the initial guess.
+ */
+ for (ln = 0; ln <= d_hierarchy->getFinestLevelNumber(); ++ln) {
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
+ hier::PatchLevel::Iterator ip(*level);
+ for ( ; ip; ip++) {
+ tbox::Pointer<hier::Patch> patch = *ip;
+ tbox::Pointer<pdat::CellData<double> > data = patch->getPatchData(
+ d_comp_soln_id);
+ data->fill(0.0);
+ }
+ }
+
+ /*
+ * Set the parameters for the Stokes equation.
+ * See classes solv::CellStokesFACSolver or
+ * solv::StokesSpecifications.
+ * (D is the diffusion coefficient.
+ * C is the source term which is not used in this example.)
+ */
+ d_stokes_fac_solver.setDConstant(1.0);
+ d_stokes_fac_solver.setCConstant(0.0);
+
+ d_stokes_fac_solver.initializeSolverState(
+ d_comp_soln_id,
+ d_rhs_id,
+ d_hierarchy,
+ 0,
+ d_hierarchy->getFinestLevelNumber());
+
+ tbox::plog << "solving..." << std::endl;
+ int solver_ret;
+ solver_ret = d_stokes_fac_solver.solveSystem(d_comp_soln_id,
+ d_rhs_id);
+ /*
+ * Present data on the solve.
+ */
+ double avg_factor, final_factor;
+ d_stokes_fac_solver.getConvergenceFactors(avg_factor, final_factor);
+ tbox::plog << "\t" << (solver_ret ? "" : "NOT ") << "converged " << "\n"
+ << " iterations: "
+ << d_stokes_fac_solver.getNumberOfIterations() << "\n"
+ << " residual: "<< d_stokes_fac_solver.getResidualNorm()
+ << "\n"
+ << " average convergence: "<< avg_factor << "\n"
+ << " final convergence: "<< final_factor << "\n"
+ << std::flush;
+
+ d_stokes_fac_solver.deallocateSolverState();
+
+ return 0;
+}
+
+#ifdef HAVE_HDF5
+/*
+ *************************************************************************
+ * Set up external plotter to plot internal data from this class. *
+ * Register variables appropriate for plotting. *
+ *************************************************************************
+ */
+int FACStokes::setupPlotter(
+ appu::VisItDataWriter& plotter) const {
+ if (d_hierarchy.isNull()) {
+ TBOX_ERROR(d_object_name << ": No hierarchy in\n"
+ << " FACStokes::setupPlotter\n"
+ << "The hierarchy must be set before calling\n"
+ << "this function.\n");
+ }
+ plotter.registerPlotQuantity("Computed solution",
+ "SCALAR",
+ d_comp_soln_id);
+ plotter.registerDerivedPlotQuantity("Error",
+ "SCALAR",
+ (appu::VisDerivedDataStrategy *)this);
+ plotter.registerPlotQuantity("Exact solution",
+ "SCALAR",
+ d_exact_id);
+ plotter.registerPlotQuantity("Stokes source",
+ "SCALAR",
+ d_rhs_id);
+
+ return 0;
+}
+#endif
+
+/*
+ *************************************************************************
+ * Write derived data to the given stream. *
+ *************************************************************************
+ */
+bool FACStokes::packDerivedDataIntoDoubleBuffer(
+ double* buffer,
+ const hier::Patch& patch,
+ const hier::Box& region,
+ const std::string& variable_name,
+ int depth_id) const
+{
+ (void)depth_id;
+
+ pdat::CellData<double>::Iterator icell(region);
+
+ if (variable_name == "Error") {
+ tbox::Pointer<pdat::CellData<double> > current_solution_ =
+ patch.getPatchData(d_comp_soln_id);
+ tbox::Pointer<pdat::CellData<double> > exact_solution_ =
+ patch.getPatchData(d_exact_id);
+ pdat::CellData<double>& current_solution = *current_solution_;
+ pdat::CellData<double>& exact_solution = *exact_solution_;
+ for ( ; icell; icell++) {
+ double diff = (current_solution(*icell) - exact_solution(*icell));
+ *buffer = diff;
+ buffer = buffer + 1;
+ }
+ } else {
+ // Did not register this name.
+ TBOX_ERROR(
+ "Unregistered variable name '" << variable_name << "' in\n"
+ <<
+ "FACStokesX::packDerivedDataIntoDoubleBuffer");
+
+ }
+ // Return true if this patch has derived data on it.
+ // False otherwise.
+ return true;
+}
+
+}
diff -r 390def91ff57 -r cd1d77fa073e FACStokes.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/FACStokes.h Thu Dec 30 20:01:24 2010 -0800
@@ -0,0 +1,211 @@
+/*************************************************************************
+ *
+ * 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: Numerical routines for example FAC Stokes solver
+ *
+ ************************************************************************/
+#ifndef included_FACStokes
+#define included_FACStokes
+
+#include "SAMRAI/solv/CellPoissonFACSolver.h"
+#include "SAMRAI/pdat/CellVariable.h"
+#include "SAMRAI/tbox/Database.h"
+#include "SAMRAI/hier/Box.h"
+#include "SAMRAI/solv/LocationIndexRobinBcCoefs.h"
+#include "SAMRAI/hier/Patch.h"
+#include "SAMRAI/hier/PatchHierarchy.h"
+#include "SAMRAI/hier/PatchLevel.h"
+#include "SAMRAI/tbox/Pointer.h"
+#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/mesh/StandardTagAndInitStrategy.h"
+#include "SAMRAI/hier/VariableContext.h"
+#include "SAMRAI/appu/VisDerivedDataStrategy.h"
+#include "SAMRAI/appu/VisItDataWriter.h"
+
+namespace SAMRAI {
+
+/*!
+ * @brief Class to solve a sample Stokes equation on a SAMR grid.
+ *
+ * This class demonstrates how use the FAC Stokes solver
+ * class to solve Stokes's equation on a single level
+ * within a hierarchy.
+ *
+ * We set up and solve the following problem:
+ *
+ * 2d: div(grad(u)) = -2 (pi^2) sin(pi x) sin(pi y)
+ *
+ * 3d: div(grad(u)) = -3 (pi^2) sin(pi x) sin(pi y) sin(pi z)
+ *
+ * which has the exact solution
+ *
+ * 2d: u = sin(pi x) sin(pi y)
+ *
+ * 3d: u = sin(pi x) sin(pi y) sin(pi z)
+ *
+ * This class inherits and implements virtual functions from
+ * - mesh::StandardTagAndInitStrategy to initialize data
+ * on the SAMR grid.
+ * - appu::VisDerivedDataStrategy to write out certain data
+ * in a vis file, such as the error of the solution.
+ *
+ * Inputs: The only input parameter for this class is
+ * "fac_stokes", the input database for the solv::CellStokesFACSolver
+ * object. See the documentation for solv::CellStokesFACSolver
+ * for its input parameters.
+ */
+class FACStokes:
+ public mesh::StandardTagAndInitStrategy,
+ public appu::VisDerivedDataStrategy
+{
+
+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 (may be NULL)
+ */
+ FACStokes(
+ const std::string& object_name,
+ const tbox::Dimension& dim,
+ tbox::Pointer<tbox::Database> database =
+ tbox::Pointer<tbox::Database>(NULL));
+
+ virtual ~FACStokes();
+
+ //@{ @name mesh::StandardTagAndInitStrategy virtuals
+
+ /*!
+ * @brief Allocate and initialize data for a new level
+ * in the patch hierarchy.
+ *
+ * This is where you implement the code for initialize data on
+ * the grid. All the information needed to initialize the grid
+ * are in the arguments.
+ *
+ * @see mesh::StandardTagAndInitStrategy::initializeLevelData()
+ */
+ virtual void
+ initializeLevelData(
+ const tbox::Pointer<hier::BasePatchHierarchy> hierarchy,
+ const int level_number,
+ const double init_data_time,
+ const bool can_be_refined,
+ const bool initial_time,
+ const tbox::Pointer<hier::BasePatchLevel> old_level,
+ const bool allocate_data);
+
+ /*!
+ * @brief Reset any internal hierarchy-dependent information.
+ */
+ virtual void
+ resetHierarchyConfiguration(
+ tbox::Pointer<hier::BasePatchHierarchy> new_hierarchy,
+ int coarsest_level,
+ int finest_level);
+
+ //@}
+
+ //@{ @name appu::VisDerivedDataStrategy virtuals
+
+ virtual bool
+ packDerivedDataIntoDoubleBuffer(
+ double* buffer,
+ const hier::Patch& patch,
+ const hier::Box& region,
+ const std::string& variable_name,
+ int depth_id) const;
+
+ //@}
+
+ /*!
+ * @brief Solve using HYPRE Stokes solver
+ *
+ * Set up the linear algebra problem and use a
+ * solv::CellStokesFACSolver object to solve it.
+ * -# Set initial guess
+ * -# Set boundary conditions
+ * -# Specify Stokes equation parameters
+ * -# Call solver
+ */
+ int
+ solveStokes();
+
+#ifdef HAVE_HDF5
+ /*!
+ * @brief Set up external plotter to plot internal
+ * data from this class.
+ *
+ * After calling this function, the external
+ * data writer may be used to write the
+ * viz file for this object.
+ *
+ * The internal hierarchy is used and must be
+ * established before calling this function.
+ * (This is commonly done by building a hierarchy
+ * with the mesh::StandardTagAndInitStrategy virtual
+ * functions implemented by this class.)
+ *
+ * @param viz_writer VisIt writer
+ */
+ int
+ setupPlotter(
+ appu::VisItDataWriter& plotter) const;
+#endif
+
+private:
+ std::string d_object_name;
+
+ const tbox::Dimension d_dim;
+
+ tbox::Pointer<hier::PatchHierarchy> d_hierarchy;
+
+ //@{
+ /*!
+ * @name Major algorithm objects.
+ */
+
+ /*!
+ * @brief FAC stokes solver.
+ */
+ solv::CellPoissonFACSolver d_stokes_fac_solver;
+
+ /*!
+ * @brief Boundary condition coefficient implementation.
+ */
+ solv::LocationIndexRobinBcCoefs d_bc_coefs;
+
+ //@}
+
+ //@{
+
+ /*!
+ * @name Private state variables for solution.
+ */
+
+ /*!
+ * @brief Context owned by this object.
+ */
+ tbox::Pointer<hier::VariableContext> d_context;
+
+ /*!
+ * @brief Descriptor indices of internal data.
+ *
+ * These are initialized in the constructor and never change.
+ */
+ int d_comp_soln_id, d_exact_id, d_rhs_id;
+
+ //@}
+
+};
+
+}
+
+#endif // included_FACStokes
diff -r 390def91ff57 -r cd1d77fa073e Makefile
--- a/Makefile Wed Dec 29 16:12:07 2010 -0800
+++ b/Makefile Thu Dec 30 20:01:24 2010 -0800
@@ -4,7 +4,7 @@
## information, see COPYRIGHT and COPYING.LESSER.
##
## Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
-## Description: makefile for SAMRAI FAC Poisson solver example
+## Description: makefile for SAMRAI FAC Stokes
##
#########################################################################
@@ -24,7 +24,7 @@ NUM_TESTS = 2
TEST_NPROCS = 0,2
-CXX_OBJS = main.o FACPoisson.o
+CXX_OBJS = main.o FACStokes.o
F_OBJS = facpoisson2d.o facpoisson3d.o
main: $(CXX_OBJS) $(F_OBJS) $(LIBSAMRAIDEPEND)
diff -r 390def91ff57 -r cd1d77fa073e Makefile.depend
--- a/Makefile.depend Wed Dec 29 16:12:07 2010 -0800
+++ b/Makefile.depend Thu Dec 30 20:01:24 2010 -0800
@@ -12,7 +12,7 @@
## This file is automatically generated by depend.pl.
-FILE_0=FACPoisson.o
+FILE_0=FACStokes.o
DEPENDS_0:=\
$(OBJECT)/include/SAMRAI/SAMRAI_config.h \
$(INCLUDE_SAM)/SAMRAI/appu/VisDerivedDataStrategy.h \
@@ -325,8 +325,8 @@ DEPENDS_0:=\
$(INCLUDE_SAM)/SAMRAI/xfer/RefineSchedule.h \
$(INCLUDE_SAM)/SAMRAI/xfer/RefineTransactionFactory.h \
$(INCLUDE_SAM)/SAMRAI/xfer/TimeInterpolateOperator.h \
- $(INCLUDE_SAM)/SAMRAI/xfer/VariableFillPattern.h FACPoisson.C \
- FACPoisson.h
+ $(INCLUDE_SAM)/SAMRAI/xfer/VariableFillPattern.h FACStokes.C \
+ FACStokes.h
ifeq (${DEPENDS_ON_TEMPLATE_IMPLEMENTATION},yes)
DEPENDS_0 +=\
@@ -737,7 +737,7 @@ DEPENDS_1:=\
$(INCLUDE_SAM)/SAMRAI/xfer/RefineSchedule.h \
$(INCLUDE_SAM)/SAMRAI/xfer/RefineTransactionFactory.h \
$(INCLUDE_SAM)/SAMRAI/xfer/TimeInterpolateOperator.h \
- $(INCLUDE_SAM)/SAMRAI/xfer/VariableFillPattern.h FACPoisson.h \
+ $(INCLUDE_SAM)/SAMRAI/xfer/VariableFillPattern.h FACStokes.h \
main.C
ifeq (${DEPENDS_ON_TEMPLATE_IMPLEMENTATION},yes)
diff -r 390def91ff57 -r cd1d77fa073e example_inputs/const_refine.2d.input
--- a/example_inputs/const_refine.2d.input Wed Dec 29 16:12:07 2010 -0800
+++ b/example_inputs/const_refine.2d.input Thu Dec 30 20:01:24 2010 -0800
@@ -4,7 +4,7 @@
* information, see COPYRIGHT and COPYING.LESSER.
*
* Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
- * Description: Input file for example FAC Poisson solver
+ * Description: Input file for example FAC Stokes solver
*
************************************************************************/
// This is the input file for the 2D FAC example
@@ -23,13 +23,13 @@ Main {
vis_writer = "Vizamrai", "VisIt"
}
-FACPoisson {
- // The FACPoisson class is the "user class" in this example.
+FACStokes {
+ // The FACStokes class is the "user class" in this example.
// It owns the solver and contains the code to set up the solver.
- // The inputs for FACPoisson is simply the inputs for the individual
- // parts owned by the FACPoisson class.
+ // The inputs for FACStokes is simply the inputs for the individual
+ // parts owned by the FACStokes class.
fac_solver {
- // This is the input for the cell-centered Poisson FAC solver
+ // This is the input for the cell-centered Stokes FAC solver
// class in the SAMRAI library.
enable_logging = TRUE // Bool flag to switch logging on/off
max_cycles = 20 // Max number of FAC cycles to use
diff -r 390def91ff57 -r cd1d77fa073e main.C
--- a/main.C Wed Dec 29 16:12:07 2010 -0800
+++ b/main.C Thu Dec 30 20:01:24 2010 -0800
@@ -4,7 +4,7 @@
* information, see COPYRIGHT and COPYING.LESSER.
*
* Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
- * Description: Main program for FAC Poisson example
+ * Description: Main program for FAC Stokes example
*
************************************************************************/
#include "SAMRAI/SAMRAI_config.h"
@@ -29,7 +29,7 @@ using namespace std;
#include "SAMRAI/tbox/Utilities.h"
#include "SAMRAI/appu/VisItDataWriter.h"
-#include "FACPoisson.h"
+#include "FACStokes.h"
using namespace SAMRAI;
@@ -37,20 +37,20 @@ using namespace SAMRAI;
************************************************************************
* *
* This is the driver program to demonstrate *
- * how to use the FAC Poisson solver. *
+ * how to use the FAC Stokes solver. *
* *
* We set up the simple problem *
* u + div(grad(u)) = sin(x)*sin(y) *
* in the domain [0:1]x[0:1], with u=0 on the *
* boundary. *
* *
- * FACPoisson is the primary object used to *
+ * FACStokes is the primary object used to *
* set up and solve the system. It maintains *
* the data for the computed solution u, the *
* exact solution, and the right hand side. *
* *
* The hierarchy created to solve this problem *
- * has only one level. (The FAC Poisson solver *
+ * has only one level. (The FAC Stokes solver *
* is a single-level solver.) *
* *
*************************************************************************
@@ -152,16 +152,16 @@ int main(
input_db->getDatabase("PatchHierarchy")));
/*
- * The FACPoisson object is the main user object specific to the
+ * The FACStokes object is the main user object specific to the
* problem being solved. It provides the implementations for setting
* up the grid and plotting data. It also wraps up the solve
* process that includes making the initial guess, specifying the
* boundary conditions and call the solver.
*/
- FACPoisson fac_poisson(base_name + "::FACPoisson",
+ FACStokes fac_stokes(base_name + "::FACStokes",
dim,
- input_db->isDatabase("FACPoisson") ?
- input_db->getDatabase("FACPoisson") :
+ input_db->isDatabase("FACStokes") ?
+ input_db->getDatabase("FACStokes") :
tbox::Pointer<tbox::Database>(NULL));
/*
@@ -172,7 +172,7 @@ int main(
new mesh::StandardTagAndInitialize(
dim,
"CellTaggingMethod",
- tbox::Pointer<mesh::StandardTagAndInitStrategy>(&fac_poisson, false),
+ tbox::Pointer<mesh::StandardTagAndInitStrategy>(&fac_stokes, false),
input_db->getDatabase("StandardTagAndInitialize")
));
tbox::Pointer<mesh::BergerRigoutsos> box_generator(
@@ -218,7 +218,7 @@ int main(
/*
* Set up the plotter for the hierarchy just created.
- * The FACPoisson object handles the data and has the
+ * The FACStokes object handles the data and has the
* function setupExternalPlotter to register its data
* with the plotter.
*/
@@ -239,7 +239,7 @@ int main(
visit_writer = new appu::VisItDataWriter(dim,
"Visit Writer",
vis_filename + ".visit");
- fac_poisson.setupPlotter(*visit_writer);
+ fac_stokes.setupPlotter(*visit_writer);
}
#endif
@@ -256,7 +256,7 @@ int main(
/*
* Solve.
*/
- fac_poisson.solvePoisson();
+ fac_stokes.solveStokes();
#ifdef HAVE_HDF5
/*
More information about the CIG-COMMITS
mailing list