diff --git a/src/PYB11/CXXTypes/CXXTypes_PYB11.py b/src/PYB11/CXXTypes/CXXTypes_PYB11.py index 6f124fc2f..42695262c 100644 --- a/src/PYB11/CXXTypes/CXXTypes_PYB11.py +++ b/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) +''') diff --git a/src/PYB11/SmoothingScale/ASPHRadialFunctor.py b/src/PYB11/SmoothingScale/ASPHRadialFunctor.py new file mode 100644 index 000000000..83b0a56d5 --- /dev/null +++ b/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" diff --git a/src/PYB11/SmoothingScale/ASPHSmoothingScale.py b/src/PYB11/SmoothingScale/ASPHSmoothingScale.py index aa3568189..ddd16b416 100644 --- a/src/PYB11/SmoothingScale/ASPHSmoothingScale.py +++ b/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", "HidealFilter", "HidealFilter", doc="Optional function to manipulate the Hideal calculation") + HidealFilter = PYB11property("std::shared_ptr", "HidealFilter", "HidealFilter", doc="Optional functor to manipulate the Hideal calculation") + RadialFunctor = PYB11property("std::shared_ptr", "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") diff --git a/src/PYB11/SmoothingScale/ASPHSmoothingScaleUserFilter.py b/src/PYB11/SmoothingScale/ASPHSmoothingScaleUserFilter.py index 1b4bf070b..0e137b32e 100644 --- a/src/PYB11/SmoothingScale/ASPHSmoothingScaleUserFilter.py +++ b/src/PYB11/SmoothingScale/ASPHSmoothingScaleUserFilter.py @@ -1,5 +1,5 @@ #------------------------------------------------------------------------------- -# ASPHSmoothingScale +# ASPHSmoothingScaleUserFilter #------------------------------------------------------------------------------- from PYB11Generator import * diff --git a/src/PYB11/SmoothingScale/SmoothingScale_PYB11.py b/src/PYB11/SmoothingScale/SmoothingScale_PYB11.py index 4a2cbcdb2..cfdb66f1f 100644 --- a/src/PYB11/SmoothingScale/SmoothingScale_PYB11.py +++ b/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" ''') diff --git a/src/SmoothingScale/ASPHRadialFunctor.hh b/src/SmoothingScale/ASPHRadialFunctor.hh new file mode 100644 index 000000000..262b6b3ba --- /dev/null +++ b/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 +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 diff --git a/src/SmoothingScale/ASPHSmoothingScale.cc b/src/SmoothingScale/ASPHSmoothingScale.cc index 4add2d84b..2b433c985 100644 --- a/src/SmoothingScale/ASPHSmoothingScale.cc +++ b/src/SmoothingScale/ASPHSmoothingScale.cc @@ -171,6 +171,7 @@ ASPHSmoothingScale(const HEvolutionType HUpdate, mCellSecondMoment(FieldStorageType::CopyFields), mRadius0(FieldStorageType::CopyFields), mHidealFilterPtr(std::make_shared>()), + mRadialFunctorPtr(std::make_shared>()), mFixShape(fixShape), mRadialOnly(radialOnly) { } @@ -203,26 +204,18 @@ registerState(DataBase& 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>(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>(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& 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); } diff --git a/src/SmoothingScale/ASPHSmoothingScale.hh b/src/SmoothingScale/ASPHSmoothingScale.hh index c22a225ad..ebac70f00 100644 --- a/src/SmoothingScale/ASPHSmoothingScale.hh +++ b/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 +#include // std::shared_ptr +#include // std::pair +#include namespace Spheral { @@ -26,6 +30,7 @@ public: using SymTensor = typename Dimension::SymTensor; using FacetedVolume = typename Dimension::FacetedVolume; using HidealFilterType = ASPHSmoothingScaleUserFilter; + using RadialFunctorType = ASPHRadialFunctor; // 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 HidealFilter() const { return mHidealFilterPtr; } void HidealFilter(std::shared_ptr functorPtr) { mHidealFilterPtr = functorPtr; } + // Optional user functor to override the radial unit normal and radius for radialOnly mode + std::shared_ptr RadialFunctor() const { return mRadialFunctorPtr; } + void RadialFunctor(std::shared_ptr functorPtr) { mRadialFunctorPtr = functorPtr; } + //**************************************************************************** // Methods required for restarting. virtual std::string label() const override { return "ASPHSmoothingScale"; } @@ -110,6 +119,7 @@ private: FieldList mSecondMoment, mCellSecondMoment; FieldList mRadius0; std::shared_ptr mHidealFilterPtr; + std::shared_ptr mRadialFunctorPtr; bool mFixShape, mRadialOnly; }; diff --git a/src/SmoothingScale/CMakeLists.txt b/src/SmoothingScale/CMakeLists.txt index 43de653d5..26c422f20 100644 --- a/src/SmoothingScale/CMakeLists.txt +++ b/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 ) diff --git a/src/SmoothingScale/IncrementASPHHtensor.cc b/src/SmoothingScale/IncrementASPHHtensor.cc index 3b6bc6fa8..0c21cb5ca 100644 --- a/src/SmoothingScale/IncrementASPHHtensor.cc +++ b/src/SmoothingScale/IncrementASPHHtensor.cc @@ -1,15 +1,16 @@ //---------------------------------Spheral++----------------------------------// // IncrementASPHHtensor // -// Specialized version of FieldUpdatePolicy for time integrating the H tensor. +// Specialized version of UpdatePolicyBase 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 "Field/FieldList.hh" #include "Hydro/HydroFieldNames.hh" +#include "Geometry/Dimension.hh" #include "Utilities/rotationMatrix.hh" #include "Utilities/GeometricUtilities.hh" #include "Utilities/DBC.hh" @@ -21,17 +22,13 @@ namespace Spheral { //------------------------------------------------------------------------------ template IncrementASPHHtensor:: -IncrementASPHHtensor(const Scalar hmin, - const Scalar hmax, - const Scalar hminratio, - const bool fixShape, - const bool radialOnly): - FieldUpdatePolicy(), - mhmin(hmin), - mhmax(hmax), - mhminratio(hminratio), +IncrementASPHHtensor(const bool fixShape, + const bool radialOnly, + std::shared_ptr radialFunctorPtr): + UpdatePolicyBase(), mFixShape(fixShape), - mRadialOnly(radialOnly) { + mRadialOnly(radialOnly), + mRadialFunctorPtr(radialFunctorPtr) { } //------------------------------------------------------------------------------ @@ -50,51 +47,59 @@ update(const KeyType& key, // Get the field name portion of the key. KeyType fieldKey, nodeListKey; StateBase::splitFieldKey(key, fieldKey, nodeListKey); - CHECK(fieldKey == HydroFieldNames::H); - - const auto hminInv = 1.0/mhmin; - const auto hmaxInv = 1.0/mhmax; + REQUIRE(fieldKey == HydroFieldNames::H and + nodeListKey == UpdatePolicyBase::wildcard()); // Get the state we're updating. - auto& H = state.field(key, SymTensor::zero); - const auto& DHDt = derivs.field(prefix() + StateBase::buildFieldKey(HydroFieldNames::H, nodeListKey), SymTensor::zero); - const auto& pos = state.field(StateBase::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(); + auto H = state.fields(HydroFieldNames::H, SymTensor::zero); + const auto DHDt = derivs.fields(prefix() + HydroFieldNames::H, SymTensor::zero); + const auto pos = state.fields(HydroFieldNames::position, Vector::zero); // Only needed if we're using radial scaling + const auto numFields = H.numFields(); + CHECK(DHDt.numFields() == numFields); + CHECK(pos.numFields() == numFields); + + // Walk the NodeLists + for (auto k = 0u; k < numFields; ++k) { + const auto& nodeList = H[k]->nodeList(); + const auto hminInv = 1.0/nodeList.hmin(); + const auto hmaxInv = 1.0/nodeList.hmax(); + const auto hminratio = nodeList.hminratio(); + const auto n = nodeList.numInternalNodes(); + + // Walk the nodes and update H (with limiting) #pragma omp parallel for - for (auto i = 0u; i < n; ++i) { + for (auto i = 0u; i < n; ++i) { - // Check for special update rules - if (mFixShape) { + // Check for special update rules + if (mFixShape) { - // Fix the shape (only volume scaling allowed) + // Fix the shape (only volume scaling allowed) + auto fi = Dimension::rootnu((H(k,i) + multiplier*DHDt(k,i)).Determinant()/H(k,i).Determinant()); + H(k,i) *= fi; - auto fi = Dimension::rootnu((H(i) + multiplier*DHDt(i)).Determinant()/H(i).Determinant()); - H(i) *= fi; + } else if (mRadialOnly) { - } else if (mRadialOnly) { + // Force only the radial component of H to be scaled + const auto nhat = mRadialFunctorPtr->radialUnitVector(k, i, pos(k,i)); + const auto T = rotationMatrix(nhat); + H(k,i).rotationalTransform(T); // Should have one eigenvector aligned with the x' axis in this frame + auto DHDti = DHDt(k,i); + DHDti.rotationalTransform(T); + H(k,i)[0] += multiplier * DHDti[0]; + H(k,i).rotationalTransform(T.Transpose()); - // 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 { - } else { + H(k,i) += multiplier * DHDt(k,i); - H(i) += multiplier * DHDt(i); + } + // Apply limiting + const auto hev = H(k,i).eigenVectors(); + const auto hminEffInv = min(hminInv, max(hmaxInv, hev.eigenValues.minElement())/hminratio); + H(k,i) = constructSymTensorWithBoundedDiagonal(hev.eigenValues, hmaxInv, hminEffInv); + H(k,i).rotationalTransform(hev.eigenVectors); } - - // 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); } } @@ -110,10 +115,7 @@ operator==(const UpdatePolicyBase& rhs) const { 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 + return (fixShape() == rhsPtr->fixShape() and radialOnly() == rhsPtr->radialOnly()); } diff --git a/src/SmoothingScale/IncrementASPHHtensor.hh b/src/SmoothingScale/IncrementASPHHtensor.hh index abb34f5b6..6197b8352 100644 --- a/src/SmoothingScale/IncrementASPHHtensor.hh +++ b/src/SmoothingScale/IncrementASPHHtensor.hh @@ -1,14 +1,17 @@ //---------------------------------Spheral++----------------------------------// // IncrementASPHHtensor // -// Specialized version of FieldUpdatePolicy for time integrating the H tensor. +// Specialized version of UpdatePolicy 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" +#include "DataBase/UpdatePolicyBase.hh" +#include "SmoothingScale/ASPHRadialFunctor.hh" + +#include // std::shared_ptr namespace Spheral { @@ -16,21 +19,20 @@ namespace Spheral { template class StateDerivatives; template -class IncrementASPHHtensor: public FieldUpdatePolicy { +class IncrementASPHHtensor: public UpdatePolicyBase { public: //--------------------------- Public Interface ---------------------------// // Useful typedefs - using KeyType = typename FieldUpdatePolicy::KeyType; + using KeyType = typename UpdatePolicyBase::KeyType; using Scalar = typename Dimension::Scalar; using Vector = typename Dimension::Vector; using SymTensor = typename Dimension::SymTensor; + using RadialFunctorType = ASPHRadialFunctor; // Constructors, destructor. - IncrementASPHHtensor(const Scalar hmin, - const Scalar hmax, - const Scalar hminratio, - const bool fixShape, - const bool radialOnly); + IncrementASPHHtensor(const bool fixShape, + const bool radialOnly, + std::shared_ptr radialFunctorPtr); virtual ~IncrementASPHHtensor() {} IncrementASPHHtensor(const IncrementASPHHtensor& rhs) = delete; IncrementASPHHtensor& operator=(const IncrementASPHHtensor& rhs) = delete; @@ -43,10 +45,7 @@ public: 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; } + // Access the internal state bool fixShape() const { return mFixShape; } bool radialOnly() const { return mRadialOnly; } @@ -57,8 +56,8 @@ public: private: //--------------------------- Private Interface ---------------------------// - Scalar mhmin, mhmax, mhminratio; bool mFixShape, mRadialOnly; + std::shared_ptr mRadialFunctorPtr; }; }