[cig-commits] commit: Make prolongation work properly with p. May not work with v.
Mercurial
hg at geodynamics.org
Fri Feb 25 14:14:43 PST 2011
changeset: 45:f4b8acb334a5
user: Walter Landry <wlandry at caltech.edu>
date: Tue Jan 11 15:56:32 2011 -0800
files: Makefile P_Refine.C P_Refine.h StokesFACOps.I StokesFACOps.h StokesFACOps/StokesFACOps.C StokesFACOps/initializeOperatorState.C StokesFACOps/smoothErrorByRedBlack.C StokesFACSolver/StokesFACSolver.C example_inputs/const_refine.2d.input main.C
description:
Make prolongation work properly with p. May not work with v.
diff -r 4c1e4cd9d6de -r f4b8acb334a5 Makefile
--- a/Makefile Mon Jan 10 13:42:56 2011 -0800
+++ b/Makefile Tue Jan 11 15:56:32 2011 -0800
@@ -30,6 +30,7 @@ CXX_OBJS = main.o FACStokes/FACStok
FACStokes/resetHierarchyConfiguration.o \
FACStokes/setupPlotter.o \
FACStokes/solveStokes.o \
+ P_Refine.o \
StokesFACOps/StokesFACOps.o \
StokesFACOps/checkInputPatchDataIndices.o \
StokesFACOps/computeCompositeResidualOnLevel.o \
diff -r 4c1e4cd9d6de -r f4b8acb334a5 P_Refine.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/P_Refine.C Tue Jan 11 15:56:32 2011 -0800
@@ -0,0 +1,138 @@
+/*************************************************************************
+ *
+ * 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: Linear refine operator for cell-centered double data on
+ * a Cartesian mesh.
+ *
+ ************************************************************************/
+
+#ifndef included_geom_P_Refine_C
+#define included_geom_P_Refine_C
+
+#include "P_Refine.h"
+
+#include <float.h>
+#include <math.h>
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/hier/Index.h"
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/pdat/CellVariable.h"
+#include "SAMRAI/tbox/Utilities.h"
+
+/*
+ *************************************************************************
+ * *
+ * External declarations for FORTRAN routines. *
+ * *
+ *************************************************************************
+ */
+
+void SAMRAI::geom::P_Refine::refine(
+ hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::BoxOverlap& fine_overlap,
+ const hier::IntVector& ratio) const
+{
+ const pdat::CellOverlap* t_overlap =
+ dynamic_cast<const pdat::CellOverlap *>(&fine_overlap);
+
+ TBOX_ASSERT(t_overlap != NULL);
+
+ const hier::BoxList& boxes = t_overlap->getDestinationBoxList();
+ for (hier::BoxList::Iterator b(boxes); b; b++) {
+ refine(fine,
+ coarse,
+ dst_component,
+ src_component,
+ b(),
+ ratio);
+ }
+}
+
+void SAMRAI::geom::P_Refine::refine(
+ hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio) const
+{
+ const tbox::Dimension& dim(getDim());
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dim, fine, coarse, fine_box, ratio);
+
+ tbox::Pointer<pdat::CellData<double> >
+ p = coarse.getPatchData(src_component);
+ tbox::Pointer<pdat::CellData<double> >
+ p_fine = fine.getPatchData(dst_component);
+#ifdef DEBUG_CHECK_ASSERTIONS
+ TBOX_ASSERT(!p.isNull());
+ TBOX_ASSERT(!p_fine.isNull());
+ TBOX_ASSERT(p->getDepth() == p_fine->getDepth());
+#endif
+
+ hier::Box coarse_box=coarse.getBox();
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = coarse.getPatchGeometry();
+
+ for(int j=fine_box.lower(1); j<=fine_box.upper(1); ++j)
+ for(int i=fine_box.lower(0); i<=fine_box.upper(0); ++i)
+ {
+ pdat::CellIndex fine(tbox::Dimension(2));
+ fine[0]=i;
+ fine[1]=j;
+
+ pdat::CellIndex ip(fine), jp(fine);
+ ip[0]=1;
+ ip[1]=0;
+ jp[0]=0;
+ jp[1]=1;
+ pdat::CellIndex center(fine/2);
+ double dp_dx,dp_dy;
+
+ /* Pressure is cell-centered, so prolongation is a
+ linear interpolation from nearby cells. */
+
+ /* This assumes that we never have coarse levels
+ that are only one grid wide. */
+
+ if(center[0]==coarse_box.lower(0)
+ && geom->getTouchesRegularBoundary(0,0))
+ {
+ dp_dx=((*p)(center+ip)-(*p)(center))/4;
+ }
+ else if(center[0]==coarse_box.upper(0)
+ && geom->getTouchesRegularBoundary(0,1))
+ {
+ dp_dx=((*p)(center)-(*p)(center-ip))/4;
+ }
+ else
+ {
+ dp_dx=((*p)(center+ip)-(*p)(center-ip))/8;
+ }
+
+ if(center[1]==coarse_box.lower(1)
+ && geom->getTouchesRegularBoundary(1,0))
+ {
+ dp_dy=((*p)(center+jp)-(*p)(center))/4;
+ }
+ else if(center[1]==coarse_box.upper(1)
+ && geom->getTouchesRegularBoundary(1,1))
+ {
+ dp_dy=((*p)(center)-(*p)(center-jp))/4;
+ }
+ else
+ {
+ dp_dy=((*p)(center+jp)-(*p)(center-jp))/8;
+ }
+
+ (*p_fine)(fine)=(*p)(center)
+ + ((i%2==0) ? (-dp_dx) : dp_dx)
+ + ((j%2==0) ? (-dp_dy) : dp_dy);
+ }
+}
+#endif
diff -r 4c1e4cd9d6de -r f4b8acb334a5 P_Refine.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/P_Refine.h Tue Jan 11 15:56:32 2011 -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: Linear refine operator for cell-centered double data on
+ * a Cartesian mesh.
+ *
+ ************************************************************************/
+
+#ifndef included_geom_P_Refine
+#define included_geom_P_Refine
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#include "SAMRAI/xfer/RefineOperator.h"
+#include "SAMRAI/hier/Box.h"
+#include "SAMRAI/hier/IntVector.h"
+#include "SAMRAI/hier/Patch.h"
+#include "SAMRAI/tbox/Pointer.h"
+#include "SAMRAI/pdat/CellVariable.h"
+
+#include <string>
+
+namespace SAMRAI {
+namespace geom {
+
+/**
+ * Class P_Refine implements linear
+ * interpolation for cell-centered double patch data defined over a Cartesian
+ * mesh. It is derived from the xfer::RefineOperator base class.
+ * CartesianCellDoubleLinearRefine does not handle the lack of real
+ * boundaries for the pressure very well, so use this instead for
+ * pressure.
+ *
+ * The findRefineOperator() operator function returns true if the input
+ * variable is cell-centered double, and the std::string is "P_REFINE".
+ *
+ * @see xfer::RefineOperator
+ */
+
+class P_Refine:
+ public xfer::RefineOperator
+{
+public:
+ /**
+ * Uninteresting default constructor.
+ */
+ explicit P_Refine(const tbox::Dimension& dim):
+ xfer::RefineOperator(dim, "P_REFINE")
+ {
+ d_name_id = "P_REFINE";
+ }
+
+
+ /**
+ * Uninteresting virtual destructor.
+ */
+ virtual ~P_Refine(){}
+
+ /**
+ * Return true if the variable and name std::string match cell-centered
+ * double linear interpolation; otherwise, return false.
+ */
+ bool findRefineOperator(const tbox::Pointer<hier::Variable>& var,
+ const std::string& op_name) const
+ {
+ TBOX_DIM_ASSERT_CHECK_ARGS2(*this, *var);
+
+ const tbox::Pointer<pdat::CellVariable<double> > cast_var(var);
+ if (!cast_var.isNull() && (op_name == d_name_id)) {
+ return true;
+ } else {
+ return false;
+ }
+ }
+ /**
+ * Return name std::string identifier of this refinement operator.
+ */
+ const std::string& getOperatorName() const
+ {
+ return d_name_id;
+ }
+
+ /**
+ * The priority of cell-centered double linear interpolation is 0.
+ * It will be performed before any user-defined interpolation operations.
+ */
+ int getOperatorPriority() const
+ {
+ return 0;
+ }
+
+ /**
+ * The stencil width of the linear interpolation operator is the vector
+ * of ones. That is, its stencil extends one cell outside the fine box.
+ */
+ hier::IntVector getStencilWidth() const
+ {
+ return hier::IntVector::getOne(getDim());
+ }
+
+ /**
+ * Refine the source component on the coarse patch to the destination
+ * component on the fine patch using the cell-centered double linear
+ * interpolation operator. Interpolation is performed on the intersection
+ * of the destination patch and the boxes contained in fine_overlap
+ * It is assumed that the coarse patch contains sufficient data for the
+ * stencil width of the refinement operator.
+ */
+ void refine(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::BoxOverlap& fine_overlap,
+ const hier::IntVector& ratio) const;
+
+ /**
+ * Refine the source component on the coarse patch to the destination
+ * component on the fine patch using the cell-centered double linear
+ * interpolation operator. Interpolation is performed on the intersection
+ * of the destination patch and the fine box. It is assumed that the
+ * coarse patch contains sufficient data for the stencil width of the
+ * refinement operator. This differs from the above refine() method
+ * only in that it operates on a single fine box instead of a BoxOverlap.
+ */
+ void refine(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio) const;
+
+private:
+ std::string d_name_id;
+
+};
+
+}
+}
+#endif
diff -r 4c1e4cd9d6de -r f4b8acb334a5 StokesFACOps.I
--- a/StokesFACOps.I Mon Jan 10 13:42:56 2011 -0800
+++ b/StokesFACOps.I Tue Jan 11 15:56:32 2011 -0800
@@ -43,7 +43,7 @@ void StokesFACOps::setPhysicalBcCoefObje
const RobinBcCoefStrategy* physical_bc_coef)
{
d_physical_bc_coef = physical_bc_coef;
- d_bc_helper.setCoefImplementation(physical_bc_coef);
+ // d_bc_helper.setCoefImplementation(physical_bc_coef);
#ifdef HAVE_HYPRE
d_hypre_solver.setPhysicalBcCoefObject(d_physical_bc_coef);
#endif
diff -r 4c1e4cd9d6de -r f4b8acb334a5 StokesFACOps.h
--- a/StokesFACOps.h Mon Jan 10 13:42:56 2011 -0800
+++ b/StokesFACOps.h Tue Jan 11 15:56:32 2011 -0800
@@ -990,7 +990,7 @@ private:
* to set and whether we are setting data with homogeneous
* boundary condition.
*/
- CartesianRobinBcHelper d_bc_helper;
+ // CartesianRobinBcHelper d_bc_helper;
//@{
/*!
diff -r 4c1e4cd9d6de -r f4b8acb334a5 StokesFACOps/StokesFACOps.C
--- a/StokesFACOps/StokesFACOps.C Mon Jan 10 13:42:56 2011 -0800
+++ b/StokesFACOps/StokesFACOps.C Tue Jan 11 15:56:32 2011 -0800
@@ -92,8 +92,8 @@ namespace SAMRAI {
),
d_cf_discretization("Ewing"),
- p_prolongation_method("CONSTANT_REFINE"),
- v_prolongation_method("CONSTANT_REFINE"),
+ p_prolongation_method("P_REFINE"),
+ v_prolongation_method("CONSERVATIVE_LINEAR_REFINE"),
d_coarse_solver_tolerance(1.e-8),
d_coarse_solver_max_iterations(10),
d_residual_tolerance_during_smoothing(-1.0),
@@ -145,8 +145,8 @@ namespace SAMRAI {
v_nocoarse_refine_operator(),
v_nocoarse_refine_algorithm(),
v_nocoarse_refine_schedules(),
- d_bc_helper(dim,
- d_object_name + "::bc helper"),
+ // d_bc_helper(dim,
+ // d_object_name + "::bc helper"),
d_enable_logging(false),
d_preconditioner(NULL),
d_hopscell(),
@@ -235,11 +235,11 @@ namespace SAMRAI {
d_cf_discretization);
p_prolongation_method =
- database->getStringWithDefault("prolongation_method",
+ database->getStringWithDefault("p_prolongation_method",
p_prolongation_method);
v_prolongation_method =
- database->getStringWithDefault("prolongation_method",
+ database->getStringWithDefault("v_prolongation_method",
v_prolongation_method);
d_enable_logging =
diff -r 4c1e4cd9d6de -r f4b8acb334a5 StokesFACOps/initializeOperatorState.C
--- a/StokesFACOps/initializeOperatorState.C Mon Jan 10 13:42:56 2011 -0800
+++ b/StokesFACOps/initializeOperatorState.C Tue Jan 11 15:56:32 2011 -0800
@@ -258,7 +258,7 @@ void SAMRAI::solv::StokesFACOps::initial
vdb->mapIndexToVariable(d_cell_scratch_id, variable);
p_ghostfill_refine_operator =
geometry->lookupRefineOperator(variable,
- "CONSERVATIVE_LINEAR_REFINE");
+ "LINEAR_REFINE");
vdb->mapIndexToVariable(d_side_scratch_id, variable);
v_ghostfill_refine_operator =
@@ -268,7 +268,7 @@ void SAMRAI::solv::StokesFACOps::initial
vdb->mapIndexToVariable(d_cell_scratch_id, variable);
p_nocoarse_refine_operator =
geometry->lookupRefineOperator(variable,
- "CONSERVATIVE_LINEAR_REFINE");
+ "LINEAR_REFINE");
vdb->mapIndexToVariable(d_side_scratch_id, variable);
v_nocoarse_refine_operator =
@@ -417,8 +417,8 @@ void SAMRAI::solv::StokesFACOps::initial
d_hierarchy->getPatchLevel(dest_ln),
tbox::Pointer<hier::PatchLevel>(),
dest_ln - 1,
- d_hierarchy,
- &d_bc_helper);
+ d_hierarchy);
+ // &d_bc_helper);
if (!p_prolongation_refine_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a refine schedule for p prolongation!\n");
@@ -429,8 +429,8 @@ void SAMRAI::solv::StokesFACOps::initial
d_hierarchy->getPatchLevel(dest_ln),
tbox::Pointer<hier::PatchLevel>(),
dest_ln - 1,
- d_hierarchy,
- &d_bc_helper);
+ d_hierarchy);
+ // &d_bc_helper);
if (!v_prolongation_refine_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a refine schedule for v prolongation!\n");
@@ -439,8 +439,8 @@ void SAMRAI::solv::StokesFACOps::initial
p_ghostfill_refine_algorithm->
createSchedule(d_hierarchy->getPatchLevel(dest_ln),
dest_ln - 1,
- d_hierarchy,
- &d_bc_helper);
+ d_hierarchy);
+ // &d_bc_helper);
if (!p_ghostfill_refine_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a refine schedule for ghost filling!\n");
@@ -449,16 +449,16 @@ void SAMRAI::solv::StokesFACOps::initial
v_ghostfill_refine_algorithm->
createSchedule(d_hierarchy->getPatchLevel(dest_ln),
dest_ln - 1,
- d_hierarchy,
- &d_bc_helper);
+ d_hierarchy);
+ // &d_bc_helper);
if (!v_ghostfill_refine_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a refine schedule for ghost filling!\n");
}
p_nocoarse_refine_schedules[dest_ln] =
p_nocoarse_refine_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- &d_bc_helper);
+ createSchedule(d_hierarchy->getPatchLevel(dest_ln));
+ // &d_bc_helper);
if (!p_nocoarse_refine_schedules[dest_ln]) {
TBOX_ERROR(
d_object_name
@@ -467,8 +467,8 @@ void SAMRAI::solv::StokesFACOps::initial
}
v_nocoarse_refine_schedules[dest_ln] =
v_nocoarse_refine_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- &d_bc_helper);
+ createSchedule(d_hierarchy->getPatchLevel(dest_ln));
+ // &d_bc_helper);
if (!v_nocoarse_refine_schedules[dest_ln]) {
TBOX_ERROR(
d_object_name
@@ -521,8 +521,8 @@ void SAMRAI::solv::StokesFACOps::initial
}
p_nocoarse_refine_schedules[d_ln_min] =
p_nocoarse_refine_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(d_ln_min),
- &d_bc_helper);
+ createSchedule(d_hierarchy->getPatchLevel(d_ln_min));
+ // &d_bc_helper);
if (!p_nocoarse_refine_schedules[d_ln_min]) {
TBOX_ERROR(
d_object_name
@@ -531,8 +531,8 @@ void SAMRAI::solv::StokesFACOps::initial
}
v_nocoarse_refine_schedules[d_ln_min] =
v_nocoarse_refine_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(d_ln_min),
- &d_bc_helper);
+ createSchedule(d_hierarchy->getPatchLevel(d_ln_min));
+ // &d_bc_helper);
if (!v_nocoarse_refine_schedules[d_ln_min]) {
TBOX_ERROR(
d_object_name
diff -r 4c1e4cd9d6de -r f4b8acb334a5 StokesFACOps/smoothErrorByRedBlack.C
--- a/StokesFACOps/smoothErrorByRedBlack.C Mon Jan 10 13:42:56 2011 -0800
+++ b/StokesFACOps/smoothErrorByRedBlack.C Tue Jan 11 15:56:32 2011 -0800
@@ -54,6 +54,11 @@ void SAMRAI::solv::StokesFACOps::smoothE
int num_sweeps,
double residual_tolerance)
{
+ tbox::plog << "Smooth "
+ << ln << " "
+ << num_sweeps << " "
+ << "\n";
+
const int p_id(solution.getComponentDescriptorIndex(0)),
p_rhs_id(residual.getComponentDescriptorIndex(0)),
v_id(solution.getComponentDescriptorIndex(1)),
@@ -139,6 +144,33 @@ void SAMRAI::solv::StokesFACOps::smoothE
--down[1];
++right[0];
--left[0];
+
+
+ if(i!=pbox.upper(0)+1 && j!=pbox.upper(1)+1)
+ tbox::plog << "relax "
+ << i << " "
+ << j << " "
+ << (*p)(center) << " "
+ << &((*p)(center)) << " "
+ << (*v)(pdat::SideIndex(center,
+ pdat::SideIndex::X,
+ pdat::SideIndex::Lower)) << " "
+
+ << (*v)(pdat::SideIndex(center,
+ pdat::SideIndex::Y,
+ pdat::SideIndex::Lower)) << " "
+ << (*p_rhs)(center) << " "
+ << (*v_rhs)(pdat::SideIndex(center,
+ pdat::SideIndex::X,
+ pdat::SideIndex::Lower)) << " "
+
+ << (*v_rhs)(pdat::SideIndex(center,
+ pdat::SideIndex::Y,
+ pdat::SideIndex::Lower)) << " "
+ << "\n";
+
+
+
/* Update p */
if(i!=pbox.upper(0)+1 && j!=pbox.upper(1)+1)
diff -r 4c1e4cd9d6de -r f4b8acb334a5 StokesFACSolver/StokesFACSolver.C
--- a/StokesFACSolver/StokesFACSolver.C Mon Jan 10 13:42:56 2011 -0800
+++ b/StokesFACSolver/StokesFACSolver.C Tue Jan 11 15:56:32 2011 -0800
@@ -137,13 +137,13 @@ namespace SAMRAI {
}
}
- /*
- * 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 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
diff -r 4c1e4cd9d6de -r f4b8acb334a5 example_inputs/const_refine.2d.input
--- a/example_inputs/const_refine.2d.input Mon Jan 10 13:42:56 2011 -0800
+++ b/example_inputs/const_refine.2d.input Tue Jan 11 15:56:32 2011 -0800
@@ -34,9 +34,14 @@ FACStokes {
enable_logging = TRUE // Bool flag to switch logging on/off
max_cycles = 20 // Max number of FAC cycles to use
residual_tol = 1e-8 // Residual tolerance to solve for
- num_pre_sweeps = 1 // Number of presmoothing sweeps to use
- num_post_sweeps = 3 // Number of postsmoothing sweeps to use
- prolongation_method = "CONSTANT_REFINE" // Type of refinement
+ num_pre_sweeps = 0 // Number of presmoothing sweeps to use
+ num_post_sweeps = 5 // Number of postsmoothing sweeps to use
+ p_prolongation_method = "P_REFINE" // Type of refinement
+ // used in prolongation.
+ // Suggested values are
+ // "LINEAR_REFINE"
+ // "CONSTANT_REFINE"
+ v_prolongation_method = "CONSERVATIVE_LINEAR_REFINE" // Type of refinement
// used in prolongation.
// Suggested values are
// "LINEAR_REFINE"
diff -r 4c1e4cd9d6de -r f4b8acb334a5 main.C
--- a/main.C Mon Jan 10 13:42:56 2011 -0800
+++ b/main.C Tue Jan 11 15:56:32 2011 -0800
@@ -28,6 +28,8 @@ using namespace std;
#include "SAMRAI/tbox/TimerManager.h"
#include "SAMRAI/tbox/Utilities.h"
#include "SAMRAI/appu/VisItDataWriter.h"
+#include "SAMRAI/xfer/RefineOperator.h"
+#include "P_Refine.h"
#include "FACStokes.h"
@@ -148,6 +150,9 @@ int main(
input_db->getDatabase("CartesianGridGeometry")));
tbox::plog << "Cartesian Geometry:" << endl;
grid_geometry->printClassData(tbox::plog);
+ grid_geometry->addSpatialRefineOperator(tbox::Pointer<SAMRAI::xfer::RefineOperator>(new SAMRAI::geom::P_Refine(dim)));
+
+
tbox::Pointer<hier::PatchHierarchy>
patch_hierarchy(new hier::PatchHierarchy
More information about the CIG-COMMITS
mailing list