OK I can get the constants I need from both frequency-dependent and g…

…ray opacs
lanl · Dec 11, 2024 · 6f45ba4 · 6f45ba4
1 parent c4b42bc
commit 6f45ba4
Show file tree

Hide file tree

Showing 4 changed files with 53 additions and 11 deletions.
diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp
@@ -37,9 +37,11 @@ namespace impl {
 // TODO(BRR) Note: It is assumed that lambda is constant for all densities and
 // temperatures
 
+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,
@@ -196,7 +198,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/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp
@@ -72,11 +72,7 @@ class MeanVariant {
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
     return mpark::visit(
-        [](auto &opac) {
-          using PC = typename std::decay_t<decltype(opac)>::PC;
-          return RuntimePhysicalConstants(PC());
-        },
-        opac_);
+        [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   PORTABLE_INLINE_FUNCTION Real

diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp
@@ -229,7 +229,8 @@ class NonCGSUnits {
 
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
-    return RuntimePhysicalConstants(PC(), time_unit_, mass_unit_, length_unit_, temp_unit_);
+    return RuntimePhysicalConstants(PC(), time_unit_, mass_unit_, length_unit_,
+                                    temp_unit_);
   }
 
  private:
@@ -242,6 +243,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,
@@ -288,6 +291,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/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp
@@ -495,12 +495,47 @@ 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::MeanOpacity>(
-          std::forward<photons::MeanOpacity>(mean_opac_host), time_unit,
-          mass_unit, length_unit, temp_unit);
+      auto mean_opac_host_base = photons::MeanOpacityBase(
+          opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT);
+      auto funny_units_host =
+          photons::MeanNonCGSUnits<photons::MeanOpacityBase>(
+              std::forward<photons::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");