Fix thermal distribution bugs

README.md

@@ -1,3 +1,4 @@
+<!-- This file was made in part with generative AI. -->
 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

plan_histories/69.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.

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 <cmath>
 
 #include <ports-of-call/portability.hpp>
@@ -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

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 <cmath>
 
 #include <ports-of-call/portability.hpp>
@@ -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));
   }

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")

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 <cmath>
 #include <iostream>
 
@@ -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<int, atomic_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

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;