Update neutrinos

lanl · Dec 11, 2024 · a5bd748 · a5bd748
1 parent 6f45ba4
commit a5bd748
Show file tree

Hide file tree

Showing 12 changed files with 130 additions and 20 deletions.
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_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 RuntimePhysicalConstants(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
@@ -369,6 +369,11 @@ class SpinerOpacity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   // TODO(JMM): Offsets probably not necessary
   PORTABLE_INLINE_FUNCTION Real toLog_(const Real x, const Real offset) const {

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/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp
@@ -104,6 +104,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 lRho = toLog_(rho);
@@ -125,7 +130,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.");
 
     lkappaPlanck_.resize(NRho, NT);
     lkappaPlanck_.setRange(0, lTMin, lTMax, NT);

diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp
@@ -222,10 +222,10 @@ class NonCGSUnits {
     return NoH * mass_unit_ / rho_unit_;
   }
 
-  template <typename T>
-  PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
-    return opac_.GetPhysicalConstants();
-  }
+  // template <typename T>
+  // PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
+  //  return opac_.GetPhysicalConstants();
+  //}
 
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {

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

diff --git a/test/test_gray_opacities.cpp b/test/test_gray_opacities.cpp
@@ -109,6 +109,38 @@ TEST_CASE("Gray neutrino opacities", "[GrayNeutrinos]") {
               neutrinos::Gray(1), time_unit, mass_unit, length_unit, temp_unit);
       auto funny_units = funny_units_host.GetOnDevice();
 
+      THEN("We can retrieve physical constants in code units") {
+        auto noncgs_rpc = funny_units.GetRuntimePhysicalConstants();
+        REQUIRE(FractionalDifference(noncgs_rpc.length, length_unit) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.time, time_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mass, mass_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.temp, temp_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.na, 6.022141e+23) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.alpha, 7.297353e-03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.h, 2.871060e-33) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.hbar, 4.569434e-34) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.kb, 2.030877e-18) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.r_gas, 4.431243e+03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.qe, 4.803205e-10) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.c, 4.673571e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.g_newt, 4.508065e-15) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.me, 1.997672e-30) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mp, 3.668030e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mn, 3.673086e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.amu, 3.641533e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.sb, 1.342760e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.ar, 1.149237e+00) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.eV, 8.538896e-17) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.Fc, 1.189435e-09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.nu_sigma0, 2.829094e-74) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.gA, -1.23) < EPS_TEST);
+      }
+
       THEN("We can convert meaningfully into and out of funny units") {
         int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS

diff --git a/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp
@@ -174,13 +174,48 @@ TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") {
       constexpr Real rho_unit =
           mass_unit / (length_unit * length_unit * length_unit);
 
+      auto mean_opac_host_base =
+          neutrinos::MeanOpacityBase(opac_host, lRhoMin, lRhoMax, NRho, lTMin,
+                                     lTMax, NT, YeMin, YeMax, NYe);
       auto funny_units_host =
-          neutrinos::MeanNonCGSUnits<neutrinos::MeanOpacity>(
-              std::forward<neutrinos::MeanOpacity>(mean_opac_host), time_unit,
-              mass_unit, length_unit, temp_unit);
+          neutrinos::MeanNonCGSUnits<neutrinos::MeanOpacityBase>(
+              std::forward<neutrinos::MeanOpacityBase>(mean_opac_host_base),
+              time_unit, mass_unit, length_unit, temp_unit);
 
       auto funny_units = funny_units_host.GetOnDevice();
 
+      THEN("We can retrieve physical constants in code units") {
+        auto noncgs_rpc = funny_units.GetRuntimePhysicalConstants();
+        REQUIRE(FractionalDifference(noncgs_rpc.length, length_unit) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.time, time_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mass, mass_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.temp, temp_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.na, 6.022141e+23) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.alpha, 7.297353e-03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.h, 2.871060e-33) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.hbar, 4.569434e-34) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.kb, 2.030877e-18) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.r_gas, 4.431243e+03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.qe, 4.803205e-10) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.c, 4.673571e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.g_newt, 4.508065e-15) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.me, 1.997672e-30) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mp, 3.668030e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mn, 3.673086e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.amu, 3.641533e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.sb, 1.342760e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.ar, 1.149237e+00) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.eV, 8.538896e-17) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.Fc, 1.189435e-09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.nu_sigma0, 2.829094e-74) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.gA, -1.23) < EPS_TEST);
+      }
+
       int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
       Kokkos::View<int, atomic_view> n_wrong_d("wrong");