Merge pull request #58 from lanl/brryan/units

Non-CGS constants
lanl · Dec 12, 2024 · b4542ed · b4542ed
2 parents 9d0a1b7 + 2679682
commit b4542ed
Show file tree

Hide file tree

Showing 18 changed files with 259 additions and 68 deletions.
diff --git a/singularity-opac/constants/constants.hpp b/singularity-opac/constants/constants.hpp
@@ -137,10 +137,6 @@ struct PhysicalConstants {
                                                  (mass * time * time);
   static constexpr Real g_newt = gravitational_constant;
 
-  static constexpr Real acceleration_from_gravity =
-      BASE::acceleration_from_gravity * length / (time * time);
-  static constexpr Real g_accel = acceleration_from_gravity;
-
   static constexpr Real electron_mass = BASE::electron_mass * mass;
   static constexpr Real me = electron_mass;
 
@@ -163,19 +159,6 @@ struct PhysicalConstants {
   static constexpr Real radiation_constant = 4.0 * sb / c;
   static constexpr Real ar = radiation_constant;
 
-  // Faraday constant (derived)
-  static constexpr Real faraday_constant = na * qe;
-  static constexpr Real faraday = faraday_constant;
-
-  // Permeability of free space
-  static constexpr Real vacuum_permeability =
-      BASE::vacuum_permeability * force / (current * current);
-  static constexpr Real mu0 = vacuum_permeability;
-
-  // Permittivity of free space (derived)
-  static constexpr Real vacuum_permittivity = 1.0 / (mu0 * c * c);
-  static constexpr Real eps0 = vacuum_permittivity;
-
   // Electron volt (derived)
   static constexpr Real electron_volt = BASE::elementary_charge * energy;
   static constexpr Real eV = electron_volt;
@@ -195,14 +178,38 @@ struct PhysicalConstants {
 };
 
 struct RuntimePhysicalConstants {
+  template <typename T>
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants(T pc)
+      : length(1.), time(1.), mass(1.), temp(1.), na(pc.na), alpha(pc.alpha),
+        h(pc.h), hbar(pc.hbar), kb(pc.kb), r_gas(pc.r_gas), qe(pc.qe), c(pc.c),
+        g_newt(pc.g_newt), me(pc.me), mp(pc.mp), mn(pc.mn), amu(pc.amu),
+        sb(pc.sb), ar(pc.ar), eV(pc.eV), Fc(pc.Fc), nu_sigma0(pc.nu_sigma0),
+        gA(pc.gA) {}
+
   template <typename T>
   PORTABLE_INLINE_FUNCTION
-  RuntimePhysicalConstants(T pc)
-      : na(pc.na), alpha(pc.alpha), h(pc.h), hbar(pc.hbar), kb(pc.kb),
-        r_gas(pc.r_gas), qe(pc.qe), c(pc.c), g_newt(pc.g_newt),
-        g_accel(pc.g_accel), me(pc.me), mp(pc.mp), mn(pc.mn), amu(pc.amu),
-        sb(pc.sb), ar(pc.ar), faraday(pc.faraday), mu0(pc.mu0), eps0(pc.eps0),
-        eV(pc.eV), Fc(pc.Fc), nu_sigma0(pc.nu_sigma0), gA(pc.gA) {}
+  RuntimePhysicalConstants(T pc, const Real time_, const Real mass_,
+                           const Real length_, const Real temp_)
+      : time(time_), mass(mass_), length(length_), temp(temp_), na(pc.na),
+        alpha(pc.alpha), h(pc.h * time / mass * std::pow(length, -2)),
+        hbar(pc.hbar * time / mass * std::pow(length, -2)),
+        kb(pc.kb * std::pow(time, 2) / mass * std::pow(length, -2) * temp),
+        r_gas(pc.r_gas * std::pow(time, 2) / mass * std::pow(length, -2)),
+        qe(pc.qe), c(pc.c * time / length),
+        g_newt(pc.g_newt * std::pow(length, -3) / mass * std::pow(time, 2)),
+        me(pc.me / mass), mp(pc.mp / mass), mn(pc.mn / mass),
+        amu(pc.amu / mass),
+        sb(pc.sb * std::pow(time, 3) / mass * std::pow(temp, 4)),
+        ar(pc.ar * std::pow(time, 2) * length / mass * std::pow(temp, 4)),
+        eV(pc.eV * std::pow(time, 2) / mass * std::pow(length, -2)),
+        Fc(pc.Fc * std::pow(time, 2) * length / mass),
+        nu_sigma0(pc.nu_sigma0 / std::pow(length, 2)), gA(pc.gA) {}
+
+  // Conversion factors to the T pc unit system from code units
+  const Real time;
+  const Real mass;
+  const Real length;
+  const Real temp;
 
   const Real na;
   const Real alpha;
@@ -213,28 +220,18 @@ struct RuntimePhysicalConstants {
   const Real qe;
   const Real c;
   const Real g_newt;
-  const Real g_accel;
   const Real me;
   const Real mp;
   const Real mn;
   const Real amu;
   const Real sb;
   const Real ar;
-  const Real faraday;
-  const Real mu0;
-  const Real eps0;
   const Real eV;
   const Real Fc;
   const Real nu_sigma0;
   const Real gA;
 };
 
-template <typename T>
-PORTABLE_INLINE_FUNCTION
-RuntimePhysicalConstants GetRuntimePhysicalConstants(T phys_constants) {
-  return RuntimePhysicalConstants(phys_constants);
-}
-
 using PhysicalConstantsUnity =
     PhysicalConstants<BaseUnity, UnitConversionDefault>;
 using PhysicalConstantsSI =

diff --git a/singularity-opac/neutrinos/brt_neutrinos.hpp b/singularity-opac/neutrinos/brt_neutrinos.hpp
@@ -173,6 +173,11 @@ class BRTOpacity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   PORTABLE_INLINE_FUNCTION
   Real GetSigmac(const RadiationType type, const Real nu) const {

diff --git a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
@@ -175,6 +175,11 @@ class GrayOpacity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   Real kappa_; // Opacity. Units of cm^2/g
   ThermalDistribution dist_;

diff --git a/singularity-opac/neutrinos/mean_neutrino_variant.hpp b/singularity-opac/neutrinos/mean_neutrino_variant.hpp
@@ -72,11 +72,7 @@ class MeanVariant {
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
     return mpark::visit(
-        [=](const auto &opac) {
-          using PC = typename std::decay_t<decltype(opac)>::PC;
-          return singularity::GetRuntimePhysicalConstants(PC());
-        },
-        opac_);
+        [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   PORTABLE_INLINE_FUNCTION Real PlanckMeanAbsorptionCoefficient(

diff --git a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
@@ -37,9 +37,11 @@ namespace impl {
 // TODO(BRR) Note: It is assumed that lambda is constant for all densities,
 // temperatures, and Ye
 
+template <typename pc = PhysicalConstantsCGS>
 class MeanOpacity {
-
  public:
+  using PC = pc;
+
   MeanOpacity() = default;
   template <typename Opacity>
   MeanOpacity(const Opacity &opac, const Real lRhoMin, const Real lRhoMax,
@@ -106,6 +108,11 @@ class MeanOpacity {
     lkappaRosseland_.finalize();
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
   PORTABLE_INLINE_FUNCTION
   Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp,
                                        const Real Ye,
@@ -133,7 +140,9 @@ class MeanOpacity {
                         const Real lTMax, const int NT, const Real YeMin,
                         const Real YeMax, const int NYe, Real lNuMin,
                         Real lNuMax, const int NNu, Real *lambda = nullptr) {
-    using PC = typename Opacity::PC;
+    static_assert(std::is_same<typename Opacity::PC, PC>::value,
+                  "Error: MeanOpacity physical constants do not match those of "
+                  "Opacity used for construction.");
 
     lkappaPlanck_.resize(NRho, NT, NYe, NEUTRINO_NTYPES);
     // index 0 is the species and is not interpolatable
@@ -220,7 +229,7 @@ class MeanOpacity {
 
 } // namespace impl
 
-using MeanOpacityBase = impl::MeanOpacity;
+using MeanOpacityBase = impl::MeanOpacity<PhysicalConstantsCGS>;
 using MeanOpacity =
     impl::MeanVariant<MeanOpacityBase, MeanNonCGSUnits<MeanOpacityBase>>;
 

diff --git a/singularity-opac/neutrinos/neutrino_variant.hpp b/singularity-opac/neutrinos/neutrino_variant.hpp
@@ -74,11 +74,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
     return mpark::visit(
-        [=](const auto &opac) {
-          using PC = typename std::decay_t<decltype(opac)>::PC;
-          return singularity::GetRuntimePhysicalConstants(PC());
-        },
-        opac_);
+        [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   // Directional absorption coefficient with units of 1/length

diff --git a/singularity-opac/neutrinos/non_cgs_neutrinos.hpp b/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
@@ -239,6 +239,12 @@ class NonCGSUnits {
     return NoH * mass_unit_ / rho_unit_;
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC(), time_unit_, mass_unit_, length_unit_,
+                                    temp_unit_);
+  }
+
  private:
   Opac opac_;
   Real time_unit_, mass_unit_, length_unit_, temp_unit_;
@@ -249,6 +255,8 @@ class NonCGSUnits {
 template <typename MeanOpac>
 class MeanNonCGSUnits {
  public:
+  using PC = typename MeanOpac::PC;
+
   MeanNonCGSUnits() = default;
   MeanNonCGSUnits(MeanOpac &&mean_opac, const Real time_unit,
                   const Real mass_unit, const Real length_unit,
@@ -298,6 +306,12 @@ class MeanNonCGSUnits {
     return alpha * length_unit_;
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PhysicalConstantsCGS(), time_unit_,
+                                    mass_unit_, length_unit_, temp_unit_);
+  }
+
  private:
   MeanOpac mean_opac_;
   Real time_unit_, mass_unit_, length_unit_, temp_unit_;

diff --git a/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp b/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
@@ -371,6 +371,11 @@ class SpinerOpacity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
 #ifdef SPINER_USE_HDF
   void LoadFromSP5_(const std::string &filename) {

diff --git a/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp b/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
@@ -191,6 +191,11 @@ class TophatEmissivity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   PORTABLE_INLINE_FUNCTION
   Real GetYeF(RadiationType type, Real Ye) const {

diff --git a/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp b/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
@@ -179,6 +179,11 @@ class EPBremsstrahlungOpacity {
     return dist_.NumberDensityFromTemperature(temp, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   Real mass_ion_;
   PlanckDistribution<pc> dist_;

diff --git a/singularity-opac/photons/gray_opacity_photons.hpp b/singularity-opac/photons/gray_opacity_photons.hpp
@@ -167,6 +167,11 @@ class GrayOpacity {
     return dist_.NumberDensityFromTemperature(temp, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   Real kappa_; // Opacity. Units of cm^2/g
   PlanckDistribution<pc> dist_;

diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp
@@ -39,9 +39,11 @@ namespace impl {
 // TODO(BRR) Note: It is assumed that lambda is constant for all densities and
 // temperatures
 
+template <typename pc = PhysicalConstantsCGS>
 class MeanOpacity {
-
  public:
+  using PC = pc;
+
   MeanOpacity() = default;
   template <typename Opacity>
   MeanOpacity(const Opacity &opac, const Real lRhoMin, const Real lRhoMax,
@@ -121,6 +123,11 @@ class MeanOpacity {
     lkappa_.finalize();
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
   PORTABLE_INLINE_FUNCTION
   Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const {
     Real lRho = toLog_(rho);
@@ -142,7 +149,9 @@ class MeanOpacity {
                         const Real lRhoMax, const int NRho, const Real lTMin,
                         const Real lTMax, const int NT, Real lNuMin,
                         Real lNuMax, const int NNu, Real *lambda = nullptr) {
-    using PC = typename Opacity::PC;
+    static_assert(std::is_same<typename Opacity::PC, PC>::value,
+                  "Error: MeanOpacity physical constants do not match those of "
+                  "Opacity used for construction.");
 
     lkappa_.resize(NRho, NT, 2);
     lkappa_.setRange(1, lTMin, lTMax, NT);
@@ -302,7 +311,7 @@ class MeanOpacity {
 
 } // namespace impl
 
-using MeanOpacityBase = impl::MeanOpacity;
+using MeanOpacityBase = impl::MeanOpacity<PhysicalConstantsCGS>;
 using MeanOpacity =
     impl::MeanVariant<MeanOpacityBase, MeanNonCGSUnits<MeanOpacityBase>>;
 

diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp
@@ -72,11 +72,7 @@ class MeanVariant {
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
     return mpark::visit(
-        [=](const auto &opac) {
-          using PC = typename std::decay_t<decltype(opac)>::PC;
-          return singularity::GetRuntimePhysicalConstants(PC());
-        },
-        opac_);
+        [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   PORTABLE_INLINE_FUNCTION Real

diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp
@@ -222,9 +222,10 @@ class NonCGSUnits {
     return NoH * mass_unit_ / rho_unit_;
   }
 
-  template <typename T>
-  PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
-    return opac_.GetPhysicalConstants();
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC(), time_unit_, mass_unit_, length_unit_,
+                                    temp_unit_);
   }
 
  private:
@@ -237,6 +238,8 @@ class NonCGSUnits {
 template <typename MeanOpac>
 class MeanNonCGSUnits {
  public:
+  using PC = typename MeanOpac::PC;
+
   MeanNonCGSUnits() = default;
   MeanNonCGSUnits(MeanOpac &&mean_opac, const Real time_unit,
                   const Real mass_unit, const Real length_unit,
@@ -283,6 +286,12 @@ class MeanNonCGSUnits {
     return alpha * length_unit_;
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PhysicalConstantsCGS(), time_unit_,
+                                    mass_unit_, length_unit_, temp_unit_);
+  }
+
  private:
   MeanOpac mean_opac_;
   Real time_unit_, mass_unit_, length_unit_, temp_unit_;

diff --git a/singularity-opac/photons/photon_variant.hpp b/singularity-opac/photons/photon_variant.hpp
@@ -74,11 +74,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
     return mpark::visit(
-        [=](const auto &opac) {
-          using PC = typename std::decay_t<decltype(opac)>::PC;
-          return singularity::GetRuntimePhysicalConstants(PC());
-        },
-        opac_);
+        [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   // Directional absorption coefficient with units of 1/length