[cig-commits] commit: Use Tackley's smoothing algorithm. Seems to be slightly faster with adapted grids.
Mercurial
hg at geodynamics.org
Fri Feb 25 14:17:22 PST 2011
changeset: 95:872d7c48020c
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 23 16:20:27 2011 -0800
files: FACStokes.h FACStokes/FACStokes.C FACStokes/initializeLevelData.C FACStokes/setupPlotter.C FACStokes/solveStokes.C StokesFACOps.I StokesFACOps.h StokesFACOps/StokesFACOps.C StokesFACOps/checkInputPatchDataIndices.C StokesFACOps/initializeOperatorState.C StokesFACOps/smoothErrorByRedBlack.C StokesFACOps/xeqScheduleGhostFillNoCoarse.C StokesFACSolver.h StokesFACSolver/initializeSolverState.C StokesFACSolver/solveSystem.C
description:
Use Tackley's smoothing algorithm. Seems to be slightly faster with adapted grids.
diff -r e9484bf65f05 -r 872d7c48020c FACStokes.h
--- a/FACStokes.h Wed Feb 23 16:19:22 2011 -0800
+++ b/FACStokes.h Wed Feb 23 16:20:27 2011 -0800
@@ -168,7 +168,7 @@ namespace SAMRAI {
*
* These are initialized in the constructor and never change.
*/
- int p_id, p_exact_id, p_rhs_id, v_id, v_rhs_id;
+ int p_id, viscosity_id, dp_id, p_exact_id, p_rhs_id, v_id, v_rhs_id;
//@}
diff -r e9484bf65f05 -r 872d7c48020c FACStokes/FACStokes.C
--- a/FACStokes/FACStokes.C Wed Feb 23 16:19:22 2011 -0800
+++ b/FACStokes/FACStokes.C Wed Feb 23 16:20:27 2011 -0800
@@ -70,6 +70,21 @@ namespace SAMRAI {
stencil widths */);
tbox::Pointer<pdat::CellVariable<double> >
+ viscosity(new pdat::CellVariable<double>(dim,
+ object_name + ":viscosity"));
+ viscosity_id = vdb->registerVariableAndContext(viscosity,d_context,
+ hier::IntVector(dim, 1)
+ /* ghost cell width is
+ 1 in case needed */);
+
+ tbox::Pointer<pdat::CellVariable<double> >
+ dp(new pdat::CellVariable<double>(dim, object_name + ":dp"));
+ dp_id = vdb->registerVariableAndContext(dp,d_context,
+ hier::IntVector(dim, 1)
+ /* ghost cell width is
+ 1 in case needed */);
+
+ tbox::Pointer<pdat::CellVariable<double> >
p_exact(new pdat::CellVariable<double>(dim, object_name + ":p exact"));
p_exact_id = vdb->registerVariableAndContext(p_exact,d_context,
hier::IntVector(dim, 1)
diff -r e9484bf65f05 -r 872d7c48020c FACStokes/initializeLevelData.C
--- a/FACStokes/initializeLevelData.C Wed Feb 23 16:19:22 2011 -0800
+++ b/FACStokes/initializeLevelData.C Wed Feb 23 16:20:27 2011 -0800
@@ -66,6 +66,8 @@ namespace SAMRAI {
if (allocate_data) {
level->allocatePatchData(p_id);
+ level->allocatePatchData(viscosity_id);
+ level->allocatePatchData(dp_id);
level->allocatePatchData(p_rhs_id);
level->allocatePatchData(p_exact_id);
level->allocatePatchData(v_id);
@@ -83,6 +85,19 @@ namespace SAMRAI {
TBOX_ERROR(d_object_name
<< ": Cannot find patch. Null patch pointer.");
}
+ tbox::Pointer<pdat::CellData<double> > viscosity_data =
+ patch->getPatchData(viscosity_id);
+
+ /* At some point this needs to do the proper interpolation for
+ lower levels */
+ viscosity_data->fill(1.0);
+
+ tbox::Pointer<pdat::CellData<double> > dp_data =
+ patch->getPatchData(dp_id);
+
+ /* This is mostly so that the outer boundaries are set properly. */
+ dp_data->fill(0.0);
+
tbox::Pointer<pdat::CellData<double> > p_rhs_data =
patch->getPatchData(p_rhs_id);
diff -r e9484bf65f05 -r 872d7c48020c FACStokes/setupPlotter.C
--- a/FACStokes/setupPlotter.C Wed Feb 23 16:19:22 2011 -0800
+++ b/FACStokes/setupPlotter.C Wed Feb 23 16:20:27 2011 -0800
@@ -37,7 +37,7 @@ namespace SAMRAI {
<< "The hierarchy must be set before calling\n"
<< "this function.\n");
}
- plotter.registerPlotQuantity("Computed solution",
+ plotter.registerPlotQuantity("Pressure",
"SCALAR",
p_id);
plotter.registerDerivedPlotQuantity("Error",
@@ -46,6 +46,9 @@ namespace SAMRAI {
plotter.registerPlotQuantity("Exact solution",
"SCALAR",
p_exact_id);
+ plotter.registerPlotQuantity("Viscosity",
+ "SCALAR",
+ viscosity_id);
plotter.registerPlotQuantity("Stokes source",
"SCALAR",
p_rhs_id);
diff -r e9484bf65f05 -r 872d7c48020c FACStokes/solveStokes.C
--- a/FACStokes/solveStokes.C Wed Feb 23 16:19:22 2011 -0800
+++ b/FACStokes/solveStokes.C Wed Feb 23 16:20:27 2011 -0800
@@ -56,12 +56,13 @@ int SAMRAI::FACStokes::solveStokes()
}
d_stokes_fac_solver.initializeSolverState
- (p_id,p_rhs_id,v_id,v_rhs_id,d_hierarchy,0,
+ (p_id,viscosity_id,dp_id,p_rhs_id,v_id,v_rhs_id,d_hierarchy,0,
d_hierarchy->getFinestLevelNumber());
tbox::plog << "solving..." << std::endl;
int solver_ret;
- solver_ret = d_stokes_fac_solver.solveSystem(p_id,p_rhs_id,v_id,v_rhs_id);
+ solver_ret = d_stokes_fac_solver.solveSystem(p_id,viscosity_id,dp_id,
+ p_rhs_id,v_id,v_rhs_id);
/*
* Present data on the solve.
*/
diff -r e9484bf65f05 -r 872d7c48020c StokesFACOps.I
--- a/StokesFACOps.I Wed Feb 23 16:19:22 2011 -0800
+++ b/StokesFACOps.I Wed Feb 23 16:20:27 2011 -0800
@@ -72,21 +72,6 @@ void StokesFACOps::enableLogging(
bool enable_logging)
{
d_enable_logging = enable_logging;
-}
-
-/*
- ********************************************************************
- * Set the patch data id for the flux. *
- ********************************************************************
- */
-
-SAMRAI_INLINE_KEYWORD
-void StokesFACOps::setFluxId(
- int flux_id) {
- d_flux_id = flux_id;
-#ifdef DEBUG_CHECK_ASSERTIONS
- checkInputPatchDataIndices();
-#endif
}
/*
diff -r e9484bf65f05 -r 872d7c48020c StokesFACOps.h
--- a/StokesFACOps.h Wed Feb 23 16:19:22 2011 -0800
+++ b/StokesFACOps.h Wed Feb 23 16:20:27 2011 -0800
@@ -296,10 +296,11 @@ public:
* flux and you would like that to be used, set flux id to the
* patch data index of that space.
*/
- void
- setFluxId(
- int flux_id);
-
+ void set_viscosity_dp_id(const int &viscosity, const int &dp)
+ {
+ viscosity_id=viscosity;
+ dp_id=dp;
+ }
//@}
/*!
@@ -869,14 +870,11 @@ private:
double d_residual_tolerance_during_smoothing;
/*!
- * @brief Id of the flux.
+ * @brief Id of viscosity and dp.
*
- * If set to -1, create and delete storage space on the fly.
- * Else, user has provided space for flux.
- *
- * @see setFluxId
+ * @see set_viscosity_dp_id.
*/
- int d_flux_id;
+ int viscosity_id, dp_id;
#ifdef HAVE_HYPRE
/*!
diff -r e9484bf65f05 -r 872d7c48020c StokesFACOps/StokesFACOps.C
--- a/StokesFACOps/StokesFACOps.C Wed Feb 23 16:19:22 2011 -0800
+++ b/StokesFACOps/StokesFACOps.C Wed Feb 23 16:20:27 2011 -0800
@@ -97,7 +97,8 @@ namespace SAMRAI {
d_coarse_solver_tolerance(1.e-8),
d_coarse_solver_max_iterations(10),
d_residual_tolerance_during_smoothing(-1.0),
- d_flux_id(-1),
+ viscosity_id(-1),
+ dp_id(-1),
#ifdef HAVE_HYPRE
d_hypre_solver(dim,
object_name + "::hypre_solver",
diff -r e9484bf65f05 -r 872d7c48020c StokesFACOps/checkInputPatchDataIndices.C
--- a/StokesFACOps/checkInputPatchDataIndices.C Wed Feb 23 16:19:22 2011 -0800
+++ b/StokesFACOps/checkInputPatchDataIndices.C Wed Feb 23 16:20:27 2011 -0800
@@ -74,14 +74,6 @@ namespace SAMRAI {
TBOX_ERROR(d_object_name << ": Bad linear term patch data index.");
}
}
-
- if (d_flux_id != -1) {
- tbox::Pointer<hier::Variable> var;
- vdb.mapIndexToVariable(d_flux_id, var);
- tbox::Pointer<pdat::SideVariable<double> > flux_var = var;
-
- TBOX_ASSERT(flux_var);
- }
}
}
diff -r e9484bf65f05 -r 872d7c48020c StokesFACOps/initializeOperatorState.C
--- a/StokesFACOps/initializeOperatorState.C Wed Feb 23 16:19:22 2011 -0800
+++ b/StokesFACOps/initializeOperatorState.C Wed Feb 23 16:20:27 2011 -0800
@@ -246,11 +246,6 @@ void SAMRAI::solv::StokesFACOps::initial
geometry->lookupCoarsenOperator(variable,
"V_COARSEN");
- vdb->mapIndexToVariable(d_oflux_scratch_id, variable);
- d_flux_coarsen_operator =
- geometry->lookupCoarsenOperator(variable,
- "CONSERVATIVE_COARSEN");
-
vdb->mapIndexToVariable(d_cell_scratch_id, variable);
p_ghostfill_refine_operator =
geometry->lookupRefineOperator(variable,
@@ -296,10 +291,6 @@ void SAMRAI::solv::StokesFACOps::initial
TBOX_ERROR(d_object_name
<< ": Cannot find v restriction coarsening operator");
}
- if (!d_flux_coarsen_operator) {
- TBOX_ERROR(d_object_name
- << ": Cannot find flux coarsening operator");
- }
if (!p_ghostfill_refine_operator) {
TBOX_ERROR(d_object_name
<< ": Cannot find ghost filling refinement operator");
@@ -318,11 +309,6 @@ void SAMRAI::solv::StokesFACOps::initial
}
#endif
- for (ln = d_ln_min + 1; ln <= d_ln_max; ++ln) {
- d_hierarchy->getPatchLevel(ln)->
- allocatePatchData(d_oflux_scratch_id);
- }
-
/*
* Make space for saving communication schedules.
* There is no need to delete the old schedules first
@@ -338,7 +324,6 @@ void SAMRAI::solv::StokesFACOps::initial
p_rrestriction_coarsen_schedules.resizeArray(d_ln_max + 1);
v_urestriction_coarsen_schedules.resizeArray(d_ln_max + 1);
v_rrestriction_coarsen_schedules.resizeArray(d_ln_max + 1);
- d_flux_coarsen_schedules.resizeArray(d_ln_max + 1);
p_prolongation_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
v_prolongation_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
@@ -346,7 +331,6 @@ void SAMRAI::solv::StokesFACOps::initial
p_rrestriction_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
v_urestriction_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
v_rrestriction_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
- d_flux_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
p_ghostfill_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
v_ghostfill_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
p_nocoarse_refine_algorithm = new xfer::RefineAlgorithm(d_dim);
@@ -388,10 +372,6 @@ void SAMRAI::solv::StokesFACOps::initial
solution.getComponentDescriptorIndex(1),
solution.getComponentDescriptorIndex(1),
v_ghostfill_refine_operator);
- d_flux_coarsen_algorithm->
- registerCoarsen(((d_flux_id != -1) ? d_flux_id : d_flux_scratch_id),
- d_oflux_scratch_id,
- d_flux_coarsen_operator);
p_nocoarse_refine_algorithm->
registerRefine(solution.getComponentDescriptorIndex(0),
solution.getComponentDescriptorIndex(0),
@@ -511,14 +491,6 @@ void SAMRAI::solv::StokesFACOps::initial
TBOX_ERROR(d_object_name
<< ": Cannot create a coarsen schedule for R v restriction!\n");
}
- d_flux_coarsen_schedules[dest_ln] =
- d_flux_coarsen_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- d_hierarchy->getPatchLevel(dest_ln + 1));
- if (!d_flux_coarsen_schedules[dest_ln]) {
- TBOX_ERROR(d_object_name
- << ": Cannot create a coarsen schedule for flux transfer!\n");
- }
}
/* Ordinary ghost fill operator on the coarsest level */
diff -r e9484bf65f05 -r 872d7c48020c StokesFACOps/smoothErrorByRedBlack.C
--- a/StokesFACOps/smoothErrorByRedBlack.C Wed Feb 23 16:19:22 2011 -0800
+++ b/StokesFACOps/smoothErrorByRedBlack.C Wed Feb 23 16:20:27 2011 -0800
@@ -91,8 +91,8 @@ void SAMRAI::solv::StokesFACOps::smoothE
}
double viscosity=1;
- double theta_momentum=1.2;
- double theta_continuity=0.3;
+ double theta_momentum=0.7;
+ double theta_continuity=1.0;
/*
* Smooth the number of sweeps specified or until
@@ -108,7 +108,10 @@ void SAMRAI::solv::StokesFACOps::smoothE
bool converged = false;
for (int sweep=0; sweep < num_sweeps*(1<<(d_ln_max-ln)) && !converged; ++sweep)
{
+
maxres=0;
+
+ /* vx sweep */
for(int rb=0;rb<2;++rb)
{
// Need to sync
@@ -139,13 +142,6 @@ void SAMRAI::solv::StokesFACOps::smoothE
hier::Box gbox=p->getGhostBox();
std::vector<bool> set_p(gbox.size(),true);
- // tbox::plog << "set_p "
- // << gbox.lower(0) << " "
- // << gbox.upper(0) << " "
- // << gbox.lower(1) << " "
- // << gbox.upper(1) << " "
- // << "\n";
-
const tbox::Array<hier::BoundaryBox >&edges
=d_cf_boundary[ln]->getEdgeBoundaries(patch->getGlobalId());
for(int mm=0; mm<edges.size(); ++mm)
@@ -153,15 +149,8 @@ void SAMRAI::solv::StokesFACOps::smoothE
j<=edges[mm].getBox().upper(1); ++j)
for(int i=edges[mm].getBox().lower(0);
i<=edges[mm].getBox().upper(0); ++i)
- {
- set_p[(i-gbox.lower(0))
- + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
- // tbox::plog << i << " "
- // << j << " "
- // << (i-gbox.lower(0))
- // + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1)) << " "
- // << "\n";
- }
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
const tbox::Array<hier::BoundaryBox >&nodes
=d_cf_boundary[ln]->getNodeBoundaries(patch->getGlobalId());
@@ -170,54 +159,27 @@ void SAMRAI::solv::StokesFACOps::smoothE
j<=nodes[mm].getBox().upper(1); ++j)
for(int i=nodes[mm].getBox().lower(0);
i<=nodes[mm].getBox().upper(0); ++i)
- {
- set_p[(i-gbox.lower(0))
- + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
- // tbox::plog << i << " "
- // << j << " "
- // << (i-gbox.lower(0))
- // + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1)) << " "
- // << "\n";
- }
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
if(geom->getTouchesRegularBoundary(0,0))
for(int j=gbox.lower(1); j<=gbox.upper(1); ++j)
- {
set_p[(gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
- // tbox::plog << gbox.lower(0) << " "
- // << j << " "
- // << "\n";
- }
if(geom->getTouchesRegularBoundary(0,1))
for(int j=gbox.lower(1); j<=gbox.upper(1); ++j)
- {
set_p[(gbox.upper(0)-gbox.lower(0))
+ (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
- // tbox::plog << gbox.upper(0) << " "
- // << j << " "
- // << "\n";
- }
if(geom->getTouchesRegularBoundary(1,0))
for(int i=gbox.lower(0); i<=gbox.upper(0); ++i)
- {
set_p[i-gbox.lower(0)]=false;
- // tbox::plog << i << " "
- // << gbox.lower(1) << " "
- // << "\n";
- }
if(geom->getTouchesRegularBoundary(1,1))
for(int i=gbox.lower(0); i<=gbox.upper(0); ++i)
- {
set_p[(i-gbox.lower(0))
+ (gbox.upper(0)-gbox.lower(0)+1)*(gbox.upper(1)-gbox.lower(1))]=
false;
- // tbox::plog << i << " "
- // << gbox.upper(1) << " "
- // << "\n";
- }
for(int j=pbox.lower(1); j<=pbox.upper(1)+1; ++j)
{
@@ -237,16 +199,214 @@ void SAMRAI::solv::StokesFACOps::smoothE
++right[0];
--left[0];
- // tbox::plog << "smooth "
- // << ln << " "
- // << sweep << " "
- // << rb << " "
- // << i << " "
- // << j << " ";
- // // << pbox.lower(0) << " "
- // // << pbox.upper(0) << " "
- // // << pbox.lower(1) << " "
- // // << pbox.upper(1) << " ";
+ /* Update v */
+ if(set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]
+ || (i==pbox.upper(0)+1 && j<pbox.upper(1)+1))
+ {
+ Update_V(0,j,pbox,geom,center,left,right,down,up,p,
+ v,v_rhs,maxres,dx,dy,viscosity,theta_momentum);
+ }
+ }
+ }
+ }
+ set_boundaries(v_id,level,true);
+ }
+
+
+ /* vy sweep */
+ for(int rb=0;rb<2;++rb)
+ {
+ // Need to sync
+ xeqScheduleGhostFillNoCoarse(p_id,v_id,ln);
+ for (hier::PatchLevel::Iterator pi(*level); pi; pi++)
+ {
+ tbox::Pointer<hier::Patch> patch = *pi;
+
+ tbox::Pointer<pdat::CellData<double> >
+ p = patch->getPatchData(p_id);
+ tbox::Pointer<pdat::CellData<double> >
+ p_rhs = patch->getPatchData(p_rhs_id);
+ tbox::Pointer<pdat::SideData<double> >
+ v = patch->getPatchData(v_id);
+ tbox::Pointer<pdat::SideData<double> >
+ v_rhs = patch->getPatchData(v_rhs_id);
+
+ hier::Box pbox=patch->getBox();
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = patch->getPatchGeometry();
+ double dx = *(geom->getDx());
+ double dy = *(geom->getDx());
+
+ /* Set an array of bools that tells me whether a point
+ should set the pressure or just let it be. This is
+ needed at coarse/fine boundaries where the pressure
+ is fixed. */
+ hier::Box gbox=p->getGhostBox();
+ std::vector<bool> set_p(gbox.size(),true);
+
+ const tbox::Array<hier::BoundaryBox >&edges
+ =d_cf_boundary[ln]->getEdgeBoundaries(patch->getGlobalId());
+ for(int mm=0; mm<edges.size(); ++mm)
+ for(int j=edges[mm].getBox().lower(1);
+ j<=edges[mm].getBox().upper(1); ++j)
+ for(int i=edges[mm].getBox().lower(0);
+ i<=edges[mm].getBox().upper(0); ++i)
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ const tbox::Array<hier::BoundaryBox >&nodes
+ =d_cf_boundary[ln]->getNodeBoundaries(patch->getGlobalId());
+ for(int mm=0; mm<nodes.size(); ++mm)
+ for(int j=nodes[mm].getBox().lower(1);
+ j<=nodes[mm].getBox().upper(1); ++j)
+ for(int i=nodes[mm].getBox().lower(0);
+ i<=nodes[mm].getBox().upper(0); ++i)
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ if(geom->getTouchesRegularBoundary(0,0))
+ for(int j=gbox.lower(1); j<=gbox.upper(1); ++j)
+ set_p[(gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ if(geom->getTouchesRegularBoundary(0,1))
+ for(int j=gbox.lower(1); j<=gbox.upper(1); ++j)
+ set_p[(gbox.upper(0)-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ if(geom->getTouchesRegularBoundary(1,0))
+ for(int i=gbox.lower(0); i<=gbox.upper(0); ++i)
+ set_p[i-gbox.lower(0)]=false;
+
+ if(geom->getTouchesRegularBoundary(1,1))
+ for(int i=gbox.lower(0); i<=gbox.upper(0); ++i)
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(gbox.upper(1)-gbox.lower(1))]=
+ false;
+
+ for(int j=pbox.lower(1); j<=pbox.upper(1)+1; ++j)
+ {
+ /* Do the red-black skip */
+ int i_min=pbox.lower(0) + (abs(pbox.lower(0) + j + rb))%2;
+ for(int i=i_min; i<=pbox.upper(0)+1; i+=2)
+ {
+ pdat::CellIndex center(tbox::Dimension(2));
+ center[0]=i;
+ center[1]=j;
+
+ pdat::CellIndex up(center), down(center), right(center),
+ left(center);
+
+ ++up[1];
+ --down[1];
+ ++right[0];
+ --left[0];
+
+ /* Update v */
+ if(set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]
+ || (i<pbox.upper(0)+1 && j==pbox.upper(1)+1))
+ {
+ Update_V(1,i,pbox,geom,center,down,up,left,right,p,
+ v,v_rhs,maxres,dy,dx,viscosity,theta_momentum);
+ }
+ }
+ }
+ }
+ set_boundaries(v_id,level,true);
+ }
+
+
+
+ /* p sweep */
+ for(int rb=0;rb<2;++rb)
+ {
+ // Need to sync
+ xeqScheduleGhostFillNoCoarse(-1,v_id,ln);
+ for (hier::PatchLevel::Iterator pi(*level); pi; pi++)
+ {
+ tbox::Pointer<hier::Patch> patch = *pi;
+
+ tbox::Pointer<pdat::CellData<double> >
+ p = patch->getPatchData(p_id);
+ tbox::Pointer<pdat::CellData<double> >
+ dp = patch->getPatchData(dp_id);
+ tbox::Pointer<pdat::CellData<double> >
+ p_rhs = patch->getPatchData(p_rhs_id);
+ tbox::Pointer<pdat::SideData<double> >
+ v = patch->getPatchData(v_id);
+ tbox::Pointer<pdat::SideData<double> >
+ v_rhs = patch->getPatchData(v_rhs_id);
+
+ hier::Box pbox=patch->getBox();
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = patch->getPatchGeometry();
+ double dx = *(geom->getDx());
+ double dy = *(geom->getDx());
+
+ /* Set an array of bools that tells me whether a point
+ should set the pressure or just let it be. This is
+ needed at coarse/fine boundaries where the pressure
+ is fixed. */
+ hier::Box gbox=p->getGhostBox();
+ std::vector<bool> set_p(gbox.size(),true);
+
+ const tbox::Array<hier::BoundaryBox >&edges
+ =d_cf_boundary[ln]->getEdgeBoundaries(patch->getGlobalId());
+ for(int mm=0; mm<edges.size(); ++mm)
+ for(int j=edges[mm].getBox().lower(1);
+ j<=edges[mm].getBox().upper(1); ++j)
+ for(int i=edges[mm].getBox().lower(0);
+ i<=edges[mm].getBox().upper(0); ++i)
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ const tbox::Array<hier::BoundaryBox >&nodes
+ =d_cf_boundary[ln]->getNodeBoundaries(patch->getGlobalId());
+ for(int mm=0; mm<nodes.size(); ++mm)
+ for(int j=nodes[mm].getBox().lower(1);
+ j<=nodes[mm].getBox().upper(1); ++j)
+ for(int i=nodes[mm].getBox().lower(0);
+ i<=nodes[mm].getBox().upper(0); ++i)
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ if(geom->getTouchesRegularBoundary(0,0))
+ for(int j=gbox.lower(1); j<=gbox.upper(1); ++j)
+ set_p[(gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ if(geom->getTouchesRegularBoundary(0,1))
+ for(int j=gbox.lower(1); j<=gbox.upper(1); ++j)
+ set_p[(gbox.upper(0)-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ if(geom->getTouchesRegularBoundary(1,0))
+ for(int i=gbox.lower(0); i<=gbox.upper(0); ++i)
+ set_p[i-gbox.lower(0)]=false;
+
+ if(geom->getTouchesRegularBoundary(1,1))
+ for(int i=gbox.lower(0); i<=gbox.upper(0); ++i)
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(gbox.upper(1)-gbox.lower(1))]=
+ false;
+
+ for(int j=pbox.lower(1); j<=pbox.upper(1)+1; ++j)
+ {
+ /* Do the red-black skip */
+ int i_min=pbox.lower(0) + (abs(pbox.lower(0) + j + rb))%2;
+ for(int i=i_min; i<=pbox.upper(0)+1; i+=2)
+ {
+ pdat::CellIndex center(tbox::Dimension(2));
+ center[0]=i;
+ center[1]=j;
+
+ pdat::CellIndex up(center), down(center), right(center),
+ left(center);
+
+ ++up[1];
+ --down[1];
+ ++right[0];
+ --left[0];
/* Update p */
if(set_p[(i-gbox.lower(0))
@@ -269,48 +429,175 @@ void SAMRAI::solv::StokesFACOps::smoothE
/* No scaling here, though there should be. */
maxres=std::max(maxres,std::fabs(delta_R_continuity));
- (*p)(center)+=
+ (*dp)(center)=
viscosity*delta_R_continuity*theta_continuity;
- // tbox::plog << "p "
- // // << (*p)(center) << " "
- // // << maxres << " "
- // // << (*p_rhs)(center) << " "
- // // << dvx_dx << " "
- // // << dvy_dy << " "
- // << delta_R_continuity << " ";
- // // << (*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::Y,
- // // pdat::SideIndex::Upper)) << " "
- // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- // // pdat::SideIndex::Lower)) << " ";
- // // << dx << " "
- // // << dy << " ";
+
+ // if(ln==2 && i==15)
+ // tbox::plog << "smooth p "
+ // << i << " "
+ // << j << " "
+ // << (*dp)(center) << " "
+ // << delta_R_continuity << " "
+ // << dvx_dx << " "
+ // << dvy_dy << " "
+ // << (*p_rhs)(center) << " "
+ // << (*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::Y,
+ // pdat::SideIndex::Upper)) << " "
+ // << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower)) << " "
+
+ // << "\n";
+
+ (*p)(center)+=(*dp)(center);
}
+ }
+ }
+ }
+ }
+
+
+ /* fix v sweep */
+ for(int rb=0;rb<2;++rb)
+ {
+ // Need to sync
+ xeqScheduleGhostFillNoCoarse(dp_id,invalid_id,ln);
+ for (hier::PatchLevel::Iterator pi(*level); pi; pi++)
+ {
+ tbox::Pointer<hier::Patch> patch = *pi;
+
+ tbox::Pointer<pdat::CellData<double> >
+ p = patch->getPatchData(p_id);
+ tbox::Pointer<pdat::CellData<double> >
+ dp = patch->getPatchData(dp_id);
+ tbox::Pointer<pdat::CellData<double> >
+ p_rhs = patch->getPatchData(p_rhs_id);
+ tbox::Pointer<pdat::SideData<double> >
+ v = patch->getPatchData(v_id);
+ tbox::Pointer<pdat::SideData<double> >
+ v_rhs = patch->getPatchData(v_rhs_id);
+
+ hier::Box pbox=patch->getBox();
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = patch->getPatchGeometry();
+ double dx = *(geom->getDx());
+ double dy = *(geom->getDx());
+
+ /* Set an array of bools that tells me whether a point
+ should set the pressure or just let it be. This is
+ needed at coarse/fine boundaries where the pressure
+ is fixed. */
+ hier::Box gbox=p->getGhostBox();
+ std::vector<bool> set_p(gbox.size(),true);
+
+ const tbox::Array<hier::BoundaryBox >&edges
+ =d_cf_boundary[ln]->getEdgeBoundaries(patch->getGlobalId());
+ for(int mm=0; mm<edges.size(); ++mm)
+ for(int j=edges[mm].getBox().lower(1);
+ j<=edges[mm].getBox().upper(1); ++j)
+ for(int i=edges[mm].getBox().lower(0);
+ i<=edges[mm].getBox().upper(0); ++i)
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ const tbox::Array<hier::BoundaryBox >&nodes
+ =d_cf_boundary[ln]->getNodeBoundaries(patch->getGlobalId());
+ for(int mm=0; mm<nodes.size(); ++mm)
+ for(int j=nodes[mm].getBox().lower(1);
+ j<=nodes[mm].getBox().upper(1); ++j)
+ for(int i=nodes[mm].getBox().lower(0);
+ i<=nodes[mm].getBox().upper(0); ++i)
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ if(geom->getTouchesRegularBoundary(0,0))
+ for(int j=gbox.lower(1); j<=gbox.upper(1); ++j)
+ set_p[(gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ if(geom->getTouchesRegularBoundary(0,1))
+ for(int j=gbox.lower(1); j<=gbox.upper(1); ++j)
+ set_p[(gbox.upper(0)-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]=false;
+
+ if(geom->getTouchesRegularBoundary(1,0))
+ for(int i=gbox.lower(0); i<=gbox.upper(0); ++i)
+ set_p[i-gbox.lower(0)]=false;
+
+ if(geom->getTouchesRegularBoundary(1,1))
+ for(int i=gbox.lower(0); i<=gbox.upper(0); ++i)
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)-gbox.lower(0)+1)*(gbox.upper(1)-gbox.lower(1))]=
+ false;
+
+ for(int j=pbox.lower(1); j<=pbox.upper(1)+1; ++j)
+ {
+ /* Do the red-black skip */
+ int i_min=pbox.lower(0) + (abs(pbox.lower(0) + j + rb))%2;
+ for(int i=i_min; i<=pbox.upper(0)+1; i+=2)
+ {
+ pdat::CellIndex center(tbox::Dimension(2));
+ center[0]=i;
+ center[1]=j;
+
+ pdat::CellIndex up(center), down(center), right(center),
+ left(center);
+
+ ++up[1];
+ --down[1];
+ ++right[0];
+ --left[0];
+
/* Update v */
if(set_p[(i-gbox.lower(0))
+ (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]
|| (i==pbox.upper(0)+1 && j<pbox.upper(1)+1))
{
- Update_V(0,j,pbox,geom,center,left,right,down,up,p,
- v,v_rhs,maxres,dx,dy,viscosity,theta_momentum);
+ if(!((center[0]==pbox.lower(0)
+ && (*v)(pdat::SideIndex(left,
+ 0,
+ pdat::SideIndex::Lower))
+ ==boundary_value)
+ || (center[0]==pbox.upper(0)+1
+ && (*v)(pdat::SideIndex
+ (right,
+ 0,
+ pdat::SideIndex::Lower))
+ ==boundary_value)))
+ (*v)(pdat::SideIndex(center,0,pdat::SideIndex::Lower))
+ -=((*dp)(center) - (*dp)(left))
+ /(dx*2*viscosity*(1/(dx*dx) + 1/(dy*dy)));
}
if(set_p[(i-gbox.lower(0))
+ (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]
|| (i<pbox.upper(0)+1 && j==pbox.upper(1)+1))
{
- Update_V(1,i,pbox,geom,center,down,up,left,right,p,
- v,v_rhs,maxres,dy,dx,viscosity,theta_momentum);
+ if(!((center[1]==pbox.lower(1)
+ && (*v)(pdat::SideIndex(down,
+ 1,
+ pdat::SideIndex::Lower))
+ ==boundary_value)
+ || (center[1]==pbox.upper(1)+1
+ && (*v)(pdat::SideIndex
+ (up,
+ 1,
+ pdat::SideIndex::Lower))
+ ==boundary_value)))
+ (*v)(pdat::SideIndex(center,1,pdat::SideIndex::Lower))
+ -=((*dp)(center) - (*dp)(down))
+ /(dy*2*viscosity*(1/(dx*dx) + 1/(dy*dy)));
}
- // tbox::plog << "\n";
}
}
}
set_boundaries(v_id,level,true);
}
+
+
+
// if (residual_tolerance >= 0.0) {
/*
* Check for early end of sweeps due to convergence
diff -r e9484bf65f05 -r 872d7c48020c StokesFACOps/xeqScheduleGhostFillNoCoarse.C
--- a/StokesFACOps/xeqScheduleGhostFillNoCoarse.C Wed Feb 23 16:19:22 2011 -0800
+++ b/StokesFACOps/xeqScheduleGhostFillNoCoarse.C Wed Feb 23 16:20:27 2011 -0800
@@ -59,6 +59,7 @@ void SAMRAI::solv::StokesFACOps::xeqSche
}
/* v */
+ if(v_id!=invalid_id)
{
if (!v_nocoarse_refine_schedules[dest_ln]) {
TBOX_ERROR("Expected side schedule not found.");
diff -r e9484bf65f05 -r 872d7c48020c StokesFACSolver.h
--- a/StokesFACSolver.h Wed Feb 23 16:19:22 2011 -0800
+++ b/StokesFACSolver.h Wed Feb 23 16:20:27 2011 -0800
@@ -180,6 +180,8 @@ public:
bool
solveSystem(
const int p,
+ const int viscosity,
+ const int dp,
const int p_rhs,
const int v,
const int v_rhs,
@@ -208,7 +210,8 @@ public:
* @see solveSystem( const int, const int, tbox::Pointer< hier::PatchHierarchy >, int, int);
*/
bool
- solveSystem(const int p, const int p_rhs, const int v, const int v_rhs);
+ solveSystem(const int p, const int viscosity, const int dp, const int p_rhs,
+ const int v, const int v_rhs);
/*!
* @brief Specify the boundary conditions that are to be used at the
@@ -505,6 +508,8 @@ public:
*/
void
initializeSolverState(const int p,
+ const int viscosity,
+ const int dp,
const int p_rhs,
const int v,
const int v_rhs,
diff -r e9484bf65f05 -r 872d7c48020c StokesFACSolver/initializeSolverState.C
--- a/StokesFACSolver/initializeSolverState.C Wed Feb 23 16:19:22 2011 -0800
+++ b/StokesFACSolver/initializeSolverState.C Wed Feb 23 16:20:27 2011 -0800
@@ -30,6 +30,8 @@
void SAMRAI::solv::StokesFACSolver::initializeSolverState
(const int p,
+ const int viscosity,
+ const int dp,
const int p_rhs,
const int v,
const int v_rhs,
@@ -101,6 +103,8 @@ void SAMRAI::solv::StokesFACSolver::init
d_fac_ops.setStokesSpecifications(d_stokes_spec);
+ d_fac_ops.set_viscosity_dp_id(viscosity,dp);
+
createVectorWrappers(p, p_rhs, v, v_rhs);
d_fac_precond.initializeSolverState(*d_uv, *d_fv);
diff -r e9484bf65f05 -r 872d7c48020c StokesFACSolver/solveSystem.C
--- a/StokesFACSolver/solveSystem.C Wed Feb 23 16:19:22 2011 -0800
+++ b/StokesFACSolver/solveSystem.C Wed Feb 23 16:20:27 2011 -0800
@@ -28,7 +28,8 @@
*************************************************************************
*/
-bool SAMRAI::solv::StokesFACSolver::solveSystem(const int p, const int p_rhs,
+bool SAMRAI::solv::StokesFACSolver::solveSystem(const int p, const int viscosity,
+ const int dp, const int p_rhs,
const int v, const int v_rhs)
{
#ifdef DEBUG_CHECK_ASSERTIONS
@@ -61,9 +62,6 @@ bool SAMRAI::solv::StokesFACSolver::solv
createVectorWrappers(p, p_rhs, v, v_rhs);
bool solver_rval;
- // d_fac_ops.relax(*d_uv, *d_fv, 0, 5, 1.0, p, p_rhs, v, v_rhs);
- // d_fac_ops.relax(*d_uv, *d_fv, 1, 5, 1.0, p, p_rhs, v, v_rhs);
- // d_fac_ops.relax(*d_uv, *d_fv, 2, 100, 1.0);
solver_rval = d_fac_precond.solveSystem(*d_uv, *d_fv);
return solver_rval;
@@ -84,6 +82,8 @@ bool SAMRAI::solv::StokesFACSolver::solv
bool SAMRAI::solv::StokesFACSolver::solveSystem
(const int p,
+ const int viscosity,
+ const int dp,
const int p_rhs,
const int v,
const int v_rhs,
@@ -116,10 +116,11 @@ bool SAMRAI::solv::StokesFACSolver::solv
<< "specified.\n");
}
#endif
- initializeSolverState(p, p_rhs, v, v_rhs, hierarchy, coarse_ln, fine_ln);
+ initializeSolverState(p, viscosity, dp, p_rhs, v, v_rhs,
+ hierarchy, coarse_ln, fine_ln);
bool solver_rval;
- solver_rval = solveSystem(p, p_rhs, v, v_rhs);
+ solver_rval = solveSystem(p, viscosity, dp, p_rhs, v, v_rhs);
deallocateSolverState();
More information about the CIG-COMMITS
mailing list