[cig-commits] commit: Adaptive somewhat implemented. Does not track boundaries correctly, so it still crashes.
Mercurial
hg at geodynamics.org
Fri Feb 25 14:15:42 PST 2011
changeset: 65:5e19536406f3
user: Walter Landry <wlandry at caltech.edu>
date: Mon Jan 24 15:49:47 2011 -0800
files: Makefile StokesFACOps/initializeOperatorState.C StokesFACOps/smoothErrorByRedBlack.C V_Boundary_Refine.C V_Boundary_Refine.h example_inputs/const_refine.2d.input main.C
description:
Adaptive somewhat implemented. Does not track boundaries correctly, so it still crashes.
diff -r c120f9b7964d -r 5e19536406f3 Makefile
--- a/Makefile Fri Jan 14 15:44:10 2011 -0800
+++ b/Makefile Mon Jan 24 15:49:47 2011 -0800
@@ -30,7 +30,7 @@ CXX_OBJS = main.o FACStokes/FACStok
FACStokes/resetHierarchyConfiguration.o \
FACStokes/setupPlotter.o \
FACStokes/solveStokes.o \
- P_Refine.o V_Refine.o V_Coarsen.o \
+ P_Refine.o V_Refine.o V_Boundary_Refine.o V_Coarsen.o \
StokesFACOps/StokesFACOps.o \
StokesFACOps/checkInputPatchDataIndices.o \
StokesFACOps/computeCompositeResidualOnLevel.o \
diff -r c120f9b7964d -r 5e19536406f3 StokesFACOps/initializeOperatorState.C
--- a/StokesFACOps/initializeOperatorState.C Fri Jan 14 15:44:10 2011 -0800
+++ b/StokesFACOps/initializeOperatorState.C Mon Jan 24 15:49:47 2011 -0800
@@ -263,7 +263,7 @@ void SAMRAI::solv::StokesFACOps::initial
vdb->mapIndexToVariable(d_side_scratch_id, variable);
v_ghostfill_refine_operator =
geometry->lookupRefineOperator(variable,
- "CONSERVATIVE_LINEAR_REFINE");
+ "V_BOUNDARY_REFINE");
vdb->mapIndexToVariable(d_cell_scratch_id, variable);
p_nocoarse_refine_operator =
@@ -273,7 +273,7 @@ void SAMRAI::solv::StokesFACOps::initial
vdb->mapIndexToVariable(d_side_scratch_id, variable);
v_nocoarse_refine_operator =
geometry->lookupRefineOperator(variable,
- "CONSERVATIVE_LINEAR_REFINE");
+ "V_REFINE");
#ifdef DEBUG_CHECK_ASSERTIONS
if (!p_prolongation_refine_operator) {
diff -r c120f9b7964d -r 5e19536406f3 StokesFACOps/smoothErrorByRedBlack.C
--- a/StokesFACOps/smoothErrorByRedBlack.C Fri Jan 14 15:44:10 2011 -0800
+++ b/StokesFACOps/smoothErrorByRedBlack.C Mon Jan 24 15:49:47 2011 -0800
@@ -77,14 +77,15 @@ void SAMRAI::solv::StokesFACOps::smoothE
#endif
tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
- // 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.
- // */
- // xeqScheduleGhostFill(data_id, ln);
- // }
+ 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);
+ }
double viscosity=1;
double theta_momentum=1.2;
diff -r c120f9b7964d -r 5e19536406f3 V_Boundary_Refine.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/V_Boundary_Refine.C Mon Jan 24 15:49:47 2011 -0800
@@ -0,0 +1,281 @@
+/*************************************************************************
+ *
+ * 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();
+
+ 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));
+ }
+ }
+ 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));
+ }
+ }
+
+ 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 c120f9b7964d -r 5e19536406f3 V_Boundary_Refine.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/V_Boundary_Refine.h Mon Jan 24 15:49:47 2011 -0800
@@ -0,0 +1,139 @@
+/*************************************************************************
+ *
+ * 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
+#define included_geom_V_Boundary_Refine
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#include "SAMRAI/xfer/RefineOperator.h"
+#include "SAMRAI/hier/Box.h"
+#include "SAMRAI/hier/IntVector.h"
+#include "SAMRAI/hier/Patch.h"
+#include "SAMRAI/tbox/Pointer.h"
+#include "SAMRAI/pdat/SideVariable.h"
+
+#include <string>
+
+namespace SAMRAI {
+namespace geom {
+
+/**
+ * Class V_Boundary_Refine implements the special interpolation needed
+ * for the boundary elements on the coarse-fine interface.
+ *
+ * The findRefineOperator() operator function returns true if the input
+ * variable is side-centered double, and the std::string is "V_BOUNDARY_REFINE".
+ *
+ * @see xfer::RefineOperator
+ */
+
+class V_Boundary_Refine:
+ public xfer::RefineOperator
+{
+public:
+ /**
+ * Uninteresting default constructor.
+ */
+ explicit V_Boundary_Refine(const tbox::Dimension& dim):
+ xfer::RefineOperator(dim, "V_BOUNDARY_REFINE")
+ {
+ d_name_id = "V_BOUNDARY_REFINE";
+ }
+
+
+ /**
+ * Uninteresting virtual destructor.
+ */
+ virtual ~V_Boundary_Refine(){}
+
+ /**
+ * Return true if the variable and name std::string match side-centered
+ * double linear interpolation; otherwise, return false.
+ */
+ bool findRefineOperator(const tbox::Pointer<hier::Variable>& var,
+ const std::string& op_name) const
+ {
+ TBOX_DIM_ASSERT_CHECK_ARGS2(*this, *var);
+
+ const tbox::Pointer<pdat::SideVariable<double> > cast_var(var);
+ if (!cast_var.isNull() && (op_name == d_name_id)) {
+ return true;
+ } else {
+ return false;
+ }
+ }
+ /**
+ * Return name std::string identifier of this refinement operator.
+ */
+ const std::string& getOperatorName() const
+ {
+ return d_name_id;
+ }
+
+ /**
+ * The priority of side-centered double linear interpolation is 0.
+ * It will be performed before any user-defined interpolation operations.
+ */
+ int getOperatorPriority() const
+ {
+ return 0;
+ }
+
+ /**
+ * The stencil width of the linear interpolation operator is the vector
+ * of ones. That is, its stencil extends one side outside the fine box.
+ */
+ hier::IntVector getStencilWidth() const
+ {
+ return hier::IntVector::getOne(getDim());
+ }
+
+ /**
+ * Refine the source component on the coarse patch to the destination
+ * component on the fine patch using the side-centered double linear
+ * interpolation operator. Interpolation is performed on the intersection
+ * of the destination patch and the boxes contained in fine_overlap
+ * It is assumed that the coarse patch contains sufficient data for the
+ * stencil width of the refinement operator.
+ */
+ void refine(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::BoxOverlap& fine_overlap,
+ const hier::IntVector& ratio) const;
+
+ /**
+ * Refine the source component on the coarse patch to the destination
+ * component on the fine patch using the side-centered double linear
+ * interpolation operator. Interpolation is performed on the intersection
+ * of the destination patch and the fine box. It is assumed that the
+ * coarse patch contains sufficient data for the stencil width of the
+ * refinement operator. This differs from the above refine() method
+ * only in that it operates on a single fine box instead of a BoxOverlap.
+ */
+ void refine(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio,
+ const int &axis) const;
+
+private:
+ std::string d_name_id;
+
+};
+
+}
+}
+#endif
diff -r c120f9b7964d -r 5e19536406f3 example_inputs/const_refine.2d.input
--- a/example_inputs/const_refine.2d.input Fri Jan 14 15:44:10 2011 -0800
+++ b/example_inputs/const_refine.2d.input Mon Jan 24 15:49:47 2011 -0800
@@ -93,9 +93,9 @@ 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_3 = [(0,0),(31,31)]
+ level_1 = [(0,0),(1,1)]
+ level_2 = [(0,0),(1,1)]
+ // level_3 = [(0,0),(1,1)]
//etc.
}
}
@@ -126,7 +126,7 @@ PatchHierarchy {
// [level 0 entry]
// etc....
// }
- max_levels = 4
+ max_levels = 3
ratio_to_coarser {
level_1 = 2, 2
level_2 = 2, 2
diff -r c120f9b7964d -r 5e19536406f3 main.C
--- a/main.C Fri Jan 14 15:44:10 2011 -0800
+++ b/main.C Mon Jan 24 15:49:47 2011 -0800
@@ -31,6 +31,7 @@ using namespace std;
#include "SAMRAI/xfer/RefineOperator.h"
#include "P_Refine.h"
#include "V_Refine.h"
+#include "V_Boundary_Refine.h"
#include "V_Coarsen.h"
#include "FACStokes.h"
@@ -152,11 +153,18 @@ int main(
input_db->getDatabase("CartesianGridGeometry")));
tbox::plog << "Cartesian Geometry:" << endl;
grid_geometry->printClassData(tbox::plog);
- grid_geometry->addSpatialRefineOperator(tbox::Pointer<SAMRAI::xfer::RefineOperator>(new SAMRAI::geom::P_Refine(dim)));
- grid_geometry->addSpatialRefineOperator(tbox::Pointer<SAMRAI::xfer::RefineOperator>(new SAMRAI::geom::V_Refine(dim)));
- grid_geometry->addSpatialCoarsenOperator(tbox::Pointer<SAMRAI::xfer::CoarsenOperator>(new SAMRAI::geom::V_Coarsen(dim)));
-
-
+ grid_geometry->addSpatialRefineOperator
+ (tbox::Pointer<SAMRAI::xfer::RefineOperator>
+ (new SAMRAI::geom::P_Refine(dim)));
+ grid_geometry->addSpatialRefineOperator
+ (tbox::Pointer<SAMRAI::xfer::RefineOperator>
+ (new SAMRAI::geom::V_Refine(dim)));
+ grid_geometry->addSpatialRefineOperator
+ (tbox::Pointer<SAMRAI::xfer::RefineOperator>
+ (new SAMRAI::geom::V_Boundary_Refine(dim)));
+ grid_geometry->addSpatialCoarsenOperator
+ (tbox::Pointer<SAMRAI::xfer::CoarsenOperator>
+ (new SAMRAI::geom::V_Coarsen(dim)));
tbox::Pointer<hier::PatchHierarchy>
patch_hierarchy(new hier::PatchHierarchy
More information about the CIG-COMMITS
mailing list