From fa652e86d22528b8856accb5a77971bbeb3a7372 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 17 Dec 2024 16:37:37 -0700 Subject: [PATCH 1/4] Add AbsorptionCoefficient and Emissivity to MeanOpacity variants. + Use gmode input to toggle Rosseland (default) to Planck. + Add PlanckDistribution function as member to MeanOpacity. + Add the new functions to generic and non-CGS MeanOpacity variants. Note: the motivation is to use one of the frequency-integrated opacity values as a grey opacity in transport, which requires using it in both absorption and emissivity, for any choice. --- .../photons/mean_opacity_photons.hpp | 21 +++++++++++++++ .../photons/mean_photon_variant.hpp | 21 +++++++++++++++ singularity-opac/photons/non_cgs_photons.hpp | 20 ++++++++++++-- test/test_mean_opacities.cpp | 27 +++++++++++++++++++ 4 files changed, 87 insertions(+), 2 deletions(-) diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp index e9bdd02..f98dbfe 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, @@ -307,6 +327,7 @@ class MeanOpacity { Rosseland = 0, Planck = 1 }; + PlanckDistribution dist_; }; #undef EPS diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp index ff860b0..d56e705 100644 --- a/singularity-opac/photons/mean_photon_variant.hpp +++ b/singularity-opac/photons/mean_photon_variant.hpp @@ -92,6 +92,27 @@ class MeanVariant { opac_); } + PORTABLE_INLINE_FUNCTION + Real AbsorptionCoefficient(const Real rho, const Real temp, + const int gmode = 0) 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 = 0, + 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..6b62cb0 100644 --- a/singularity-opac/photons/non_cgs_photons.hpp +++ b/singularity-opac/photons/non_cgs_photons.hpp @@ -246,7 +246,8 @@ 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_ / mass_unit_) {} auto GetOnDevice() { return MeanNonCGSUnits(mean_opac_.GetOnDevice(), time_unit_, @@ -286,6 +287,21 @@ class MeanNonCGSUnits { return alpha * length_unit_; } + PORTABLE_INLINE_FUNCTION + Real AbsorptionCoefficient(const Real rho, const Real temp, + const int gmode = 0) 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 = 0, + Real *lambda = nullptr) const { + const Real J = mean_opac_.Emissivity(rho, temp, gmode); + return J * inv_emiss_unit_ * time_unit_; + } + PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants GetRuntimePhysicalConstants() const { return RuntimePhysicalConstants(PhysicalConstantsCGS(), time_unit_, @@ -295,7 +311,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 From 2357ffb98679cb4d9fe3e8c60717710cb2801703 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Thu, 19 Dec 2024 16:05:01 -0700 Subject: [PATCH 2/4] Use Rosseland/Planck enum in photon MeanVariant. + Move enum from MeanOpacity class to mean_photon_variant.hpp. + Use Rosseland instead of 0 in variant template function defaults. --- singularity-opac/photons/mean_opacity_photons.hpp | 4 ---- singularity-opac/photons/mean_photon_variant.hpp | 10 ++++++++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp index f98dbfe..9e97557 100644 --- a/singularity-opac/photons/mean_opacity_photons.hpp +++ b/singularity-opac/photons/mean_opacity_photons.hpp @@ -323,10 +323,6 @@ class MeanOpacity { } Spiner::DataBox lkappa_; const char *filename_; - enum OpacityAveraging { - Rosseland = 0, - Planck = 1 - }; PlanckDistribution dist_; }; diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp index d56e705..63f7ec4 100644 --- a/singularity-opac/photons/mean_photon_variant.hpp +++ b/singularity-opac/photons/mean_photon_variant.hpp @@ -28,6 +28,12 @@ namespace singularity { namespace photons { namespace impl { +// mean-opacity mode +enum OpacityAveraging { + Rosseland = 0, + Planck = 1 +}; + template class MeanVariant { private: @@ -94,7 +100,7 @@ class MeanVariant { PORTABLE_INLINE_FUNCTION Real AbsorptionCoefficient(const Real rho, const Real temp, - const int gmode = 0) const { + const int gmode = Rosseland) const { return mpark::visit( [=](const auto &opac) { return opac.AbsorptionCoefficient(rho, temp, gmode); @@ -104,7 +110,7 @@ class MeanVariant { PORTABLE_INLINE_FUNCTION Real Emissivity(const Real rho, const Real temp, - const int gmode = 0, + const int gmode = Rosseland, Real *lambda = nullptr) const { return mpark::visit( [=](const auto &opac) { From 08dfff7e794685acfa8f44c67b8dc02bac4ffb37 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 7 Jan 2025 18:31:34 -0700 Subject: [PATCH 3/4] Move mean opacity enum to its own header: mean_photon_types.hpp. --- .../photons/mean_photon_types.hpp | 31 +++++++++++++++++++ .../photons/mean_photon_variant.hpp | 7 +---- singularity-opac/photons/non_cgs_photons.hpp | 5 +-- 3 files changed, 35 insertions(+), 8 deletions(-) create mode 100644 singularity-opac/photons/mean_photon_types.hpp 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 63f7ec4..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 @@ -28,12 +29,6 @@ namespace singularity { namespace photons { namespace impl { -// mean-opacity mode -enum OpacityAveraging { - Rosseland = 0, - Planck = 1 -}; - template class MeanVariant { private: diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp index 6b62cb0..4d29001 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 { @@ -289,14 +290,14 @@ class MeanNonCGSUnits { PORTABLE_INLINE_FUNCTION Real AbsorptionCoefficient(const Real rho, const Real temp, - const int gmode = 0) const { + 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 = 0, + const int gmode = Rosseland, Real *lambda = nullptr) const { const Real J = mean_opac_.Emissivity(rho, temp, gmode); return J * inv_emiss_unit_ * time_unit_; From df1dd30f7ddac253438dfca53f8bc18701b5eee8 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 7 Jan 2025 18:42:07 -0700 Subject: [PATCH 4/4] Roll time unit into mean non-cgs inverse emissivity unit. --- singularity-opac/photons/non_cgs_photons.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp index 4d29001..36aff7e 100644 --- a/singularity-opac/photons/non_cgs_photons.hpp +++ b/singularity-opac/photons/non_cgs_photons.hpp @@ -248,7 +248,8 @@ class MeanNonCGSUnits { : 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_)), - inv_emiss_unit_(length_unit_ * time_unit_ * time_unit_ / mass_unit_) {} + inv_emiss_unit_(length_unit_ * time_unit_ * time_unit_ * time_unit_ / + mass_unit_) {} auto GetOnDevice() { return MeanNonCGSUnits(mean_opac_.GetOnDevice(), time_unit_, @@ -300,7 +301,7 @@ class MeanNonCGSUnits { const int gmode = Rosseland, Real *lambda = nullptr) const { const Real J = mean_opac_.Emissivity(rho, temp, gmode); - return J * inv_emiss_unit_ * time_unit_; + return J * inv_emiss_unit_; } PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants