Return physical constants object from opacities

singularity-opac/neutrinos/gray_opacity_neutrinos.hpp

@@ -173,6 +173,9 @@ 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_;

singularity-opac/neutrinos/mean_neutrino_variant.hpp

@@ -28,7 +28,7 @@ namespace singularity {
 namespace neutrinos {
 namespace impl {
 
-template <typename... Opacs>
+template <typename CONSTANTS, typename... Opacs>
 class MeanVariant {
  private:
   opac_variant<Opacs...> opac_;
@@ -88,6 +88,16 @@ class MeanVariant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION CONSTANTS GetPhysicalConstants() const {
+    return mpark::visit([](auto &opac) { return opac.GetPhysicalConstants(); },
+                        opac_);
+  }
+  // static CONSTANTS GetPhysicalConstants() {
+  //  return mpark::visit(
+  //      [](auto &opac) { return decltype(opac)::GetPhysicalConstants(); },
+  //     opac_);
+  // }
+
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
   }

singularity-opac/neutrinos/mean_opacity_neutrinos.hpp

@@ -102,6 +102,10 @@ class MeanOpacity {
     return other;
   }
 
+  PORTABLE_INLINE_FUNCTION
+  pc GetPhysicalConstants() const { return pc(); }
+  // static pc GetPhysicalConstants() { return pc(); }
+
   void Finalize() {
     lkappaPlanck_.finalize();
     lkappaRosseland_.finalize();

singularity-opac/photons/gray_opacity_photons.hpp

@@ -165,6 +165,9 @@ class GrayOpacity {
     return dist_.NumberDensityFromTemperature(temp, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION
+  pc GetPhysicalConstants() const { return pc(); }
+
  private:
   Real kappa_; // Opacity. Units of cm^2/g
   PlanckDistribution<pc> dist_;

singularity-opac/photons/mean_opacity_photons.hpp

@@ -119,6 +119,9 @@ 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,

singularity-opac/photons/mean_photon_variant.hpp

@@ -28,7 +28,7 @@ namespace singularity {
 namespace photons {
 namespace impl {
 
-template <typename... Opacs>
+template <typename CONSTANTS, typename... Opacs>
 class MeanVariant {
  private:
   opac_variant<Opacs...> opac_;
@@ -86,6 +86,11 @@ class MeanVariant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION CONSTANTS GetPhysicalConstants() const {
+    return mpark::visit([](auto &opac) { return opac.GetPhysicalConstants(); },
+                        opac_);
+  }
+
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
   }

singularity-opac/photons/photon_variant.hpp

@@ -30,7 +30,7 @@ namespace impl {
 template <typename... Ts>
 using opac_variant = mpark::variant<Ts...>;
 
-template <typename... Opacs>
+template <typename CONSTANTS, typename... Opacs>
 class Variant {
  private:
   opac_variant<Opacs...> opac_;
@@ -275,6 +275,11 @@ class Variant {
                         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION CONSTANTS GetPhysicalConstants() const {
+    return mpark::visit([](auto &opac) { return opac.GetPhysicalConstants(); },
+                        opac_);
+  }
+
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
   }

test/test_gray_opacities.cpp

@@ -57,6 +57,9 @@ TEST_CASE("Gray neutrino opacities", "[GrayNeutrinos]") {
     neutrinos::Gray opac_host(1);
     neutrinos::Opacity opac = opac_host.GetOnDevice();
 
+    auto constants = opac_host.GetPhysicalConstants();
+    printf("cl: %e\n", constants.c);
+
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
       int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -283,8 +286,7 @@ TEST_CASE("Gray photon opacities", "[GrayPhotons]") {
 
       Real *nu_bins = (Real *)PORTABLE_MALLOC(nbins * sizeof(Real));
       Real *temp_bins = (Real *)PORTABLE_MALLOC(ntemps * sizeof(Real));
-      DataBox loglin_bins(Spiner::AllocationTarget::Device, ntemps,
-                                  nbins);
+      DataBox loglin_bins(Spiner::AllocationTarget::Device, ntemps, nbins);
 
       portableFor(
           "set nu bins", 0, nbins, PORTABLE_LAMBDA(const int &i) {

test/test_mean_opacities.cpp

@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -46,7 +46,7 @@ using atomic_view = Kokkos::MemoryTraits<Kokkos::Atomic>;
 
 template <typename T>
 PORTABLE_INLINE_FUNCTION T FractionalDifference(const T &a, const T &b) {
-  return 2 * std::abs(b - a) / (std::abs(a) + std::abs(b) + 1e-20);
+  return 2 * std::abs(b - a) / (std::abs(a) + std::abs(b) + 1e-100);
 }
 constexpr Real EPS_TEST = 1e-3;
 template <typename T>
@@ -89,9 +89,24 @@ TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") {
 
     neutrinos::MeanOpacityCGS mean_opac_host(
         opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin, YeMax, NYe);
-    // neutrinos::MeanOpacity mean_opac = mean_opac_host.GetOnDevice();
     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
@@ -611,7 +626,8 @@ TEST_CASE("Mean photon scattering opacities", "[MeanPhotonS]") {
         int n_wrong = 0;
         portableReduce(
             "rebuilt table vs gray", 0, NRho, 0, NT, 0, 0,
-            PORTABLE_LAMBDA(const int iRho, const int iT, const int igarbage, int &accumulate) {
+            PORTABLE_LAMBDA(const int iRho, const int iT, const int igarbage,
+                            int &accumulate) {
               const Real lRho =
                   lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho;
               const Real rho = std::pow(10, lRho);