diff --git a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
index fc9a57e..b1f527b 100644
--- a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
+++ b/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_;
diff --git a/singularity-opac/neutrinos/mean_neutrino_variant.hpp b/singularity-opac/neutrinos/mean_neutrino_variant.hpp
index fd50081..079e463 100644
--- 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... 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_);
}
diff --git a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
index 61a6396..4124319 100644
--- a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
+++ b/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();
diff --git a/singularity-opac/photons/gray_opacity_photons.hpp b/singularity-opac/photons/gray_opacity_photons.hpp
index 9bc0358..434df97 100644
--- a/singularity-opac/photons/gray_opacity_photons.hpp
+++ b/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_;
diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp
index 8a58859..4e0b410 100644
--- a/singularity-opac/photons/mean_opacity_photons.hpp
+++ b/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,
diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp
index af292ea..604d061 100644
--- 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... 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_);
}
diff --git a/singularity-opac/photons/photon_variant.hpp b/singularity-opac/photons/photon_variant.hpp
index 440e5e8..b0c9979 100644
--- a/singularity-opac/photons/photon_variant.hpp
+++ b/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_);
}
diff --git a/test/test_gray_opacities.cpp b/test/test_gray_opacities.cpp
index 464dbd1..a26ad73 100644
--- a/test/test_gray_opacities.cpp
+++ b/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) {
diff --git a/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp
index b6ad91f..467f583 100644
--- a/test/test_mean_opacities.cpp
+++ b/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);