Constants now work, runtime only

lanl · Sep 27, 2024 · abbef3c · abbef3c
1 parent 6d5323c
commit abbef3c
Show file tree

Hide file tree

Showing 13 changed files with 107 additions and 26 deletions.
diff --git a/singularity-opac/constants/constants.hpp b/singularity-opac/constants/constants.hpp
@@ -194,6 +194,45 @@ struct PhysicalConstants {
   static constexpr Real gA = axial_vector_coupling_constant;
 };
 
+struct RuntimePhysicalConstants {
+  template <typename T>
+  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) {}
+
+  const Real na;
+  const Real alpha;
+  const Real h;
+  const Real hbar;
+  const Real kb;
+  const Real r_gas;
+  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>
+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
@@ -32,6 +32,8 @@ namespace neutrinos {
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class BRTOpacity {
  public:
+  using PC = pc;
+
   BRTOpacity() = default;
   BRTOpacity(const ThermalDistribution &dist) : dist_(dist) {}
   BRTOpacity GetOnDevice() { return *this; }

diff --git a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
@@ -30,6 +30,8 @@ namespace neutrinos {
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class GrayOpacity {
  public:
+  using PC = pc;
+
   GrayOpacity() = default;
   GrayOpacity(const Real kappa) : kappa_(kappa) {}
   GrayOpacity(const ThermalDistribution &dist, const Real kappa)
@@ -42,6 +44,12 @@ class GrayOpacity {
   void PrintParams() const noexcept {
     printf("Gray opacity. kappa = %g\n", kappa_);
   }
+
+  PORTABLE_INLINE_FUNCTION
+  RuntimePhysicalConstants GetPhysicalConstants() const {
+    return GetRuntimePhysicalConstants(pc());
+  }
+
   inline void Finalize() noexcept {}
 
   PORTABLE_INLINE_FUNCTION
@@ -173,8 +181,8 @@ class GrayOpacity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
-  PORTABLE_INLINE_FUNCTION
-  pc GetPhysicalConstants() const { return pc(); }
+  // PORTABLE_INLINE_FUNCTION
+  // pc GetPhysicalConstants() const { return pc(); }
 
  private:
   Real kappa_; // Opacity. Units of cm^2/g

diff --git a/singularity-opac/neutrinos/neutrino_variant.hpp b/singularity-opac/neutrinos/neutrino_variant.hpp
@@ -71,6 +71,16 @@ class Variant {
         opac_);
   }
 
+  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_);
+  }
+
   // Directional absorption coefficient with units of 1/length
   // Signature should be at least
   // rho, temp, Ye, type, nu, lambda
@@ -290,11 +300,12 @@ class Variant {
                         opac_);
   }
 
-  template <typename T>
-  PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
-    return mpark::visit([](auto &opac) { return opac.GetPhysicalConstants(); },
-                        opac_);
-  }
+  // template <typename T>
+  // PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
+  //  return mpark::visit([](auto &opac) { return opac.GetPhysicalConstants();
+  //  },
+  //                      opac_);
+  //}
 
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);

diff --git a/singularity-opac/neutrinos/non_cgs_neutrinos.hpp b/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
@@ -29,6 +29,8 @@ namespace neutrinos {
 template <typename Opac>
 class NonCGSUnits {
  public:
+  using PC = typename Opac::PC;
+
   NonCGSUnits() = default;
   NonCGSUnits(Opac &&opac, const Real time_unit, const Real mass_unit,
               const Real length_unit, const Real temp_unit)
@@ -48,11 +50,6 @@ class NonCGSUnits {
                              length_unit_, temp_unit_);
   }
 
-  template <typename T>
-  PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
-    return opac_.GetPhysicalConstants();
-  }
-
   inline void Finalize() noexcept { opac_.Finalize(); }
 
   PORTABLE_INLINE_FUNCTION

diff --git a/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp b/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
@@ -61,6 +61,7 @@ enum class DataStatus { Deallocated, OnDevice, OnHost };
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class SpinerOpacity {
  public:
+  using PC = pc;
   using DataBox = Spiner::DataBox<Real>;
   static constexpr Real EPS = 10.0 * std::numeric_limits<Real>::min();
   static constexpr Real Hz2MeV = pc::h / (1e6 * pc::eV);
@@ -110,7 +111,8 @@ class SpinerOpacity {
             Real J = std::max(opac.Emissivity(rho, T * MeV2K, Ye, type), 0.0);
             Real lJ = toLog_(J);
             lJ_(iRho, iT, iYe, idx) = lJ;
-            Real JYe = std::max(opac.NumberEmissivity(rho, T * MeV2K, Ye, type), 0.0);
+            Real JYe =
+                std::max(opac.NumberEmissivity(rho, T * MeV2K, Ye, type), 0.0);
             lJYe_(iRho, iT, iYe, idx) = toLog_(JYe);
             for (int ie = 0; ie < Ne; ++ie) {
               Real lE = lalphanu_.range(0).x(ie);
@@ -119,8 +121,8 @@ class SpinerOpacity {
               Real alpha = std::max(
                   opac.AbsorptionCoefficient(rho, T, Ye, type, nu), 0.0);
               lalphanu_(iRho, iT, iYe, idx, ie) = toLog_(alpha);
-              Real j = std::max(opac.EmissivityPerNuOmega(rho, T * MeV2K, Ye, type, nu),
-                                0.0);
+              Real j = std::max(
+                  opac.EmissivityPerNuOmega(rho, T * MeV2K, Ye, type, nu), 0.0);
               ljnu_(iRho, iT, iYe, idx, ie) = toLog_(j);
             }
           }
@@ -131,8 +133,8 @@ class SpinerOpacity {
 
   // DataBox constructor. Note that this constructor *shallow* copies
   // the databoxes, so they must be managed externally.
-  SpinerOpacity(const DataBox &lalphanu, const DataBox ljnu,
-                const DataBox lJ, const DataBox lJYe)
+  SpinerOpacity(const DataBox &lalphanu, const DataBox ljnu, const DataBox lJ,
+                const DataBox lJYe)
       : memoryStatus_(impl::DataStatus::OnHost), lalphanu_(lalphanu),
         ljnu_(ljnu), lJ_(lJ), lJYe_(lJYe) {}
 

diff --git a/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp b/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
@@ -32,6 +32,8 @@ namespace neutrinos {
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class TophatEmissivity {
  public:
+  using PC = pc;
+
   TophatEmissivity(const Real C, const Real numin, const Real numax)
       : C_(C), numin_(numin), numax_(numax) {}
   TophatEmissivity(const ThermalDistribution &dist, const Real C,

diff --git a/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp b/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
@@ -32,6 +32,8 @@ namespace photons {
 template <typename pc = PhysicalConstantsCGS>
 class EPBremsstrahlungOpacity {
  public:
+  using PC = pc;
+
   EPBremsstrahlungOpacity() = default;
   EPBremsstrahlungOpacity(const PlanckDistribution<pc> &dist) : dist_(dist) {}
 

diff --git a/singularity-opac/photons/gray_opacity_photons.hpp b/singularity-opac/photons/gray_opacity_photons.hpp
@@ -30,6 +30,8 @@ namespace photons {
 template <typename pc = PhysicalConstantsCGS>
 class GrayOpacity {
  public:
+  using PC = pc;
+
   GrayOpacity() = default;
   GrayOpacity(const Real kappa) : kappa_(kappa) {}
   GrayOpacity(const PlanckDistribution<pc> &dist, const Real kappa)

diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp
@@ -29,6 +29,8 @@ namespace photons {
 template <typename Opac>
 class NonCGSUnits {
  public:
+  using PC = typename Opac::PC;
+
   NonCGSUnits() = default;
   NonCGSUnits(Opac &&opac, const Real time_unit, const Real mass_unit,
               const Real length_unit, const Real temp_unit)

diff --git a/singularity-opac/photons/photon_variant.hpp b/singularity-opac/photons/photon_variant.hpp
@@ -71,6 +71,16 @@ class Variant {
         opac_);
   }
 
+  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_);
+  }
+
   // Directional absorption coefficient with units of 1/length
   // Signature should be at least
   // rho, temp, nu, lambda
@@ -275,11 +285,12 @@ class Variant {
                         opac_);
   }
 
-  template <typename T>
-  PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
-    return mpark::visit([](auto &opac) { return opac.GetPhysicalConstants(); },
-                        opac_);
-  }
+  // template <typename T>
+  // PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
+  //  return mpark::visit([](auto &opac) { return opac.GetPhysicalConstants();
+  //  },
+  //                      opac_);
+  //}
 
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);

diff --git a/singularity-opac/photons/powerlaw_opacity_photons.hpp b/singularity-opac/photons/powerlaw_opacity_photons.hpp
@@ -30,6 +30,8 @@ namespace photons {
 template <typename pc = PhysicalConstantsCGS>
 class PowerLawOpacity {
  public:
+  using PC = pc;
+
   PowerLawOpacity() = default;
   PowerLawOpacity(const Real kappa0, const Real rho_exp, const Real temp_exp)
       : kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp) {}

diff --git a/test/test_gray_opacities.cpp b/test/test_gray_opacities.cpp
@@ -54,12 +54,13 @@ TEST_CASE("Gray neutrino opacities", "[GrayNeutrinos]") {
     constexpr RadiationType type = RadiationType::NU_ELECTRON;
     constexpr Real nu = 1.25 * MeV2Hz; // 1 MeV
 
-    neutrinos::Gray opac_host(1);
+    // neutrinos::Gray opac_host(1);
+    neutrinos::Opacity opac_host = neutrinos::Gray(1);
     neutrinos::Opacity opac = opac_host.GetOnDevice();
 
-    // Check constants from mean opacity
+    // Check constants from opacity
     THEN("Check constants from mean opacity for consistency") {
-      auto constants = opac_host.GetPhysicalConstants();
+      auto constants = opac_host.GetRuntimePhysicalConstants();
       int n_wrong = 0;
       if (FractionalDifference(pc::eV, constants.eV) > EPS_TEST) {
         n_wrong += 1;
@@ -222,7 +223,7 @@ TEST_CASE("Gray photon opacities", "[GrayPhotons]") {
 
     // Check constants from mean opacity
     THEN("Check constants from mean opacity for consistency") {
-      auto constants = opac_host.GetPhysicalConstants();
+      auto constants = opac_host.GetRuntimePhysicalConstants();
       int n_wrong = 0;
       if (FractionalDifference(pc::eV, constants.eV) > EPS_TEST) {
         n_wrong += 1;