Refactor volume calculation

src/ArtificialConduction/ArtificialConduction.cc

@@ -300,9 +300,10 @@ dt(const DataBase<Dimension>& dataBase,
   const FieldList<Dimension, Scalar> eps = state.fields(HydroFieldNames::specificThermalEnergy, Scalar());
   const FieldList<Dimension, Scalar> soundSpeed = state.fields(HydroFieldNames::soundSpeed, 0.0);
   const FieldList<Dimension, Scalar> vsigMax = state.fields("Maximum Artificial Cond Signal Speed", 0.0);
-  const ConnectivityMap<Dimension>& connectivityMap = dataBase.connectivityMap(this->requireGhostConnectivity(),
-                                                                               this->requireOverlapConnectivity(),
-                                                                               this->requireIntersectionConnectivity());
+  const auto& [conn, ghost, overlap, intersect] = this->requireConnectivity();
+  const ConnectivityMap<Dimension>& connectivityMap = dataBase.connectivityMap(ghost,
+                                                                               overlap,
+                                                                               intersect);
   //const int numNodeLists = connectivityMap.nodeLists().size();
 
   // Initialize the return value to some impossibly high value.

src/ArtificialViscosity/TensorCRKSPHViscosity.hh

@@ -20,6 +20,7 @@ public:
   using Tensor = typename Dimension::Tensor;
   using SymTensor = typename Dimension::SymTensor;
   using ThirdRankTensor = typename Dimension::ThirdRankTensor;
+  using RKRequirements = typename Physics<Dimension>::RKRequirements;
 
   // Constructors, destructor
   TensorCRKSPHViscosity(const Scalar Clinear,
@@ -42,7 +43,7 @@ public:
                                       const StateDerivatives<Dimension>& derivs) override;
 
   // Override the method listing our RK requirements
-  virtual std::set<RKOrder> requireReproducingKernels() const override { return std::set<RKOrder>({mOrder}); }
+  virtual RKRequirements requireReproducingKernels() const override { return {{mOrder}, {}, false}; }
 
   RKOrder order()                                       const          { return mOrder; }
   void order(const RKOrder x)                                          { mOrder = x; }

src/Boundary/Boundary.hh

@@ -176,6 +176,10 @@ public:
   // Provide a hook to note such cases.
   virtual bool meshGhostNodes() const                                            { return true; };
 
+  // Controls the ordering of boundaries when applied.
+  // Lower priority boundaries are applied first.
+  virtual int priority() const { return 0; };
+
   // Allow read access to the map of NodeList->BoundaryNodes.
   const std::map<NodeList<Dimension>*, BoundaryNodes>& boundaryNodeMap() const;
 

src/Boundary/ConstantBoundary.hh

@@ -67,6 +67,9 @@ public:
 
   // Apply the boundary condition to the violation node values in the given Field.
   virtual void enforceBoundary(FieldBase<Dimension>& fieldBase) const override;
+
+  // ConstantBoundary ghosts must precede distributed ghosts
+  virtual int priority() const override { return -4; }
   //**********************************************************************
 
   // After physics have been initialized we take a snapshot of the node state.

src/Boundary/InflowOutflowBoundary.hh

@@ -75,6 +75,9 @@ public:
   // We need to not cull ghost nodes, since they might need to cross the boundary
   // and become new inflow nodes.
   virtual bool allowGhostCulling() const override { return false; }
+
+  // InflowOutflow ghosts must precede distributed ghosts
+  virtual int priority() const override { return -2; }
   //**********************************************************************
 
   //**********************************************************************

src/CRKSPH/CRKSPHBase.cc

@@ -348,10 +348,10 @@ enforceBoundaries(State<Dimension>& state,
 // Return the RK orders we want to use
 //------------------------------------------------------------------------------
 template<typename Dimension>
-std::set<RKOrder>
+typename CRKSPHBase<Dimension>::RKRequirements
 CRKSPHBase<Dimension>::
 requireReproducingKernels() const {
-  return std::set<RKOrder>({RKOrder::ZerothOrder, mOrder});
+  return {{RKOrder::ZerothOrder, mOrder}, {}, false};
 }
 
 //------------------------------------------------------------------------------

src/CRKSPH/CRKSPHBase.hh

@@ -41,6 +41,7 @@ public:
   using FacetedVolume = typename Dimension::FacetedVolume;
 
   using ConstBoundaryIterator = typename Physics<Dimension>::ConstBoundaryIterator;
+  using RKRequirements = typename Physics<Dimension>::RKRequirements;
 
   // Constructors.
   CRKSPHBase(DataBase<Dimension>& dataBase,
@@ -103,7 +104,7 @@ public:
                          StateDerivatives<Dimension>& derivs) override;
 
   // We require RK corrections
-  virtual std::set<RKOrder> requireReproducingKernels() const override;
+  virtual RKRequirements requireReproducingKernels() const override;
 
   // The spatial order
   RKOrder correctionOrder()                                             const { return mOrder; }

src/CRKSPH/CRKSPHRZ.cc

@@ -28,6 +28,7 @@
 #include "Field/NodeIterators.hh"
 #include "Boundary/Boundary.hh"
 #include "Neighbor/ConnectivityMap.hh"
+#include "VoronoiCells/GeometryScaling.hh"
 #include "Neighbor/PairwiseField.hh"
 #include "Utilities/safeInv.hh"
 #include "Utilities/range.hh"
@@ -90,28 +91,11 @@ initializeProblemStartupDependencies(DataBase<Dim<2>>& dataBase,
                                      State<Dim<2>>& state,
                                      StateDerivatives<Dim<2>>& derivs) {
 
-  // Correct the mass to mass/r.
   auto mass = dataBase.fluidMass();
   const auto pos = dataBase.fluidPosition();
-  const unsigned numNodeLists = mass.numFields();
-  for (unsigned nodeListi = 0; nodeListi != numNodeLists; ++nodeListi) {
-    const unsigned n = mass[nodeListi]->numElements();
-    for (unsigned i = 0; i != n; ++i) {
-      const Scalar circi = 2.0*M_PI*abs(pos(nodeListi, i).y());
-      mass(nodeListi, i) /= circi;
-    }
-  }
-
-  // Do general initializations.
-  CRKSPHBase<Dim<2>>::initializeProblemStartupDependencies(dataBase, state, derivs);
-
-  // Convert back to mass.
-  for (unsigned nodeListi = 0; nodeListi != numNodeLists; ++nodeListi) {
-    const unsigned n = mass[nodeListi]->numElements();
-    for (unsigned i = 0; i != n; ++i) {
-      const Scalar circi = 2.0*M_PI*abs(pos(nodeListi, i).y());
-      mass(nodeListi, i) *= circi;
-    }
+  {
+    auto guard = unscaledRegion(pos, mass);
+    CRKSPHBase<Dim<2>>::initializeProblemStartupDependencies(dataBase, state, derivs);
   }
 }
 
@@ -179,29 +163,11 @@ preStepInitialize(const DataBase<Dim<2>>& dataBase,
                   State<Dim<2>>& state,
                   StateDerivatives<Dim<2>>& derivs) {
 
-  // Convert the mass to mass per unit length first.
   auto mass = state.fields(HydroFieldNames::mass, 0.0);
   const auto pos = state.fields(HydroFieldNames::position, Vector::zero());
-  const unsigned numNodeLists = mass.numFields();
-  for (unsigned nodeListi = 0; nodeListi != numNodeLists; ++nodeListi) {
-    const unsigned n = mass[nodeListi]->numElements();
-    for (unsigned i = 0; i != n; ++i) {
-      const auto circi = 2.0*M_PI*abs(pos(nodeListi, i).y());
-      mass(nodeListi, i) /= circi;
-    }
-  }
-
-  // Base class does most of the work.
-  CRKSPHBase<Dimension>::preStepInitialize(dataBase, state, derivs);
-
-  // Now convert back to true masses.
-  for (unsigned nodeListi = 0; nodeListi != numNodeLists; ++nodeListi) {
-    const unsigned n = mass[nodeListi]->numElements();
-    for (unsigned i = 0; i != n; ++i) {
-      const auto& xi = pos(nodeListi, i);
-      const auto circi = 2.0*M_PI*abs(xi.y());
-      mass(nodeListi, i) *= circi;
-    }
+  {
+    auto guard = unscaledRegion(pos, mass);
+    CRKSPHBase<Dimension>::preStepInitialize(dataBase, state, derivs);
   }
 }
 
@@ -514,31 +480,12 @@ CRKSPHRZ::
 applyGhostBoundaries(State<Dim<2>>& state,
                      StateDerivatives<Dim<2>>& derivs) {
 
-  // Convert the mass to mass/length before BCs are applied.
-  FieldList<Dimension, Scalar> mass = state.fields(HydroFieldNames::mass, 0.0);
-  const FieldList<Dimension, Vector> pos = state.fields(HydroFieldNames::position, Vector::zero());
-  const unsigned numNodeLists = mass.numFields();
-  for (unsigned nodeListi = 0; nodeListi != numNodeLists; ++nodeListi) {
-    const unsigned n = mass[nodeListi]->numElements();
-    for (unsigned i = 0; i != n; ++i) {
-      const Scalar circi = 2.0*M_PI*abs(pos(nodeListi, i).y());
-      CHECK(circi > 0.0);
-      mass(nodeListi, i) /= circi;
-    }
-  }
-
-  // Apply ordinary CRKSPH BCs.
-  CRKSPHBase<Dim<2>>::applyGhostBoundaries(state, derivs);
-  for (auto boundaryPtr: range(this->boundaryBegin(), this->boundaryEnd())) boundaryPtr->finalizeGhostBoundary();
-
-  // Scale back to mass.
-  for (unsigned nodeListi = 0; nodeListi != numNodeLists; ++nodeListi) {
-    const unsigned n = mass[nodeListi]->numElements();
-    for (unsigned i = 0; i != n; ++i) {
-      const Scalar circi = 2.0*M_PI*abs(pos(nodeListi, i).y());
-      CHECK(circi > 0.0);
-      mass(nodeListi, i) *= circi;
-    }
+  auto mass = state.fields(HydroFieldNames::mass, 0.0);
+  const auto pos = state.fields(HydroFieldNames::position, Vector::zero());
+  {
+    auto guard = unscaledRegion(pos, mass);
+    CRKSPHBase<Dim<2>>::applyGhostBoundaries(state, derivs);
+    for (auto boundaryPtr: range(this->boundaryBegin(), this->boundaryEnd())) boundaryPtr->finalizeGhostBoundary();
   }
 }
 
@@ -550,33 +497,11 @@ CRKSPHRZ::
 enforceBoundaries(State<Dim<2>>& state,
                   StateDerivatives<Dim<2>>& derivs) {
 
-  // Convert the mass to mass/length before BCs are applied.
-  FieldList<Dimension, Scalar> mass = state.fields(HydroFieldNames::mass, 0.0);
-  FieldList<Dimension, Vector> pos = state.fields(HydroFieldNames::position, Vector::zero());
-  const unsigned numNodeLists = mass.numFields();
-  for (unsigned nodeListi = 0; nodeListi != numNodeLists; ++nodeListi) {
-    const unsigned n = mass[nodeListi]->numInternalElements();
-    for (unsigned i = 0; i != n; ++i) {
-      const Scalar circi = 2.0*M_PI*abs(pos(nodeListi, i).y());
-      CHECK(circi > 0.0);
-      mass(nodeListi, i) /= circi;
-    }
-  }
-
-  // Apply ordinary CRKSPH BCs.
-  CRKSPHBase<Dim<2>>::enforceBoundaries(state, derivs);
-
-  // Scale back to mass.
-  // We also ensure no point approaches the z-axis too closely.
-  FieldList<Dimension, SymTensor> H = state.fields(HydroFieldNames::H, SymTensor::zero());
-  for (unsigned nodeListi = 0; nodeListi != numNodeLists; ++nodeListi) {
-    const unsigned n = mass[nodeListi]->numInternalElements();
-    //const Scalar nPerh = mass[nodeListi]->nodeList().nodesPerSmoothingScale();
-    for (unsigned i = 0; i != n; ++i) {
-      Vector& posi = pos(nodeListi, i);
-      const Scalar circi = 2.0*M_PI*abs(posi.y());
-      mass(nodeListi, i) *= circi;
-    }
+  auto mass = state.fields(HydroFieldNames::mass, 0.0);
+  const auto pos = state.fields(HydroFieldNames::position, Vector::zero());
+  {
+    auto guard = unscaledInternalRegion(pos, mass);
+    CRKSPHBase<Dim<2>>::enforceBoundaries(state, derivs);
   }
 }