From 11ce39c3d461a87ff632f0b7fed3a4564e3ef4b4 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Thu, 9 Apr 2026 13:29:51 -0400 Subject: [PATCH 1/2] Fix thermal distribution bugs --- README.md | 13 ++- plan_histories/62.md | 18 ++++ .../thermal_distributions_neutrinos.hpp | 15 ++-- .../photons/thermal_distributions_photons.hpp | 16 ++-- test/CMakeLists.txt | 4 +- test/test_mean_opacities.cpp | 61 ++++++------- test/test_spiner_opac_neutrinos.cpp | 4 +- test/test_thermal_dist.cpp | 86 ++++++++++++++++++- 8 files changed, 170 insertions(+), 47 deletions(-) create mode 100644 plan_histories/62.md diff --git a/README.md b/README.md index 88c137d..72ade99 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,4 @@ + singularity-opac == @@ -21,7 +22,7 @@ For frequency-dependent absorption opacities, the following functions are provid | NumberEmissivity | $\int \frac{1}{h \nu} j_{\nu} d\Omega d\nu$ | Total number emissivity | ${\rm cm}^{-3}{\rm s}^{-1}$ | | ThermalDistributionOfTNu | $B_{\nu} = \frac{dE}{dA dt d\Omega d\nu}$ | Specific intensity of thermal distribution | ${\rm erg}{\rm cm}^{-2}{\rm s}^{-1}{\rm Sr}^{-1}{\rm Hz}^{-1}$ | | DThermalDistributionOfTNuDT | $dB_{\nu}/dT$ | Temperature derivative of specific intensity of thermal distribution | ${\rm erg}{\rm cm}^{-2}{\rm s}^{-1}{\rm Sr}^{-1}{\rm Hz}^{-1}{\rm K}^{-1}$ | -| ThermalDistributionOfT | $B = \int B_{\nu} d\Omega d\nu$ | Frequency- and angle-integrated intensity of thermal distribution | ${\rm erg}{\rm cm}^{-2}{\rm s}^{-1}$ | +| ThermalDistributionOfT | $4 \pi \int B_{\nu} d\nu = \int B_{\nu} d\Omega d\nu$ | Frequency- and angle-integrated intensity of thermal distribution | ${\rm erg}{\rm cm}^{-2}{\rm s}^{-1}$ | | ThermalNumberDistributionOfT | $B = \int \frac{1}{h \nu} B_{\nu} d\Omega d\nu$ | Frequency- and angle-integrated intensity of thermal distribution | ${\rm erg}{\rm cm}^{-2}{\rm s}^{-1}$ | | EnergyDensityFromTemperature | $E_{\rm R}$ | Radiation energy density | ${\rm erg}{\rm cm}^{-3}$ | | TemperatureFromEnergyDensity | $T_{\rm R}$ | Radiation temperature | ${\rm K}$ | @@ -82,7 +83,13 @@ with the following function signatures: PlanckMeanScatteringCoefficient(density, temperature) RosselandMeanScatteringCoefficient(density, temperature) -Note that the thermal radiation energy density `u = 1/c ThermalDistributionOfT` and the thermal radiation number density `n = 1/c ThermalNumberDistributionOfT`. +Note that `ThermalDistributionOfTNu` is the per-steradian Planck function `B_\nu`, so +`\int B_\nu d\nu = (c / 4 \pi) a T^4`, while +`ThermalDistributionOfT = 4 \pi \int B_\nu d\nu = c a T^4`. Therefore the +thermal radiation energy density is +`u = (4 \pi / c) \int B_\nu d\nu = (1 / c) ThermalDistributionOfT = a T^4`, and the +thermal radiation number density is +`n = (1 / c) ThermalNumberDistributionOfT`. Opacity variant constructors are specific to the opacity model being requested; consult the source code for individual opacities. @@ -150,7 +157,7 @@ max bounds the header. ## Copyright -© 2021. Triad National Security, LLC. All rights reserved. This +© 2021-2026. 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 diff --git a/plan_histories/62.md b/plan_histories/62.md new file mode 100644 index 0000000..397f88e --- /dev/null +++ b/plan_histories/62.md @@ -0,0 +1,18 @@ +## Goal + +Correct the photon and neutrino `dB_\nu/dT` implementations used for Rosseland +weighting and validate the fix. + +## Functional Plan + +- Correct the photon `dB_\nu/dT` formula to match the derivative of the Planck spectrum. +- Correct the analogous neutrino `dB_\nu/dT` formula for the Fermi-Dirac spectrum. +- Add focused thermal-distribution regression coverage: + - verify `4 \pi \int B_\nu d\nu = c a T^4` + - verify the expected `T^4` scaling + - verify `dB_\nu/dT` against a finite-difference reference + +## Scope Boundaries + +- Always produce comments in terms of a T^4 to prevent ambiguity surrounding + `ThermalDistributionFromT` definition. diff --git a/singularity-opac/neutrinos/thermal_distributions_neutrinos.hpp b/singularity-opac/neutrinos/thermal_distributions_neutrinos.hpp index c31b046..b66bc1e 100644 --- a/singularity-opac/neutrinos/thermal_distributions_neutrinos.hpp +++ b/singularity-opac/neutrinos/thermal_distributions_neutrinos.hpp @@ -1,5 +1,5 @@ // ====================================================================== -// © 2021. Triad National Security, LLC. All rights reserved. This +// © 2021-2026. 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. @@ -16,6 +16,8 @@ #ifndef SINGULARITY_OPAC_NEUTRINOS_THERMAL_DISTRIBUTIONS_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_THERMAL_DISTRIBUTIONS_NEUTRINOS_ +// This file was made in part with generative AI. + #include #include @@ -42,11 +44,12 @@ struct FermiDiracDistributionNoMu { Real DThermalDistributionOfTNuDT(const Real temp, const RadiationType type, const Real nu, Real *lambda = nullptr) const { - Real x = pc::h * nu / (pc::kb * temp); - Real dBnudT = NSPECIES * - (2. * pc::h * pc::h * nu * nu * nu * nu / - (temp * temp * pc::c * pc::c * pc::kb)) * - 1. / (std::exp(x) + 1.); + const Real x = pc::h * nu / (pc::kb * temp); + const Real expmx = std::exp(-x); + const Real denom = 1. + expmx; + const Real prefactor = NSPECIES * (2. * pc::h * pc::h * nu * nu * nu * nu / + (temp * temp * pc::c * pc::c * pc::kb)); + const Real dBnudT = prefactor * expmx / (denom * denom); return dBnudT; } PORTABLE_INLINE_FUNCTION diff --git a/singularity-opac/photons/thermal_distributions_photons.hpp b/singularity-opac/photons/thermal_distributions_photons.hpp index 1df7a4f..05d2565 100644 --- a/singularity-opac/photons/thermal_distributions_photons.hpp +++ b/singularity-opac/photons/thermal_distributions_photons.hpp @@ -1,5 +1,5 @@ // ====================================================================== -// © 2021. Triad National Security, LLC. All rights reserved. This +// © 2021-2026. 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. @@ -16,6 +16,8 @@ #ifndef SINGULARITY_OPAC_PHOTONS_THERMAL_DISTRIBUTIONS_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_THERMAL_DISTRIBUTIONS_PHOTONS_ +// This file was made in part with generative AI. + #include #include @@ -38,14 +40,18 @@ struct PlanckDistribution { PORTABLE_INLINE_FUNCTION Real DThermalDistributionOfTNuDT(const Real temp, const Real nu, Real *lambda = nullptr) const { - Real x = pc::h * nu / (pc::kb * temp); - Real dBnudT = (2. * pc::h * pc::h * nu * nu * nu * nu / - (temp * temp * pc::c * pc::c * pc::kb)) * - 1. / std::expm1(x); + const Real x = pc::h * nu / (pc::kb * temp); + const Real expmx = std::exp(-x); + const Real one_minus_expmx = -std::expm1(-x); + const Real prefactor = 2. * pc::h * pc::h * nu * nu * nu * nu / + (temp * temp * pc::c * pc::c * pc::kb); + const Real dBnudT = prefactor * expmx / (one_minus_expmx * one_minus_expmx); return dBnudT; } PORTABLE_INLINE_FUNCTION Real ThermalDistributionOfT(const Real temp, Real *lambda = nullptr) const { + // This is the 4 pi angle-integrated Planck intensity, not the + // per-steradian quantity (c / 4 pi) a T^4. return 8. * std::pow(M_PI, 5) * std::pow(pc::kb, 4) * std::pow(temp, 4) / (15. * std::pow(pc::c, 2) * std::pow(pc::h, 3)); } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 22d3860..27c4eee 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,4 +1,4 @@ -# © 2021. Triad National Security, LLC. All rights reserved. This +# © 2021-2026. 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 @@ -10,6 +10,8 @@ # prepare derivative works, distribute copies to the public, perform # publicly and display publicly, and to permit others to do so. +# This file was made in part with generative AI. + if(NOT Catch2_FOUND) # Get catch2 message(STATUS "Fetching Catch2 as needed") diff --git a/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp index 0e036c7..ab62cc3 100644 --- a/test/test_mean_opacities.cpp +++ b/test/test_mean_opacities.cpp @@ -1,5 +1,5 @@ // ====================================================================== -// © 2022-2024. Triad National Security, LLC. All rights reserved. This +// © 2022-2026. 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. @@ -13,6 +13,8 @@ // publicly, and to permit others to do so. // ====================================================================== +// This file was made in part with generative AI. + #include #include @@ -72,11 +74,11 @@ TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") { constexpr Real nu = 1.25 * MeV2Hz; // 1 MeV constexpr int nT = 10; - constexpr Real lRhoMin = std::log10(0.1 * rho); - constexpr Real lRhoMax = std::log10(10. * rho); + const Real lRhoMin = std::log10(0.1 * rho); + const Real lRhoMax = std::log10(10. * rho); constexpr int NRho = 2; - constexpr Real lTMin = std::log10(0.1 * temp); - constexpr Real lTMax = std::log10(10. * temp); + const Real lTMin = std::log10(0.1 * temp); + const Real lTMax = std::log10(10. * temp); constexpr int NT = 10; constexpr Real YeMin = 0.1; constexpr Real YeMax = 0.5; @@ -268,11 +270,11 @@ TEST_CASE("Mean neutrino scattering opacities", "[MeanNeutrinosS]") { constexpr Real nu = 1.25 * MeV2Hz; // 1 MeV constexpr int nT = 10; - constexpr Real lRhoMin = std::log10(0.1 * rho); - constexpr Real lRhoMax = std::log10(10. * rho); + const Real lRhoMin = std::log10(0.1 * rho); + const Real lRhoMax = std::log10(10. * rho); constexpr int NRho = 2; - constexpr Real lTMin = std::log10(0.1 * temp); - constexpr Real lTMax = std::log10(10. * temp); + const Real lTMin = std::log10(0.1 * temp); + const Real lTMax = std::log10(10. * temp); constexpr int NT = 10; constexpr Real YeMin = 0.1; constexpr Real YeMax = 0.5; @@ -289,7 +291,7 @@ TEST_CASE("Mean neutrino scattering opacities", "[MeanNeutrinosS]") { opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin, YeMax, NYe); auto mean_opac = mean_opac_host.GetOnDevice(); - THEN("The emissivity per nu omega is consistent with the emissity per nu") { + THEN("The gray Planck and Rosseland means reduce to rho kappa") { int n_wrong_h = 0; #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::View n_wrong_d("wrong"); @@ -435,11 +437,11 @@ TEST_CASE("Mean photon opacities", "[MeanPhotons]") { constexpr Real nu = 1.e15; // Hz constexpr int nT = 10; - constexpr Real lRhoMin = std::log10(0.1 * rho); - constexpr Real lRhoMax = std::log10(10. * rho); + const Real lRhoMin = std::log10(0.1 * rho); + const Real lRhoMax = std::log10(10. * rho); constexpr int NRho = 2; - constexpr Real lTMin = std::log10(0.1 * temp); - constexpr Real lTMax = std::log10(10. * temp); + const Real lTMin = std::log10(0.1 * temp); + const Real lTMax = std::log10(10. * temp); constexpr int NT = 10; constexpr Real kappa = 1.e-20; @@ -464,11 +466,11 @@ TEST_CASE("Mean photon opacities", "[MeanPhotons]") { Real alphaPlanck = mean_opac.PlanckMeanAbsorptionCoefficient(rho, temp); Real alphaRosseland = - mean_opac.PlanckMeanAbsorptionCoefficient(rho, temp); - if (FractionalDifference(kappa * rho, alphaPlanck) > EPS_TEST) { + mean_opac.RosselandMeanAbsorptionCoefficient(rho, temp); + if (FractionalDifference(kappa * rho, alphaPlanck) > 1.e-12) { n_wrong_d() += 1; } - if (FractionalDifference(kappa * rho, alphaRosseland) > EPS_TEST) { + if (FractionalDifference(kappa * rho, alphaRosseland) > 1.e-12) { n_wrong_d() += 1; } }); @@ -621,11 +623,11 @@ TEST_CASE("Mean photon scattering opacities", "[MeanPhotonS]") { constexpr Real nu = 1.e15; // Hz constexpr int nT = 10; - constexpr Real lRhoMin = std::log10(0.1 * rho); - constexpr Real lRhoMax = std::log10(10. * rho); + const Real lRhoMin = std::log10(0.1 * rho); + const Real lRhoMax = std::log10(10. * rho); constexpr int NRho = 2; - constexpr Real lTMin = std::log10(0.1 * temp); - constexpr Real lTMax = std::log10(10. * temp); + const Real lTMin = std::log10(0.1 * temp); + const Real lTMax = std::log10(10. * temp); constexpr int NT = 10; constexpr Real sigma = 1.e-20; @@ -771,14 +773,14 @@ TEST_CASE("ASCII-parsed Mean photon opacities", "[MeanPhotons]") { WHEN("We initialize a mean photon opacity from an ASCII table") { - constexpr Real rho_min = 1e-14; // g/cc. - constexpr Real temp_min = 1.0; // Kelvin. - constexpr Real ross_at_min = 0.001; // cm^2/g - constexpr Real plnk_at_min = 0.1; // cm^2/g - constexpr Real rho_max = 0.7943282347241912; // g/cc. + constexpr Real rho_min = 1e-14; // g/cc. + constexpr Real temp_min = 1.0; // Kelvin. + constexpr Real ross_at_min = 0.001; // cm^2/g + constexpr Real plnk_at_min = 0.1; // cm^2/g + constexpr Real rho_max = 0.7943282347241912; // g/cc. constexpr Real temp_max = 7943282.347242886; // Kelvin. - constexpr Real ross_at_max = 0.001; // cm^2/g - constexpr Real plnk_at_max = 0.1; // cm^2/g + constexpr Real ross_at_max = 0.001; // cm^2/g + constexpr Real plnk_at_max = 0.1; // cm^2/g photons::MeanOpacity mean_opac_host = photons::MeanOpacity(grayname); auto mean_opac = mean_opac_host.GetOnDevice(); @@ -792,7 +794,8 @@ TEST_CASE("ASCII-parsed Mean photon opacities", "[MeanPhotons]") { #endif // get a test value at min rho-T point - Real mross = mean_opac.RosselandMeanAbsorptionCoefficient(rho_min, temp_min); + Real mross = + mean_opac.RosselandMeanAbsorptionCoefficient(rho_min, temp_min); Real mplnk = mean_opac.PlanckMeanAbsorptionCoefficient(rho_min, temp_min); // check min rho-T point diff --git a/test/test_spiner_opac_neutrinos.cpp b/test/test_spiner_opac_neutrinos.cpp index a927459..13f442a 100644 --- a/test/test_spiner_opac_neutrinos.cpp +++ b/test/test_spiner_opac_neutrinos.cpp @@ -60,8 +60,8 @@ TEST_CASE("Spiner opacities, filled with gray data", constexpr Real lRhoMin = 8; constexpr Real lRhoMax = 12; constexpr int NRho = 32; - constexpr Real lTMin = -2 + std::log10(MeV2K); - constexpr Real lTMax = 2 + std::log10(MeV2K); + const Real lTMin = -2 + std::log10(MeV2K); + const Real lTMax = 2 + std::log10(MeV2K); constexpr int NT = 32; constexpr Real YeMin = 0.1; constexpr Real YeMax = 0.5; diff --git a/test/test_thermal_dist.cpp b/test/test_thermal_dist.cpp index af80b00..32b6cca 100644 --- a/test/test_thermal_dist.cpp +++ b/test/test_thermal_dist.cpp @@ -1,5 +1,5 @@ // ====================================================================== -// © 2022. Triad National Security, LLC. All rights reserved. This +// © 2022-2026. 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. @@ -13,6 +13,9 @@ // publicly, and to permit others to do so. // ====================================================================== +// This file was made in part with generative AI. + +#include #include #include @@ -41,6 +44,21 @@ template PORTABLE_INLINE_FUNCTION T FractionalDifference(const T &a, const T &b) { return 2 * std::abs(b - a) / (std::abs(a) + std::abs(b) + 1e-20); } +template +Real IntegrateOverLogNu(const Function &f, const Real nu_min, const Real nu_max, + const int nnu) { + const Real lnu_min = std::log(nu_min); + const Real lnu_max = std::log(nu_max); + const Real dlnu = (lnu_max - lnu_min) / (nnu - 1); + Real integral = 0.; + for (int inu = 0; inu < nnu; ++inu) { + const Real weight = (inu == 0 || inu == nnu - 1) ? 0.5 : 1.; + const Real lnu = lnu_min + inu * dlnu; + const Real nu = std::exp(lnu); + integral += weight * f(nu) * nu * dlnu; + } + return integral; +} constexpr Real EPS_TEST = 1e-3; TEST_CASE("Neutrino thermal distribution", "[NeutrinoThermalDistribution]") { @@ -122,3 +140,69 @@ TEST_CASE("Photon thermal distribution", "[PhotonThermalDistribution]") { opac.Finalize(); } } + +TEST_CASE("Photon thermal distribution integral scales as T^4", + "[PhotonThermalDistribution]") { + constexpr int nnu = 2049; + const std::array temps = {1.e4, 3.e5, 2.e6}; + photons::PlanckDistribution dist; + + std::array integrated_B{}; + for (int i = 0; i < static_cast(temps.size()); ++i) { + const Real temp = temps[i]; + const Real nu_scale = pc::kb * temp / pc::h; + // In x = h nu / (k_B T), these bounds are x in [1e-6, 1e2], which is + // effectively [0, infinity): the low-frequency tail is negligible as + // x -> 0, and the high-frequency tail is exponentially suppressed. + const Real nu_min = 1.e-6 * nu_scale; + const Real nu_max = 1.e2 * nu_scale; + + const Real Bnu_integral = IntegrateOverLogNu( + [&](const Real nu) { return dist.ThermalDistributionOfTNu(temp, nu); }, + nu_min, nu_max, nnu); + integrated_B[i] = 4. * M_PI * Bnu_integral; + + // singularity-opac uses the angle-integrated quantity + // + // ThermalDistributionOfT = 4 pi \int B_nu dnu. + // + // The radiation energy density E_r = a T^4. Since + // E_r = (1 / c) * ThermalDistributionOfT, the corresponding analytic + // value here is + // + // ThermalDistributionOfT = c a T^4. + // + // The per-steradian integral is \int B_nu dnu = (c / 4 pi) a T^4. + const Real B_from_dist = dist.ThermalDistributionOfT(temp); + const Real B_analytic = pc::ar * pc::c * std::pow(temp, 4); + const Real B_analytic_per_sr = + pc::ar * pc::c * std::pow(temp, 4) / (4. * M_PI); + REQUIRE(FractionalDifference(integrated_B[i], B_from_dist) < 1.e-8); + REQUIRE(FractionalDifference(integrated_B[i], B_analytic) < 1.e-8); + REQUIRE(FractionalDifference(Bnu_integral, B_analytic_per_sr) < 1.e-8); + } + + const Real expected_ratio = std::pow(temps[2] / temps[0], 4); + REQUIRE(FractionalDifference(integrated_B[2] / integrated_B[0], + expected_ratio) < 1.e-8); +} + +TEST_CASE("Photon dBnu/dT matches a finite difference derivative", + "[PhotonThermalDistribution]") { + constexpr Real temp = 2.e6; + constexpr Real rel_step = 1.e-6; + const std::array x_values = {1.e-3, 1.e-1, 1., 3.}; + photons::PlanckDistribution dist; + + for (const Real x : x_values) { + const Real nu = x * pc::kb * temp / pc::h; + const Real temp_lo = temp * (1. - rel_step); + const Real temp_hi = temp * (1. + rel_step); + const Real dBnudT_fd = (dist.ThermalDistributionOfTNu(temp_hi, nu) - + dist.ThermalDistributionOfTNu(temp_lo, nu)) / + (temp_hi - temp_lo); + const Real dBnudT = dist.DThermalDistributionOfTNuDT(temp, nu); + + REQUIRE(FractionalDifference(dBnudT, dBnudT_fd) < 1.e-6); + } +} From b080ef47e8939429e2e0e582d85041457f1e4768 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Thu, 9 Apr 2026 13:31:38 -0400 Subject: [PATCH 2/2] Correct PR number --- plan_histories/{62.md => 69.md} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename plan_histories/{62.md => 69.md} (100%) diff --git a/plan_histories/62.md b/plan_histories/69.md similarity index 100% rename from plan_histories/62.md rename to plan_histories/69.md