diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp index e9bdd02..9e97557 100644 --- a/singularity-opac/photons/mean_opacity_photons.hpp +++ b/singularity-opac/photons/mean_opacity_photons.hpp @@ -29,6 +29,7 @@ #include #include +#include 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 void MeanOpacityImpl_(const Opacity &opac, const Real lRhoMin, @@ -303,10 +323,7 @@ class MeanOpacity { } Spiner::DataBox lkappa_; const char *filename_; - enum OpacityAveraging { - Rosseland = 0, - Planck = 1 - }; + PlanckDistribution dist_; }; #undef EPS diff --git a/singularity-opac/photons/mean_photon_types.hpp b/singularity-opac/photons/mean_photon_types.hpp new file mode 100644 index 0000000..763a3af --- /dev/null +++ b/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_ diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp index ff860b0..a307d0f 100644 --- a/singularity-opac/photons/mean_photon_variant.hpp +++ b/singularity-opac/photons/mean_photon_variant.hpp @@ -21,6 +21,7 @@ #include #include #include +#include #include #include @@ -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_); } diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp index 36857b0..36aff7e 100644 --- a/singularity-opac/photons/non_cgs_photons.hpp +++ b/singularity-opac/photons/non_cgs_photons.hpp @@ -22,6 +22,7 @@ #include #include +#include namespace singularity { namespace photons { @@ -246,7 +247,9 @@ class MeanNonCGSUnits { const Real temp_unit) : mean_opac_(std::forward(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(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 diff --git a/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp index 5205f4c..ba1d371 100644 --- a/test/test_mean_opacities.cpp +++ b/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 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