Skip to content

Commit

Permalink
Improving radialOnly option for ASPH. Added expected analytic evoluti…
Browse files Browse the repository at this point in the history 
…on in the

non-radial direction to account for changing r, and forced special
mode for iterateH which switches to fixShape for the iteration.
  • Loading branch information
@jmikeowen
jmikeowen committed Oct 9, 2024
1 parent 19022d1 commit 131410a
Show file tree
Hide file tree
Showing 6 changed files with 105 additions and 17 deletions.
9 changes: 9 additions & 0 deletions src/PYB11/SmoothingScale/ASPHSmoothingScale.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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<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")
49 changes: 45 additions & 4 deletions src/SmoothingScale/ASPHSmoothingScale.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}

Expand All @@ -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;
Expand All @@ -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;
}
Expand All @@ -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;
}
Expand All @@ -154,6 +169,7 @@ ASPHSmoothingScale(const HEvolutionType HUpdate,
mZerothMoment(FieldStorageType::CopyFields),
mSecondMoment(FieldStorageType::CopyFields),
mCellSecondMoment(FieldStorageType::CopyFields),
mRadius0(FieldStorageType::CopyFields),
mHidealFilterPtr(std::make_shared<ASPHSmoothingScaleUserFilter<Dimension>>()),
mFixShape(fixShape),
mRadialOnly(radialOnly) {
Expand All @@ -171,6 +187,7 @@ initializeProblemStartup(DataBase<Dimension>& 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);
}

//------------------------------------------------------------------------------
Expand Down Expand Up @@ -266,6 +283,30 @@ evaluateDerivatives(const typename Dimension::Scalar time,
TIME_END("ASPHSmoothingScaleDerivs");
}

//------------------------------------------------------------------------------
// Initialize at the beginning of a timestep.
//------------------------------------------------------------------------------
template<typename Dimension>
void
ASPHSmoothingScale<Dimension>::
preStepInitialize(const DataBase<Dimension>& dataBase,
State<Dimension>& state,
StateDerivatives<Dimension>& 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
Expand Down Expand Up @@ -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());

}

Expand Down
9 changes: 8 additions & 1 deletion src/SmoothingScale/ASPHSmoothingScale.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,12 @@ public:
const State<Dimension>& state,
StateDerivatives<Dimension>& 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<Dimension>& dataBase,
State<Dimension>& state,
StateDerivatives<Dimension>& 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,
Expand All @@ -81,6 +86,7 @@ public:
const FieldList<Dimension, Scalar>& zerothMoment() const { return mZerothMoment; }
const FieldList<Dimension, SymTensor>& secondMoment() const { return mSecondMoment; }
const FieldList<Dimension, SymTensor>& cellSecondMoment() const { return mCellSecondMoment; }
const FieldList<Dimension, Scalar>& radius0() const { return mRadius0; }

// Special evolution flags
bool fixShape() const { return mFixShape; }
Expand All @@ -102,6 +108,7 @@ private:
const TableKernel<Dimension>& mWT;
FieldList<Dimension, Scalar> mZerothMoment;
FieldList<Dimension, SymTensor> mSecondMoment, mCellSecondMoment;
FieldList<Dimension, Scalar> mRadius0;
std::shared_ptr<HidealFilterType> mHidealFilterPtr;
bool mFixShape, mRadialOnly;
};
Expand Down
21 changes: 21 additions & 0 deletions src/Utilities/iterateIdealH.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include "DataBase/IncrementBoundedState.hh"
#include "DataBase/ReplaceBoundedState.hh"
#include "Geometry/GeometryRegistrar.hh"
#include "SmoothingScale/ASPHSmoothingScale.hh"

#include <ctime>
using std::vector;
Expand Down Expand Up @@ -101,6 +102,19 @@ iterateIdealH(DataBase<Dimension>& dataBase,
}
}

// Check if we're using ASPH and radialOnly. If so we'll switch to fixShape for the iteration.
auto radialOnly = false;
ASPHSmoothingScale<Dimension>* asphPkg = nullptr;
for (auto* pkg: packages) {
asphPkg = dynamic_cast<ASPHSmoothingScale<Dimension>*>(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");

Expand Down Expand Up @@ -251,6 +265,13 @@ iterateIdealH(DataBase<Dimension>& 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)
Expand Down
23 changes: 14 additions & 9 deletions tests/functional/Hydro/Noh/Noh-cylindrical-2d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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]

Expand Down
11 changes: 8 additions & 3 deletions tests/functional/Hydro/Noh/Noh-spherical-3d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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,
Expand Down Expand Up @@ -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]

Expand Down

0 comments on commit 131410a

Please sign in to comment.