[cig-commits] commit: Adaptive runs, but pressure boundary conditions are not set correctly
Mercurial
hg at geodynamics.org
Fri Feb 25 14:16:01 PST 2011
changeset: 71:589e7f172c2b
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 02 15:40:08 2011 -0800
files: Boundary.h FACStokes/initializeLevelData.C Makefile P_Refine.C StokesFACOps.I StokesFACOps.h StokesFACOps/StokesFACOps.C StokesFACOps/Update_V.C StokesFACOps/computeCompositeResidualOnLevel.C StokesFACOps/initializeOperatorState.C StokesFACOps/prolongErrorAndCorrect.C StokesFACOps/restrictResidual.C StokesFACOps/restrictSolution.C StokesFACOps/set_boundaries.C StokesFACOps/smoothErrorByRedBlack.C StokesFACOps/xeqScheduleGhostFill.C StokesFACOps/xeqScheduleGhostFillNoCoarse.C StokesFACSolver/setBcObject.C StokesFACSolver/setBoundaries.C V_Boundary_Refine/Update_V.C V_Boundary_Refine/refine.C V_Coarsen.C V_Coarsen_Patch_Strategy.C V_Coarsen_Patch_Strategy.h V_Refine.C V_Refine_Patch_Strategy.C V_Refine_Patch_Strategy.h example_inputs/const_refine.2d.input
description:
Adaptive runs, but pressure boundary conditions are not set correctly
diff -r cf514c7ba5e6 -r 589e7f172c2b Boundary.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Boundary.h Wed Feb 02 15:40:08 2011 -0800
@@ -0,0 +1,4 @@
+#ifndef GAMR_BOUNDARY_H
+#define GAMR_BOUNDARY_H
+const double boundary_value=1e100;
+#endif
diff -r cf514c7ba5e6 -r 589e7f172c2b FACStokes/initializeLevelData.C
--- a/FACStokes/initializeLevelData.C Thu Jan 27 14:33:28 2011 -0800
+++ b/FACStokes/initializeLevelData.C Wed Feb 02 15:40:08 2011 -0800
@@ -88,6 +88,17 @@ namespace SAMRAI {
p_rhs_data->fill(0.0);
+ // for(pdat::CellIterator ci(p_rhs_data->getBox()); ci; ci++)
+ // {
+ // pdat::CellIndex cc=ci();
+ // tbox::plog << "p_rhs "
+ // << cc[0] << " "
+ // << cc[1] << " "
+ // << (*p_rhs_data)(cc) << " "
+ // << (&(*p_rhs_data)(cc)) << " "
+ // << "\n";
+ // }
+
tbox::Pointer<pdat::SideData<double> > v_rhs_data =
patch->getPatchData(v_rhs_id);
diff -r cf514c7ba5e6 -r 589e7f172c2b Makefile
--- a/Makefile Thu Jan 27 14:33:28 2011 -0800
+++ b/Makefile Wed Feb 02 15:40:08 2011 -0800
@@ -33,6 +33,8 @@ CXX_OBJS = main.o FACStokes/FACStok
P_Refine.o V_Refine.o V_Coarsen.o \
V_Boundary_Refine/refine.o \
V_Boundary_Refine/Update_V.o \
+ V_Refine_Patch_Strategy.o \
+ V_Coarsen_Patch_Strategy.o \
StokesFACOps/StokesFACOps.o \
StokesFACOps/checkInputPatchDataIndices.o \
StokesFACOps/computeCompositeResidualOnLevel.o \
diff -r cf514c7ba5e6 -r 589e7f172c2b P_Refine.C
--- a/P_Refine.C Thu Jan 27 14:33:28 2011 -0800
+++ b/P_Refine.C Wed Feb 02 15:40:08 2011 -0800
@@ -89,8 +89,11 @@ void SAMRAI::geom::P_Refine::refine(
/* 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. */
+ /* This assumes that the levels are always properly nested,
+ so that we always have an extra grid space for
+ interpolation. So we only have to have a special case for
+ physical boundaries, where we do not have an extra grid
+ space. */
if(center[0]==coarse_box.lower(0)
&& geom->getTouchesRegularBoundary(0,0))
@@ -125,6 +128,16 @@ void SAMRAI::geom::P_Refine::refine(
(*p_fine)(fine)=(*p)(center)
+ ((i%2==0) ? (-dp_dx) : dp_dx)
+ ((j%2==0) ? (-dp_dy) : dp_dy);
+
+ // tbox::plog << "P_Refine "
+ // << i << " "
+ // << j << " "
+ // << (*p_fine)(fine) << " "
+ // << (*p)(center) << " "
+ // << dp_dx << " "
+ // << dp_dy << " "
+ // << "\n";
+
}
}
#endif
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps.I
--- a/StokesFACOps.I Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps.I Wed Feb 02 15:40:08 2011 -0800
@@ -38,16 +38,16 @@ void StokesFACOps::setUseSMG(
********************************************************************
*/
-SAMRAI_INLINE_KEYWORD
-void StokesFACOps::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
-}
+// SAMRAI_INLINE_KEYWORD
+// void StokesFACOps::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
+// }
/*
********************************************************************
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps.h
--- a/StokesFACOps.h Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps.h Wed Feb 02 15:40:08 2011 -0800
@@ -12,10 +12,8 @@
#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 "StokesHypreSolver.h"
#include "SAMRAI/solv/SAMRAIVectorReal.h"
#include "StokesSpecifications.h"
@@ -44,6 +42,8 @@
#include "SAMRAI/tbox/Database.h"
#include "SAMRAI/tbox/Pointer.h"
#include "SAMRAI/tbox/Timer.h"
+#include "V_Refine_Patch_Strategy.h"
+#include "V_Coarsen_Patch_Strategy.h"
#include <string>
@@ -320,9 +320,9 @@ public:
* @param physical_bc_coef tbox::Pointer to an object that can
* set the Robin bc coefficients.
*/
- void
- setPhysicalBcCoefObject(
- const RobinBcCoefStrategy* physical_bc_coef);
+ // void
+ // setPhysicalBcCoefObject(
+ // const RobinBcCoefStrategy* physical_bc_coef);
//@{
@@ -483,6 +483,11 @@ public:
const SAMRAIVectorReal<double>& current_soln,
const SAMRAIVectorReal<double>& residual);
+ void set_boundaries(const int &v_id, const int &l)
+ {
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(l);
+ set_boundaries(v_id,level);
+ }
void set_boundaries(const int &v_id, tbox::Pointer<hier::PatchLevel> &level);
//@}
@@ -492,6 +497,8 @@ private:
/*!
* @name Private workhorse functions.
*/
+
+ const int invalid_id;
/*!
* @brief Red-black Gauss-Seidel error smoothing on a level.
@@ -872,7 +879,7 @@ private:
*
* see setPhysicalBcCoefObject()
*/
- const RobinBcCoefStrategy* d_physical_bc_coef;
+ // const RobinBcCoefStrategy* d_physical_bc_coef;
//@}
@@ -1008,11 +1015,12 @@ private:
* 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
+ * we must tell the patch strategies the patch data index we want
* to set and whether we are setting data with homogeneous
* boundary condition.
*/
- // CartesianRobinBcHelper d_bc_helper;
+ V_Refine_Patch_Strategy v_refine_patch_strategy;
+ V_Coarsen_Patch_Strategy v_coarsen_patch_strategy;
//@{
/*!
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/StokesFACOps.C
--- a/StokesFACOps/StokesFACOps.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/StokesFACOps.C Wed Feb 02 15:40:08 2011 -0800
@@ -105,13 +105,14 @@ namespace SAMRAI {
database->getDatabase("hypre_solver"):
tbox::Pointer<tbox::Database>(NULL)),
#endif
- d_physical_bc_coef(NULL),
+ // d_physical_bc_coef(NULL),
d_context(hier::VariableDatabase::getDatabase()
->getContext(object_name + "::PRIVATE_CONTEXT")),
d_cell_scratch_id(-1),
d_side_scratch_id(-1),
d_flux_scratch_id(-1),
d_oflux_scratch_id(-1),
+ invalid_id(-1),
p_prolongation_refine_operator(),
p_prolongation_refine_algorithm(),
p_prolongation_refine_schedules(),
@@ -145,8 +146,10 @@ namespace SAMRAI {
v_nocoarse_refine_operator(),
v_nocoarse_refine_algorithm(),
v_nocoarse_refine_schedules(),
- // d_bc_helper(dim,
- // d_object_name + "::bc helper"),
+ v_refine_patch_strategy(dim,
+ d_object_name + "::refine patch strategy"),
+ v_coarsen_patch_strategy(dim,
+ d_object_name + "::coarsen patch strategy"),
d_enable_logging(false),
d_preconditioner(NULL),
d_hopscell(),
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/Update_V.C
--- a/StokesFACOps/Update_V.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/Update_V.C Wed Feb 02 15:40:08 2011 -0800
@@ -40,6 +40,7 @@
#include "SAMRAI/xfer/RefineSchedule.h"
#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
+#include "Boundary.h"
/*
********************************************************************
* Updates one component of the velocity during a red-black *
@@ -73,13 +74,13 @@ void SAMRAI::solv::StokesFACOps::Update_
&& (*v)(pdat::SideIndex(left,
axis,
pdat::SideIndex::Lower))
- ==std::numeric_limits<double>::max())
+ ==boundary_value)
|| (center[axis]==pbox.upper(axis)+1
&& (*v)(pdat::SideIndex
(right,
axis,
pdat::SideIndex::Lower))
- ==std::numeric_limits<double>::max())))
+ ==boundary_value)))
{
double dp_dx, d2vx_dxx, d2vx_dyy, C_vx;
/* If y==0 */
@@ -140,6 +141,13 @@ void SAMRAI::solv::StokesFACOps::Update_
pdat::SideIndex::Lower))
- viscosity*(d2vx_dxx + d2vx_dyy) + dp_dx;
+
+ tbox::plog << "v " << axis << " "
+ << (*v_rhs)(pdat::SideIndex(center,
+ axis,
+ pdat::SideIndex::Lower))
+ << " ";
+
// tbox::plog << "Update "
// << axis << " "
// << off_axis << " "
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/computeCompositeResidualOnLevel.C
--- a/StokesFACOps/computeCompositeResidualOnLevel.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Wed Feb 02 15:40:08 2011 -0800
@@ -40,6 +40,7 @@
#include "SAMRAI/xfer/RefineSchedule.h"
#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
+#include "Boundary.h"
/*
********************************************************************
* FACOperatorStrategy virtual *
@@ -73,8 +74,8 @@ void SAMRAI::solv::StokesFACOps::compute
*/
const int p_id = solution.getComponentDescriptorIndex(0);
const int v_id = solution.getComponentDescriptorIndex(1);
- // d_bc_helper.setTargetDataId(soln_id);
- // d_bc_helper.setHomogeneousBc(error_equation_indicator);
+ v_refine_patch_strategy.setTargetDataId(v_id);
+ // v_refine_patch_strategy.setHomogeneousBc(error_equation_indicator);
/*
* Assumptions:
@@ -114,6 +115,8 @@ void SAMRAI::solv::StokesFACOps::compute
// }
/* S1. Fill solution ghost data. */
+
+ set_boundaries(v_id,ln);
if (ln > d_ln_min) {
/* Fill from current, next coarser level and physical boundary */
xeqScheduleGhostFill(p_id, v_id, ln);
@@ -214,9 +217,15 @@ void SAMRAI::solv::StokesFACOps::compute
{
/* If x==0 */
if((center[0]==pbox.lower(0)
- && geom->getTouchesRegularBoundary(0,0))
+ && (*v)(pdat::SideIndex(left,
+ pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ ==boundary_value)
|| (center[0]==pbox.upper(0)+1
- && geom->getTouchesRegularBoundary(0,1)))
+ && (*v)(pdat::SideIndex(right,
+ pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ ==boundary_value))
{
(*v_resid)(pdat::SideIndex(center,pdat::SideIndex::X,
pdat::SideIndex::Lower))=0;
@@ -224,44 +233,17 @@ void SAMRAI::solv::StokesFACOps::compute
else
{
double dp_dx, d2vx_dxx, d2vx_dyy, C_vx;
- /* If y==0 */
- if(center[1]==pbox.lower(1)
- && geom->getTouchesRegularBoundary(1,0))
- {
- d2vx_dyy=
- ((*v)(pdat::SideIndex(up,pdat::SideIndex::X,
+ d2vx_dyy=
+ ((*v)(pdat::SideIndex(up,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ - 2*(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
pdat::SideIndex::Lower))
- - (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- /(dy*dy);
- C_vx=-viscosity*(2/(dx*dx) + 1/(dy*dy));
- }
- /* If y==max_y */
- else if(center[1]==pbox.upper(1)
- && geom->getTouchesRegularBoundary(1,1))
- {
- d2vx_dyy=
- (-(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex(down,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- /(dy*dy);
- C_vx=-viscosity*(2/(dx*dx) + 1/(dy*dy));
- // tbox::plog << "vx y1 boundary ";
- }
- else
- {
- d2vx_dyy=
- ((*v)(pdat::SideIndex(up,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- - 2*(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex(down,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- /(dy*dy);
+ + (*v)(pdat::SideIndex(down,pdat::SideIndex::X,
+ pdat::SideIndex::Lower)))
+ /(dy*dy);
- C_vx=-2*viscosity*(1/(dx*dx) + 1/(dy*dy));
- }
+ C_vx=-2*viscosity*(1/(dx*dx) + 1/(dy*dy));
+
d2vx_dxx=((*v)(pdat::SideIndex(left,pdat::SideIndex::X,
pdat::SideIndex::Lower))
- 2*(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
@@ -278,6 +260,30 @@ void SAMRAI::solv::StokesFACOps::compute
pdat::SideIndex::Lower))
- viscosity*(d2vx_dxx + d2vx_dyy) + dp_dx;
}
+
+
+
+ tbox::plog << "resid "
+ << i << " "
+ << j << " "
+ // << (*p_resid)(center) << " "
+ // << (*p_rhs)(center) << " "
+ // << dvx_dx << " "
+ // << dvy_dy << " "
+ << "vx "
+ << (*v_resid)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ << " "
+ // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Upper))
+ // << " "
+ // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))
+ // << " "
+ // << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Upper)))
+ // << " "
+ << "\n";
}
/* vy */
@@ -285,9 +291,15 @@ void SAMRAI::solv::StokesFACOps::compute
{
/* If y==0 */
if((center[1]==pbox.lower(1)
- && geom->getTouchesRegularBoundary(1,0))
+ && (*v)(pdat::SideIndex(down,
+ pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ ==boundary_value)
|| (center[1]==pbox.upper(1)+1
- && geom->getTouchesRegularBoundary(1,1)))
+ && (*v)(pdat::SideIndex(up,
+ pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ ==boundary_value))
{
(*v_resid)(pdat::SideIndex(center,pdat::SideIndex::Y,
pdat::SideIndex::Lower))=0;
@@ -295,43 +307,16 @@ void SAMRAI::solv::StokesFACOps::compute
else
{
double dp_dy, d2vy_dxx, d2vy_dyy, C_vy;
- /* If x==0 */
- if(center[0]==pbox.lower(0)
- && geom->getTouchesRegularBoundary(0,0))
- {
- d2vy_dxx=
- ((*v)(pdat::SideIndex(right,pdat::SideIndex::Y,
+ d2vy_dxx=
+ ((*v)(pdat::SideIndex(left,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ - 2*(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
pdat::SideIndex::Lower))
- - (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- /(dx*dx);
- C_vy=-viscosity*(1/(dx*dx) + 2/(dy*dy));
- }
- /* If x==max_x */
- else if(center[0]==pbox.upper(0)
- && geom->getTouchesRegularBoundary(0,1))
- {
- d2vy_dxx=
- ((*v)(pdat::SideIndex(left,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- - (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- /(dx*dx);
- C_vy=-viscosity*(1/(dx*dx) + 2/(dy*dy));
- }
- else
- {
- d2vy_dxx=
- ((*v)(pdat::SideIndex(left,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- - 2*(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex(right,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- /(dx*dx);
+ + (*v)(pdat::SideIndex(right,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower)))
+ /(dx*dx);
- C_vy=-2*viscosity*(1/(dx*dx) + 1/(dy*dy));
- }
+ C_vy=-2*viscosity*(1/(dx*dx) + 1/(dy*dy));
d2vy_dyy=((*v)(pdat::SideIndex(up,pdat::SideIndex::Y,
pdat::SideIndex::Lower))
- 2*(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
@@ -369,6 +354,13 @@ void SAMRAI::solv::StokesFACOps::compute
// }
}
+ /* We also need to set the boundaries of the rhs so that coarsening
+ works correctly. */
+ const int v_rhs_id = rhs.getComponentDescriptorIndex(1);
+ set_boundaries(v_rhs_id,ln);
+ xeqScheduleGhostFillNoCoarse(invalid_id, v_rhs_id, ln);
+
+
// if (deallocate_flux_data_when_done) {
// level->deallocatePatchData(flux_id);
// }
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/initializeOperatorState.C
--- a/StokesFACOps/initializeOperatorState.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/initializeOperatorState.C Wed Feb 02 15:40:08 2011 -0800
@@ -69,22 +69,22 @@ void SAMRAI::solv::StokesFACOps::initial
#ifdef DEBUG_CHECK_ASSERTIONS
- if (d_physical_bc_coef == NULL) {
- /*
- * It's an error not to have bc object set.
- * Note that the bc object cannot be passed in through
- * the argument because the interface is inherited.
- */
- TBOX_ERROR(
- d_object_name << ": No physical bc object in\n"
- <<
- "StokesFACOps::initializeOperatorState\n"
- << "You must use "
- <<
- "StokesFACOps::setPhysicalBcCoefObject\n"
- <<
- "to set one before calling initializeOperatorState\n");
- }
+ // if (d_physical_bc_coef == NULL) {
+ // /*
+ // * It's an error not to have bc object set.
+ // * Note that the bc object cannot be passed in through
+ // * the argument because the interface is inherited.
+ // */
+ // TBOX_ERROR(
+ // d_object_name << ": No physical bc object in\n"
+ // <<
+ // "StokesFACOps::initializeOperatorState\n"
+ // << "You must use "
+ // <<
+ // "StokesFACOps::setPhysicalBcCoefObject\n"
+ // <<
+ // "to set one before calling initializeOperatorState\n");
+ // }
if (solution.getNumberOfComponents() != 1) {
TBOX_WARNING(d_object_name
@@ -201,18 +201,18 @@ void SAMRAI::solv::StokesFACOps::initial
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
+// #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.
@@ -401,6 +401,7 @@ void SAMRAI::solv::StokesFACOps::initial
solution.getComponentDescriptorIndex(1),
v_nocoarse_refine_operator);
+ /* Refinement and ghost fill operators */
for (int dest_ln = d_ln_min + 1; dest_ln <= d_ln_max; ++dest_ln) {
tbox::Pointer<xfer::PatchLevelFullFillPattern>
@@ -412,7 +413,7 @@ void SAMRAI::solv::StokesFACOps::initial
tbox::Pointer<hier::PatchLevel>(),
dest_ln - 1,
d_hierarchy);
- // &d_bc_helper);
+ // &v_refine_patch_strategy);
if (!p_prolongation_refine_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a refine schedule for p prolongation!\n");
@@ -424,7 +425,7 @@ void SAMRAI::solv::StokesFACOps::initial
tbox::Pointer<hier::PatchLevel>(),
dest_ln - 1,
d_hierarchy);
- // &d_bc_helper);
+ // &v_refine_patch_strategy);
if (!v_prolongation_refine_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a refine schedule for v prolongation!\n");
@@ -434,7 +435,7 @@ void SAMRAI::solv::StokesFACOps::initial
createSchedule(d_hierarchy->getPatchLevel(dest_ln),
dest_ln - 1,
d_hierarchy);
- // &d_bc_helper);
+ // &v_refine_patch_strategy);
if (!p_ghostfill_refine_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a refine schedule for ghost filling!\n");
@@ -443,8 +444,9 @@ 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_hierarchy,
+ &v_refine_patch_strategy);
if (!v_ghostfill_refine_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a refine schedule for ghost filling!\n");
@@ -452,7 +454,7 @@ void SAMRAI::solv::StokesFACOps::initial
p_nocoarse_refine_schedules[dest_ln] =
p_nocoarse_refine_algorithm->
createSchedule(d_hierarchy->getPatchLevel(dest_ln));
- // &d_bc_helper);
+ // &v_refine_patch_strategy);
if (!p_nocoarse_refine_schedules[dest_ln]) {
TBOX_ERROR(
d_object_name
@@ -462,7 +464,7 @@ void SAMRAI::solv::StokesFACOps::initial
v_nocoarse_refine_schedules[dest_ln] =
v_nocoarse_refine_algorithm->
createSchedule(d_hierarchy->getPatchLevel(dest_ln));
- // &d_bc_helper);
+ // &v_refine_patch_strategy);
if (!v_nocoarse_refine_schedules[dest_ln]) {
TBOX_ERROR(
d_object_name
@@ -470,6 +472,8 @@ void SAMRAI::solv::StokesFACOps::initial
": Cannot create a refine schedule for ghost filling on bottom level!\n");
}
}
+
+ /* Coarsening operators */
for (int dest_ln = d_ln_min; dest_ln < d_ln_max; ++dest_ln) {
p_urestriction_coarsen_schedules[dest_ln] =
p_urestriction_coarsen_algorithm->
@@ -491,7 +495,8 @@ void SAMRAI::solv::StokesFACOps::initial
v_urestriction_coarsen_schedules[dest_ln] =
v_urestriction_coarsen_algorithm->
createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- d_hierarchy->getPatchLevel(dest_ln + 1));
+ d_hierarchy->getPatchLevel(dest_ln + 1),
+ &v_coarsen_patch_strategy);
if (!v_urestriction_coarsen_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a coarsen schedule for U v restriction!\n");
@@ -499,7 +504,8 @@ void SAMRAI::solv::StokesFACOps::initial
v_rrestriction_coarsen_schedules[dest_ln] =
v_rrestriction_coarsen_algorithm->
createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- d_hierarchy->getPatchLevel(dest_ln + 1));
+ d_hierarchy->getPatchLevel(dest_ln + 1),
+ &v_coarsen_patch_strategy);
if (!v_rrestriction_coarsen_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a coarsen schedule for R v restriction!\n");
@@ -513,10 +519,12 @@ void SAMRAI::solv::StokesFACOps::initial
<< ": Cannot create a coarsen schedule for flux transfer!\n");
}
}
+
+ /* Ordinary ghost fill operator on the coarsest level */
p_nocoarse_refine_schedules[d_ln_min] =
p_nocoarse_refine_algorithm->
createSchedule(d_hierarchy->getPatchLevel(d_ln_min));
- // &d_bc_helper);
+ // &v_refine_patch_strategy);
if (!p_nocoarse_refine_schedules[d_ln_min]) {
TBOX_ERROR(
d_object_name
@@ -526,7 +534,7 @@ 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);
+ // &v_refine_patch_strategy);
if (!v_nocoarse_refine_schedules[d_ln_min]) {
TBOX_ERROR(
d_object_name
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/prolongErrorAndCorrect.C
--- a/StokesFACOps/prolongErrorAndCorrect.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/prolongErrorAndCorrect.C Wed Feb 02 15:40:08 2011 -0800
@@ -76,13 +76,19 @@ void SAMRAI::solv::StokesFACOps::prolong
fine_level->allocatePatchData(d_cell_scratch_id);
fine_level->allocatePatchData(d_side_scratch_id);
+ // int p_src(s.getComponentDescriptorIndex(0)),
+ // v_src(s.getComponentDescriptorIndex(1)),
+ // p_dst(d.getComponentDescriptorIndex(0)),
+ // v_dst(d.getComponentDescriptorIndex(1));
+ // xeqScheduleGhostFillNoCoarse(invalid_id,v_src,dest_ln+1);
+
/*
* 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);
+ v_refine_patch_strategy.setTargetDataId(d_side_scratch_id);
+ // v_refine_patch_strategy.setHomogeneousBc(true);
xeqScheduleProlongation(d_cell_scratch_id,
s.getComponentDescriptorIndex(0),
d_cell_scratch_id,
@@ -91,6 +97,7 @@ void SAMRAI::solv::StokesFACOps::prolong
d_side_scratch_id,
dest_ln);
+ set_boundaries(s.getComponentDescriptorIndex(1),fine_level);
/*
* Add the refined error in the scratch space to the error currently
* residing in the destination level.
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/restrictResidual.C
--- a/StokesFACOps/restrictResidual.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/restrictResidual.C Wed Feb 02 15:40:08 2011 -0800
@@ -59,7 +59,8 @@ void SAMRAI::solv::StokesFACOps::restric
v_dst(d.getComponentDescriptorIndex(1));
/* Need to do a sync because the coarsening for v uses ghost zones */
- xeqScheduleGhostFillNoCoarse(-1,v_src,dest_ln+1);
+ v_coarsen_patch_strategy.setSourceDataId(v_src);
+ xeqScheduleGhostFillNoCoarse(invalid_id,v_src,dest_ln+1);
xeqScheduleRRestriction(p_dst,p_src,v_dst,v_src,dest_ln);
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/restrictSolution.C
--- a/StokesFACOps/restrictSolution.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/restrictSolution.C Wed Feb 02 15:40:08 2011 -0800
@@ -61,12 +61,15 @@ void SAMRAI::solv::StokesFACOps::restric
v_dst(d.getComponentDescriptorIndex(1));
/* Need to do a sync because the coarsening for v uses ghost zones. */
- xeqScheduleGhostFillNoCoarse(-1,v_src,dest_ln+1);
+ v_coarsen_patch_strategy.setSourceDataId(v_src);
+ xeqScheduleGhostFillNoCoarse(invalid_id,v_src,dest_ln+1);
xeqScheduleURestriction(p_dst,p_src,v_dst,v_src,dest_ln);
- // d_bc_helper.setHomogeneousBc(false);
- // d_bc_helper.setTargetDataId(d.getComponentDescriptorIndex(0));
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(dest_ln);
+ set_boundaries(v_dst,level);
+ // v_refine_patch_strategy.setHomogeneousBc(false);
+ v_refine_patch_strategy.setTargetDataId(d.getComponentDescriptorIndex(1));
if (dest_ln == d_ln_min) {
xeqScheduleGhostFillNoCoarse(p_dst,v_dst,dest_ln);
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/set_boundaries.C
--- a/StokesFACOps/set_boundaries.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/set_boundaries.C Wed Feb 02 15:40:08 2011 -0800
@@ -40,6 +40,7 @@
#include "SAMRAI/xfer/RefineSchedule.h"
#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
+#include "Boundary.h"
/*
********************************************************************
* Workhorse function to smooth error using red-black *
@@ -89,7 +90,7 @@ void SAMRAI::solv::StokesFACOps::set_bou
{
(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
pdat::SideIndex::Lower))=
- std::numeric_limits<double>::max();
+ boundary_value;
}
/* Set the value so the derivative=0 */
else if(j<pbox.lower(0))
@@ -106,6 +107,17 @@ void SAMRAI::solv::StokesFACOps::set_bou
(*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
pdat::SideIndex::Lower));
}
+ tbox::plog << "set boundary vx "
+ << level->getLevelNumber() << " "
+ << i << " "
+ << j << " "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ << " "
+ << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower)))
+ << " "
+ << "\n";
}
/* vy */
if(i<=gbox.upper(0))
@@ -114,7 +126,7 @@ void SAMRAI::solv::StokesFACOps::set_bou
{
(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
pdat::SideIndex::Lower))=
- std::numeric_limits<double>::max();
+ boundary_value;
}
else if(i<pbox.lower(0))
{
@@ -130,6 +142,17 @@ void SAMRAI::solv::StokesFACOps::set_bou
(*v)(pdat::SideIndex(center-ip,pdat::SideIndex::Y,
pdat::SideIndex::Lower));
}
+ tbox::plog << "set boundary vy "
+ << level->getLevelNumber() << " "
+ << i << " "
+ << j << " "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ << " "
+ << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower)))
+ << " "
+ << "\n";
}
}
}
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/smoothErrorByRedBlack.C
--- a/StokesFACOps/smoothErrorByRedBlack.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/smoothErrorByRedBlack.C Wed Feb 02 15:40:08 2011 -0800
@@ -40,6 +40,7 @@
#include "SAMRAI/xfer/RefineSchedule.h"
#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
+#include "Boundary.h"
/*
********************************************************************
* Workhorse function to smooth error using red-black *
@@ -59,12 +60,6 @@ void SAMRAI::solv::StokesFACOps::smoothE
v_id(solution.getComponentDescriptorIndex(1)),
v_rhs_id(residual.getComponentDescriptorIndex(1));
- /* Only need to sync the rhs once. This sync is needed because
- calculating a new pressure update requires computing in the ghost
- region so that the update for the velocity inside the box will be
- correct. */
- xeqScheduleGhostFillNoCoarse(p_rhs_id,v_rhs_id,ln);
-
checkInputPatchDataIndices();
#ifdef DEBUG_CHECK_ASSERTIONS
@@ -77,6 +72,30 @@ void SAMRAI::solv::StokesFACOps::smoothE
#endif
tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
+ {
+ hier::PatchLevel::Iterator pi(*level);
+ for (pi.initialize(*level); pi; pi++) {
+ tbox::Pointer<hier::Patch> patch = *pi;
+ tbox::Pointer<pdat::CellData<double> > p_rhs_data =
+ patch->getPatchData(p_rhs_id);
+ for(pdat::CellIterator ci(p_rhs_data->getBox()); ci; ci++)
+ {
+ pdat::CellIndex cc=ci();
+ tbox::plog << "p_rhs "
+ << cc[0] << " "
+ << cc[1] << " "
+ << (*p_rhs_data)(cc) << " "
+ // << (&(*p_rhs_data)(cc)) << " "
+ << "\n";
+ }
+ }
+ }
+ /* Only need to sync the rhs once. This sync is needed because
+ calculating a new pressure update requires computing in the ghost
+ region so that the update for the velocity inside the box will be
+ correct. */
+ v_refine_patch_strategy.setTargetDataId(v_id);
+ xeqScheduleGhostFillNoCoarse(p_rhs_id,v_rhs_id,ln);
set_boundaries(v_id,level);
if (ln > d_ln_min) {
@@ -148,46 +167,71 @@ void SAMRAI::solv::StokesFACOps::smoothE
++right[0];
--left[0];
+ tbox::plog << "smooth "
+ << ln << " "
+ << i << " "
+ << j << " "
+ << pbox.lower(0) << " "
+ << pbox.upper(0) << " "
+ << pbox.lower(1) << " "
+ << pbox.upper(1) << " ";
+
/* Update p */
- double dvx_dx=
- ((*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Upper))
- - (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))/dx;
- double dvy_dy=
- ((*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Upper))
- - (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))/dy;
+ if(!((*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ ==boundary_value
+ || (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Upper))
+ ==boundary_value
+ || (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ ==boundary_value
+ || (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Upper))
+ ==boundary_value))
+ {
+ double dvx_dx=
+ ((*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Upper))
+ - (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower)))/dx;
+ double dvy_dy=
+ ((*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Upper))
+ - (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower)))/dy;
- double delta_R_continuity=
- (*p_rhs)(center) - dvx_dx - dvy_dy;
+ double delta_R_continuity=
+ (*p_rhs)(center) - dvx_dx - dvy_dy;
- /* No scaling here, though there should be. */
- maxres=std::max(maxres,delta_R_continuity);
+ /* No scaling here, though there should be. */
+ maxres=std::max(maxres,delta_R_continuity);
- (*p)(center)+=
- viscosity*delta_R_continuity*theta_continuity;
+ (*p)(center)+=
+ viscosity*delta_R_continuity*theta_continuity;
+ tbox::plog << "p "
+ << (*p)(center) << " "
+ << (*p_rhs)(center) << " "
+ << dvx_dx << " "
+ << dvy_dy << " "
+ << delta_R_continuity << " ";
+ }
/* Update v */
Update_V(0,j,pbox,center,left,right,down,up,p,v,v_rhs,
maxres,dx,dy,viscosity,theta_momentum);
Update_V(1,i,pbox,center,down,up,left,right,p,v,v_rhs,
maxres,dy,dx,viscosity,theta_momentum);
- // tbox::plog << "smooth "
- // << i << " "
- // << j << " "
- // << (*p)(center) << " "
- // << (*v)(pdat::SideIndex
- // (center,
- // pdat::SideIndex::X,
- // pdat::SideIndex::Lower)) << " "
- // << (*v)(pdat::SideIndex
- // (center,
- // pdat::SideIndex::Y,
- // pdat::SideIndex::Lower)) << " "
- // << "\n";
+ // << (*v)(pdat::SideIndex
+ // (center,
+ // pdat::SideIndex::X,
+ // pdat::SideIndex::Lower)) << " "
+ // << (*v)(pdat::SideIndex
+ // (center,
+ // pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower)) << " "
+ tbox::plog << "\n";
}
}
}
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/xeqScheduleGhostFill.C
--- a/StokesFACOps/xeqScheduleGhostFill.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/xeqScheduleGhostFill.C Wed Feb 02 15:40:08 2011 -0800
@@ -61,6 +61,7 @@ void SAMRAI::solv::StokesFACOps::xeqSche
if (!v_ghostfill_refine_schedules[dest_ln]) {
TBOX_ERROR("Expected schedule not found.");
}
+ set_boundaries(v_id,dest_ln-1);
xfer::RefineAlgorithm refiner(d_dim);
refiner.registerRefine(v_id,v_id,v_id,v_ghostfill_refine_operator);
refiner.resetSchedule(v_ghostfill_refine_schedules[dest_ln]);
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACOps/xeqScheduleGhostFillNoCoarse.C
--- a/StokesFACOps/xeqScheduleGhostFillNoCoarse.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACOps/xeqScheduleGhostFillNoCoarse.C Wed Feb 02 15:40:08 2011 -0800
@@ -45,7 +45,7 @@ void SAMRAI::solv::StokesFACOps::xeqSche
int dest_ln)
{
/* p */
- if(p_id!=-1)
+ if(p_id!=invalid_id)
{
if (!p_nocoarse_refine_schedules[dest_ln]) {
TBOX_ERROR("Expected cell schedule not found.");
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACSolver/setBcObject.C
--- a/StokesFACSolver/setBcObject.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACSolver/setBcObject.C Wed Feb 02 15:40:08 2011 -0800
@@ -29,7 +29,7 @@ namespace SAMRAI {
}
#endif
d_bc_object = bc_object;
- d_fac_ops.setPhysicalBcCoefObject(d_bc_object);
+ // d_fac_ops.setPhysicalBcCoefObject(d_bc_object);
}
}
diff -r cf514c7ba5e6 -r 589e7f172c2b StokesFACSolver/setBoundaries.C
--- a/StokesFACSolver/setBoundaries.C Thu Jan 27 14:33:28 2011 -0800
+++ b/StokesFACSolver/setBoundaries.C Wed Feb 02 15:40:08 2011 -0800
@@ -37,7 +37,7 @@ namespace SAMRAI {
flags,
bdry_types);
d_bc_object = &d_simple_bc;
- d_fac_ops.setPhysicalBcCoefObject(d_bc_object);
+ // d_fac_ops.setPhysicalBcCoefObject(d_bc_object);
}
diff -r cf514c7ba5e6 -r 589e7f172c2b V_Boundary_Refine/Update_V.C
--- a/V_Boundary_Refine/Update_V.C Thu Jan 27 14:33:28 2011 -0800
+++ b/V_Boundary_Refine/Update_V.C Wed Feb 02 15:40:08 2011 -0800
@@ -20,6 +20,7 @@
#include "SAMRAI/tbox/Utilities.h"
#include "SAMRAI/pdat/CellData.h"
+#include "Boundary.h"
/* This is written from the perspective of axis==x. For axis==y, we
switch i and j and everything works out. */
@@ -48,8 +49,34 @@ void SAMRAI::geom::V_Boundary_Refine::Up
double v_minus=(*v)(center)
- (5.0/32)*dv_minus + (3.0/32)*dv_plus;
- (*v_fine)(fine)=v_minus*(2*(*v)(center))/(v_plus + v_minus);
- (*v_fine)(fine+jp)=v_plus*(2*(*v)(center))/(v_plus + v_minus);
+ const double epsilon(1.0e-8);
+ if(std::abs(v_plus+v_minus)<=epsilon*std::abs((*v)(center)))
+ {
+ (*v_fine)(fine)=(*v_fine)(fine+jp)=(*v)(center);
+ }
+ else
+ {
+ (*v_fine)(fine)=v_minus*(2*(*v)(center))/(v_plus + v_minus);
+ (*v_fine)(fine+jp)=v_plus*(2*(*v)(center))/(v_plus + v_minus);
+ }
+
+ tbox::plog << "Update V "
+ << axis << " "
+ << fine[0] << " "
+ << fine[1] << " "
+ << center[0] << " "
+ << center[1] << " "
+ << (*v_fine)(fine) << " "
+ << (*v_fine)(fine+jp) << " "
+ << (*v)(center) << " "
+ << (*v)(center+jp) << " "
+ << (*v)(center-jp) << " "
+ << (&(*v)(center-jp)) << " "
+ << v_plus << " "
+ << v_minus << " "
+ << dv_plus << " "
+ << dv_minus << " "
+ << "\n";
/* Set the point outside of the boundary to be max_double. This
give us a marker for whether the current value is a boundary
@@ -60,7 +87,7 @@ void SAMRAI::geom::V_Boundary_Refine::Up
offset[axis]=-1;
}
(*v_fine)(fine+offset)=(*v_fine)(fine+offset+jp)=
- std::numeric_limits<double>::max();
+ boundary_value;
++j;
}
else
diff -r cf514c7ba5e6 -r 589e7f172c2b V_Boundary_Refine/refine.C
--- a/V_Boundary_Refine/refine.C Thu Jan 27 14:33:28 2011 -0800
+++ b/V_Boundary_Refine/refine.C Wed Feb 02 15:40:08 2011 -0800
@@ -66,24 +66,6 @@ void SAMRAI::geom::V_Boundary_Refine::re
hier::Box coarse_box=coarse.getBox();
hier::Box fine_box=fine.getBox();
-
- std::cout << "VBR "
- << fine.getPatchLevelNumber() << " "
- << axis << " "
- << coarse_box.lower(0) << " "
- << coarse_box.upper(0) << " "
- << coarse_box.lower(1) << " "
- << coarse_box.upper(1) << " "
- << fine_box.lower(0) << " "
- << fine_box.upper(0) << " "
- << fine_box.lower(1) << " "
- << fine_box.upper(1) << " "
-
- << overlap_box.lower(0) << " "
- << overlap_box.upper(0) << " "
- << overlap_box.lower(1) << " "
- << overlap_box.upper(1) << " "
- << "\n";
/* We have to infer where the boundary is from the boxes */
int boundary_direction;
@@ -162,6 +144,31 @@ void SAMRAI::geom::V_Boundary_Refine::re
}
}
+ tbox::plog << "VBR "
+ << fine.getPatchLevelNumber() << " "
+ << axis << " "
+ << boundary_direction << " "
+ << std::boolalpha
+ << boundary_positive << " "
+ << coarse_box.lower(0) << " "
+ << coarse_box.upper(0) << " "
+ << coarse_box.lower(1) << " "
+ << coarse_box.upper(1) << " "
+ << fine_box.lower(0) << " "
+ << fine_box.upper(0) << " "
+ << fine_box.lower(1) << " "
+ << fine_box.upper(1) << " "
+
+ << overlap_box.lower(0) << " "
+ << overlap_box.upper(0) << " "
+ << overlap_box.lower(1) << " "
+ << overlap_box.upper(1) << " "
+ << i_min << " "
+ << i_max << " "
+ << j_min << " "
+ << j_max << " "
+ << "\n";
+
for(int j=j_min; j<=j_max; ++j)
for(int i=i_min; i<=i_max; ++i)
{
diff -r cf514c7ba5e6 -r 589e7f172c2b V_Coarsen.C
--- a/V_Coarsen.C Thu Jan 27 14:33:28 2011 -0800
+++ b/V_Coarsen.C Wed Feb 02 15:40:08 2011 -0800
@@ -21,6 +21,7 @@
#include "SAMRAI/pdat/SideData.h"
#include "SAMRAI/pdat/SideVariable.h"
#include "SAMRAI/tbox/Utilities.h"
+#include "Boundary.h"
void SAMRAI::geom::V_Coarsen::coarsen(hier::Patch& coarse,
@@ -91,19 +92,21 @@ void SAMRAI::geom::V_Coarsen::coarsen(hi
for(int j=coarse_box.lower(1); j<=coarse_box.upper(1)+1; ++j)
for(int i=coarse_box.lower(0); i<=coarse_box.upper(0)+1; ++i)
{
+ tbox::plog << "v coarsen "
+ << coarse.getPatchLevelNumber() << " "
+ << i << " "
+ << j << " ";
+
hier::Index ip(1,0), jp(0,1);
if(directions(0) && j!=coarse_box.upper(1)+1)
{
pdat::SideIndex coarse(hier::Index(i,j),0,
pdat::SideIndex::Lower);
pdat::SideIndex center(coarse*2);
- if(i==coarse_box.lower(0)
- && cgeom->getTouchesRegularBoundary(0,0))
- {
- (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+jp))/2;
- }
- else if(i==coarse_box.upper(0)+1
- && cgeom->getTouchesRegularBoundary(0,1))
+ if((i==coarse_box.lower(0)
+ && (*v_fine)(center-ip)==boundary_value)
+ || (i==coarse_box.upper(0)+1
+ && (*v_fine)(center+ip)==boundary_value))
{
(*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+jp))/2;
}
@@ -113,19 +116,21 @@ void SAMRAI::geom::V_Coarsen::coarsen(hi
+ ((*v_fine)(center-ip) + (*v_fine)(center-ip+jp)
+ (*v_fine)(center+ip) + (*v_fine)(center+ip+jp))/8;
}
+ tbox::plog << "vx "
+ << (*v)(coarse) << " "
+ << coarse_box.lower(0) << " "
+ << (*v_fine)(center-ip) << " "
+ << &((*v_fine)(center-ip)) << " ";
}
if(directions(1) && i!=coarse_box.upper(0)+1)
{
pdat::SideIndex coarse(hier::Index(i,j),1,
pdat::SideIndex::Lower);
pdat::SideIndex center(coarse*2);
- if(j==coarse_box.lower(1)
- && cgeom->getTouchesRegularBoundary(1,0))
- {
- (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+ip))/2;
- }
- else if(j==coarse_box.upper(1)+1
- && cgeom->getTouchesRegularBoundary(1,1))
+ if((j==coarse_box.lower(1)
+ && (*v_fine)(center-jp)==boundary_value)
+ || (j==coarse_box.upper(1)+1
+ && (*v_fine)(center+jp)==boundary_value))
{
(*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+ip))/2;
}
@@ -135,7 +140,10 @@ void SAMRAI::geom::V_Coarsen::coarsen(hi
+ ((*v_fine)(center-jp) + (*v_fine)(center-jp+ip)
+ (*v_fine)(center+jp) + (*v_fine)(center+ip+jp))/8;
}
+ tbox::plog << "vy "
+ << (*v)(coarse) << " ";
}
+ tbox::plog << "\n";
}
}
diff -r cf514c7ba5e6 -r 589e7f172c2b V_Coarsen_Patch_Strategy.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/V_Coarsen_Patch_Strategy.C Wed Feb 02 15:40:08 2011 -0800
@@ -0,0 +1,121 @@
+#include "V_Coarsen_Patch_Strategy.h"
+
+void
+SAMRAI::solv::V_Coarsen_Patch_Strategy::preprocessCoarsen
+(hier::Patch& coarse,
+ const hier::Patch& fine,
+ const hier::Box& coarse_box,
+ const hier::IntVector& ratio)
+{
+ tbox::Pointer<pdat::SideData<double> >
+ v = fine.getPatchData(v_id);
+
+ hier::Box pbox=fine.getBox();
+
+ hier::Box gbox=v->getGhostBox();
+
+ // tbox::plog << "boundary "
+ // << gbox.lower(0) << " "
+ // << gbox.upper(0) << " "
+ // << gbox.lower(1) << " "
+ // << gbox.upper(1) << " "
+ // << pbox.lower(0) << " "
+ // << pbox.upper(0) << " "
+ // << pbox.lower(1) << " "
+ // << pbox.upper(1) << " "
+ // << "\n";
+ for(int j=gbox.lower(1); j<=gbox.upper(1)+1; ++j)
+ for(int i=gbox.lower(0); i<=gbox.upper(0)+1; ++i)
+ {
+ pdat::CellIndex center(tbox::Dimension(2));
+ center[0]=i;
+ center[1]=j;
+ hier::Index ip(1,0), jp(0,1);
+
+ /* vx */
+ if(j<=gbox.upper(1))
+ {
+ /* Set a sentinel value */
+ if(i<pbox.lower(0) || i>pbox.upper(0)+1)
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=
+ boundary_value;
+ }
+ /* Set the value so the derivative=0 */
+ else if(j<pbox.lower(0))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
+ pdat::SideIndex::Lower));
+ tbox::plog << "V Coarsen Patch vx lower "
+ << fine.getPatchLevelNumber() << " "
+ << i << " "
+ << j << " "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ << " "
+ << (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ << " "
+ << gbox.lower(0) << " "
+ << pbox.lower(0) << " "
+ << "\n";
+ }
+ else if(j>pbox.upper(0))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
+ pdat::SideIndex::Lower));
+ }
+ tbox::plog << "V Coarsen Patch vx "
+ << fine.getPatchLevelNumber() << " "
+ << i << " "
+ << j << " "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ << " "
+ << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower)))
+ << " "
+ << "\n";
+ }
+ /* vy */
+ if(i<=gbox.upper(0))
+ {
+ if(j<pbox.lower(1) || j>pbox.upper(1)+1)
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=
+ boundary_value;
+ }
+ else if(i<pbox.lower(0))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center+ip,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower));
+ }
+ else if(i>pbox.upper(0))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center-ip,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower));
+ }
+ tbox::plog << "V Coarsen Patch vy "
+ << fine.getPatchLevelNumber() << " "
+ << i << " "
+ << j << " "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ << " "
+ << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower)))
+ << " "
+ << "\n";
+ }
+ }
+}
diff -r cf514c7ba5e6 -r 589e7f172c2b V_Coarsen_Patch_Strategy.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/V_Coarsen_Patch_Strategy.h Wed Feb 02 15:40:08 2011 -0800
@@ -0,0 +1,394 @@
+/*************************************************************************
+ *
+ * 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: Robin boundary condition support on cartesian grids.
+ *
+ ************************************************************************/
+#ifndef included_solv_V_Coarsen_Patch_Strategy
+#define included_solv_V_Coarsen_Patch_Strategy
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#include "SAMRAI/xfer/CoarsenPatchStrategy.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/pdat/CellIndex.h"
+#include "Boundary.h"
+
+namespace SAMRAI {
+namespace solv {
+
+/*!
+ * @brief Helper utility for setting boundary conditions on V.
+ *
+ * This class inherits and implements virtual functions from
+ * xfer::CoarsenPatchStrategy so it may be used to help create
+ * communication schedules if desired.
+ */
+class V_Coarsen_Patch_Strategy:
+ public xfer::CoarsenPatchStrategy
+{
+
+public:
+ /*!
+ * @brief Constructor using.
+ *
+ * @param object_name Name of the object, for general referencing.
+ * @param coef_strategy Coefficients strategy being helped.
+ */
+ V_Coarsen_Patch_Strategy(
+ const tbox::Dimension& dim,
+ std::string object_name = std::string()):
+ xfer::CoarsenPatchStrategy(dim), d_dim(dim),
+ d_object_name(object_name) {}
+
+ /*!
+ * @brief Destructor.
+ */
+ virtual ~V_Coarsen_Patch_Strategy(void) {}
+
+ //@{ @name xfer::CoarsenPatchStrategy virtuals
+
+ virtual hier::IntVector
+ getCoarsenOpStencilWidth() const
+ { return hier::IntVector::getOne(d_dim); }
+ // virtual void
+ // preprocessCoarsenBoxes(
+ // hier::Patch& fine,
+ // const hier::Patch& coarse,
+ // const hier::BoxList& fine_boxes,
+ // const hier::IntVector& ratio) {}
+ virtual void
+ preprocessCoarsen(hier::Patch& coarse,
+ const hier::Patch& fine,
+ const hier::Box& coarse_box,
+ const hier::IntVector& ratio);
+
+ virtual void
+ postprocessCoarsen(
+ hier::Patch& coarse,
+ const hier::Patch& fine,
+ const hier::Box& coarse_box,
+ const hier::IntVector& ratio) {}
+
+ //@}
+
+ //@{
+
+ /*!
+ * @name Functions to set boundary condition values
+ */
+
+ /*!
+ * @brief Set the physical boundary condition by setting the
+ * value of the first ghost cells.
+ *
+ * This function has an interface similar to the virtual function
+ * xfer::CoarsenPatchStrategy::setPhysicalBoundaryConditions(),
+ * and it may be used to help implement that function,
+ * but it does not serve the same purpose. The primary
+ * differences are:
+ * -# It specializes to cell-centered variables.
+ * -# Only one ghost cell width is filled. Setting a Robin
+ * boundary condition for cell-centered quantities requires
+ * only one ghost cell to be set.
+ * (More ghost cells can be filled by continuing the linear
+ * distribution of data beyond the first cell, but that is
+ * not implemented at this time.)
+ * -# User must specify the index of the data whose ghost
+ * cells need to be filled. This index is used to determine
+ * the variable for which to set the boundary coefficients
+ * and to get the data to be set.
+ *
+ * This function calls RobinBcStrategy::setBcCoefs() to
+ * get the coefficients, then it sets the values in the first
+ * ghost cell on the boundary.
+ *
+ * To determine the value for the ghost cell,
+ * a @em linear approximation in the direction normal to the
+ * boundary is assumed. We write the following discrete
+ * approximations:
+ * @f[ u_b = \frac{ u_i + u_o }{2} @f]
+ * @f[ [u_n]_b = \frac{ u_o - u_i }{h} @f]
+ * where the subscript b stands for the the boundary,
+ * i stands for the first cell inside the boundary and
+ * o stands for the first cell outside the boundary
+ * and h is the grid spacing normal to the boundary.
+ * Applying this to the Robin formula gives
+ * @f[ u_o = \frac{ h\gamma + u_i( \beta - \frac{h}{2} \alpha ) }
+ * { \beta + \frac{h}{2} \alpha } @f] or equivalently
+ * @f[ u_o = \frac{ hg + u_i (1-a(1+\frac{h}{2})) }{ 1-a(1-\frac{h}{2}) } @f]
+ *
+ * After setting the edge (face in 3D) boundary conditions,
+ * linear approximations are used to set the boundary conditions
+ * of higher boundary types (nodes in 2D, edges and nodes in 3D).
+ *
+ * In some cases, the calling function wants to set the
+ * boundary condition homogeneously, with g=0.
+ * This is useful in problems where the the solution of the
+ * homogeneous problem is required in solving the inhomogeneous
+ * problem. This function respects such requests specified
+ * through the argument @c homogeneous_bc.
+ *
+ * @internal To be more general to other data types,
+ * there could be versions for other data types also,
+ * such as ...InNodes, ...InFaces, etc. However, I'm not
+ * sure exactly how to implement the Robin boundary condition
+ * on the faces and nodes when m != 1. Should the boundary
+ * value be set or should the first ghost value be set?
+ *
+ * @internal I have not addressed possibility of differences
+ * in chosing the discrete formulation with which to compute
+ * the boundary value. The above formulation is obviously
+ * one specific approximation, but there could be others.
+ * If anoter approximation is required, there should be
+ * another class like this or this class can offer a choice
+ * to be set by the user. I favor another class.
+ *
+ * @internal Since the data alignment can be found through
+ * the target_data_id, these types of functions may be changed
+ * to just plain setBoundaryValues or setBoundaryValuesInBoundaryBoxes
+ * since it does assume boundary boxes. This may have to be
+ * expanded to later include coarse-fine boundary boxes
+ * for more generality.
+ *
+ * @param patch hier::Patch on which to set boundary condition
+ * @param fill_time Solution time corresponding to filling
+ * @param ghost_width_to_fill Max ghost width requiring fill
+ * @param target_data_id hier::Patch data index of data to be set.
+ * This data must be a cell-centered double.
+ * @param homogeneous_bc Set a homogeneous boundary condition.
+ * This means g=0 for the boundary.
+ */
+ // void
+ // setBoundaryValuesInCells(
+ // hier::Patch& patch,
+ // const double fill_time,
+ // const hier::IntVector& ghost_width_to_fill,
+ // int target_data_id,
+ // bool homogeneous_bc = false) const;
+
+ /*!
+ * @brief Set ghost cells for an entire level.
+ *
+ * Loop through all patches on the given level and call
+ * setBoundaryValuesInCells(hier::Patch &patch,
+ * const double fill_time ,
+ * const hier::IntVector &ghost_width_to_fill ,
+ * int target_data_id ,
+ * bool homogeneous_bc=false )
+ * for each.
+ *
+ * @param level PatchLevel on which to set boundary condition
+ * @param fill_time Solution time corresponding to filling
+ * @param ghost_width_to_fill Max ghost width requiring fill
+ * @param target_data_id hier::Patch data index of data to be set.
+ * This data must be a cell-centered double.
+ * @param homogeneous_bc Set a homogeneous boundary condition.
+ * This means g=0 for the boundary.
+ */
+ // void
+ // setBoundaryValuesInCells(
+ // hier::PatchLevel& level,
+ // const double fill_time,
+ // const hier::IntVector& ghost_width_to_fill,
+ // int target_data_id,
+ // bool homogeneous_bc = false) const;
+
+ /*!
+ * @brief Set the physical boundary condition by setting the
+ * value of the boundary nodes.
+ *
+ * This function is not yet implemented!
+ *
+ * There are some decisions that must be made before
+ * the implementation can be written.
+ * -# Do we set the values on the boundary or one cell
+ * away from the boundary?
+ * -# What is the discrete formulation we should use
+ * to compute the value to be set?
+ *
+ * This function has an interface similar to the virtual function
+ * xfer::CoarsenPatchStrategy::setPhysicalBoundaryConditions(),
+ * and it may be used to help implement that function,
+ * but it does not serve the same purpose. The primary
+ * differences are:
+ * -# It specializes to node-centered variables.
+ * -# User must specify the index of the data whose ghost
+ * cells need to be filled. This index is used to determine
+ * the variable for which to set the boundary coefficients
+ * and to get the data to be set.
+ *
+ * This function calls RobinBcStrategy::setBcCoefs() to get the
+ * coefficients, then it sets the values at the boundary nodes.
+ *
+ * In some cases, the calling function wants to set the
+ * boundary condition homogeneously, with g=0.
+ * This is useful in problems where the the solution of the
+ * homogeneous problem is required to solving the inhomogeneous
+ * problem. This function respects such requests specified
+ * through the argument @c homogeneous_bc.
+ *
+ * @param patch hier::Patch on which to set boundary condition
+ * @param fill_time Solution time corresponding to filling
+ * @param target_data_id hier::Patch data index of data to be set.
+ * @param homogeneous_bc Set a homogeneous boundary condition.
+ * This means g=0 for the boundary.
+ */
+ // void
+ // setBoundaryValuesAtNodes(
+ // hier::Patch& patch,
+ // const double fill_time,
+ // int target_data_id,
+ // bool homogeneous_bc = false) const;
+
+ //@}
+
+ //@{
+ /*!
+ * @name Ways to provide the Robin bc coefficients
+ */
+
+ /*!
+ * @brief Provide an implementation of the RobinBcCoefStrategy
+ * for determining the boundary coefficients.
+ *
+ * Provide the implementation that can be used to set the
+ * Robin bc coefficients.
+ *
+ * @param coef_strategy tbox::Pointer to a concrete inmplementation of
+ * the coefficient strategy.
+ */
+ // void
+ // setCoefImplementation(
+ // const RobinBcCoefStrategy* coef_strategy);
+
+ /*!
+ * @brief Set the data id that should be filled when setting
+ * physical boundary conditions.
+ *
+ * When setPhysicalBoundaryConditions is called, the data
+ * specified will be set. This information is required because
+ * the it is not passed in through the argument list of
+ * setPhysicalBounaryConditions.
+ */
+ void setSourceDataId(int id)
+ {
+ v_id=id;
+ }
+
+ /*!
+ * @brief Set whether boundary filling should assume homogeneous
+ * conditions.
+ *
+ * In certain circumstances, only the value of a is needed, while
+ * the value of g is temporarily not required and taken to be zero.
+ * (An example is in setting the boundary condition for error
+ * value in an iterative method.) In such cases, use this function
+ * to set a flag that will cause a null pointer to be given to
+ * setBcCoefs() to indicate that fact.
+ */
+ // void
+ // setHomogeneousBc(
+ // bool homogeneous_bc);
+
+ //@}
+
+private:
+ /*!
+ * @brief Trim a boundary box so that it does not stick out
+ * past a given box.
+ *
+ * Certain boundary-related operations occur on patch data that
+ * do not or cannot extend past the edgr or corner of a patch.
+ * This function is used to trim down the parts of the boundary box
+ * that extend past those points so that a suitable index range
+ * is achieved.
+ *
+ * The boundary box trimmed must be of type 1 or 2.
+ *
+ * @param boundary_box Boundary box to be trimmed.
+ * @param limit_box hier::Box to not stick past
+ *
+ * @return New trimmed boundary box.
+ */
+ // hier::BoundaryBox
+ // trimBoundaryBox(
+ // const hier::BoundaryBox& boundary_box,
+ // const hier::Box& limit_box) const;
+
+ /*!
+ * @brief Return box describing the index space of boundary nodes
+ * defined by a boundary box.
+ *
+ * Define a box describing the indices of the nodes corresponding
+ * to the input boundary box. These nodes lie on the boundary
+ * itself.
+ *
+ * The input boundary_box must be of type 1
+ * (see hier::BoundaryBox::getBoundaryType()).
+ *
+ * @param boundary_box input boundary box
+ * @return a box to define the node indices corresponding to
+ * boundary_box
+ */
+ // hier::Box
+ // makeNodeBoundaryBox(
+ // const hier::BoundaryBox& boundary_box) const;
+
+ /*!
+ * @brief Return box describing the index space of faces
+ * defined by a boundary box.
+ *
+ * Define a box describing the indices of the codimension 1
+ * surface corresponding to the input boundary box.
+ *
+ * The input boundary_box must be of type 1
+ * (see hier::BoundaryBox::getBoundaryType()).
+ *
+ * This is a utility function for working with the
+ * indices coresponding to a boundary box but coincide
+ * with the patch boundary.
+ *
+ * @param boundary_box input boundary box
+ * @return a box to define the face indices corresponding to
+ * boundary_box
+ */
+ // hier::Box
+ // makeFaceBoundaryBox(
+ // const hier::BoundaryBox& boundary_box) const;
+
+ const tbox::Dimension d_dim;
+
+ std::string d_object_name;
+
+ /*!
+ * @brief Coefficient strategy giving a way to get to
+ * Robin bc coefficients.
+ */
+ // const RobinBcCoefStrategy* d_coef_strategy;
+
+ /*!
+ * @brief hier::Index of target patch data when filling ghosts.
+ */
+ int v_id;
+
+ /*!
+ * @brief Whether to assumg g=0 when filling ghosts.
+ */
+ // bool d_homogeneous_bc;
+
+ /*!
+ * @brief Timers for performance measurement.
+ */
+ // tbox::Pointer<tbox::Timer> t_set_boundary_values_in_cells;
+ // tbox::Pointer<tbox::Timer> t_use_set_bc_coefs;
+};
+
+}
+}
+
+#endif // included_solv_V_Coarsen_Patch_Strategy
diff -r cf514c7ba5e6 -r 589e7f172c2b V_Refine.C
--- a/V_Refine.C Thu Jan 27 14:33:28 2011 -0800
+++ b/V_Refine.C Wed Feb 02 15:40:08 2011 -0800
@@ -81,12 +81,22 @@ void SAMRAI::geom::V_Refine::refine(hier
center[0]/=2;
center[1]/=2;
+ /* This assumes that the levels are always properly nested,
+ so that we always have an extra grid space for
+ interpolation. So we only have to have a special case for
+ physical boundaries, where we do not have an extra grid
+ space. */
+
if(axis==0)
{
double dvx_dy;
if(i%2==0)
{
+ /* Maybe this has to be fixed when dvx/dy != 0 on the
+ outer boundary because the approximation to the
+ derivative is not accurate enough? */
+
if(center[1]==coarse_box.lower(1)
&& geom->getTouchesRegularBoundary(1,0))
{
diff -r cf514c7ba5e6 -r 589e7f172c2b V_Refine_Patch_Strategy.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/V_Refine_Patch_Strategy.C Wed Feb 02 15:40:08 2011 -0800
@@ -0,0 +1,107 @@
+#include "V_Refine_Patch_Strategy.h"
+
+void
+SAMRAI::solv::V_Refine_Patch_Strategy::preprocessRefine
+(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio)
+{
+ tbox::Pointer<pdat::SideData<double> >
+ v = coarse.getPatchData(v_id);
+
+ hier::Box pbox=coarse.getBox();
+ hier::Box gbox=v->getGhostBox();
+
+ // tbox::plog << "boundary "
+ // << gbox.lower(0) << " "
+ // << gbox.upper(0) << " "
+ // << gbox.lower(1) << " "
+ // << gbox.upper(1) << " "
+ // << pbox.lower(0) << " "
+ // << pbox.upper(0) << " "
+ // << pbox.lower(1) << " "
+ // << pbox.upper(1) << " "
+ // << "\n";
+ for(int j=gbox.lower(1); j<=gbox.upper(1)+1; ++j)
+ for(int i=gbox.lower(0); i<=gbox.upper(0)+1; ++i)
+ {
+ pdat::CellIndex center(tbox::Dimension(2));
+ center[0]=i;
+ center[1]=j;
+ hier::Index ip(1,0), jp(0,1);
+
+ /* vx */
+ if(j<=gbox.upper(1))
+ {
+ /* Set a sentinel value */
+ if(i<pbox.lower(0) || i>pbox.upper(0)+1)
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=
+ boundary_value;
+ }
+ /* Set the value so the derivative=0 */
+ else if(j<pbox.lower(0))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
+ pdat::SideIndex::Lower));
+ }
+ else if(j>pbox.upper(0))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
+ pdat::SideIndex::Lower));
+ }
+ tbox::plog << "V Refine Patch vx "
+ << coarse.getPatchLevelNumber() << " "
+ << i << " "
+ << j << " "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ << " "
+ << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower)))
+ << " "
+ << "\n";
+ }
+ /* vy */
+ if(i<=gbox.upper(0))
+ {
+ if(j<pbox.lower(1) || j>pbox.upper(1)+1)
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=
+ boundary_value;
+ }
+ else if(i<pbox.lower(0))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center+ip,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower));
+ }
+ else if(i>pbox.upper(0))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center-ip,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower));
+ }
+ tbox::plog << "V Refine Patch vy "
+ << coarse.getPatchLevelNumber() << " "
+ << i << " "
+ << j << " "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ << " "
+ << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower)))
+ << " "
+ << "\n";
+ }
+ }
+}
diff -r cf514c7ba5e6 -r 589e7f172c2b V_Refine_Patch_Strategy.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/V_Refine_Patch_Strategy.h Wed Feb 02 15:40:08 2011 -0800
@@ -0,0 +1,404 @@
+/*************************************************************************
+ *
+ * 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: Robin boundary condition support on cartesian grids.
+ *
+ ************************************************************************/
+#ifndef included_solv_V_Refine_Patch_Strategy
+#define included_solv_V_Refine_Patch_Strategy
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#include "SAMRAI/xfer/RefinePatchStrategy.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/pdat/CellIndex.h"
+#include "Boundary.h"
+
+namespace SAMRAI {
+namespace solv {
+
+/*!
+ * @brief Helper utility for setting boundary conditions on V.
+ *
+ * This class inherits and implements virtual functions from
+ * xfer::RefinePatchStrategy so it may be used to help create
+ * communication schedules if desired.
+ */
+class V_Refine_Patch_Strategy:
+ public xfer::RefinePatchStrategy
+{
+
+public:
+ /*!
+ * @brief Constructor using.
+ *
+ * @param object_name Name of the object, for general referencing.
+ * @param coef_strategy Coefficients strategy being helped.
+ */
+ V_Refine_Patch_Strategy(
+ const tbox::Dimension& dim,
+ std::string object_name = std::string()):
+ xfer::RefinePatchStrategy(dim), d_dim(dim),d_object_name(object_name) {}
+
+ /*!
+ * @brief Destructor.
+ */
+ virtual ~V_Refine_Patch_Strategy(void) {}
+
+ //@{ @name xfer::RefinePatchStrategy virtuals
+
+ virtual void
+ setPhysicalBoundaryConditions(
+ hier::Patch& patch,
+ const double fill_time,
+ const hier::IntVector& ghost_width_to_fill) {}
+ hier::IntVector
+ getRefineOpStencilWidth() const
+ { return hier::IntVector::getOne(d_dim); }
+ // virtual void
+ // preprocessRefineBoxes(
+ // hier::Patch& fine,
+ // const hier::Patch& coarse,
+ // const hier::BoxList& fine_boxes,
+ // const hier::IntVector& ratio) {}
+ virtual void
+ preprocessRefine(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio);
+
+ virtual void
+ postprocessRefineBoxes(
+ hier::Patch& fine,
+ const hier::Patch& coarse,
+ const hier::BoxList& fine_boxes,
+ const hier::IntVector& ratio) {}
+ virtual void
+ postprocessRefine(
+ hier::Patch& fine,
+ const hier::Patch& coarse,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio) {}
+
+ //@}
+
+ //@{
+
+ /*!
+ * @name Functions to set boundary condition values
+ */
+
+ /*!
+ * @brief Set the physical boundary condition by setting the
+ * value of the first ghost cells.
+ *
+ * This function has an interface similar to the virtual function
+ * xfer::RefinePatchStrategy::setPhysicalBoundaryConditions(),
+ * and it may be used to help implement that function,
+ * but it does not serve the same purpose. The primary
+ * differences are:
+ * -# It specializes to cell-centered variables.
+ * -# Only one ghost cell width is filled. Setting a Robin
+ * boundary condition for cell-centered quantities requires
+ * only one ghost cell to be set.
+ * (More ghost cells can be filled by continuing the linear
+ * distribution of data beyond the first cell, but that is
+ * not implemented at this time.)
+ * -# User must specify the index of the data whose ghost
+ * cells need to be filled. This index is used to determine
+ * the variable for which to set the boundary coefficients
+ * and to get the data to be set.
+ *
+ * This function calls RobinBcStrategy::setBcCoefs() to
+ * get the coefficients, then it sets the values in the first
+ * ghost cell on the boundary.
+ *
+ * To determine the value for the ghost cell,
+ * a @em linear approximation in the direction normal to the
+ * boundary is assumed. We write the following discrete
+ * approximations:
+ * @f[ u_b = \frac{ u_i + u_o }{2} @f]
+ * @f[ [u_n]_b = \frac{ u_o - u_i }{h} @f]
+ * where the subscript b stands for the the boundary,
+ * i stands for the first cell inside the boundary and
+ * o stands for the first cell outside the boundary
+ * and h is the grid spacing normal to the boundary.
+ * Applying this to the Robin formula gives
+ * @f[ u_o = \frac{ h\gamma + u_i( \beta - \frac{h}{2} \alpha ) }
+ * { \beta + \frac{h}{2} \alpha } @f] or equivalently
+ * @f[ u_o = \frac{ hg + u_i (1-a(1+\frac{h}{2})) }{ 1-a(1-\frac{h}{2}) } @f]
+ *
+ * After setting the edge (face in 3D) boundary conditions,
+ * linear approximations are used to set the boundary conditions
+ * of higher boundary types (nodes in 2D, edges and nodes in 3D).
+ *
+ * In some cases, the calling function wants to set the
+ * boundary condition homogeneously, with g=0.
+ * This is useful in problems where the the solution of the
+ * homogeneous problem is required in solving the inhomogeneous
+ * problem. This function respects such requests specified
+ * through the argument @c homogeneous_bc.
+ *
+ * @internal To be more general to other data types,
+ * there could be versions for other data types also,
+ * such as ...InNodes, ...InFaces, etc. However, I'm not
+ * sure exactly how to implement the Robin boundary condition
+ * on the faces and nodes when m != 1. Should the boundary
+ * value be set or should the first ghost value be set?
+ *
+ * @internal I have not addressed possibility of differences
+ * in chosing the discrete formulation with which to compute
+ * the boundary value. The above formulation is obviously
+ * one specific approximation, but there could be others.
+ * If anoter approximation is required, there should be
+ * another class like this or this class can offer a choice
+ * to be set by the user. I favor another class.
+ *
+ * @internal Since the data alignment can be found through
+ * the target_data_id, these types of functions may be changed
+ * to just plain setBoundaryValues or setBoundaryValuesInBoundaryBoxes
+ * since it does assume boundary boxes. This may have to be
+ * expanded to later include coarse-fine boundary boxes
+ * for more generality.
+ *
+ * @param patch hier::Patch on which to set boundary condition
+ * @param fill_time Solution time corresponding to filling
+ * @param ghost_width_to_fill Max ghost width requiring fill
+ * @param target_data_id hier::Patch data index of data to be set.
+ * This data must be a cell-centered double.
+ * @param homogeneous_bc Set a homogeneous boundary condition.
+ * This means g=0 for the boundary.
+ */
+ // void
+ // setBoundaryValuesInCells(
+ // hier::Patch& patch,
+ // const double fill_time,
+ // const hier::IntVector& ghost_width_to_fill,
+ // int target_data_id,
+ // bool homogeneous_bc = false) const;
+
+ /*!
+ * @brief Set ghost cells for an entire level.
+ *
+ * Loop through all patches on the given level and call
+ * setBoundaryValuesInCells(hier::Patch &patch,
+ * const double fill_time ,
+ * const hier::IntVector &ghost_width_to_fill ,
+ * int target_data_id ,
+ * bool homogeneous_bc=false )
+ * for each.
+ *
+ * @param level PatchLevel on which to set boundary condition
+ * @param fill_time Solution time corresponding to filling
+ * @param ghost_width_to_fill Max ghost width requiring fill
+ * @param target_data_id hier::Patch data index of data to be set.
+ * This data must be a cell-centered double.
+ * @param homogeneous_bc Set a homogeneous boundary condition.
+ * This means g=0 for the boundary.
+ */
+ // void
+ // setBoundaryValuesInCells(
+ // hier::PatchLevel& level,
+ // const double fill_time,
+ // const hier::IntVector& ghost_width_to_fill,
+ // int target_data_id,
+ // bool homogeneous_bc = false) const;
+
+ /*!
+ * @brief Set the physical boundary condition by setting the
+ * value of the boundary nodes.
+ *
+ * This function is not yet implemented!
+ *
+ * There are some decisions that must be made before
+ * the implementation can be written.
+ * -# Do we set the values on the boundary or one cell
+ * away from the boundary?
+ * -# What is the discrete formulation we should use
+ * to compute the value to be set?
+ *
+ * This function has an interface similar to the virtual function
+ * xfer::RefinePatchStrategy::setPhysicalBoundaryConditions(),
+ * and it may be used to help implement that function,
+ * but it does not serve the same purpose. The primary
+ * differences are:
+ * -# It specializes to node-centered variables.
+ * -# User must specify the index of the data whose ghost
+ * cells need to be filled. This index is used to determine
+ * the variable for which to set the boundary coefficients
+ * and to get the data to be set.
+ *
+ * This function calls RobinBcStrategy::setBcCoefs() to get the
+ * coefficients, then it sets the values at the boundary nodes.
+ *
+ * In some cases, the calling function wants to set the
+ * boundary condition homogeneously, with g=0.
+ * This is useful in problems where the the solution of the
+ * homogeneous problem is required to solving the inhomogeneous
+ * problem. This function respects such requests specified
+ * through the argument @c homogeneous_bc.
+ *
+ * @param patch hier::Patch on which to set boundary condition
+ * @param fill_time Solution time corresponding to filling
+ * @param target_data_id hier::Patch data index of data to be set.
+ * @param homogeneous_bc Set a homogeneous boundary condition.
+ * This means g=0 for the boundary.
+ */
+ // void
+ // setBoundaryValuesAtNodes(
+ // hier::Patch& patch,
+ // const double fill_time,
+ // int target_data_id,
+ // bool homogeneous_bc = false) const;
+
+ //@}
+
+ //@{
+ /*!
+ * @name Ways to provide the Robin bc coefficients
+ */
+
+ /*!
+ * @brief Provide an implementation of the RobinBcCoefStrategy
+ * for determining the boundary coefficients.
+ *
+ * Provide the implementation that can be used to set the
+ * Robin bc coefficients.
+ *
+ * @param coef_strategy tbox::Pointer to a concrete inmplementation of
+ * the coefficient strategy.
+ */
+ // void
+ // setCoefImplementation(
+ // const RobinBcCoefStrategy* coef_strategy);
+
+ /*!
+ * @brief Set the data id that should be filled when setting
+ * physical boundary conditions.
+ *
+ * When setPhysicalBoundaryConditions is called, the data
+ * specified will be set. This information is required because
+ * the it is not passed in through the argument list of
+ * setPhysicalBounaryConditions.
+ */
+ void setTargetDataId(int id)
+ {
+ v_id=id;
+ }
+
+ /*!
+ * @brief Set whether boundary filling should assume homogeneous
+ * conditions.
+ *
+ * In certain circumstances, only the value of a is needed, while
+ * the value of g is temporarily not required and taken to be zero.
+ * (An example is in setting the boundary condition for error
+ * value in an iterative method.) In such cases, use this function
+ * to set a flag that will cause a null pointer to be given to
+ * setBcCoefs() to indicate that fact.
+ */
+ // void
+ // setHomogeneousBc(
+ // bool homogeneous_bc);
+
+ //@}
+
+private:
+ /*!
+ * @brief Trim a boundary box so that it does not stick out
+ * past a given box.
+ *
+ * Certain boundary-related operations occur on patch data that
+ * do not or cannot extend past the edgr or corner of a patch.
+ * This function is used to trim down the parts of the boundary box
+ * that extend past those points so that a suitable index range
+ * is achieved.
+ *
+ * The boundary box trimmed must be of type 1 or 2.
+ *
+ * @param boundary_box Boundary box to be trimmed.
+ * @param limit_box hier::Box to not stick past
+ *
+ * @return New trimmed boundary box.
+ */
+ // hier::BoundaryBox
+ // trimBoundaryBox(
+ // const hier::BoundaryBox& boundary_box,
+ // const hier::Box& limit_box) const;
+
+ /*!
+ * @brief Return box describing the index space of boundary nodes
+ * defined by a boundary box.
+ *
+ * Define a box describing the indices of the nodes corresponding
+ * to the input boundary box. These nodes lie on the boundary
+ * itself.
+ *
+ * The input boundary_box must be of type 1
+ * (see hier::BoundaryBox::getBoundaryType()).
+ *
+ * @param boundary_box input boundary box
+ * @return a box to define the node indices corresponding to
+ * boundary_box
+ */
+ // hier::Box
+ // makeNodeBoundaryBox(
+ // const hier::BoundaryBox& boundary_box) const;
+
+ /*!
+ * @brief Return box describing the index space of faces
+ * defined by a boundary box.
+ *
+ * Define a box describing the indices of the codimension 1
+ * surface corresponding to the input boundary box.
+ *
+ * The input boundary_box must be of type 1
+ * (see hier::BoundaryBox::getBoundaryType()).
+ *
+ * This is a utility function for working with the
+ * indices coresponding to a boundary box but coincide
+ * with the patch boundary.
+ *
+ * @param boundary_box input boundary box
+ * @return a box to define the face indices corresponding to
+ * boundary_box
+ */
+ // hier::Box
+ // makeFaceBoundaryBox(
+ // const hier::BoundaryBox& boundary_box) const;
+
+ const tbox::Dimension d_dim;
+
+ std::string d_object_name;
+
+ /*!
+ * @brief Coefficient strategy giving a way to get to
+ * Robin bc coefficients.
+ */
+ // const RobinBcCoefStrategy* d_coef_strategy;
+
+ /*!
+ * @brief hier::Index of target patch data when filling ghosts.
+ */
+ int v_id;
+
+ /*!
+ * @brief Whether to assumg g=0 when filling ghosts.
+ */
+ // bool d_homogeneous_bc;
+
+ /*!
+ * @brief Timers for performance measurement.
+ */
+ // tbox::Pointer<tbox::Timer> t_set_boundary_values_in_cells;
+ // tbox::Pointer<tbox::Timer> t_use_set_bc_coefs;
+};
+
+}
+}
+
+#endif // included_solv_V_Refine_Patch_Strategy
diff -r cf514c7ba5e6 -r 589e7f172c2b example_inputs/const_refine.2d.input
--- a/example_inputs/const_refine.2d.input Thu Jan 27 14:33:28 2011 -0800
+++ b/example_inputs/const_refine.2d.input Wed Feb 02 15:40:08 2011 -0800
@@ -32,7 +32,7 @@ FACStokes {
// 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
+ max_cycles = 10 // Max number of FAC cycles to use
residual_tol = 1e-8 // Residual tolerance to solve for
num_pre_sweeps = 0 // Number of presmoothing sweeps to use
num_post_sweeps = 0 // Number of postsmoothing sweeps to use
@@ -95,10 +95,10 @@ StandardTagAndInitialize {
tagging_method = "REFINE_BOXES"
RefineBoxes {
level_0 = [(0,0),(3,3)]
- level_1 = [(0,0),(7,7)]
- level_2 = [(0,0),(15,15)]
- // level_1 = [(0,0),(1,1)]
- // level_2 = [(0,0),(1,1)]
+ // level_1 = [(0,0),(7,7)]
+ // level_2 = [(0,0),(15,15)]
+ level_1 = [(0,0),(1,1)]
+ level_2 = [(0,0),(1,1)]
// level_3 = [(0,0),(1,1)]
//etc.
}
More information about the CIG-COMMITS
mailing list