Adding specialized ASPH H advance options (fixShape and radialOnly),

and creating a special H time advance policy.  Note when using these
options the ASPH idealH does not require the Voronoi (big time savings).

src/PYB11/SmoothingScale/ASPHSmoothingScale.py

@@ -21,7 +21,9 @@ class ASPHSmoothingScale(SmoothingScaleBase):
     # Constructors
     def pyinit(self,
                HUpdate = "HEvolutionType",
-               W = "const TableKernel<%(Dimension)s>&"):
+               W = "const TableKernel<%(Dimension)s>&",
+               fixShape = ("const bool", "false"),
+               radialOnly = ("const bool", "false")):
         "ASPHSmoothingScale constructor"
 
     #...........................................................................
@@ -95,3 +97,5 @@ def label(self):
     secondMoment = PYB11property("const FieldList<%(Dimension)s, SymTensor>&", "secondMoment", doc="The second moment storage FieldList")
     cellSecondMoment = PYB11property("const FieldList<%(Dimension)s, SymTensor>&", "cellSecondMoment", doc="The second moment of the Voronoi cells")
     HidealFilter = PYB11property("std::shared_ptr<HidealFilterType>", "HidealFilter", "HidealFilter", doc="Optional function to manipulate the Hideal calculation")
+    fixShape = PYB11property("bool", "fixShape", "fixShape", doc="Force the H tensor shape to be fixed -- only adjust volume")
+    radialOnly = PYB11property("bool", "radialOnly", "radialOnly", doc="Force the H tensor to evolve solely in the radial direction")

src/SmoothingScale/ASPHSmoothingScale.hh

@@ -29,7 +29,9 @@ public:
 
   // Constructors, destructor.
   ASPHSmoothingScale(const HEvolutionType HUpdate,
-                     const TableKernel<Dimension>& W);
+                     const TableKernel<Dimension>& W,
+                     const bool fixShape = false,
+                     const bool radialOnly = false);
   ASPHSmoothingScale() = delete;
   virtual ~ASPHSmoothingScale() {}
 
@@ -71,14 +73,21 @@ public:
                                     StateDerivatives<Dimension>& derivs) override;
 
   // We require the Voronoi-like cells per point
-  virtual bool requireVoronoiCells() const override                               { return this->HEvolution() == HEvolutionType::IdealH; }
+  virtual bool requireVoronoiCells() const override                               { return (this->HEvolution() == HEvolutionType::IdealH and
+                                                                                            not (mFixShape or mRadialOnly)); }
 
   // Access our internal data
   const TableKernel<Dimension>&                          WT()            const    { return mWT; }
   const FieldList<Dimension, Scalar>&                    zerothMoment()  const    { return mZerothMoment; }
   const FieldList<Dimension, SymTensor>&                 secondMoment()  const    { return mSecondMoment; }
   const FieldList<Dimension, SymTensor>&                 cellSecondMoment() const { return mCellSecondMoment; }
 
+  // Special evolution flags
+  bool fixShape() const                                                           { return mFixShape; }
+  bool radialOnly() const                                                         { return mRadialOnly; }
+  void fixShape(const bool x)                                                     { mFixShape = x; }
+  void radialOnly(const bool x)                                                   { mRadialOnly = x; }
+
   // Optional user hook providing a functor to manipulate the ideal H vote
   std::shared_ptr<HidealFilterType> HidealFilter() const                          { return mHidealFilterPtr; }
   void HidealFilter(std::shared_ptr<HidealFilterType> functorPtr)                 { mHidealFilterPtr = functorPtr; }
@@ -94,6 +103,7 @@ private:
   FieldList<Dimension, Scalar> mZerothMoment;
   FieldList<Dimension, SymTensor> mSecondMoment, mCellSecondMoment;
   std::shared_ptr<HidealFilterType> mHidealFilterPtr;
+  bool mFixShape, mRadialOnly;
 };
 
 }

src/SmoothingScale/CMakeLists.txt

@@ -3,6 +3,7 @@ set(SmoothingScale_inst
   SmoothingScaleBase
   SPHSmoothingScale
   ASPHSmoothingScale
+  IncrementASPHHtensor
 )
 
 set(SmoothingScale_sources

src/SmoothingScale/IncrementASPHHtensor.cc

@@ -0,0 +1,120 @@
+//---------------------------------Spheral++----------------------------------//
+// IncrementASPHHtensor
+//
+// Specialized version of FieldUpdatePolicy for time integrating the H tensor.
+//
+// Created by JMO, Mon Oct  7 13:31:02 PDT 2024
+//----------------------------------------------------------------------------//
+#include "IncrementASPHHtensor.hh"
+#include "DataBase/State.hh"
+#include "DataBase/StateDerivatives.hh"
+#include "Field/Field.hh"
+#include "Hydro/HydroFieldNames.hh"
+#include "Utilities/rotationMatrix.hh"
+#include "Utilities/GeometricUtilities.hh"
+#include "Utilities/DBC.hh"
+
+namespace Spheral {
+
+//------------------------------------------------------------------------------
+// Constructor
+//------------------------------------------------------------------------------
+template<typename Dimension>
+IncrementASPHHtensor<Dimension>::
+IncrementASPHHtensor(const Scalar hmin,
+                     const Scalar hmax,
+                     const Scalar hminratio,
+                     const bool fixShape,
+                     const bool radialOnly):
+  FieldUpdatePolicy<Dimension>(),
+  mhmin(hmin),
+  mhmax(hmax),
+  mhminratio(hminratio),
+  mFixShape(fixShape),
+  mRadialOnly(radialOnly) {
+}
+
+//------------------------------------------------------------------------------
+// Update the field.
+//------------------------------------------------------------------------------
+template<typename Dimension>
+void
+IncrementASPHHtensor<Dimension>::
+update(const KeyType& key,
+       State<Dimension>& state,
+       StateDerivatives<Dimension>& derivs,
+       const double multiplier,
+       const double t,
+       const double dt) {
+
+  // Get the field name portion of the key.
+  KeyType fieldKey, nodeListKey;
+  StateBase<Dimension>::splitFieldKey(key, fieldKey, nodeListKey);
+  CHECK(fieldKey == HydroFieldNames::H);
+
+  const auto hminInv = 1.0/mhmin;
+  const auto hmaxInv = 1.0/mhmax;
+
+  // Get the state we're updating.
+  auto&       H = state.field(key, SymTensor::zero);
+  const auto& DHDt = derivs.field(prefix() + StateBase<Dimension>::buildFieldKey(HydroFieldNames::H, nodeListKey), SymTensor::zero);
+  const auto& pos = state.field(StateBase<Dimension>::buildFieldKey(HydroFieldNames::position, nodeListKey), Vector::zero);  // Only needed if we're using radial scaling
+
+  // Walk the nodes and update H (with limiting)
+  const auto n = H.numInternalElements();
+#pragma omp parallel for
+  for (auto i = 0u; i < n; ++i) {
+
+    // Check for special update rules
+    if (mFixShape) {
+
+      // Fix the shape (only volume scaling allowed)
+
+      auto fi = Dimension::rootnu((H(i) + multiplier*DHDt(i)).Determinant()/H(i).Determinant());
+      H(i) *= fi;
+
+    } else if (mRadialOnly) {
+
+      // Force only the radial component of H to be scaled
+      const auto nhat = pos(i).unitVector();
+      const auto T = rotationMatrix(nhat);
+      H(i).rotationalTransform(T);          // Should have one eigenvector aligned with the x' axis in this frame
+      auto DHDti = DHDt(i);
+      DHDti.rotationalTransform(T);
+      H(i)[0] += multiplier * DHDti[0];
+      H(i).rotationalTransform(T.Transpose());
+
+    } else {
+
+      H(i) += multiplier * DHDt(i);
+
+    }
+
+    // Apply limiting
+    const auto hev = H(i).eigenVectors();
+    const auto hminEffInv = min(hminInv, max(hmaxInv, hev.eigenValues.minElement())/mhminratio);
+    H(i) = constructSymTensorWithBoundedDiagonal(hev.eigenValues, hmaxInv, hminEffInv);
+    H(i).rotationalTransform(hev.eigenVectors);
+  }
+}
+
+//------------------------------------------------------------------------------
+// Equivalence operator.
+//------------------------------------------------------------------------------
+template<typename Dimension>
+inline
+bool
+IncrementASPHHtensor<Dimension>::
+operator==(const UpdatePolicyBase<Dimension>& rhs) const {
+  const auto rhsPtr = dynamic_cast<const IncrementASPHHtensor<Dimension>*>(&rhs);
+  if (rhsPtr == nullptr) return false;
+
+  // Ok, now do we agree on min & max?
+  return (hmin() == rhsPtr->hmin() and
+          hmax() == rhsPtr->hmax() and
+          hminratio() == rhsPtr->hminratio() and
+          fixShape() == rhsPtr->fixShape() and
+          radialOnly() == rhsPtr->radialOnly());
+}
+
+}

src/SmoothingScale/IncrementASPHHtensor.hh

@@ -0,0 +1,66 @@
+//---------------------------------Spheral++----------------------------------//
+// IncrementASPHHtensor
+//
+// Specialized version of FieldUpdatePolicy for time integrating the H tensor.
+//
+// Created by JMO, Mon Oct  7 13:31:02 PDT 2024
+//----------------------------------------------------------------------------//
+#ifndef __Spheral_IncrementASPHHtensor_hh__
+#define __Spheral_IncrementASPHHtensor_hh__
+
+#include "DataBase/FieldUpdatePolicy.hh"
+
+namespace Spheral {
+
+// Forward declarations.
+template<typename Dimension> class StateDerivatives;
+
+template<typename Dimension>
+class IncrementASPHHtensor: public FieldUpdatePolicy<Dimension> {
+public:
+  //--------------------------- Public Interface ---------------------------//
+  // Useful typedefs
+  using KeyType = typename FieldUpdatePolicy<Dimension>::KeyType;
+  using Scalar = typename Dimension::Scalar;
+  using Vector = typename Dimension::Vector;
+  using SymTensor = typename Dimension::SymTensor;
+
+  // Constructors, destructor.
+  IncrementASPHHtensor(const Scalar hmin,
+                       const Scalar hmax,
+                       const Scalar hminratio,
+                       const bool fixShape,
+                       const bool radialOnly);
+  virtual ~IncrementASPHHtensor()   {}
+  IncrementASPHHtensor(const IncrementASPHHtensor& rhs) = delete;
+  IncrementASPHHtensor& operator=(const IncrementASPHHtensor& rhs) = delete;
+
+  // Overload the methods describing how to update Fields.
+  virtual void update(const KeyType& key,
+                      State<Dimension>& state,
+                      StateDerivatives<Dimension>& derivs,
+                      const double multiplier,
+                      const double t,
+                      const double dt) override;
+
+  // Access the min and max's.
+  Scalar hmin() const               { return mhmin; }
+  Scalar hmax() const               { return mhmax; }
+  Scalar hminratio() const          { return mhminratio; }
+  bool fixShape() const             { return mFixShape; }
+  bool radialOnly() const           { return mRadialOnly; }
+
+  // Equivalence.
+  virtual bool operator==(const UpdatePolicyBase<Dimension>& rhs) const override;
+
+  static const std::string prefix() { return "delta "; }
+
+private:
+  //--------------------------- Private Interface ---------------------------//
+  Scalar mhmin, mhmax, mhminratio;
+  bool mFixShape, mRadialOnly;
+};
+
+}
+
+#endif

src/SmoothingScale/IncrementASPHHtensorInst.cc.py

@@ -0,0 +1,11 @@
+text = """
+//------------------------------------------------------------------------------
+// Explicit instantiation.
+//------------------------------------------------------------------------------
+#include "SmoothingScale/IncrementASPHHtensor.cc"
+#include "Geometry/Dimension.hh"
+
+namespace Spheral {
+  template class IncrementASPHHtensor<Dim<%(ndim)s>>;
+}
+"""

tests/functional/RK/RKInterpolation.py

@@ -232,7 +232,7 @@
 # Randomize nodes
 #-------------------------------------------------------------------------------
 import random
-seed = 2
+seed = 459297849234
 rangen = random.Random()
 rangen.seed(seed)