diff --git a/src/PYB11/SmoothingScale/ASPHSmoothingScale.py b/src/PYB11/SmoothingScale/ASPHSmoothingScale.py index be9ae7d38..aa3568189 100644 --- a/src/PYB11/SmoothingScale/ASPHSmoothingScale.py +++ b/src/PYB11/SmoothingScale/ASPHSmoothingScale.py @@ -62,6 +62,14 @@ def evaluateDerivatives(self, "Increment the derivatives." return "void" + @PYB11virtual + def preStepInitialize(self, + dataBase = "const DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&", + derivs = "StateDerivatives<%(Dimension)s>&"): + "Optional hook to be called at the beginning of a time step." + return "void" + @PYB11virtual def finalize(self, time = "const Scalar", @@ -96,6 +104,7 @@ def label(self): zerothMoment = PYB11property("const FieldList<%(Dimension)s, Scalar>&", "zerothMoment", doc="The zeroth moment storage FieldList") 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") 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/SmoothingScale/ASPHSmoothingScale.cc b/src/SmoothingScale/ASPHSmoothingScale.cc index 2b140d9c1..4add2d84b 100644 --- a/src/SmoothingScale/ASPHSmoothingScale.cc +++ b/src/SmoothingScale/ASPHSmoothingScale.cc @@ -98,7 +98,9 @@ inline Dim<1>::SymTensor radialEvolution(const Dim<1>::SymTensor& Hi, const Dim<1>::Vector& nhat, - const Dim<1>::Scalar s) { + const Dim<1>::Scalar s, + const Dim<1>::Scalar r0, + const Dim<1>::Scalar r1) { return Hi / s; } @@ -107,7 +109,9 @@ inline Dim<2>::SymTensor radialEvolution(const Dim<2>::SymTensor& Hi, const Dim<2>::Vector& nhat, - const Dim<2>::Scalar s) { + const Dim<2>::Scalar s, + const Dim<2>::Scalar r0, + const Dim<2>::Scalar r1) { const auto T = rotationMatrix(nhat).Transpose(); const auto hev0 = Hi.eigenVectors(); Dim<2>::SymTensor result; @@ -118,7 +122,10 @@ radialEvolution(const Dim<2>::SymTensor& Hi, result(0,0) = hev0.eigenValues(1); result(1,1) = hev0.eigenValues(0); } + const auto fr = r1*safeInvVar(r0); + CHECK(fr > 0.0); result(0,0) /= s; + result(1,1) /= fr; result.rotationalTransform(T); return result; } @@ -128,12 +135,20 @@ inline Dim<3>::SymTensor radialEvolution(const Dim<3>::SymTensor& Hi, const Dim<3>::Vector& nhat, - const Dim<3>::Scalar s) { + const Dim<3>::Scalar s, + const Dim<3>::Scalar r0, + const Dim<3>::Scalar r1) { const auto Tprinciple = rotationMatrix(nhat); const auto Tlab = Tprinciple.Transpose(); auto result = Hi; result.rotationalTransform(Tprinciple); + const auto fr = r1*safeInvVar(r0); + CHECK(fr > 0.0); result(0,0) /= s; + result(1,1) /= fr; + result(1,2) /= fr; + result(2,1) /= fr; + result(2,2) /= fr; result.rotationalTransform(Tlab); return result; } @@ -154,6 +169,7 @@ ASPHSmoothingScale(const HEvolutionType HUpdate, mZerothMoment(FieldStorageType::CopyFields), mSecondMoment(FieldStorageType::CopyFields), mCellSecondMoment(FieldStorageType::CopyFields), + mRadius0(FieldStorageType::CopyFields), mHidealFilterPtr(std::make_shared>()), mFixShape(fixShape), mRadialOnly(radialOnly) { @@ -171,6 +187,7 @@ initializeProblemStartup(DataBase& dataBase) { dataBase.resizeFluidFieldList(mZerothMoment, 0.0, HydroFieldNames::massZerothMoment, false); dataBase.resizeFluidFieldList(mSecondMoment, SymTensor::zero, HydroFieldNames::massSecondMoment, false); dataBase.resizeFluidFieldList(mCellSecondMoment, SymTensor::zero, HydroFieldNames::massSecondMoment + " cells", false); + if (mRadialOnly) dataBase.resizeFluidFieldList(mRadius0, 0.0, "Start of step radius", false); } //------------------------------------------------------------------------------ @@ -266,6 +283,30 @@ evaluateDerivatives(const typename Dimension::Scalar time, TIME_END("ASPHSmoothingScaleDerivs"); } +//------------------------------------------------------------------------------ +// Initialize at the beginning of a timestep. +//------------------------------------------------------------------------------ +template +void +ASPHSmoothingScale:: +preStepInitialize(const DataBase& dataBase, + State& state, + StateDerivatives& derivs) { + // If we're using the radial H scaling, take a snapshot of the initial radius of + // each point. + if (mRadialOnly) { + const auto pos = state.fields(HydroFieldNames::position, Vector::zero); + const auto numFields = pos.numFields(); + for (auto k = 0u; k < numFields; ++k) { + const auto n = pos[k]->numInternalElements(); +#pragma omp parallel for + for (auto i = 0u; i < n; ++i) { + mRadius0(k,i) = pos(k,i).magnitude(); + } + } + } +} + //------------------------------------------------------------------------------ // Finalize at the end of the step. // This is where we compute the Voronoi cell geometry and use it to set our @@ -564,7 +605,7 @@ 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); + Hideali = radialEvolution(Hi, nhat, 1.0 - a + a*s, mRadius0(k,i), pos(k,i).magnitude()); } diff --git a/src/SmoothingScale/ASPHSmoothingScale.hh b/src/SmoothingScale/ASPHSmoothingScale.hh index 3f9c24b10..c22a225ad 100644 --- a/src/SmoothingScale/ASPHSmoothingScale.hh +++ b/src/SmoothingScale/ASPHSmoothingScale.hh @@ -60,7 +60,12 @@ public: const State& state, StateDerivatives& derivatives) const override; - // Similarly packages might want a hook to do some post-step finalizations. + // Optional hook to be called at the beginning of a time step. + virtual void preStepInitialize(const DataBase& dataBase, + State& state, + StateDerivatives& derivs) override; + + // Similarly packages might want a hook to do some post-step finalizations. // Really we should rename this post-step finalize. virtual void finalize(const Scalar time, const Scalar dt, @@ -81,6 +86,7 @@ public: const FieldList& zerothMoment() const { return mZerothMoment; } const FieldList& secondMoment() const { return mSecondMoment; } const FieldList& cellSecondMoment() const { return mCellSecondMoment; } + const FieldList& radius0() const { return mRadius0; } // Special evolution flags bool fixShape() const { return mFixShape; } @@ -102,6 +108,7 @@ private: const TableKernel& mWT; FieldList mZerothMoment; FieldList mSecondMoment, mCellSecondMoment; + FieldList mRadius0; std::shared_ptr mHidealFilterPtr; bool mFixShape, mRadialOnly; }; diff --git a/src/Utilities/iterateIdealH.cc b/src/Utilities/iterateIdealH.cc index 6e7797e27..9b4a1a66c 100644 --- a/src/Utilities/iterateIdealH.cc +++ b/src/Utilities/iterateIdealH.cc @@ -13,6 +13,7 @@ #include "DataBase/IncrementBoundedState.hh" #include "DataBase/ReplaceBoundedState.hh" #include "Geometry/GeometryRegistrar.hh" +#include "SmoothingScale/ASPHSmoothingScale.hh" #include using std::vector; @@ -101,6 +102,19 @@ iterateIdealH(DataBase& dataBase, } } + // Check if we're using ASPH and radialOnly. If so we'll switch to fixShape for the iteration. + auto radialOnly = false; + ASPHSmoothingScale* asphPkg = nullptr; + for (auto* pkg: packages) { + asphPkg = dynamic_cast*>(pkg); + if (asphPkg != nullptr and asphPkg->radialOnly()) { + radialOnly = true; + asphPkg->radialOnly(false); + asphPkg->fixShape(true); + break; + } + } + // Build a list of flags to indicate which nodes have been completed. auto flagNodeDone = dataBase.newFluidFieldList(0, "node completed"); @@ -251,6 +265,13 @@ iterateIdealH(DataBase& dataBase, for (auto* boundaryPtr: range(boundaries.begin(), boundaries.end())) boundaryPtr->applyFieldListGhostBoundary(m); for (auto* boundaryPtr: range(boundaries.begin(), boundaries.end())) boundaryPtr->finalizeGhostBoundary(); + // Restore ASPH radialOnly choice if necessary + if (radialOnly) { + CHECK(asphPkg != nullptr); + asphPkg->radialOnly(true); + asphPkg->fixShape(false); + } + // Report the final timing. const auto t1 = clock(); if (Process::getRank() == 0 && maxIterations > 1) diff --git a/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py b/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py index e367c5653..4cbb98ea8 100644 --- a/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py +++ b/tests/functional/Hydro/Noh/Noh-cylindrical-2d.py @@ -78,14 +78,15 @@ # hydro type (only one!) svph = False, - crksph = False, # high order conservative formulation of SPH - psph = False, # pressure-based formulation of SPH - fsisph = False, # formulation for multimaterial problems - gsph = False, # godunov SPH - mfm = False, # moving finite mass of Hopkins 2015 - mfv=False, # moving finite volume of Hopkins 2015 - asph = False, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. - solid = False, # If true, use the fluid limit of the solid hydro option + crksph = False, # high order conservative formulation of SPH + psph = False, # pressure-based formulation of SPH + fsisph = False, # formulation for multimaterial problems + gsph = False, # godunov SPH + mfm = False, # moving finite mass of Hopkins 2015 + mfv=False, # moving finite volume of Hopkins 2015 + asph = False, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. + solid = False, # If true, use the fluid limit of the solid hydro option + radialOnly = False, # Force ASPH tensors to be aligned and evolve radially # general hydro options densityUpdate = RigorousSumDensity, # (IntegrateDensity) @@ -462,9 +463,13 @@ output("hydro.cfl") output("hydro.compatibleEnergyEvolution") output("hydro.densityUpdate") -#output("hydro._smoothingScaleMethod.HEvolution") +output("hydro._smoothingScaleMethod.HEvolution") if crksph: output("hydro.correctionOrder") +if radialOnly: + assert asph + hydro._smoothingScaleMethod.radialOnly = True + output("hydro._smoothingScaleMethod.radialOnly") packages = [hydro] diff --git a/tests/functional/Hydro/Noh/Noh-spherical-3d.py b/tests/functional/Hydro/Noh/Noh-spherical-3d.py index 43226af1a..66aa574fa 100644 --- a/tests/functional/Hydro/Noh/Noh-spherical-3d.py +++ b/tests/functional/Hydro/Noh/Noh-spherical-3d.py @@ -47,7 +47,7 @@ gamma = 5.0/3.0, mu = 1.0, - solid = False, # If true, use the fluid limit of the solid hydro option + solid = False, # If true, use the fluid limit of the solid hydro option svph = False, crksph = False, @@ -57,7 +57,8 @@ mfm = False, mfv = False, - asph = False, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. + asph = False, # This just chooses the H algorithm -- you can use this with CRKSPH for instance. + radialOnly = False, # Force ASPH tensors to be aligned and evolve radially boolReduceViscosity = False, HopkinsConductivity = False, # For PSPH nhQ = 5.0, @@ -379,11 +380,15 @@ output("hydro.kernel") output("hydro.cfl") output("hydro.compatibleEnergyEvolution") -output("hydro.HEvolution") if not (gsph or mfm or mfv or fsisph): output("hydro.PiKernel") if not fsisph: output("hydro.densityUpdate") +output("hydro._smoothingScaleMethod.HEvolution") +if radialOnly: + assert asph + hydro._smoothingScaleMethod.radialOnly = True + output("hydro._smoothingScaleMethod.radialOnly") packages = [hydro]