Merge pull request #59 from RyanWollaeger/meanopac_grey_api

Add AbsorptionCoefficient and Emissivity to MeanOpacity variants.

singularity-opac/photons/mean_opacity_photons.hpp

@@ -29,6 +29,7 @@
 
 #include <singularity-opac/photons/mean_photon_variant.hpp>
 #include <singularity-opac/photons/non_cgs_photons.hpp>
+#include <singularity-opac/photons/thermal_distributions_photons.hpp>
 
 namespace singularity {
 namespace photons {
@@ -143,6 +144,25 @@ class MeanOpacity {
     return rho * fromLog_(lkappa_.interpToReal(lRho, lT, Rosseland));
   }
 
+  // standard grey opacity functions
+  PORTABLE_INLINE_FUNCTION
+  Real AbsorptionCoefficient(const Real rho, const Real temp,
+                            const int gmode = Rosseland) const {
+    Real lRho = toLog_(rho);
+    Real lT = toLog_(temp);
+    return rho * fromLog_(lkappa_.interpToReal(lRho, lT, gmode));
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real Emissivity(const Real rho, const Real temp,
+                  const int gmode = Rosseland,
+                  Real *lambda = nullptr) const {
+    Real B = dist_.ThermalDistributionOfT(temp, lambda);
+    Real lRho = toLog_(rho);
+    Real lT = toLog_(temp);
+    return rho * fromLog_(lkappa_.interpToReal(lRho, lT, gmode)) * B;
+  }
+
  private:
   template <typename Opacity, bool AUTOFREQ>
   void MeanOpacityImpl_(const Opacity &opac, const Real lRhoMin,
@@ -303,10 +323,7 @@ class MeanOpacity {
   }
   Spiner::DataBox<Real> lkappa_;
   const char *filename_;
-  enum OpacityAveraging {
-    Rosseland = 0,
-    Planck = 1
-  };
+  PlanckDistribution<PC> dist_;
 };
 
 #undef EPS

singularity-opac/photons/mean_photon_types.hpp

@@ -0,0 +1,31 @@
+// ======================================================================
+// © 2025. 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.
+// Department of Energy/National Nuclear Security Administration. All
+// rights in the program are reserved by Triad National Security, LLC,
+// and the U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others
+// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works,
+// distribute copies to the public, perform publicly and display
+// publicly, and to permit others to do so.
+// ======================================================================
+
+#ifndef SINGULARITY_OPAC_PHOTONS_MEAN_PHOTON_TYPES_
+#define SINGULARITY_OPAC_PHOTONS_MEAN_PHOTON_TYPES_
+
+namespace singularity {
+namespace photons {
+
+// mean-opacity mode
+enum OpacityAveraging {
+  Rosseland = 0,
+  Planck = 1
+};
+
+} // namespace photons
+} // namespace singularity
+
+#endif // SINGULARITY_OPAC_PHOTONS_MEAN_PHOTON_TYPES_

singularity-opac/photons/mean_photon_variant.hpp

@@ -21,6 +21,7 @@
 #include <ports-of-call/portability.hpp>
 #include <singularity-opac/base/opac_error.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
+#include <singularity-opac/photons/mean_photon_types.hpp>
 #include <singularity-opac/photons/photon_variant.hpp>
 #include <variant/include/mpark/variant.hpp>
 
@@ -92,6 +93,27 @@ class MeanVariant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION
+  Real AbsorptionCoefficient(const Real rho, const Real temp,
+                            const int gmode = Rosseland) const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          return opac.AbsorptionCoefficient(rho, temp, gmode);
+        },
+        opac_);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real Emissivity(const Real rho, const Real temp,
+                  const int gmode = Rosseland,
+                  Real *lambda = nullptr) const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          return opac.Emissivity(rho, temp, gmode, lambda);
+        },
+        opac_);
+  }
+
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
   }

singularity-opac/photons/non_cgs_photons.hpp

@@ -22,6 +22,7 @@
 
 #include <ports-of-call/portability.hpp>
 #include <singularity-opac/base/opac_error.hpp>
+#include <singularity-opac/photons/mean_photon_types.hpp>
 
 namespace singularity {
 namespace photons {
@@ -246,7 +247,9 @@ class MeanNonCGSUnits {
                   const Real temp_unit)
       : mean_opac_(std::forward<MeanOpac>(mean_opac)), time_unit_(time_unit),
         mass_unit_(mass_unit), length_unit_(length_unit), temp_unit_(temp_unit),
-        rho_unit_(mass_unit_ / (length_unit_ * length_unit_ * length_unit_)) {}
+        rho_unit_(mass_unit_ / (length_unit_ * length_unit_ * length_unit_)),
+        inv_emiss_unit_(length_unit_ * time_unit_ * time_unit_ * time_unit_ /
+          mass_unit_) {}
 
   auto GetOnDevice() {
     return MeanNonCGSUnits<MeanOpac>(mean_opac_.GetOnDevice(), time_unit_,
@@ -286,6 +289,21 @@ class MeanNonCGSUnits {
     return alpha * length_unit_;
   }
 
+  PORTABLE_INLINE_FUNCTION
+  Real AbsorptionCoefficient(const Real rho, const Real temp,
+                            const int gmode = Rosseland) const {
+    const Real alpha = mean_opac_.AbsorptionCoefficient(rho, temp, gmode);
+    return alpha * length_unit_;
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real Emissivity(const Real rho, const Real temp,
+                  const int gmode = Rosseland,
+                  Real *lambda = nullptr) const {
+    const Real J = mean_opac_.Emissivity(rho, temp, gmode);
+    return J * inv_emiss_unit_;
+  }
+
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
     return RuntimePhysicalConstants(PhysicalConstantsCGS(), time_unit_,
@@ -295,7 +313,7 @@ class MeanNonCGSUnits {
  private:
   MeanOpac mean_opac_;
   Real time_unit_, mass_unit_, length_unit_, temp_unit_;
-  Real rho_unit_;
+  Real rho_unit_, inv_emiss_unit_;
 };
 
 } // namespace photons

test/test_mean_opacities.cpp

@@ -815,6 +815,33 @@ TEST_CASE("ASCII-parsed Mean photon opacities", "[MeanPhotons]") {
         n_wrong_d() += 1;
       }
 
+      // test absorption and emission in Rossland and Planck modes
+      Real mabs_ross = mean_opac.AbsorptionCoefficient(rho_max, temp_max, 0);
+      Real mabs_plnk = mean_opac.AbsorptionCoefficient(rho_max, temp_max, 1);
+
+      // compare to Rossland and Planck
+      if (FractionalDifference(mross, mabs_ross) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+      if (FractionalDifference(mplnk, mabs_plnk) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+
+      Real *lambda = nullptr;
+      singularity::photons::PlanckDistribution<PhysicalConstantsCGS> dist;
+      Real B = dist.ThermalDistributionOfT(temp_max, lambda);
+
+      Real memiss_ross = mean_opac.Emissivity(rho_max, temp_max, 0);
+      Real memiss_plnk = mean_opac.Emissivity(rho_max, temp_max, 0);
+
+      // compare to Rossland and Planck
+      if (FractionalDifference(mross * B, memiss_ross) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+      if (FractionalDifference(mplnk * B, memiss_plnk) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+
 #ifdef PORTABILITY_STRATEGY_KOKKOS
       Kokkos::deep_copy(n_wrong_h, n_wrong_d);
 #endif