Move default implementations of MeanAtomicMass/MeanAtomicNumber to macro

singularity-eos/eos/eos_base.hpp

@@ -82,9 +82,17 @@ char *StrCat(char *destination, const char *source) {
   using EosBase<EOSDERIVED>::EntropyFromDensityTemperature;                              \
   using EosBase<EOSDERIVED>::EntropyFromDensityInternalEnergy;                           \
   using EosBase<EOSDERIVED>::GibbsFreeEnergyFromDensityTemperature;                      \
-  using EosBase<EOSDERIVED>::GibbsFreeEnergyFromDensityInternalEnergy;                   \
-  using EosBase<EOSDERIVED>::MeanAtomicMass;                                             \
-  using EosBase<EOSDERIVED>::MeanAtomicNumber;
+  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,
@@ -103,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;
@@ -135,8 +147,9 @@ struct Transform {
 };
 
 /*
-  This is a utility struct used for analytic equations of state,
-  allowing them to store/report mean atomic mass/number easily.
+  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;
@@ -149,6 +162,15 @@ struct MeanAtomicProperties {
   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);
+  }
 };
 
 /*
@@ -758,26 +780,14 @@ class EosBase {
   PORTABLE_INLINE_FUNCTION
   Real RhoPmin(const Real temp) const { return 0.0; }
 
-  // Defaults. Likely not accurate. You should almost always
-  // overwrite. For a given atom, always an integer, but this is a
-  // mean.
   // 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.
-  PORTABLE_INLINE_FUNCTION
-  Real MeanAtomicMass() const {
-    PORTABLE_THROW_OR_ABORT("Mean atomic mass not implemented!");
-    return 1.0;
-  }
-  PORTABLE_INLINE_FUNCTION
-  Real MeanAtomicNumber() const {
-    PORTABLE_THROW_OR_ABORT("Mean atomic number not implemented!");
-    return 1.0;
-  }
-  // TODO(JMM): Should we provide vector implementations if we depend on rho, T, etc?
+  // 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,

singularity-eos/eos/eos_carnahan_starling.hpp

@@ -74,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(
@@ -224,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)
@@ -238,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
@@ -250,11 +254,6 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
   static std::string EosType() { return std::string("CarnahanStarling"); }
   static std::string EosPyType() { return EosType(); }
 
-  PORTABLE_INLINE_FUNCTION
-  Real MeanAtomicMass() const { return _AZbar.Abar; }
-  PORTABLE_INLINE_FUNCTION
-  Real MeanAtomicNumber() const { return _AZbar.Zbar; }
-
  private:
   Real _Cv, _gm1, _bb, _qq;
   // optional reference state variables

singularity-eos/eos/eos_davis.hpp

@@ -46,6 +46,7 @@ class DavisReactants : public EosBase<DavisReactants> {
   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(
@@ -147,10 +148,6 @@ class DavisReactants : public EosBase<DavisReactants> {
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     return Gamma(rho);
   }
-  PORTABLE_INLINE_FUNCTION
-  Real MeanAtomicMass() const { return _AZbar.Abar; }
-  PORTABLE_INLINE_FUNCTION
-  Real MeanAtomicNumber() const { return _AZbar.Zbar; }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
@@ -165,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)
@@ -180,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"); }
@@ -209,6 +209,7 @@ class DavisProducts : public EosBase<DavisProducts> {
       : _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; }
@@ -292,10 +293,6 @@ class DavisProducts : public EosBase<DavisProducts> {
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     return Gamma(rho);
   }
-  PORTABLE_INLINE_FUNCTION
-  Real MeanAtomicMass() const { return _AZbar.Abar; }
-  PORTABLE_INLINE_FUNCTION
-  Real MeanAtomicNumber() const { return _AZbar.Zbar; }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   DensityEnergyFromPressureTemperature(const Real press, const Real temp,
@@ -310,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)
@@ -324,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"); }

singularity-eos/eos/eos_eospac.hpp

@@ -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; }
 
@@ -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
@@ -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_; }
@@ -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 = {
@@ -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_)};

singularity-eos/eos/eos_gruneisen.hpp

@@ -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);
@@ -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
@@ -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;
@@ -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; }
@@ -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
@@ -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 {

singularity-eos/eos/eos_helmholtz.hpp

@@ -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,

singularity-eos/eos/eos_ideal.hpp

@@ -71,6 +71,7 @@ class IdealGas : public EosBase<IdealGas> {
                             "Entropy reference temperature must be positive");
     PORTABLE_ALWAYS_REQUIRE(_EntropyRho0 >= 0,
                             "Entropy reference density must be positive");
+    _AZbar.CheckParams();
     return;
   }
   template <typename Indexer_t = Real *>
@@ -153,11 +154,6 @@ class IdealGas : public EosBase<IdealGas> {
           const unsigned long output,
           Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
 
-  PORTABLE_INLINE_FUNCTION
-  Real MeanAtomicMass() const { return _AZbar.Abar; }
-  PORTABLE_INLINE_FUNCTION
-  Real MeanAtomicNumber() const { return _AZbar.Zbar; }
-
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
@@ -175,6 +171,7 @@ class IdealGas : public EosBase<IdealGas> {
   }
   // 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(IdealGas)
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return 0; }
@@ -185,6 +182,7 @@ class IdealGas : public EosBase<IdealGas> {
   static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; }
   PORTABLE_INLINE_FUNCTION void PrintParams() const {
     printf("Ideal Gas Parameters:\nGamma = %g\nCv    = %g\n", _gm1 + 1.0, _Cv);
+    _AZbar.PrintParams();
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void