Adding ability to override radial geometry for radialOnly option in ASPH

src/PYB11/CXXTypes/CXXTypes_PYB11.py

@@ -74,7 +74,8 @@ def pyinit(self,
 # RKCoefficients
 
 for ndim in dims:
-    exec('''
-vector_of_RKCoefficients%(ndim)id = PYB11_bind_vector("Spheral::RKCoefficients<%(Dimension)s>", opaque=True, local=False)
-''' % {"ndim"      : ndim,
-       "Dimension" : "Spheral::Dim<" + str(ndim) + ">"})
+    Dimension = f"Spheral::Dim<{ndim}>"
+    Vector = f"{Dimension}::Vector"
+    exec(f'''
+vector_of_RKCoefficients{ndim}d = PYB11_bind_vector("Spheral::RKCoefficients<{Dimension}>", opaque=True, local=False)
+''')

src/PYB11/SmoothingScale/ASPHRadialFunctor.py

@@ -0,0 +1,38 @@
+#-------------------------------------------------------------------------------
+# ASPHRadialFunctor
+#-------------------------------------------------------------------------------
+from PYB11Generator import *
+
+@PYB11template("Dimension")
+@PYB11holder("std::shared_ptr")
+class ASPHRadialFunctor:
+
+    PYB11typedefs = """
+    using Scalar = typename %(Dimension)s::Scalar;
+    using Vector = typename %(Dimension)s::Vector;
+    using Tensor = typename %(Dimension)s::Tensor;
+    using SymTensor = typename %(Dimension)s::SymTensor;
+"""
+
+    #...........................................................................
+    # Constructors
+    def pyinit(self):
+        "ASPHRadialFunctor constructor"
+
+    #...........................................................................
+    # Virtual methods
+    @PYB11virtual
+    @PYB11const
+    def radialUnitVector(self,
+                         nodeListi = "const size_t",
+                         i = "const size_t",
+                         posi = "const Vector&"):
+        return "Vector"
+
+    @PYB11virtual
+    @PYB11const
+    def radialCoordinate(self,
+                         nodeListi = "const size_t",
+                         i = "const size_t",
+                         posi = "const Vector&"):
+        return "Scalar"

src/PYB11/SmoothingScale/ASPHSmoothingScale.py

@@ -15,6 +15,7 @@ class ASPHSmoothingScale(SmoothingScaleBase):
     using ThirdRankTensor = typename %(Dimension)s::ThirdRankTensor;
     using TimeStepType = typename Physics<%(Dimension)s>::TimeStepType;
     using HidealFilterType = typename ASPHSmoothingScale<%(Dimension)s>::HidealFilterType;
+    using RadialFunctorType = typename ASPHSmoothingScale<%(Dimension)s>::RadialFunctorType;
 """
 
     #...........................................................................
@@ -105,6 +106,7 @@ 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")
     radius0 = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "radius0", doc="The radius of a point at the beginning of a timestep (if using radialOnly=True)")
-    HidealFilter = PYB11property("std::shared_ptr<HidealFilterType>", "HidealFilter", "HidealFilter", doc="Optional function to manipulate the Hideal calculation")
+    HidealFilter = PYB11property("std::shared_ptr<HidealFilterType>", "HidealFilter", "HidealFilter", doc="Optional functor to manipulate the Hideal calculation")
+    RadialFunctor = PYB11property("std::shared_ptr<RadialFunctorType>", "RadialFunctor", "RadialFunctor", doc="Optional functor to manipulate the radial normal and radius are computed when using radialOnly")
     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/PYB11/SmoothingScale/ASPHSmoothingScaleUserFilter.py

@@ -1,5 +1,5 @@
 #-------------------------------------------------------------------------------
-# ASPHSmoothingScale
+# ASPHSmoothingScaleUserFilter
 #-------------------------------------------------------------------------------
 from PYB11Generator import *
 

src/PYB11/SmoothingScale/SmoothingScale_PYB11.py

@@ -17,6 +17,7 @@
                   '"SmoothingScale/SPHSmoothingScale.hh"',
                   '"SmoothingScale/ASPHSmoothingScale.hh"',
                   '"SmoothingScale/ASPHSmoothingScaleUserFilter.hh"',
+                  '"SmoothingScale/ASPHRadialFunctor.hh"',
                   '"SmoothingScale/polySecondMoment.hh"',
                   '"Kernel/TableKernel.hh"',
                   '"Neighbor/ConnectivityMap.hh"',
@@ -41,20 +42,24 @@
 from SPHSmoothingScale import SPHSmoothingScale
 from ASPHSmoothingScale import ASPHSmoothingScale
 from ASPHSmoothingScaleUserFilter import ASPHSmoothingScaleUserFilter
+from ASPHRadialFunctor import ASPHRadialFunctor
 
 for ndim in dims:
+    Dimension = f"Dim<{ndim}>"
+    Vector = f"{Dimension}::Vector"
     exec(f'''
-SmoothingScaleBase{ndim}d = PYB11TemplateClass(SmoothingScaleBase, template_parameters="Dim<{ndim}>")
-FixedSmoothingScale{ndim}d = PYB11TemplateClass(FixedSmoothingScale, template_parameters="Dim<{ndim}>")
-SPHSmoothingScale{ndim}d = PYB11TemplateClass(SPHSmoothingScale, template_parameters="Dim<{ndim}>")
-ASPHSmoothingScale{ndim}d = PYB11TemplateClass(ASPHSmoothingScale, template_parameters="Dim<{ndim}>")
-ASPHSmoothingScaleUserFilter{ndim}d = PYB11TemplateClass(ASPHSmoothingScaleUserFilter, template_parameters="Dim<{ndim}>")
+SmoothingScaleBase{ndim}d = PYB11TemplateClass(SmoothingScaleBase, template_parameters="{Dimension}")
+FixedSmoothingScale{ndim}d = PYB11TemplateClass(FixedSmoothingScale, template_parameters="{Dimension}")
+SPHSmoothingScale{ndim}d = PYB11TemplateClass(SPHSmoothingScale, template_parameters="{Dimension}")
+ASPHSmoothingScale{ndim}d = PYB11TemplateClass(ASPHSmoothingScale, template_parameters="{Dimension}")
+ASPHSmoothingScaleUserFilter{ndim}d = PYB11TemplateClass(ASPHSmoothingScaleUserFilter, template_parameters="{Dimension}")
+ASPHRadialFunctor{ndim}d = PYB11TemplateClass(ASPHRadialFunctor, template_parameters="{Dimension}")
 
 @PYB11cppname("polySecondMoment")
-def polySecondMoment{ndim}d(poly = "const Dim<{ndim}>::FacetedVolume&",
-                            center = "const Dim<{ndim}>::Vector&"):
+def polySecondMoment{ndim}d(poly = "const {Dimension}::FacetedVolume&",
+                            center = "const {Dimension}::Vector&"):
     "Return the second moment of a convex polytope"
-    return "Dim<{ndim}>::SymTensor"
+    return "{Dimension}::SymTensor"
 ''')
 
 

src/SmoothingScale/ASPHRadialFunctor.hh

@@ -0,0 +1,39 @@
+//---------------------------------Spheral++----------------------------------//
+// ASPHRadialFunctor
+//
+// Provides user-overridable hooks to modify how the ASPH object computes
+// radial normals and magnitude
+//
+// Created by JMO, Thu Oct 10 14:12:37 PDT 2024
+//----------------------------------------------------------------------------//
+#ifndef __Spheral_ASPHRadialFunctor__
+#define __Spheral_ASPHRadialFunctor__
+
+namespace Spheral {
+
+template<typename Dimension>
+class ASPHRadialFunctor {
+
+public:
+  //--------------------------- Public Interface ---------------------------//
+  using Scalar = typename Dimension::Scalar;
+  using Vector = typename Dimension::Vector;
+
+  // Constructors, destructor.
+  ASPHRadialFunctor()           {}
+  virtual ~ASPHRadialFunctor()  {}
+
+  // Compute the outward pointing radial unit vector
+  virtual Vector radialUnitVector(const size_t nodeListi,
+                                  const size_t i,
+                                  const Vector& posi) const { return posi.unitVector(); }
+
+  // Compute the radial coordinate
+  virtual Scalar radialCoordinate(const size_t nodeListi,
+                                  const size_t i,
+                                  const Vector& posi) const { return posi.magnitude(); }
+};
+
+}
+
+#endif

src/SmoothingScale/ASPHSmoothingScale.cc

@@ -171,6 +171,7 @@ ASPHSmoothingScale(const HEvolutionType HUpdate,
   mCellSecondMoment(FieldStorageType::CopyFields),
   mRadius0(FieldStorageType::CopyFields),
   mHidealFilterPtr(std::make_shared<ASPHSmoothingScaleUserFilter<Dimension>>()),
+  mRadialFunctorPtr(std::make_shared<ASPHRadialFunctor<Dimension>>()),
   mFixShape(fixShape),
   mRadialOnly(radialOnly) {
 }
@@ -203,26 +204,18 @@ registerState(DataBase<Dimension>& dataBase,
 
   const auto Hupdate = this->HEvolution();
   auto Hfields = dataBase.fluidHfield();
-  const auto numFields = Hfields.numFields();
-  for (auto k = 0u; k < numFields; ++k) {
-    auto& Hfield = *Hfields[k];
-    const auto& nodeList = Hfield.nodeList();
-    const auto hmin = nodeList.hmin();
-    const auto hmax = nodeList.hmax();
-    const auto hminratio = nodeList.hminratio();
-    switch (Hupdate) {
-      case HEvolutionType::IntegrateH:
-      case HEvolutionType::IdealH:
-        state.enroll(Hfield, make_policy<IncrementASPHHtensor<Dimension>>(hmin, hmax, hminratio, mFixShape, mRadialOnly));
-        break;
-
-      case HEvolutionType::FixedH:
-        state.enroll(Hfield);
-        break;
-
-       default:
-         VERIFY2(false, "ASPHSmoothingScale ERROR: Unknown Hevolution option ");
-    }
+  switch (Hupdate) {
+  case HEvolutionType::IntegrateH:
+  case HEvolutionType::IdealH:
+    state.enroll(Hfields, make_policy<IncrementASPHHtensor<Dimension>>(mFixShape, mRadialOnly, mRadialFunctorPtr));
+    break;
+
+  case HEvolutionType::FixedH:
+    state.enroll(Hfields);
+    break;
+
+  default:
+    VERIFY2(false, "ASPHSmoothingScale ERROR: Unknown Hevolution option ");
   }
 }
 
@@ -301,7 +294,7 @@ preStepInitialize(const DataBase<Dimension>& dataBase,
       const auto n = pos[k]->numInternalElements();
 #pragma omp parallel for
       for (auto i = 0u; i < n; ++i) {
-        mRadius0(k,i) = pos(k,i).magnitude();
+        mRadius0(k,i) = mRadialFunctorPtr->radialCoordinate(k, i, pos(k,i));
       }
     }
   }
@@ -604,8 +597,9 @@ finalize(const Scalar time,
 
           // We scale H in the radial direction only (also force H to be aligned radially).
           CHECK(mRadialOnly);
-          const auto nhat = pos(k, i).unitVector();
-          Hideali = radialEvolution(Hi, nhat, 1.0 - a + a*s, mRadius0(k,i), pos(k,i).magnitude());
+          const auto nhat = mRadialFunctorPtr->radialUnitVector(k, i, pos(k,i));
+          const auto r1 = mRadialFunctorPtr->radialCoordinate(k, i, pos(k,i));
+          Hideali = radialEvolution(Hi, nhat, 1.0 - a + a*s, mRadius0(k,i), r1);
 
         }
 

src/SmoothingScale/ASPHSmoothingScale.hh

@@ -10,8 +10,12 @@
 
 #include "SmoothingScale/SmoothingScaleBase.hh"
 #include "SmoothingScale/ASPHSmoothingScaleUserFilter.hh"
+#include "SmoothingScale/ASPHRadialFunctor.hh"
+#include "Utilities/Functors.hh"
 
-#include <memory>
+#include <memory>   // std::shared_ptr
+#include <utility>  // std::pair
+#include <tuple>
 
 namespace Spheral {
 
@@ -26,6 +30,7 @@ public:
   using SymTensor = typename Dimension::SymTensor;
   using FacetedVolume = typename Dimension::FacetedVolume;
   using HidealFilterType = ASPHSmoothingScaleUserFilter<Dimension>;
+  using RadialFunctorType = ASPHRadialFunctor<Dimension>;
 
   // Constructors, destructor.
   ASPHSmoothingScale(const HEvolutionType HUpdate,
@@ -94,10 +99,14 @@ public:
   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
+  // Optional user functor to manipulate the final ideal H vote
   std::shared_ptr<HidealFilterType> HidealFilter() const                          { return mHidealFilterPtr; }
   void HidealFilter(std::shared_ptr<HidealFilterType> functorPtr)                 { mHidealFilterPtr = functorPtr; }
 
+  // Optional user functor to override the radial unit normal and radius for radialOnly mode
+  std::shared_ptr<RadialFunctorType> RadialFunctor() const                        { return mRadialFunctorPtr; }
+  void RadialFunctor(std::shared_ptr<RadialFunctorType> functorPtr)               { mRadialFunctorPtr = functorPtr; }
+
   //****************************************************************************
   // Methods required for restarting.
   virtual std::string label() const override { return "ASPHSmoothingScale"; }
@@ -110,6 +119,7 @@ private:
   FieldList<Dimension, SymTensor> mSecondMoment, mCellSecondMoment;
   FieldList<Dimension, Scalar> mRadius0;
   std::shared_ptr<HidealFilterType> mHidealFilterPtr;
+  std::shared_ptr<RadialFunctorType> mRadialFunctorPtr;
   bool mFixShape, mRadialOnly;
 };
 

src/SmoothingScale/CMakeLists.txt

@@ -18,7 +18,9 @@ set(SmoothingScale_headers
   FixedSmoothingScale.hh
   SPHSmoothingScale.hh
   ASPHSmoothingScale.hh
+  IncrementASPHHtensor.hh
   ASPHSmoothingScaleUserFilter.hh
+  ASPHRadialFunctor.hh
   polySecondMoment.hh
 )