Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add mean atomic mass and number #439

Open
wants to merge 11 commits into
base: main
Choose a base branch
from
Open

Add mean atomic mass and number #439

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
bb4b797
start adding mean atomic mass and number
jonahm-LANL Nov 29, 2024
3e9a1f2
Rela -> Real
jonahm-LANL Nov 30, 2024
f37996c
dont shadow T parameter
jonahm-LANL Nov 30, 2024
6b2c05c
Add MeanAtomicProperties for analytic EOS's to use and implement for …
jonahm-LANL Dec 2, 2024
2b7d0b5
test
jonahm-LANL Dec 2, 2024
ade7cc3
formatting. Add carnahan starling AZbar
jonahm-LANL Dec 2, 2024
5edeed4
eos davis
jonahm-LANL Dec 2, 2024
91865e5
Merge branch 'main' into jmm/mean-atomic-mass
jonahm-LANL Dec 3, 2024
388a596
Move default implementations of MeanAtomicMass/MeanAtomicNumber to macro
jonahm-LANL Dec 3, 2024
44b9fc8
spiner already had Abar and Zbar in metadata. Just add it to class.
jonahm-LANL Dec 3, 2024
738b6f5
temperature -> T
jonahm-LANL Dec 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 65 additions & 1 deletion singularity-eos/eos/eos_base.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,16 @@ char *StrCat(char *destination, const char *source) {
using EosBase<EOSDERIVED>::GibbsFreeEnergyFromDensityTemperature; \
using EosBase<EOSDERIVED>::GibbsFreeEnergyFromDensityInternalEnergy;

// This macro adds these methods to a derived class. Due to scope,
// these can't be implemented in the base class, unless we make
// _AZbar public. Not all EOS's may want these default functions
// TODO(JMM): Should we go the alternate route and make _AZbar public?
#define SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) \
PORTABLE_INLINE_FUNCTION \
Real MeanAtomicMass() const { return _AZbar.Abar; } \
PORTABLE_INLINE_FUNCTION \
Real MeanAtomicNumber() const { return _AZbar.Zbar; }

// This macro adds several methods that most modifiers will
// want. Not ALL modifiers will want these methods as written here,
// so use this macro with care.
Expand All @@ -101,7 +111,11 @@ char *StrCat(char *destination, const char *source) {
} \
constexpr bool AllDynamicMemoryIsShareable() const { \
return t_.AllDynamicMemoryIsShareable(); \
}
} \
PORTABLE_INLINE_FUNCTION \
Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } \
PORTABLE_INLINE_FUNCTION \
Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); }

class Factor {
Real value_ = 1.0;
Expand Down Expand Up @@ -132,6 +146,33 @@ struct Transform {
Factor x, y, f;
};

/*
This is a utility struct used to bundle mean atomic
mass/number. Used in the default implementations of MeanAtomicMass
and MeanAtomicNumber provided by the base class.
*/
struct MeanAtomicProperties {
Real Abar, Zbar;

// default is hydrogen
static constexpr Real DEFAULT_ABAR = 1.0;
static constexpr Real DEFAULT_ZBAR = 1.0;

PORTABLE_INLINE_FUNCTION
MeanAtomicProperties(Real Abar_, Real Zbar_) : Abar(Abar_), Zbar(Zbar_) {}
PORTABLE_INLINE_FUNCTION
MeanAtomicProperties() : Abar(DEFAULT_ABAR), Zbar(DEFAULT_ZBAR) {}
PORTABLE_INLINE_FUNCTION
void CheckParams() const {
PORTABLE_ALWAYS_REQUIRE(Abar > 0, "Positive mean atomic mass");
PORTABLE_ALWAYS_REQUIRE(Zbar > 0, "Positive mean atomic number");
}
void PrintParams() const {
printf(" Abar = %g\n", Abar);
printf(" Zbar = %g\n", Zbar);
}
};

/*
This is a CRTP that allows for static inheritance so that default behavior for
various member functions can be defined.
Expand Down Expand Up @@ -739,6 +780,29 @@ class EosBase {
PORTABLE_INLINE_FUNCTION
Real RhoPmin(const Real temp) const { return 0.0; }

// JMM: EOS's which encapsulate a mix or reactions may wish to vary
// this. For example, Helmholtz and StellarCollapse. This isn't the
// default, so by default the base class provides a specialization.
// for models where density and temperature are required, the EOS
// developer is in charge of either throwing an error or choosing
// reasonable defaults.
// TODO(JMM): Should we provide vector implementations if we depend
// on rho, T, etc?
template <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
const Real rho, const Real T,
Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
CRTP copy = *(static_cast<CRTP const *>(this));
return copy.MeanAtomicMass();
}
template <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
const Real rho, const Real T,
Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
CRTP copy = *(static_cast<CRTP const *>(this));
return copy.MeanAtomicNumber();
}

// Default entropy behavior is to cause an error
PORTABLE_FORCEINLINE_FUNCTION
void EntropyIsNotEnabled(const char *eosname) const {
Expand Down
18 changes: 13 additions & 5 deletions singularity-eos/eos/eos_carnahan_starling.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,26 +38,29 @@ using namespace eos_base;
class CarnahanStarling : public EosBase<CarnahanStarling> {
public:
CarnahanStarling() = default;
PORTABLE_INLINE_FUNCTION CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq)
PORTABLE_INLINE_FUNCTION
CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq,
const MeanAtomicProperties &AZbar = MeanAtomicProperties())
: _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(ROOM_TEMPERATURE),
_P0(ATMOSPHERIC_PRESSURE), _qp(0.0),
_rho0(DensityFromPressureTemperature(_P0, _T0)), _vol0(robust::ratio(1.0, _rho0)),
_sie0(Cv * _T0 + qq),
_bmod0(_rho0 * Cv * _T0 *
(PartialRhoZedFromDensity(_rho0) +
ZedFromDensity(_rho0) * ZedFromDensity(_rho0) * gm1)),
_dpde0(_rho0 * ZedFromDensity(_rho0) * gm1) {
_dpde0(_rho0 * ZedFromDensity(_rho0) * gm1), _AZbar(AZbar) {
CheckParams();
}
PORTABLE_INLINE_FUNCTION CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq, Real qp,
Real T0, Real P0)
PORTABLE_INLINE_FUNCTION
CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq, Real qp, Real T0, Real P0,
const MeanAtomicProperties &AZbar = MeanAtomicProperties())
: _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(T0), _P0(P0), _qp(qp),
_rho0(DensityFromPressureTemperature(P0, T0)), _vol0(robust::ratio(1.0, _rho0)),
_sie0(Cv * T0 + qq),
_bmod0(_rho0 * Cv * T0 *
(PartialRhoZedFromDensity(_rho0) +
ZedFromDensity(_rho0) * ZedFromDensity(_rho0) * gm1)),
_dpde0(_rho0 * ZedFromDensity(_rho0) * gm1) {
_dpde0(_rho0 * ZedFromDensity(_rho0) * gm1), _AZbar(AZbar) {
CheckParams();
}
CarnahanStarling GetOnDevice() { return *this; }
Expand All @@ -71,6 +74,7 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive");
PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive");
PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive");
_AZbar.CheckParams();
}
template <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION Real ZedFromDensity(
Expand Down Expand Up @@ -221,6 +225,8 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
bmod = _bmod0;
dpde = _dpde0;
}
// Default accessors for mean atomic mass/number
SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
// Generic functions provided by the base class. These contain e.g. the vector
// overloads that use the scalar versions declared here
SG_ADD_BASE_CLASS_USINGS(CarnahanStarling)
Expand All @@ -235,6 +241,7 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
printf("Carnahan-Starling Parameters:\nGamma = %g\nCv = %g\nb = %g\nq = "
"%g\n",
_gm1 + 1.0, _Cv, _bb, _qq);
_AZbar.PrintParams();
}
template <typename Indexer_t>
PORTABLE_INLINE_FUNCTION void
Expand All @@ -253,6 +260,7 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
Real _T0, _P0, _qp;
// reference values
Real _rho0, _vol0, _sie0, _bmod0, _dpde0;
MeanAtomicProperties _AZbar;
// static constexpr const Real _T0 = ROOM_TEMPERATURE;
// static constexpr const Real _P0 = ATMOSPHERIC_PRESSURE;
static constexpr const unsigned long _preferred_input =
Expand Down
20 changes: 16 additions & 4 deletions singularity-eos/eos/eos_davis.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,16 +35,18 @@ class DavisReactants : public EosBase<DavisReactants> {
PORTABLE_INLINE_FUNCTION
DavisReactants(const Real rho0, const Real e0, const Real P0, const Real T0,
const Real A, const Real B, const Real C, const Real G0, const Real Z,
const Real alpha, const Real Cv0)
const Real alpha, const Real Cv0,
const MeanAtomicProperties &AZbar = MeanAtomicProperties())
: _rho0(rho0), _e0(e0), _P0(P0), _T0(T0), _A(A), _B(B), _C(C), _G0(G0), _Z(Z),
_alpha(alpha), _Cv0(Cv0) {
_alpha(alpha), _Cv0(Cv0), _AZbar(AZbar) {
CheckParams();
}
DavisReactants GetOnDevice() { return *this; }
PORTABLE_INLINE_FUNCTION
void CheckParams() const {
PORTABLE_REQUIRE(_rho0 >= 0, "Density must be positive");
PORTABLE_REQUIRE(_T0 >= 0, "Temperature must be positive");
_AZbar.CheckParams();
}
template <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
Expand Down Expand Up @@ -160,6 +162,8 @@ class DavisReactants : public EosBase<DavisReactants> {
PORTABLE_INLINE_FUNCTION void
DensityEnergyFromPressureTemperature(const Real press, const Real temp,
Indexer_t &&lambda, Real &rho, Real &sie) const;
// Default accessors for mean atomic mass/number
SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
// Generic functions provided by the base class. These contain e.g. the vector
// overloads that use the scalar versions declared here
SG_ADD_BASE_CLASS_USINGS(DavisReactants)
Expand All @@ -175,6 +179,7 @@ class DavisReactants : public EosBase<DavisReactants> {
printf("%srho0:%e e0:%e P0:%e\nT0:%e A:%e B:%e\nC:%e G0:%e Z:%e\nalpha:%e "
"Cv0:%e\n",
s1, _rho0, _e0, _P0, _T0, _A, _B, _C, _G0, _Z, _alpha, _Cv0);
_AZbar.PrintParams();
}
void inline Finalize() {}
static std::string EosType() { return std::string("DavisReactants"); }
Expand All @@ -183,6 +188,7 @@ class DavisReactants : public EosBase<DavisReactants> {
private:
static constexpr Real onethird = 1.0 / 3.0;
Real _rho0, _e0, _P0, _T0, _A, _B, _C, _G0, _Z, _alpha, _Cv0;
MeanAtomicProperties _AZbar;
// static constexpr const char _eos_type[] = "DavisReactants";
static constexpr unsigned long _preferred_input =
thermalqs::density | thermalqs::specific_internal_energy;
Expand All @@ -198,10 +204,12 @@ class DavisProducts : public EosBase<DavisProducts> {
DavisProducts() = default;
PORTABLE_INLINE_FUNCTION
DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc,
const Real pc, const Real Cv)
: _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv) {}
const Real pc, const Real Cv,
const MeanAtomicProperties &AZbar = MeanAtomicProperties())
: _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv), _AZbar(AZbar) {}
PORTABLE_INLINE_FUNCTION
void CheckParams() const {
_AZbar.CheckParams();
// TODO(JMM): Stub.
}
DavisProducts GetOnDevice() { return *this; }
Expand Down Expand Up @@ -299,6 +307,8 @@ class DavisProducts : public EosBase<DavisProducts> {
ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
Real &bmod, Real &dpde, Real &dvdt,
Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
// Default accessors for mean atomic mass/number
SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
// Generic functions provided by the base class. These contain e.g. the vector
// overloads that use the scalar versions declared here
SG_ADD_BASE_CLASS_USINGS(DavisProducts)
Expand All @@ -313,6 +323,7 @@ class DavisProducts : public EosBase<DavisProducts> {
static constexpr char s1[]{"DavisProducts Params: "};
printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e \n", s1, _a, _b, _k, _n, _vc, _pc,
_Cv);
_AZbar.PrintParams();
}
inline void Finalize() {}
static std::string EosType() { return std::string("DavisProducts"); }
Expand All @@ -321,6 +332,7 @@ class DavisProducts : public EosBase<DavisProducts> {
private:
static constexpr Real onethird = 1.0 / 3.0;
Real _a, _b, _k, _n, _vc, _pc, _Cv;
MeanAtomicProperties _AZbar;
// static constexpr const char _eos_type[] = "DavisProducts";
static constexpr const unsigned long _preferred_input =
thermalqs::density | thermalqs::specific_internal_energy;
Expand Down
8 changes: 8 additions & 0 deletions singularity-eos/eos/eos_eospac.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,7 @@ class EOSPAC : public EosBase<EOSPAC> {
// TODO(JMM): More validation checks?
PORTABLE_ALWAYS_REQUIRE(rho_min_ >= 0, "Non-negative minimum density");
PORTABLE_ALWAYS_REQUIRE(temp_min_ >= 0, "Non-negative minimum temperature");
AZbar_.CheckParams();
}
inline EOSPAC GetOnDevice() { return *this; }

Expand Down Expand Up @@ -457,6 +458,8 @@ class EOSPAC : public EosBase<EOSPAC> {
using EosBase<EOSPAC>::FillEos;
using EosBase<EOSPAC>::EntropyIsNotEnabled;

SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_)

// EOSPAC vector implementations
template <typename LambdaIndexer>
inline void
Expand Down Expand Up @@ -1097,6 +1100,7 @@ class EOSPAC : public EosBase<EOSPAC> {
static std::string EosPyType() { return EosType(); }
PORTABLE_INLINE_FUNCTION void PrintParams() const {
printf("EOSPAC parameters:\nmatid = %i\n", matid_);
_AZbar.PrintParams();
}
PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rho_min_; }
PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return temp_min_; }
Expand Down Expand Up @@ -1127,6 +1131,7 @@ class EOSPAC : public EosBase<EOSPAC> {
// TODO(JMM): Is the fact that EOS_INTEGER isn't a size_t a
// problem? Could it ever realistically overflow?
EOS_INTEGER shared_size_, packed_size_;
MeanAtomicProperties AZbar_;

static inline std::map<std::string, unsigned int> &scratch_nbuffers() {
static std::map<std::string, unsigned int> nbuffers = {
Expand Down Expand Up @@ -1194,6 +1199,9 @@ inline EOSPAC::EOSPAC(const int matid, bool invert_at_setup, Real insert_data,
rho_ref_ = m.normalDensity;
rho_min_ = m.rhoMin;
temp_min_ = m.TMin;
// use std::max to hydrogen, in case of bad table
AZbar_.Abar = std::max(1.0, m.meanAtomicMass);
AZbar_.Zbar = std::max(1.0, m.meanAtomicNumber);

EOS_REAL R[1] = {rho_ref_};
EOS_REAL T[1] = {temperatureToSesame(temp_ref_)};
Expand Down
15 changes: 11 additions & 4 deletions singularity-eos/eos/eos_gruneisen.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,9 +39,10 @@ class Gruneisen : public EosBase<Gruneisen> {
PORTABLE_INLINE_FUNCTION
Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0,
const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv,
const Real rho_max)
const Real rho_max,
const MeanAtomicProperties &AZbar = MeanAtomicProperties())
: _C0(C0), _s1(s1), _s2(s2), _s3(s3), _G0(G0), _b(b), _rho0(rho0), _T0(T0), _P0(P0),
_Cv(Cv), _rho_max(rho_max) {
_Cv(Cv), _rho_max(rho_max), _AZbar(AZbar) {
// Warn user when provided rho_max is greater than the computed rho_max
#ifndef NDEBUG
const Real computed_rho_max = ComputeRhoMax(s1, s2, s3, rho0);
Expand All @@ -58,9 +59,11 @@ class Gruneisen : public EosBase<Gruneisen> {
// Constructor when rho_max isn't specified automatically determines _rho_max
PORTABLE_INLINE_FUNCTION
Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0,
const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv)
const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv,
const MeanAtomicProperties &AZbar = MeanAtomicProperties())
: _C0(C0), _s1(s1), _s2(s2), _s3(s3), _G0(G0), _b(b), _rho0(rho0), _T0(T0), _P0(P0),
_Cv(Cv), _rho_max(RHOMAX_SAFETY * ComputeRhoMax(s1, s2, s3, rho0)) {}
_Cv(Cv), _rho_max(RHOMAX_SAFETY * ComputeRhoMax(s1, s2, s3, rho0)),
_AZbar(AZbar) {}
static PORTABLE_INLINE_FUNCTION Real ComputeRhoMax(const Real s1, const Real s2,
const Real s3, const Real rho0);
PORTABLE_INLINE_FUNCTION
Expand All @@ -71,6 +74,7 @@ class Gruneisen : public EosBase<Gruneisen> {
PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Non-negative heat capacity required");
PORTABLE_ALWAYS_REQUIRE(_rho_max > _rho0,
"Maximum density must be greater than reference");
_AZbar.CheckParams();
}
PORTABLE_INLINE_FUNCTION Real
MaxStableDensityAtTemperature(const Real temperature) const;
Expand Down Expand Up @@ -150,6 +154,7 @@ class Gruneisen : public EosBase<Gruneisen> {
Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
// Generic functions provided by the base class. These contain e.g. the vector
// overloads that use the scalar versions declared here
SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar)
SG_ADD_BASE_CLASS_USINGS(Gruneisen)
PORTABLE_INLINE_FUNCTION
int nlambda() const noexcept { return 0; }
Expand All @@ -163,6 +168,7 @@ class Gruneisen : public EosBase<Gruneisen> {
printf("%s C0:%e s1:%e s2:%e s3:%e\n G0:%e b:%e rho0:%e T0:%e\n P0:%eCv:%e "
"rho_max:%e\n",
s1, _C0, _s1, _s2, _s3, _G0, _b, _rho0, _T0, _P0, _Cv, _rho_max);
_AZbar.PrintParams();
}
template <typename Indexer_t>
PORTABLE_INLINE_FUNCTION void
Expand All @@ -174,6 +180,7 @@ class Gruneisen : public EosBase<Gruneisen> {

private:
Real _C0, _s1, _s2, _s3, _G0, _b, _rho0, _T0, _P0, _Cv, _rho_max;
MeanAtomicProperties _AZbar;
// static constexpr const char _eos_type[] = {"Gruneisen"};
PORTABLE_INLINE_FUNCTION
Real Gamma(const Real rho_in) const {
Expand Down
25 changes: 25 additions & 0 deletions singularity-eos/eos/eos_helmholtz.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -640,6 +640,31 @@ class Helmholtz : public EosBase<Helmholtz> {
return gamma3 - 1.0;
}

PORTABLE_INLINE_FUNCTION
Real MeanAtomicMass() const {
PORTABLE_THROW_OR_ABORT("For Helmholtz EOS, mean atomic mass is an input!\n");
return 1.0;
}
PORTABLE_INLINE_FUNCTION
Real MeanAtomicNumber() const {
PORTABLE_THROW_OR_ABORT("For Helmholtz EOS, mean atomic number is an input!\n");
return 1.0;
}
template <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature(
const Real rho, const Real T,
Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
using namespace HelmUtils;
return lambda[Lambda::Abar];
}
template <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
const Real rho, const Real T,
Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
using namespace HelmUtils;
return lambda[Lambda::Zbar];
}

template <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION void
FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
Expand Down
Loading
Loading