From fa652e86d22528b8856accb5a77971bbeb3a7372 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 17 Dec 2024 16:37:37 -0700 Subject: [PATCH] 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