Constants returned from opacities but not mean opacities; developer b…

…eware
lanl · Sep 28, 2024 · bfaeb1d · bfaeb1d
1 parent 2a26ada
commit bfaeb1d
Show file tree

Hide file tree

Showing 7 changed files with 49 additions and 65 deletions.
diff --git a/singularity-opac/neutrinos/brt_neutrinos.hpp b/singularity-opac/neutrinos/brt_neutrinos.hpp
@@ -44,7 +44,6 @@ class BRTOpacity {
   void PrintParams() const noexcept {
     printf("Burrows-Reddy-Thompson analytic neutrino opacity.\n");
   }
-
   inline void Finalize() noexcept {}
 
   PORTABLE_INLINE_FUNCTION

diff --git a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
@@ -44,7 +44,6 @@ class GrayOpacity {
   void PrintParams() const noexcept {
     printf("Gray opacity. kappa = %g\n", kappa_);
   }
-
   inline void Finalize() noexcept {}
 
   PORTABLE_INLINE_FUNCTION
@@ -176,9 +175,6 @@ class GrayOpacity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
-  // PORTABLE_INLINE_FUNCTION
-  // pc GetPhysicalConstants() const { return 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
@@ -28,7 +28,7 @@ namespace singularity {
 namespace neutrinos {
 namespace impl {
 
-template <typename PC, typename... Opacs>
+template <typename... Opacs>
 class MeanVariant {
  private:
   opac_variant<Opacs...> opac_;
@@ -69,6 +69,16 @@ class MeanVariant {
         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_);
+  }
+
   PORTABLE_INLINE_FUNCTION Real PlanckMeanAbsorptionCoefficient(
       const Real rho, const Real temp, const Real Ye,
       const RadiationType type) const {
@@ -88,11 +98,6 @@ class MeanVariant {
         opac_);
   }
 
-  PORTABLE_INLINE_FUNCTION PC 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/mean_opacity_neutrinos.hpp b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
@@ -41,6 +41,8 @@ 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,
@@ -102,9 +104,6 @@ class MeanOpacity {
     return other;
   }
 
-  PORTABLE_INLINE_FUNCTION
-  pc GetPhysicalConstants() const { return pc(); }
-
   void Finalize() {
     lkappaPlanck_.finalize();
     lkappaRosseland_.finalize();
@@ -137,6 +136,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) {
+    static_assert(std::is_same<PC, typename Opacity::PC>::value,
+                  "Mean opacity constants must match opacity constants");
+
     lkappaPlanck_.resize(NRho, NT, NYe, NEUTRINO_NTYPES);
     // index 0 is the species and is not interpolatable
     lkappaPlanck_.setRange(1, YeMin, YeMax, NYe);
@@ -162,8 +164,8 @@ class MeanOpacity {
             // Choose default temperature-specific frequency grid if frequency
             // grid not specified
             if (AUTOFREQ) {
-              lNuMin = toLog_(1.e-3 * pc::kb * fromLog_(lTMin) / pc::h);
-              lNuMax = toLog_(1.e3 * pc::kb * fromLog_(lTMax) / pc::h);
+              lNuMin = toLog_(1.e-3 * PC::kb * fromLog_(lTMin) / PC::h);
+              lNuMax = toLog_(1.e3 * PC::kb * fromLog_(lTMax) / PC::h);
             }
             const Real dlnu = (lNuMax - lNuMin) / (NNu - 1);
             // Integrate over frequency

diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp
@@ -41,7 +41,8 @@ template <typename pc = PhysicalConstantsCGS>
 class MeanOpacity {
 
  public:
-  using DataBox = Spiner::DataBox<Real>;
+  using PC = pc;
+
   MeanOpacity() = default;
   template <typename Opacity>
   MeanOpacity(const Opacity &opac, const Real lRhoMin, const Real lRhoMax,
@@ -119,15 +120,15 @@ class MeanOpacity {
     return rho * fromLog_(lkappaRosseland_.interpToReal(lRho, lT));
   }
 
-  PORTABLE_INLINE_FUNCTION
-  pc GetPhysicalConstants() const { return pc(); }
-
  private:
   template <typename Opacity, bool AUTOFREQ>
   void MeanOpacityImpl_(const Opacity &opac, const Real lRhoMin,
                         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) {
+    static_assert(std::is_same<PC, typename Opacity::PC>::value,
+                  "Mean opacity constants must match opacity constants");
+
     lkappaPlanck_.resize(NRho, NT);
     lkappaPlanck_.setRange(0, lTMin, lTMax, NT);
     lkappaPlanck_.setRange(1, lRhoMin, lRhoMax, NRho);
@@ -145,8 +146,8 @@ class MeanOpacity {
         Real kappaRosselandNum = 0.;
         Real kappaRosselandDenom = 0.;
         if (AUTOFREQ) {
-          lNuMin = toLog_(1.e-3 * pc::kb * fromLog_(lTMin) / pc::h);
-          lNuMax = toLog_(1.e3 * pc::kb * fromLog_(lTMax) / pc::h);
+          lNuMin = toLog_(1.e-3 * PC::kb * fromLog_(lTMin) / PC::h);
+          lNuMax = toLog_(1.e3 * PC::kb * fromLog_(lTMax) / PC::h);
         }
         const Real dlnu = (lNuMax - lNuMin) / (NNu - 1);
         // Integrate over frequency
@@ -190,8 +191,8 @@ class MeanOpacity {
   PORTABLE_INLINE_FUNCTION Real fromLog_(const Real lx) const {
     return std::pow(10., lx);
   }
-  DataBox lkappaPlanck_;
-  DataBox lkappaRosseland_;
+  Spiner::DataBox<Real> lkappaPlanck_;
+  Spiner::DataBox<Real> lkappaRosseland_;
   const char *filename_;
 };
 

diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp
@@ -28,7 +28,7 @@ namespace singularity {
 namespace photons {
 namespace impl {
 
-template <typename PC, typename... Opacs>
+template <typename... Opacs>
 class MeanVariant {
  private:
   opac_variant<Opacs...> opac_;
@@ -86,14 +86,19 @@ class MeanVariant {
         opac_);
   }
 
-  PORTABLE_INLINE_FUNCTION PC GetPhysicalConstants() const {
-    return mpark::visit([](auto &opac) { return opac.GetPhysicalConstants(); },
-                        opac_);
-  }
-
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, 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_);
+  }
 };
 
 } // namespace impl

diff --git a/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp
@@ -87,26 +87,10 @@ TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") {
     neutrinos::Gray opac_host(kappa);
     neutrinos::Opacity opac = opac_host.GetOnDevice();
 
-    neutrinos::MeanOpacityCGS mean_opac_host(
+    neutrinos::MeanOpacity mean_opac_host = neutrinos::MeanOpacityCGS(
         opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin, YeMax, NYe);
     auto mean_opac = mean_opac_host.GetOnDevice();
 
-    // Check constants from mean opacity
-    THEN("Check constants from mean opacity for consistency") {
-      auto constants = mean_opac_host.GetPhysicalConstants();
-      int n_wrong = 0;
-      if (FractionalDifference(pc::eV, constants.eV) > EPS_TEST) {
-        n_wrong += 1;
-      }
-      if (FractionalDifference(pc::kb, constants.kb) > EPS_TEST) {
-        n_wrong += 1;
-      }
-      if (FractionalDifference(pc::h, constants.h) > EPS_TEST) {
-        n_wrong += 1;
-      }
-      REQUIRE(n_wrong == 0);
-    }
-
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
       int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -138,7 +122,7 @@ TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") {
 #ifdef SPINER_USE_HDF
     THEN("We can save to disk and reload") {
       mean_opac.Save(grayname);
-      neutrinos::MeanOpacityCGS mean_opac_host_load(grayname);
+      neutrinos::MeanOpacity mean_opac_host_load(grayname);
       AND_THEN("The reloaded table matches the gray opacities") {
 
         auto mean_opac_load = mean_opac_host_load.GetOnDevice();
@@ -428,26 +412,14 @@ TEST_CASE("Mean photon opacities", "[MeanPhotons]") {
     photons::Gray opac_host(kappa);
     photons::Opacity opac = opac_host.GetOnDevice();
 
+    // photons::MeanOpacity mean_opac_host = photons::MeanOpacityCGS(
+    //    opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT);
+    // auto mean_opac = mean_opac_host.GetOnDevice();
+
     photons::MeanOpacityCGS mean_opac_host(opac_host, lRhoMin, lRhoMax, NRho,
                                            lTMin, lTMax, NT);
     auto mean_opac = mean_opac_host.GetOnDevice();
 
-    // Check constants from mean opacity
-    THEN("Check constants from mean opacity for consistency") {
-      auto constants = mean_opac_host.GetPhysicalConstants();
-      int n_wrong = 0;
-      if (FractionalDifference(pc::eV, constants.eV) > EPS_TEST) {
-        n_wrong += 1;
-      }
-      if (FractionalDifference(pc::kb, constants.kb) > EPS_TEST) {
-        n_wrong += 1;
-      }
-      if (FractionalDifference(pc::h, constants.h) > EPS_TEST) {
-        n_wrong += 1;
-      }
-      REQUIRE(n_wrong == 0);
-    }
-
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
       int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -479,7 +451,7 @@ TEST_CASE("Mean photon opacities", "[MeanPhotons]") {
 #ifdef SPINER_USE_HDF
     THEN("We can save to disk and reload") {
       mean_opac.Save(grayname);
-      photons::MeanOpacityCGS mean_opac_host_load(grayname);
+      photons::MeanOpacity mean_opac_host_load(grayname);
       AND_THEN("The reloaded table matches the gray opacities") {
 
         auto mean_opac_load = mean_opac_host_load.GetOnDevice();
@@ -526,6 +498,10 @@ TEST_CASE("Mean photon opacities", "[MeanPhotons]") {
       constexpr Real rho_unit =
           mass_unit / (length_unit * length_unit * length_unit);
 
+      // auto funny_units_host =
+      // photons::MeanNonCGSUnits<photons::MeanOpacityCGS>(
+      //    std::forward<photons::MeanOpacityCGS>(mean_opac_host), time_unit,
+      //    mean_opac_host, time_unit, mass_unit, length_unit, temp_unit);
       auto funny_units_host = photons::MeanNonCGSUnits<photons::MeanOpacity>(
           std::forward<photons::MeanOpacity>(mean_opac_host), time_unit,
           mass_unit, length_unit, temp_unit);