[cig-commits] commit: Add MDPI. Does not work all that much better.
Mercurial
hg at geodynamics.org
Sat Mar 12 05:37:11 PST 2011
changeset: 115:f4f4ca0e12df
user: Walter Landry <wlandry at caltech.edu>
date: Fri Mar 11 12:06:38 2011 -0800
files: P_MDPI_Refine.C P_MDPI_Refine.h dRc_dp.h input/shear_corner.input main.C wscript
description:
Add MDPI. Does not work all that much better.
diff -r 61f1ba62bedf -r f4f4ca0e12df P_MDPI_Refine.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/P_MDPI_Refine.C Fri Mar 11 12:06:38 2011 -0800
@@ -0,0 +1,145 @@
+/*************************************************************************
+ *
+ * This file is part of the SAMRAI distribution. For full copyright
+ * information, see COPYRIGHT and COPYING.LESSER.
+ *
+ * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
+ * Description: Linear refine operator for cell-centered double data on
+ * a Cartesian mesh.
+ *
+ ************************************************************************/
+
+#ifndef included_geom_P_MDPI_Refine_C
+#define included_geom_P_MDPI_Refine_C
+
+#include "P_MDPI_Refine.h"
+
+#include <float.h>
+#include <math.h>
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/hier/Index.h"
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/pdat/NodeData.h"
+#include "SAMRAI/tbox/Utilities.h"
+#include "dRc_dp.h"
+
+void SAMRAI::geom::P_MDPI_Refine::refine(
+ hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::BoxOverlap& fine_overlap,
+ const hier::IntVector& ratio) const
+{
+ const pdat::CellOverlap* t_overlap =
+ dynamic_cast<const pdat::CellOverlap *>(&fine_overlap);
+
+ TBOX_ASSERT(t_overlap != NULL);
+
+ const hier::BoxList& boxes = t_overlap->getDestinationBoxList();
+ for (hier::BoxList::Iterator b(boxes); b; b++) {
+ refine(fine,
+ coarse,
+ dst_component,
+ src_component,
+ b(),
+ ratio);
+ }
+}
+
+void SAMRAI::geom::P_MDPI_Refine::refine(
+ hier::Patch& fine_patch,
+ const hier::Patch& coarse_patch,
+ const int dst_component,
+ const int src_component,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio) const
+{
+ const tbox::Dimension& dim(getDim());
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dim, fine_patch, coarse_patch,
+ fine_box, ratio);
+
+ tbox::Pointer<pdat::CellData<double> >
+ p_ptr = coarse_patch.getPatchData(src_component);
+ pdat::CellData<double> &p(*p_ptr);
+ tbox::Pointer<pdat::CellData<double> >
+ p_fine_ptr = fine_patch.getPatchData(dst_component);
+ pdat::CellData<double> &p_fine(*p_fine_ptr);
+ tbox::Pointer<pdat::SideData<double> >
+ v_ptr = fine_patch.getPatchData(v_id);
+ pdat::SideData<double> &v(*v_ptr);
+ tbox::Pointer<pdat::CellData<double> >
+ cell_viscosity_ptr = fine_patch.getPatchData(cell_viscosity_id);
+ pdat::CellData<double> &cell_viscosity(*cell_viscosity_ptr);
+ tbox::Pointer<pdat::NodeData<double> >
+ edge_viscosity_ptr = fine_patch.getPatchData(edge_viscosity_id);
+ pdat::NodeData<double> &edge_viscosity(*edge_viscosity_ptr);
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+ TBOX_ASSERT(!p_ptr.isNull());
+ TBOX_ASSERT(!p_fine_ptr.isNull());
+ TBOX_ASSERT(p.getDepth() == p_fine.getDepth());
+#endif
+
+ hier::Box coarse_box=coarse_patch.getBox();
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = coarse_patch.getPatchGeometry();
+
+ const double dx = geom->getDx()[0];
+ const double dy = geom->getDx()[1];
+
+ for(int j=fine_box.lower(1); j<=fine_box.upper(1); ++j)
+ for(int i=fine_box.lower(0); i<=fine_box.upper(0); ++i)
+ {
+ pdat::CellIndex fine(tbox::Dimension(2));
+ fine[0]=i;
+ fine[1]=j;
+
+ pdat::CellIndex ip(fine), jp(fine);
+ ip[0]=1;
+ ip[1]=0;
+ jp[0]=0;
+ jp[1]=1;
+ pdat::CellIndex center(hier::Index::coarsen(fine,hier::Index(2,2)));
+
+ double dRc_dp_total(0), dRc_dp_fine(0);
+ /* This is horribly inefficient */
+ for(int xx=0;xx<2;++xx)
+ for(int yy=0;yy<2;++yy)
+ {
+ pdat::CellIndex c_fine(center*2);
+ c_fine[0]+=xx;
+ c_fine[1]+=yy;
+
+ pdat::CellIndex left(c_fine-ip),right(c_fine+ip),
+ down(c_fine-jp), up(c_fine+jp);
+ pdat::SideIndex left_x(left,0,pdat::SideIndex::Lower),
+ right_x(right,0,pdat::SideIndex::Lower),
+ down_y(down,1,pdat::SideIndex::Lower),
+ up_y(up,1,pdat::SideIndex::Lower);
+
+ double dRc_dp_weight=dRc_dp(fine_box,c_fine,left,right,down,up,
+ left_x,right_x,down_y,up_y,
+ cell_viscosity,edge_viscosity,v,
+ dx,dy);
+ if(c_fine==fine)
+ dRc_dp_fine=dRc_dp_weight;
+ dRc_dp_total+=dRc_dp_weight;
+ }
+
+ p_fine(fine)=p(center)*dRc_dp_total/(4*dRc_dp_fine);
+
+ // tbox::plog << "P_MDPI_Refine "
+ // << fine_patch.getPatchLevelNumber() << " "
+ // << i << " "
+ // << j << " "
+ // << (*p_fine)(fine) << " "
+ // << (*p)(center) << " "
+ // << dp_dx << " "
+ // << dp_dy << " "
+ // << "\n";
+
+ }
+}
+#endif
diff -r 61f1ba62bedf -r f4f4ca0e12df P_MDPI_Refine.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/P_MDPI_Refine.h Fri Mar 11 12:06:38 2011 -0800
@@ -0,0 +1,144 @@
+/*************************************************************************
+ *
+ * This file is part of the SAMRAI distribution. For full copyright
+ * information, see COPYRIGHT and COPYING.LESSER.
+ *
+ * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
+ * Description: Linear refine operator for cell-centered double data on
+ * a Cartesian mesh.
+ *
+ ************************************************************************/
+
+#ifndef included_geom_P_MDPI_Refine
+#define included_geom_P_MDPI_Refine
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#include "SAMRAI/xfer/RefineOperator.h"
+#include "SAMRAI/hier/Box.h"
+#include "SAMRAI/hier/IntVector.h"
+#include "SAMRAI/hier/Patch.h"
+#include "SAMRAI/tbox/Pointer.h"
+#include "SAMRAI/pdat/CellVariable.h"
+
+#include <string>
+
+namespace SAMRAI {
+namespace geom {
+
+/**
+ * Class P_MDPI_Refine implements matrix dependent
+ * interpolation for cell-centered double patch data defined over a Cartesian
+ * mesh. It should work better for variable viscosity.
+ *
+ * The findRefineOperator() operator function returns true if the input
+ * variable is cell-centered double, and the std::string is "P_MDPI_REFINE".
+ *
+ * @see xfer::RefineOperator
+ */
+
+class P_MDPI_Refine:
+ public xfer::RefineOperator
+{
+public:
+ /**
+ * Uninteresting default constructor.
+ */
+ explicit P_MDPI_Refine(const tbox::Dimension& dim, const int &v,
+ const int &cell_viscosity,
+ const int &edge_viscosity):
+ xfer::RefineOperator(dim, "P_MDPI_REFINE"),
+ v_id(v),
+ cell_viscosity_id(cell_viscosity),
+ edge_viscosity_id(edge_viscosity)
+ {
+ d_name_id = "P_MDPI_REFINE";
+ }
+
+
+ /**
+ * Uninteresting virtual destructor.
+ */
+ virtual ~P_MDPI_Refine(){}
+
+ /**
+ * Return true if the variable and name std::string match cell-centered
+ * double linear interpolation; otherwise, return false.
+ */
+ bool findRefineOperator(const tbox::Pointer<hier::Variable>& var,
+ const std::string& op_name) const
+ {
+ TBOX_DIM_ASSERT_CHECK_ARGS2(*this, *var);
+
+ const tbox::Pointer<pdat::CellVariable<double> > cast_var(var);
+ if (!cast_var.isNull() && (op_name == d_name_id)) {
+ return true;
+ } else {
+ return false;
+ }
+ }
+ /**
+ * Return name std::string identifier of this refinement operator.
+ */
+ const std::string& getOperatorName() const
+ {
+ return d_name_id;
+ }
+
+ /**
+ * The priority of cell-centered double linear interpolation is 0.
+ * It will be performed before any user-defined interpolation operations.
+ */
+ int getOperatorPriority() const
+ {
+ return 0;
+ }
+
+ /**
+ * The stencil width of the linear interpolation operator is the vector
+ * of ones. That is, its stencil extends one cell outside the fine box.
+ */
+ hier::IntVector getStencilWidth() const
+ {
+ return hier::IntVector::getOne(getDim());
+ }
+
+ /**
+ * Refine the source component on the coarse patch to the destination
+ * component on the fine patch using the cell-centered double linear
+ * interpolation operator. Interpolation is performed on the intersection
+ * of the destination patch and the boxes contained in fine_overlap
+ * It is assumed that the coarse patch contains sufficient data for the
+ * stencil width of the refinement operator.
+ */
+ void refine(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::BoxOverlap& fine_overlap,
+ const hier::IntVector& ratio) const;
+
+ /**
+ * Refine the source component on the coarse patch to the destination
+ * component on the fine patch using the cell-centered double linear
+ * interpolation operator. Interpolation is performed on the intersection
+ * of the destination patch and the fine box. It is assumed that the
+ * coarse patch contains sufficient data for the stencil width of the
+ * refinement operator. This differs from the above refine() method
+ * only in that it operates on a single fine box instead of a BoxOverlap.
+ */
+ void refine(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio) const;
+
+private:
+ std::string d_name_id;
+ int v_id, cell_viscosity_id, edge_viscosity_id;
+};
+
+}
+}
+#endif
diff -r 61f1ba62bedf -r f4f4ca0e12df dRc_dp.h
--- a/dRc_dp.h Fri Mar 11 10:59:45 2011 -0800
+++ b/dRc_dp.h Fri Mar 11 12:06:38 2011 -0800
@@ -5,6 +5,7 @@
#include "SAMRAI/pdat/NodeData.h"
#include "SAMRAI/pdat/SideData.h"
#include "dRm_dv.h"
+#include "Boundary.h"
/* The derivative of the continuity equation with respect to
pressure. Note that pressure does not appear in the continuity
diff -r 61f1ba62bedf -r f4f4ca0e12df input/shear_corner.input
--- a/input/shear_corner.input Fri Mar 11 10:59:45 2011 -0800
+++ b/input/shear_corner.input Fri Mar 11 12:06:38 2011 -0800
@@ -42,17 +42,9 @@ FACStokes {
// coarse_solver_choice = "Gerya"
coarse_solver_max_iterations = 25000
coarse_solver_tolerance = 1e-8
- p_prolongation_method = "P_REFINE" // Type of refinement
- // used in prolongation.
- // Suggested values are
- // "LINEAR_REFINE"
- // "CONSTANT_REFINE"
- v_prolongation_method = "V_REFINE" // Type of refinement
- // used in prolongation.
- // Suggested values are
- // "LINEAR_REFINE"
- // "CONSTANT_REFINE"
- use_smg = TRUE // Whether to use HYPRE's SMG instead of PFMG.
+ // p_prolongation_method = "P_MDPI_REFINE"
+ p_prolongation_method = "P_REFINE"
+ v_prolongation_method = "V_REFINE"
}
bc_coefs {
// These are the boundary condition specifications. The number
diff -r 61f1ba62bedf -r f4f4ca0e12df main.C
--- a/main.C Fri Mar 11 10:59:45 2011 -0800
+++ b/main.C Fri Mar 11 12:06:38 2011 -0800
@@ -36,6 +36,7 @@ using namespace std;
#include "V_Coarsen.h"
#include "Cell_Viscosity_Coarsen.h"
#include "Edge_Viscosity_Coarsen.h"
+#include "P_MDPI_Refine.h"
#include "FACStokes.h"
@@ -198,6 +199,11 @@ int main(
(tbox::Pointer<SAMRAI::xfer::CoarsenOperator>
(new SAMRAI::geom::Edge_Viscosity_Coarsen
(dim,fac_stokes.cell_viscosity_id)));
+ grid_geometry->addSpatialRefineOperator
+ (tbox::Pointer<SAMRAI::xfer::RefineOperator>
+ (new SAMRAI::geom::P_MDPI_Refine(dim,fac_stokes.v_id,
+ fac_stokes.cell_viscosity_id,
+ fac_stokes.edge_viscosity_id)));
/*
* Create the tag-and-initializer, box-generator and load-balancer
diff -r 61f1ba62bedf -r f4f4ca0e12df wscript
--- a/wscript Fri Mar 11 10:59:45 2011 -0800
+++ b/wscript Fri Mar 11 12:06:38 2011 -0800
@@ -19,6 +19,7 @@ def build(bld):
'FACStokes/setupPlotter.C',
'FACStokes/solveStokes.C',
'P_Refine.C',
+ 'P_MDPI_Refine.C',
'V_Refine.C',
'V_Coarsen.C',
'P_Boundary_Refine.C',
More information about the CIG-COMMITS
mailing list