[cig-commits] commit: Use sentinel values for boundaries
Mercurial
hg at geodynamics.org
Fri Feb 25 14:15:58 PST 2011
changeset: 70:cf514c7ba5e6
user: Walter Landry <wlandry at caltech.edu>
date: Thu Jan 27 14:33:28 2011 -0800
files: FACStokes/solveStokes.C Makefile StokesFACOps.h StokesFACOps/Update_V.C StokesFACOps/set_boundaries.C StokesFACOps/smoothErrorByRedBlack.C StokesFACSolver.h V_Boundary_Refine.C V_Boundary_Refine.h V_Boundary_Refine/Update_V.C V_Boundary_Refine/refine.C example_inputs/const_refine.2d.input
description:
Use sentinel values for boundaries
diff -r 817017ebaf83 -r cf514c7ba5e6 FACStokes/solveStokes.C
--- a/FACStokes/solveStokes.C Wed Jan 26 12:02:22 2011 -0800
+++ b/FACStokes/solveStokes.C Thu Jan 27 14:33:28 2011 -0800
@@ -46,24 +46,14 @@ int SAMRAI::FACStokes::solveStokes()
for ( ; ip; ip++) {
tbox::Pointer<hier::Patch> patch = *ip;
tbox::Pointer<pdat::CellData<double> >
- pdata = patch->getPatchData(p_id);
- pdata->fill(0.0);
+ p = patch->getPatchData(p_id);
+ p->fill(0.0);
tbox::Pointer<pdat::SideData<double> >
- vdata = patch->getPatchData(v_id);
- vdata->fill(0.0);
- /* This implicitly sets boundary conditions for v as well */
+ v = patch->getPatchData(v_id);
+ v->fill(0.0);
}
+ d_stokes_fac_solver.set_boundaries(v_id,level);
}
-
- /*
- * Set the parameters for the Stokes equation.
- * See classes solv::StokesFACSolver or
- * solv::StokesSpecifications.
- * (D is the diffusion coefficient.
- * C is the source term which is not used in this example.)
- */
- d_stokes_fac_solver.setDConstant(1.0);
- d_stokes_fac_solver.setCConstant(0.0);
d_stokes_fac_solver.initializeSolverState
(p_id,p_rhs_id,v_id,v_rhs_id,d_hierarchy,0,
diff -r 817017ebaf83 -r cf514c7ba5e6 Makefile
--- a/Makefile Wed Jan 26 12:02:22 2011 -0800
+++ b/Makefile Thu Jan 27 14:33:28 2011 -0800
@@ -30,7 +30,9 @@ CXX_OBJS = main.o FACStokes/FACStok
FACStokes/resetHierarchyConfiguration.o \
FACStokes/setupPlotter.o \
FACStokes/solveStokes.o \
- P_Refine.o V_Refine.o V_Boundary_Refine.o V_Coarsen.o \
+ P_Refine.o V_Refine.o V_Coarsen.o \
+ V_Boundary_Refine/refine.o \
+ V_Boundary_Refine/Update_V.o \
StokesFACOps/StokesFACOps.o \
StokesFACOps/checkInputPatchDataIndices.o \
StokesFACOps/computeCompositeResidualOnLevel.o \
@@ -47,8 +49,10 @@ CXX_OBJS = main.o FACStokes/FACStok
StokesFACOps/restrictSolution.o \
StokesFACOps/smoothError.o \
StokesFACOps/smoothErrorByRedBlack.o \
+ StokesFACOps/set_boundaries.o \
StokesFACOps/solveCoarsestLevel.o \
StokesFACOps/solveCoarsestLevel_HYPRE.o \
+ StokesFACOps/Update_V.o \
StokesFACOps/xeqScheduleFluxCoarsen.o \
StokesFACOps/xeqScheduleGhostFill.o \
StokesFACOps/xeqScheduleGhostFillNoCoarse.o \
diff -r 817017ebaf83 -r cf514c7ba5e6 StokesFACOps.h
--- a/StokesFACOps.h Wed Jan 26 12:02:22 2011 -0800
+++ b/StokesFACOps.h Thu Jan 27 14:33:28 2011 -0800
@@ -483,6 +483,8 @@ public:
const SAMRAIVectorReal<double>& current_soln,
const SAMRAIVectorReal<double>& residual);
+ void set_boundaries(const int &v_id, tbox::Pointer<hier::PatchLevel> &level);
+
//@}
private:
@@ -510,6 +512,22 @@ private:
int ln,
int num_sweeps,
double residual_tolerance = -1.0);
+
+ void Update_V(const int &axis, const int j,
+ const hier::Box &pbox,
+ const pdat::CellIndex ¢er,
+ const pdat::CellIndex &left,
+ const pdat::CellIndex &right,
+ const pdat::CellIndex &down,
+ const pdat::CellIndex &up,
+ tbox::Pointer<pdat::CellData<double> > &p,
+ tbox::Pointer<pdat::SideData<double> > &v,
+ tbox::Pointer<pdat::SideData<double> > &v_rhs,
+ double &maxres,
+ const double &dx,
+ const double &dy,
+ const double &viscosity,
+ const double &theta_momentum);
/*!
* @brief Solve the coarsest level using HYPRE
diff -r 817017ebaf83 -r cf514c7ba5e6 StokesFACOps/Update_V.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/Update_V.C Thu Jan 27 14:33:28 2011 -0800
@@ -0,0 +1,186 @@
+/*************************************************************************
+ *
+ * This file is part of the SAMRAI distribution. For full copyright
+ * information, see COPYRIGHT and COPYING.LESSER.
+ *
+ * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
+ * Description: Operator class for cell-centered scalar Stokes using FAC
+ *
+ ************************************************************************/
+#include "StokesFACOps.h"
+
+#include IOMANIP_HEADER_FILE
+
+#include "SAMRAI/hier/BoundaryBoxUtils.h"
+#include "SAMRAI/geom/CartesianGridGeometry.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/hier/Index.h"
+#include "SAMRAI/hier/Variable.h"
+#include "SAMRAI/hier/VariableDatabase.h"
+#include "SAMRAI/pdat/CellDoubleConstantRefine.h"
+#include "SAMRAI/pdat/CellVariable.h"
+#include "SAMRAI/pdat/OutersideData.h"
+#include "SAMRAI/pdat/OutersideVariable.h"
+#include "SAMRAI/hier/PatchData.h"
+#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/solv/FACPreconditioner.h"
+#include "StokesHypreSolver.h"
+#include "SAMRAI/tbox/Array.h"
+#include "SAMRAI/tbox/MathUtilities.h"
+#include "SAMRAI/tbox/StartupShutdownManager.h"
+#include "SAMRAI/tbox/Timer.h"
+#include "SAMRAI/tbox/TimerManager.h"
+#include "SAMRAI/tbox/Utilities.h"
+#include "SAMRAI/tbox/MathUtilities.h"
+#include "SAMRAI/xfer/CoarsenAlgorithm.h"
+#include "SAMRAI/xfer/CoarsenOperator.h"
+#include "SAMRAI/xfer/CoarsenSchedule.h"
+#include "SAMRAI/xfer/RefineAlgorithm.h"
+#include "SAMRAI/xfer/RefineOperator.h"
+#include "SAMRAI/xfer/RefineSchedule.h"
+#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
+
+/*
+********************************************************************
+* Updates one component of the velocity during a red-black *
+* Gauss-Seidel iteration. *
+********************************************************************
+*/
+void SAMRAI::solv::StokesFACOps::Update_V
+(const int &axis,
+ const int j,
+ const hier::Box &pbox,
+ const pdat::CellIndex ¢er,
+ const pdat::CellIndex &left,
+ const pdat::CellIndex &right,
+ const pdat::CellIndex &down,
+ const pdat::CellIndex &up,
+ tbox::Pointer<pdat::CellData<double> > &p,
+ tbox::Pointer<pdat::SideData<double> > &v,
+ tbox::Pointer<pdat::SideData<double> > &v_rhs,
+ double &maxres,
+ const double &dx,
+ const double &dy,
+ const double &viscosity,
+ const double &theta_momentum)
+{
+ const int off_axis=(axis==0) ? 1 : 0;
+ /* Update vx */
+ if(j<pbox.upper(off_axis)+1)
+ {
+ /* If at the 'x' boundaries, leave vx as is */
+ if(!((center[axis]==pbox.lower(axis)
+ && (*v)(pdat::SideIndex(left,
+ axis,
+ pdat::SideIndex::Lower))
+ ==std::numeric_limits<double>::max())
+ || (center[axis]==pbox.upper(axis)+1
+ && (*v)(pdat::SideIndex
+ (right,
+ axis,
+ pdat::SideIndex::Lower))
+ ==std::numeric_limits<double>::max())))
+ {
+ double dp_dx, d2vx_dxx, d2vx_dyy, C_vx;
+ /* If y==0 */
+ hier::Index offset(0,0);
+ bool set_boundary(false);
+ if(center[off_axis]==pbox.lower(off_axis))
+ {
+ offset[off_axis]=-1;
+ set_boundary=true;
+ }
+ else if(center[off_axis]==pbox.upper(off_axis))
+ {
+ offset[off_axis]=1;
+ set_boundary=true;
+ }
+
+ double dv;
+ if(set_boundary)
+ {
+ dv=(*v)(pdat::SideIndex
+ (center,
+ axis,
+ pdat::SideIndex::Lower))
+ - (*v)(pdat::SideIndex
+ (center+offset,axis,
+ pdat::SideIndex::Lower));
+ }
+
+ d2vx_dyy=
+ ((*v)(pdat::SideIndex(up,axis,
+ pdat::SideIndex::Lower))
+ - 2*(*v)(pdat::SideIndex
+ (center,axis,
+ pdat::SideIndex::Lower))
+ + (*v)(pdat::SideIndex
+ (down,axis,
+ pdat::SideIndex::Lower)))
+ /(dy*dy);
+
+ C_vx=-2*viscosity*(1/(dx*dx) + 1/(dy*dy));
+
+ d2vx_dxx=((*v)(pdat::SideIndex
+ (left,axis,
+ pdat::SideIndex::Lower))
+ - 2*(*v)(pdat::SideIndex
+ (center,axis,
+ pdat::SideIndex::Lower))
+ + (*v)(pdat::SideIndex
+ (right,axis,
+ pdat::SideIndex::Lower)))
+ /(dx*dx);
+
+ dp_dx=((*p)(center)-(*p)(left))/dx;
+
+ double delta_Rx=
+ (*v_rhs)(pdat::SideIndex(center,
+ axis,
+ pdat::SideIndex::Lower))
+ - viscosity*(d2vx_dxx + d2vx_dyy) + dp_dx;
+
+ // tbox::plog << "Update "
+ // << axis << " "
+ // << off_axis << " "
+ // << j << " "
+ // << center[axis] << " "
+ // << pbox.lower(axis) << " "
+ // << pbox.upper(axis) << " "
+ // << pbox.lower(off_axis) << " "
+ // << pbox.upper(off_axis) << " "
+ // << (*v_rhs)(pdat::SideIndex(center,
+ // axis,
+ // pdat::SideIndex::Lower)) << " "
+ // << right[0] << " "
+ // << right[1] << " "
+ // << (*v)(pdat::SideIndex
+ // (right,
+ // axis,
+ // pdat::SideIndex::Lower)) << " "
+ // << delta_Rx << " "
+ // << (theta_momentum/C_vx) << " "
+ // << "\n";
+
+
+ /* No scaling here, though there should be. */
+ maxres=std::max(maxres,delta_Rx);
+
+ (*v)(pdat::SideIndex(center,axis,
+ pdat::SideIndex::Lower))+=
+ delta_Rx*theta_momentum/C_vx;
+
+ /* Set the boundary elements so that the
+ derivative is zero. */
+ if(set_boundary)
+ {
+ (*v)(pdat::SideIndex
+ (center+offset,axis,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center,axis,
+ pdat::SideIndex::Lower))
+ -dv;
+ }
+ }
+ }
+}
diff -r 817017ebaf83 -r cf514c7ba5e6 StokesFACOps/set_boundaries.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/StokesFACOps/set_boundaries.C Thu Jan 27 14:33:28 2011 -0800
@@ -0,0 +1,137 @@
+/*************************************************************************
+ *
+ * This file is part of the SAMRAI distribution. For full copyright
+ * information, see COPYRIGHT and COPYING.LESSER.
+ *
+ * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
+ * Description: Operator class for cell-centered scalar Stokes using FAC
+ *
+ ************************************************************************/
+#include "StokesFACOps.h"
+
+#include IOMANIP_HEADER_FILE
+
+#include "SAMRAI/hier/BoundaryBoxUtils.h"
+#include "SAMRAI/geom/CartesianGridGeometry.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/hier/Index.h"
+#include "SAMRAI/hier/Variable.h"
+#include "SAMRAI/hier/VariableDatabase.h"
+#include "SAMRAI/pdat/CellDoubleConstantRefine.h"
+#include "SAMRAI/pdat/CellVariable.h"
+#include "SAMRAI/pdat/OutersideData.h"
+#include "SAMRAI/pdat/OutersideVariable.h"
+#include "SAMRAI/hier/PatchData.h"
+#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/solv/FACPreconditioner.h"
+#include "StokesHypreSolver.h"
+#include "SAMRAI/tbox/Array.h"
+#include "SAMRAI/tbox/MathUtilities.h"
+#include "SAMRAI/tbox/StartupShutdownManager.h"
+#include "SAMRAI/tbox/Timer.h"
+#include "SAMRAI/tbox/TimerManager.h"
+#include "SAMRAI/tbox/Utilities.h"
+#include "SAMRAI/tbox/MathUtilities.h"
+#include "SAMRAI/xfer/CoarsenAlgorithm.h"
+#include "SAMRAI/xfer/CoarsenOperator.h"
+#include "SAMRAI/xfer/CoarsenSchedule.h"
+#include "SAMRAI/xfer/RefineAlgorithm.h"
+#include "SAMRAI/xfer/RefineOperator.h"
+#include "SAMRAI/xfer/RefineSchedule.h"
+#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
+
+/*
+********************************************************************
+* Workhorse function to smooth error using red-black *
+* Gauss-Seidel iterations. *
+********************************************************************
+*/
+
+void SAMRAI::solv::StokesFACOps::set_boundaries
+(const int &v_id, tbox::Pointer<hier::PatchLevel> &level)
+{
+ for (hier::PatchLevel::Iterator pi(*level); pi; pi++)
+ {
+ tbox::Pointer<hier::Patch> patch = *pi;
+
+ tbox::Pointer<pdat::SideData<double> >
+ v = patch->getPatchData(v_id);
+
+ hier::Box pbox=patch->getBox();
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = patch->getPatchGeometry();
+
+ 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))=
+ std::numeric_limits<double>::max();
+ }
+ /* 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));
+ }
+ }
+ /* 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))=
+ std::numeric_limits<double>::max();
+ }
+ 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));
+ }
+ }
+ }
+ }
+
+}
diff -r 817017ebaf83 -r cf514c7ba5e6 StokesFACOps/smoothErrorByRedBlack.C
--- a/StokesFACOps/smoothErrorByRedBlack.C Wed Jan 26 12:02:22 2011 -0800
+++ b/StokesFACOps/smoothErrorByRedBlack.C Thu Jan 27 14:33:28 2011 -0800
@@ -77,13 +77,14 @@ void SAMRAI::solv::StokesFACOps::smoothE
#endif
tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
+ set_boundaries(v_id,level);
+
if (ln > d_ln_min) {
/*
* Perform a one-time transfer of data from coarser level,
* to fill ghost boundaries that will not change through
* the smoothing loop.
*/
- std::cout << "smooth\n";
xeqScheduleGhostFill(p_id, v_id, ln);
}
@@ -129,10 +130,6 @@ void SAMRAI::solv::StokesFACOps::smoothE
double dx = *(geom->getDx());
double dy = *(geom->getDx());
- hier::Box pgbox=p->getGhostBox();
- hier::Box prgbox=p_rhs->getGhostBox();
- hier::Box vgbox=v->getGhostBox();
-
for(int j=pbox.lower(1); j<=pbox.upper(1)+1; ++j)
{
/* Do the red-black skip */
@@ -152,199 +149,45 @@ void SAMRAI::solv::StokesFACOps::smoothE
--left[0];
/* Update p */
- if((!(i==pbox.upper(0)+1
- && geom->getTouchesRegularBoundary(0,1))
- && !(j==pbox.upper(1)+1
- && geom->getTouchesRegularBoundary(1,1))))
- {
- 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 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;
- /* Update vx */
- if(j!=pbox.upper(1)+1)
- {
- /* If x==0 */
- if((center[0]==pbox.lower(0)
- && geom->getTouchesRegularBoundary(0,0))
- || (center[0]==pbox.upper(0)+1
- && geom->getTouchesRegularBoundary(0,1)))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=0;
- }
- 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,
- 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));
- }
- 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);
+ /* 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);
- 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,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex
- (right,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- /(dx*dx);
-
- dp_dx=((*p)(center)-(*p)(left))/dx;
-
- double delta_Rx=
- (*v_rhs)(pdat::SideIndex(center,
- pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- - viscosity*(d2vx_dxx + d2vx_dyy) + dp_dx;
-
- /* No scaling here, though there should be. */
- maxres=std::max(maxres,delta_Rx);
-
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))+=
- delta_Rx*theta_momentum/C_vx;
- }
- }
-
- /* Update vy */
- if(i!=pbox.upper(0)+1)
- {
- /* If y==0 */
- if((center[1]==pbox.lower(1)
- && geom->getTouchesRegularBoundary(1,0))
- || (center[1]==pbox.upper(1)+1
- && geom->getTouchesRegularBoundary(1,1)))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=0;
- }
- 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,
- 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);
-
- 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,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex
- (down,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- /(dy*dy);
-
- dp_dy=((*p)(center)-(*p)(down))/dy;
-
- double delta_Ry=
- (*v_rhs)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- - viscosity*(d2vy_dxx + d2vy_dyy) + dp_dy;
-
- /* No scaling here, though there should be. */
- maxres=std::max(maxres,delta_Ry);
-
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))+=
- delta_Ry*theta_momentum/C_vy;
- }
- }
+ // 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";
}
}
}
diff -r 817017ebaf83 -r cf514c7ba5e6 StokesFACSolver.h
--- a/StokesFACSolver.h Wed Jan 26 12:02:22 2011 -0800
+++ b/StokesFACSolver.h Thu Jan 27 14:33:28 2011 -0800
@@ -558,6 +558,12 @@ public:
double
getResidualNorm() const;
+ void set_boundaries(const int &v_id, tbox::Pointer<hier::PatchLevel> &level)
+ {
+ d_fac_ops.set_boundaries(v_id,level);
+ }
+
+
//@}
private:
diff -r 817017ebaf83 -r cf514c7ba5e6 V_Boundary_Refine.C
--- a/V_Boundary_Refine.C Wed Jan 26 12:02:22 2011 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,297 +0,0 @@
-/*************************************************************************
- *
- * This file is part of the SAMRAI distribution. For full copyright
- * information, see COPYRIGHT and COPYING.LESSER.
- *
- * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
- * Description: Linear refine operator for side-centered double data on
- * a Cartesian mesh.
- *
- ************************************************************************/
-
-#ifndef included_geom_V_Boundary_Refine_C
-#define included_geom_V_Boundary_Refine_C
-
-#include "V_Boundary_Refine.h"
-
-#include <float.h>
-#include <math.h>
-#include "SAMRAI/geom/CartesianPatchGeometry.h"
-#include "SAMRAI/hier/Index.h"
-#include "SAMRAI/pdat/SideData.h"
-#include "SAMRAI/pdat/SideVariable.h"
-#include "SAMRAI/tbox/Utilities.h"
-
-void SAMRAI::geom::V_Boundary_Refine::refine(
- hier::Patch& fine,
- const hier::Patch& coarse,
- const int dst_component,
- const int src_component,
- const hier::BoxOverlap& fine_overlap,
- const hier::IntVector& ratio) const
-{
- const pdat::SideOverlap* t_overlap =
- dynamic_cast<const pdat::SideOverlap *>(&fine_overlap);
-
- TBOX_ASSERT(t_overlap != NULL);
-
- for(int axis=0; axis<getDim().getValue(); ++axis)
- {
- const hier::BoxList& boxes = t_overlap->getDestinationBoxList(axis);
- for (hier::BoxList::Iterator b(boxes); b; b++)
- {
- refine(fine,coarse,dst_component,src_component,b(),ratio,axis);
- }
- }
-}
-
-void SAMRAI::geom::V_Boundary_Refine::refine(hier::Patch& fine,
- const hier::Patch& coarse,
- const int dst_component,
- const int src_component,
- const hier::Box& overlap_box,
- const hier::IntVector& ratio,
- const int &axis) const
-{
- const tbox::Dimension& dim(getDim());
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dim, fine, coarse, overlap_box, ratio);
-
- tbox::Pointer<pdat::SideData<double> >
- v = coarse.getPatchData(src_component);
- tbox::Pointer<pdat::SideData<double> >
- v_fine = fine.getPatchData(dst_component);
-#ifdef DEBUG_CHECK_ASSERTIONS
- TBOX_ASSERT(!v.isNull());
- TBOX_ASSERT(!v_fine.isNull());
- TBOX_ASSERT(v->getDepth() == v_fine->getDepth());
- TBOX_ASSERT(v->getDepth() == 1);
-#endif
-
- hier::Box coarse_box=coarse.getBox();
- hier::Box fine_box=fine.getBox();
- tbox::Pointer<geom::CartesianPatchGeometry>
- geom = coarse.getPatchGeometry();
-
- double dx = geom->getDx()[0];
- double dy = geom->getDx()[1];
-
- 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;
- bool boundary_positive(false);
- if(std::abs(overlap_box.lower(0)-overlap_box.upper(0))==(axis==0 ? 1 : 0))
- {
- boundary_direction=0;
- if(fine_box.upper(0)<overlap_box.lower(0))
- boundary_positive=true;
- else if(fine_box.lower(0)>overlap_box.upper(0))
- boundary_positive=false;
- else
- abort();
- }
- else if(std::abs(overlap_box.lower(1)-overlap_box.upper(1))==
- (axis==1 ? 1 : 0))
- {
- boundary_direction=1;
- if(fine_box.upper(1)<overlap_box.lower(1))
- boundary_positive=true;
- else if(fine_box.lower(1)>overlap_box.upper(1))
- boundary_positive=false;
- else
- abort();
- }
- else
- {
- abort();
- }
- // hier::Index lower_offset(0,0), upper_offset(0,0);
- // if(axis==boundary_direction)
- // {
- // if(boundary_sign==1)
- // upper_offset[axis]=-1;
- // else
- // lower_offset[axis]=1;
- // }
-
- int i_min(overlap_box.lower(0)), i_max(overlap_box.upper(0)),
- j_min(overlap_box.lower(1)), j_max(overlap_box.upper(1));
- if(axis==0)
- {
- if(boundary_direction==0)
- {
- if(boundary_positive)
- {
- i_max=i_min;
- }
- else
- {
- i_min=i_max;
- }
- j_min=std::max(j_min,fine_box.lower(1));
- j_max=std::min(j_max,fine_box.upper(1));
- }
- /* We need to shrink the box because we do not want the edges.
- Those are points that are either covered by other patches or
- are ghost points that we do not care about. */
- else
- {
- --i_max;
- ++i_min;
- }
- }
- else if(axis==1)
- {
- if(boundary_direction==1)
- {
- if(boundary_positive)
- {
- j_max=j_min;
- }
- else
- {
- j_min=j_max;
- }
- i_min=std::max(i_min,fine_box.lower(0));
- i_max=std::min(i_max,fine_box.upper(0));
- }
- else
- {
- --j_max;
- ++j_min;
- }
- }
-
- for(int j=j_min; j<=j_max; ++j)
- for(int i=i_min; i<=i_max; ++i)
- {
- pdat::SideIndex fine(hier::Index(i,j),axis,pdat::SideIndex::Lower);
-
- hier::Index ip(1,0), jp(0,1);
- pdat::SideIndex center(fine);
- center[0]/=2;
- center[1]/=2;
-
- /* Note that at boundaries that are not in the same direction
- as the axis, we store the derivative in the variable */
- if(axis==0)
- {
- if(boundary_direction==0)
- {
- /* Interpolate in the y direction */
- double dv_plus, dv_minus;
-
- if(center[1]==coarse_box.lower(1)
- && geom->getTouchesRegularBoundary(1,0))
- {
- dv_plus=(*v)(center+jp)-(*v)(center);
- dv_minus=(*v)(center-jp)*dy;
- }
- else if(center[1]==coarse_box.upper(1)
- && geom->getTouchesRegularBoundary(1,1))
- {
- dv_plus=(*v)(center+jp)*dy;
- dv_minus=(*v)(center)-(*v)(center-jp);
- }
- else
- {
- dv_plus=(*v)(center+jp)-(*v)(center);
- dv_minus=(*v)(center)-(*v)(center-jp);
- }
-
- double v_plus=(*v)(center)
- + (5.0/32)*dv_plus - (3.0/32)*dv_minus;
- 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_minus*(2*(*v)(center))/(v_plus + v_minus);
- ++j;
- }
- else
- {
- /* We are computing derivatives here */
- if(i%2==0)
- {
- (*v_fine)(fine)=((*v)(center+jp) - (*v)(center))/dy;
- }
- else
- {
- (*v_fine)(fine)=
- 0.5*((*v)(center+jp) - (*v)(center)
- + (*v)(center+jp+ip) - (*v)(center+ip))/dy;
- }
- }
- }
- else if(axis==1)
- {
- if(boundary_direction==1)
- {
- /* Interpolate in the x direction */
- double dv_plus, dv_minus;
-
- if(center[0]==coarse_box.lower(0)
- && geom->getTouchesRegularBoundary(0,0))
- {
- dv_plus=(*v)(center+ip)-(*v)(center);
- dv_minus=(*v)(center-ip)*dx;
- }
- else if(center[0]==coarse_box.upper(0)
- && geom->getTouchesRegularBoundary(0,1))
- {
- dv_plus=(*v)(center+ip)*dx;
- dv_minus=(*v)(center)-(*v)(center-ip);
- }
- else
- {
- dv_plus=(*v)(center+ip)-(*v)(center);
- dv_minus=(*v)(center)-(*v)(center-ip);
- }
-
- double v_plus=(*v)(center)
- + (5.0/32)*dv_plus - (3.0/32)*dv_minus;
- 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+ip)=v_minus*(2*(*v)(center))/(v_plus + v_minus);
- ++i;
- }
- else
- {
- /* We are computing derivatives here */
- if(j%2==0)
- {
- (*v_fine)(fine)=((*v)(center+ip) - (*v)(center))/dx;
- }
- else
- {
- (*v_fine)(fine)=
- 0.5*((*v)(center+ip) - (*v)(center)
- + (*v)(center+ip+jp) - (*v)(center+jp))/dx;
- }
- }
- }
- else
- {
- abort();
- }
- }
-}
-#endif
diff -r 817017ebaf83 -r cf514c7ba5e6 V_Boundary_Refine.h
--- a/V_Boundary_Refine.h Wed Jan 26 12:02:22 2011 -0800
+++ b/V_Boundary_Refine.h Thu Jan 27 14:33:28 2011 -0800
@@ -20,6 +20,7 @@
#include "SAMRAI/hier/Patch.h"
#include "SAMRAI/tbox/Pointer.h"
#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
#include <string>
@@ -132,6 +133,16 @@ private:
private:
std::string d_name_id;
+ void Update_V
+ (const int &axis,
+ const int &boundary_direction,
+ const bool &boundary_positive,
+ const pdat::SideIndex &fine,
+ const hier::Index &ip, const hier::Index &jp,
+ int &j,
+ tbox::Pointer<pdat::SideData<double> > &v,
+ tbox::Pointer<pdat::SideData<double> > &v_fine) const;
+
};
}
diff -r 817017ebaf83 -r cf514c7ba5e6 V_Boundary_Refine/Update_V.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/V_Boundary_Refine/Update_V.C Thu Jan 27 14:33:28 2011 -0800
@@ -0,0 +1,86 @@
+/*************************************************************************
+ *
+ * This file is part of the SAMRAI distribution. For full copyright
+ * information, see COPYRIGHT and COPYING.LESSER.
+ *
+ * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
+ * Description: Linear refine operator for side-centered double data on
+ * a Cartesian mesh.
+ *
+ ************************************************************************/
+
+#include "V_Boundary_Refine.h"
+
+#include <float.h>
+#include <math.h>
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/hier/Index.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/tbox/Utilities.h"
+#include "SAMRAI/pdat/CellData.h"
+
+
+/* This is written from the perspective of axis==x. For axis==y, we
+ switch i and j and everything works out. */
+void SAMRAI::geom::V_Boundary_Refine::Update_V
+(const int &axis,
+ const int &boundary_direction,
+ const bool &boundary_positive,
+ const pdat::SideIndex &fine,
+ const hier::Index &ip, const hier::Index &jp,
+ int &j,
+ tbox::Pointer<SAMRAI::pdat::SideData<double> > &v,
+ tbox::Pointer<SAMRAI::pdat::SideData<double> > &v_fine) const
+{
+ pdat::SideIndex center(fine);
+ center/=2;
+ const int off_axis= (axis==0) ? 1 : 0;
+ if(boundary_direction==axis)
+ {
+ /* Interpolate in the y direction */
+ const double dv_plus=(*v)(center+jp)-(*v)(center);
+ const double dv_minus=(*v)(center)-(*v)(center-jp);
+
+ /* Quadratic interpolation */
+ double v_plus=(*v)(center)
+ + (5.0/32)*dv_plus - (3.0/32)*dv_minus;
+ 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);
+
+ /* Set the point outside of the boundary to be max_double. This
+ give us a marker for whether the current value is a boundary
+ condition. */
+ hier::Index offset(ip);
+ if(!boundary_positive)
+ {
+ offset[axis]=-1;
+ }
+ (*v_fine)(fine+offset)=(*v_fine)(fine+offset+jp)=
+ std::numeric_limits<double>::max();
+ ++j;
+ }
+ else
+ {
+ double dv;
+ /* Compute derivatives and use that to set the new vx */
+ if(fine[axis]%2==0)
+ {
+ dv=((*v)(center+jp) - (*v)(center))/2;
+ }
+ else
+ {
+ dv=((*v)(center+jp) - (*v)(center)
+ + (*v)(center+jp+ip) - (*v)(center+ip))/4;
+ }
+ hier::Index offset(jp);
+ if(!boundary_positive)
+ {
+ offset[off_axis]=-1;
+ }
+ (*v_fine)(fine)=(*v_fine)(fine-offset)+dv;
+ }
+}
diff -r 817017ebaf83 -r cf514c7ba5e6 V_Boundary_Refine/refine.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/V_Boundary_Refine/refine.C Thu Jan 27 14:33:28 2011 -0800
@@ -0,0 +1,182 @@
+/*************************************************************************
+ *
+ * This file is part of the SAMRAI distribution. For full copyright
+ * information, see COPYRIGHT and COPYING.LESSER.
+ *
+ * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
+ * Description: Linear refine operator for side-centered double data on
+ * a Cartesian mesh.
+ *
+ ************************************************************************/
+
+#include "V_Boundary_Refine.h"
+
+#include <float.h>
+#include <math.h>
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/hier/Index.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/tbox/Utilities.h"
+
+void SAMRAI::geom::V_Boundary_Refine::refine(
+ hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::BoxOverlap& fine_overlap,
+ const hier::IntVector& ratio) const
+{
+ const pdat::SideOverlap* t_overlap =
+ dynamic_cast<const pdat::SideOverlap *>(&fine_overlap);
+
+ TBOX_ASSERT(t_overlap != NULL);
+
+ for(int axis=0; axis<getDim().getValue(); ++axis)
+ {
+ const hier::BoxList& boxes = t_overlap->getDestinationBoxList(axis);
+ for (hier::BoxList::Iterator b(boxes); b; b++)
+ {
+ refine(fine,coarse,dst_component,src_component,b(),ratio,axis);
+ }
+ }
+}
+
+void SAMRAI::geom::V_Boundary_Refine::refine(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::Box& overlap_box,
+ const hier::IntVector& ratio,
+ const int &axis) const
+{
+ const tbox::Dimension& dim(getDim());
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dim, fine, coarse, overlap_box, ratio);
+
+ tbox::Pointer<pdat::SideData<double> >
+ v = coarse.getPatchData(src_component);
+ tbox::Pointer<pdat::SideData<double> >
+ v_fine = fine.getPatchData(dst_component);
+#ifdef DEBUG_CHECK_ASSERTIONS
+ TBOX_ASSERT(!v.isNull());
+ TBOX_ASSERT(!v_fine.isNull());
+ TBOX_ASSERT(v->getDepth() == v_fine->getDepth());
+ TBOX_ASSERT(v->getDepth() == 1);
+#endif
+
+ 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;
+ bool boundary_positive(false);
+ if(std::abs(overlap_box.lower(0)-overlap_box.upper(0))==(axis==0 ? 1 : 0))
+ {
+ boundary_direction=0;
+ if(fine_box.upper(0)<overlap_box.lower(0))
+ boundary_positive=true;
+ else if(fine_box.lower(0)>overlap_box.upper(0))
+ boundary_positive=false;
+ else
+ abort();
+ }
+ else if(std::abs(overlap_box.lower(1)-overlap_box.upper(1))==
+ (axis==1 ? 1 : 0))
+ {
+ boundary_direction=1;
+ if(fine_box.upper(1)<overlap_box.lower(1))
+ boundary_positive=true;
+ else if(fine_box.lower(1)>overlap_box.upper(1))
+ boundary_positive=false;
+ else
+ abort();
+ }
+ else
+ {
+ abort();
+ }
+
+ int i_min(overlap_box.lower(0)), i_max(overlap_box.upper(0)),
+ j_min(overlap_box.lower(1)), j_max(overlap_box.upper(1));
+ if(axis==0)
+ {
+ if(boundary_direction==0)
+ {
+ if(boundary_positive)
+ {
+ i_max=i_min;
+ }
+ else
+ {
+ i_min=i_max;
+ }
+ j_min=std::max(j_min,fine_box.lower(1));
+ j_max=std::min(j_max,fine_box.upper(1));
+ }
+ /* We need to shrink the box because we do not want the edges.
+ Those are points that are either covered by other patches or
+ are ghost points that we do not care about. */
+ else
+ {
+ --i_max;
+ ++i_min;
+ }
+ }
+ else if(axis==1)
+ {
+ if(boundary_direction==1)
+ {
+ if(boundary_positive)
+ {
+ j_max=j_min;
+ }
+ else
+ {
+ j_min=j_max;
+ }
+ i_min=std::max(i_min,fine_box.lower(0));
+ i_max=std::min(i_max,fine_box.upper(0));
+ }
+ else
+ {
+ --j_max;
+ ++j_min;
+ }
+ }
+
+ for(int j=j_min; j<=j_max; ++j)
+ for(int i=i_min; i<=i_max; ++i)
+ {
+ pdat::SideIndex fine(hier::Index(i,j),axis,pdat::SideIndex::Lower);
+ hier::Index ip(1,0), jp(0,1);
+
+ if(axis==0)
+ {
+ Update_V(axis,boundary_direction,boundary_positive,fine,
+ ip,jp,j,v,v_fine);
+ }
+ else if(axis==1)
+ {
+ Update_V(axis,boundary_direction,boundary_positive,fine,
+ jp,ip,i,v,v_fine);
+ }
+ }
+}
diff -r 817017ebaf83 -r cf514c7ba5e6 example_inputs/const_refine.2d.input
--- a/example_inputs/const_refine.2d.input Wed Jan 26 12:02:22 2011 -0800
+++ b/example_inputs/const_refine.2d.input Thu Jan 27 14:33:28 2011 -0800
@@ -34,8 +34,10 @@ FACStokes {
enable_logging = TRUE // Bool flag to switch logging on/off
max_cycles = 20 // Max number of FAC cycles to use
residual_tol = 1e-8 // Residual tolerance to solve for
- num_pre_sweeps = 5 // Number of presmoothing sweeps to use
- num_post_sweeps = 5 // Number of postsmoothing sweeps to use
+ num_pre_sweeps = 0 // Number of presmoothing sweeps to use
+ num_post_sweeps = 0 // Number of postsmoothing sweeps to use
+ coarse_solver_max_iterations = 200
+ coarse_solver_tolerance = 1e-10
p_prolongation_method = "P_REFINE" // Type of refinement
// used in prolongation.
// Suggested values are
@@ -93,8 +95,10 @@ StandardTagAndInitialize {
tagging_method = "REFINE_BOXES"
RefineBoxes {
level_0 = [(0,0),(3,3)]
- 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.
}
@@ -132,8 +136,8 @@ PatchHierarchy {
level_2 = 2, 2
}
largest_patch_size {
- //level_0 = 32, 32
- level_0 = 8, 8
+ level_0 = 32, 32
+ // level_0 = 8, 8
// all finer levels will use same values as level_0...
}
}
More information about the CIG-COMMITS
mailing list