diff --git a/README.md b/README.md index 72ade99..7e091e4 100644 --- a/README.md +++ b/README.md @@ -44,20 +44,56 @@ with the following function signatures: TemperatureFromEnergyDensity(radiation energy density) NumberDensityFromTemperature(temperature) -For mean absorption opacities, the following functions are provided: +For mean absorption opacities, the following functions are provided. +Mean opacities compute frequency-averaged (Planck or Rosseland) absorption +coefficients over one or more energy groups. When `ngroups=1` with group bounds +`[0, ∞)`, mean opacities reduce to traditional gray opacities. For `ngroups>1`, +they provide multigroup radiation transport capabilities. + | Function | Expression | Description | Units | | --------------------- | ---------- | --------------------- | ------- | -| PlankMeanAbsorptionCoefficient | $n \sigma$ | Absorption coefficient | ${\rm cm}^{-1}$ | -| RosselandMeanAbsorptionCoefficient | $n \sigma$ | Absorption coefficient | ${\rm cm}^{-1}$ | -| AbsorptionCoefficient | $n \sigma$ | Absorption coefficient | ${\rm cm}^{-1}$ | +| PlankMeanAbsorptionCoefficient | $n \sigma$ | Planck mean absorption coefficient (gray, ngroups=1) | ${\rm cm}^{-1}$ | +| RosselandMeanAbsorptionCoefficient | $n \sigma$ | Rosseland mean absorption coefficient (gray, ngroups=1) | ${\rm cm}^{-1}$ | +| PlanckGroupAbsorptionCoefficient | $n \sigma_g$ | Planck-weighted absorption coefficient in a frequency group | ${\rm cm}^{-1}$ | +| RosselandGroupAbsorptionCoefficient | $n \sigma_g$ | Rosseland-weighted absorption coefficient in a frequency group | ${\rm cm}^{-1}$ | +| AbsorptionCoefficient | $n \sigma$ or $n \sigma_g$ | Absorption coefficient (gray or group) | ${\rm cm}^{-1}$ | | Emissivity | $\int j_{\nu} d\nu d\Omega$ | Total emissivity | ${\rm erg}{\rm cm}^{-3}{\rm s}^{-1}$ | +| GroupOfNu | $g(\nu)$ | Group index containing a given frequency | dimensionless | +| PlanckGroupAbsorptionCoefficientFromNu | $n \sigma_{g(\nu)}$ | Planck-weighted group coefficient selected by frequency | ${\rm cm}^{-1}$ | +| RosselandGroupAbsorptionCoefficientFromNu | $n \sigma_{g(\nu)}$ | Rosseland-weighted group coefficient selected by frequency | ${\rm cm}^{-1}$ | +| AbsorptionCoefficientFromNu | $n \sigma_{g(\nu)}$ | Group coefficient selected by frequency | ${\rm cm}^{-1}$ | with the following function signatures: + ngroups() PlanckMeanAbsorptionCoefficient(density, temperature) RosselandMeanAbsorptionCoefficient(density, temperature) + PlanckGroupAbsorptionCoefficient(density, temperature, group index) + RosselandGroupAbsorptionCoefficient(density, temperature, group index) AbsorptionCoefficient(density, temperature, gmode [Planck, Rosseland]) + AbsorptionCoefficient(density, temperature, group index, gmode [Planck, Rosseland]) Emissivity(density, temperature) + GroupOfNu(frequency) + PlanckGroupAbsorptionCoefficientFromNu(density, temperature, frequency) + RosselandGroupAbsorptionCoefficientFromNu(density, temperature, frequency) + AbsorptionCoefficientFromNu(density, temperature, frequency, gmode [Planck, Rosseland]) + +Mean absorption opacities can either be constructed from a frequency-dependent +opacity model plus user-provided group bounds, or loaded directly from +precomputed Spiner tables of group Planck and Rosseland opacities. The direct +table-backed path does not require singularity-opac to recompute opacity +integrals, but it must include explicit group bounds. Mean opacities +always carry group bounds, and the convention is half-open groups +`[nu_g, nu_{g+1})`, with the final upper bound included in the last group. +If you want singularity-opac to interpret a group as extending to infinity, +then the final entry of the `group bounds` array must literally be IEEE +positive infinity, not just a very large finite cutoff. In C++, that means +using `std::numeric_limits::infinity()`, while in Python that would +typically be `float("inf")` or `numpy.inf`. In an HDF/Spiner table, that same +IEEE `+infinity` value must appear in the final element of the `"group bounds"` +dataset. Likewise, a lower tail group `[0, nu_1)` is represented by setting +the first bound to `0.`. A very large finite number is still interpreted as a +finite bound. For frequency-dependent scattering opacities, the following functions are provided | Function | Expression | Description | Units | @@ -141,19 +177,6 @@ A number of options are avaialable for compiling: | SINGULARITY_USE_HDF5 | ON | Enables HDF5. Required for Spiner opacities. | | SINGULARITY_KOKKOS_IN_TREE | OFF | Force cmake to use Kokkos source included in tree. | -### Loading ASCII Data - -Currently, the MeanOpacity class, defined in singularity-opac/photons/mean_opacity_photons.hpp, supports -loading grey Rosseland and Planck opacity data in an ASCII format. An example of this format is -provided by singularity-opac/photons/example_ascii/kap_plaw.txt. The 1st row of the header has the -number of density points, NRho, then the number of temperature points, NT. The 2nd (3rd) row of the header -has min and max density (temperature) bounds. These bounds are inclusive, so the opacity data in the file -should have evaluations at these min and max values. The rest of the ASCII file is a two-column table, where -the 1st (2nd) column is Rosseland (Planck) opacity. The number of rows in each column is NRhoxNT, where -density is the slow index and temperature is the fast index (thus the row index = temperature index -+ NT x (density index), indexing from 0). Each opacity is assumed to be evaluated on log-spaced -density and temperature grids, where these grids are defined by NRho, NT, and the (again inclusive) min and -max bounds the header. ## Copyright diff --git a/plan_histories/69.md b/plan_histories/69.md index 397f88e..b019a3d 100644 --- a/plan_histories/69.md +++ b/plan_histories/69.md @@ -1,3 +1,5 @@ + + ## Goal Correct the photon and neutrino `dB_\nu/dT` implementations used for Rosseland diff --git a/plan_histories/70.md b/plan_histories/70.md index e760790..ae95d84 100644 --- a/plan_histories/70.md +++ b/plan_histories/70.md @@ -1,3 +1,5 @@ + + ## Goal Extend photon `PowerLaw` opacities to support an optional frequency power law diff --git a/plan_histories/71.md b/plan_histories/71.md new file mode 100644 index 0000000..fc604d3 --- /dev/null +++ b/plan_histories/71.md @@ -0,0 +1,38 @@ + + +## Goal + +Add a photon multigroup absorption API that supports group-index lookups, +frequency-to-group lookup, direct table-backed loading, and non-CGS wrapping. + +## Functional Plan + +- Add a photon multigroup opacity type that stores Planck and Rosseland group + absorption coefficients on `(rho, T, group)`. +- Support two construction paths: + - build group means from a monochromatic photon opacity plus group bounds + - load pretabulated group Planck and Rosseland opacities from Spiner data +- Keep group bounds as first-class multigroup data: + - preserve half-open group semantics `[nu_g, nu_{g+1})` + - allow a lower tail group by setting the first bound to `0.` + - allow an upper tail group by setting the final bound to IEEE `+infinity` +- Add explicit frequency-entry APIs: + - `GroupOfNu` + - `PlanckGroupAbsorptionCoefficientFromNu` + - `RosselandGroupAbsorptionCoefficientFromNu` + - `AbsorptionCoefficientFromNu` +- Keep the `NonCGSUnits` wrapper working for multigroup photon opacities. +- Add focused regression coverage for: + - gray exactness across groups + - autogenerated logarithmic groups with tail groups + - direct table-backed construction + - SP5/HDF round-trip + - half-open frequency-bound behavior + - non-CGS wrapper behavior +- Unify multigroup into existing mean opacity classes rather than creating separate types. + +## Scope Boundaries + +- Do not add neutrino multigroup support in this PR (coming in later PR). +- Do not add integrated multigroup emissivity in this PR (add in later PR, as downstream + codes may independently handle their emissivity integrals). diff --git a/singularity-opac/base/indexers.hpp b/singularity-opac/base/indexers.hpp index c642668..f173bdb 100644 --- a/singularity-opac/base/indexers.hpp +++ b/singularity-opac/base/indexers.hpp @@ -12,13 +12,12 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_BASE_INDEXERS_ #define SINGULARITY_OPAC_BASE_INDEXERS_ #include -#include #include +#include #include @@ -59,8 +58,8 @@ class Linear { SetRange_(numin, numax, N); } - PORTABLE_INLINE_FUNCTION Linear(const DataBox &data, Real numin, - Real numax, int N) + PORTABLE_INLINE_FUNCTION Linear(const DataBox &data, Real numin, Real numax, + int N) : data_(data) { SetRange_(numin, numax, N); } @@ -99,8 +98,7 @@ class LogLinear { } PORTABLE_INLINE_FUNCTION - LogLinear(const DataBox &data, Real numin, Real numax, int N) - : data_(data) { + LogLinear(const DataBox &data, Real numin, Real numax, int N) : data_(data) { SetRange_(numin, numax, N); } diff --git a/singularity-opac/base/opac_error.hpp b/singularity-opac/base/opac_error.hpp index 0d320a6..4b92624 100644 --- a/singularity-opac/base/opac_error.hpp +++ b/singularity-opac/base/opac_error.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_BASE_OPAC_ERROR_ #define SINGULARITY_OPAC_BASE_OPAC_ERROR_ diff --git a/singularity-opac/base/radiation_types.hpp b/singularity-opac/base/radiation_types.hpp index 9ce1add..6feb27a 100644 --- a/singularity-opac/base/radiation_types.hpp +++ b/singularity-opac/base/radiation_types.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_BASE_RADIATION_TYPES_ #define SINGULARITY_OPAC_BASE_RADIATION_TYPES_ diff --git a/singularity-opac/base/robust_utils.hpp b/singularity-opac/base/robust_utils.hpp index 470e544..7d8889e 100644 --- a/singularity-opac/base/robust_utils.hpp +++ b/singularity-opac/base/robust_utils.hpp @@ -11,7 +11,6 @@ // prepare derivative works, distribute copies to the public, perform // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ - #ifndef SINGULARITY_OPAC_BASE_ROBUST_UTILS_HPP_ #define SINGULARITY_OPAC_BASE_ROBUST_UTILS_HPP_ @@ -19,39 +18,39 @@ #include namespace singularity_opac { - namespace robust { - - template - PORTABLE_FORCEINLINE_FUNCTION constexpr auto SMALL() { - return 10 * std::numeric_limits::min(); - } - - template - PORTABLE_FORCEINLINE_FUNCTION constexpr auto EPS() { - return 10 * std::numeric_limits::epsilon(); - } - - template - PORTABLE_FORCEINLINE_FUNCTION auto make_positive(const T val) { - return std::max(val, EPS()); - } - - PORTABLE_FORCEINLINE_FUNCTION - Real make_bounded(const Real val, const Real vmin, const Real vmax) { - return std::min(std::max(val, vmin + EPS()), vmax * (1.0 - EPS())); - } - - template - PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val) { - return (T(0) <= val) - (val < T(0)); - } - - template - PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) { - return a / (b + sgn(b) * SMALL()); - } - - } // namespace robust +namespace robust { + +template +PORTABLE_FORCEINLINE_FUNCTION constexpr auto SMALL() { + return 10 * std::numeric_limits::min(); +} + +template +PORTABLE_FORCEINLINE_FUNCTION constexpr auto EPS() { + return 10 * std::numeric_limits::epsilon(); +} + +template +PORTABLE_FORCEINLINE_FUNCTION auto make_positive(const T val) { + return std::max(val, EPS()); +} + +PORTABLE_FORCEINLINE_FUNCTION +Real make_bounded(const Real val, const Real vmin, const Real vmax) { + return std::min(std::max(val, vmin + EPS()), vmax * (1.0 - EPS())); +} + +template +PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val) { + return (T(0) <= val) - (val < T(0)); +} + +template +PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) { + return a / (b + sgn(b) * SMALL()); +} + +} // namespace robust } // namespace singularity_opac #endif // SINGULARITY_OPAC_BASE_ROBUST_UTILS_HPP_ diff --git a/singularity-opac/base/sp5.hpp b/singularity-opac/base/sp5.hpp index eea9d63..540ae04 100644 --- a/singularity-opac/base/sp5.hpp +++ b/singularity-opac/base/sp5.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. @@ -12,9 +12,9 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_BASE_SP5_ #define SINGULARITY_OPAC_BASE_SP5_ +// This file was made in part with generative AI. namespace SP5 { @@ -38,6 +38,18 @@ constexpr char PlanckMeanSOpacity[] = "Planck mean scattering opacity"; constexpr char RosselandMeanSOpacity[] = "Rosseland mean scattering opacity"; } // namespace MeanSOpac +namespace MultigroupOpac { +constexpr char PlanckGroupOpacity[] = "Planck group opacity"; +constexpr char RosselandGroupOpacity[] = "Rosseland group opacity"; +constexpr char GroupBounds[] = "group bounds"; +} // namespace MultigroupOpac + +namespace MultigroupSOpac { +constexpr char PlanckGroupSOpacity[] = "Planck group scattering opacity"; +constexpr char RosselandGroupSOpacity[] = "Rosseland group scattering opacity"; +constexpr char GroupBounds[] = "group bounds"; +} // namespace MultigroupSOpac + } // namespace SP5 #endif // SINGULARITY_OPAC_BASE_SP5_ diff --git a/singularity-opac/chebyshev/chebyshev.hpp b/singularity-opac/chebyshev/chebyshev.hpp index 3045033..0aeeeea 100644 --- a/singularity-opac/chebyshev/chebyshev.hpp +++ b/singularity-opac/chebyshev/chebyshev.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_CHEBYSHEV_CHEBYSHEV_ #define SINGULARITY_OPAC_CHEBYSHEV_CHEBYSHEV_ diff --git a/singularity-opac/chebyshev/vandermonde.hpp b/singularity-opac/chebyshev/vandermonde.hpp index 4d61569..e93050f 100644 --- a/singularity-opac/chebyshev/vandermonde.hpp +++ b/singularity-opac/chebyshev/vandermonde.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_CHEBYSHEV_VANDERMONDE_ #define SINGULARITY_OPAC_CHEBYSHEV_VANDERMONDE_ @@ -800,44 +799,41 @@ constexpr Real Vandermonde21[21][21] = { -0.02119246990060137183703833947588521518727195892995666241, 0.00711715177013564326580378530807301460700464334926828038}}; - PORTABLE_INLINE_FUNCTION void get_vmbox(Real *r) - { - constexpr Real Vandermonde9[9][9] = { - {0.11111111111111, 0.11111111111111, 0.11111111111111, 0.11111111111111, - 0.11111111111111, 0.11111111111111, 0.11111111111111, 0.11111111111111, - 0.11111111111111}, - {-0.21884616733605, -0.19245008972988, -0.14284169104145, - -0.076004476294593, 0, 0.076004476294593, 0.14284169104145, - 0.19245008972988, 0.21884616733605}, - {0.20882058239687, 0.11111111111111, -0.038588483925985, -0.17023209847088, - -0.22222222222222, -0.17023209847088, -0.038588483925985, 0.11111111111111, - 0.20882058239687}, - {-0.19245008972988, 0, 0.19245008972988, 0.19245008972988, 0, - -0.19245008972988, -0.19245008972988, 0, 0.19245008972988}, - {0.17023209847088, -0.11111111111111, -0.20882058239687, 0.038588483925985, - 0.22222222222222, 0.038588483925985, -0.20882058239687, -0.11111111111111, - 0.17023209847088}, - {-0.14284169104145, 0.19245008972988, 0.076004476294593, -0.21884616733605, - 0, 0.21884616733605, -0.076004476294593, -0.19245008972988, - 0.14284169104145}, - {0.11111111111111, -0.22222222222222, 0.11111111111111, 0.11111111111111, - -0.22222222222222, 0.11111111111111, 0.11111111111111, -0.22222222222222, - 0.11111111111111}, - {-0.076004476294593, 0.19245008972988, -0.21884616733605, 0.14284169104145, - 0, -0.14284169104145, 0.21884616733605, -0.19245008972988, - 0.076004476294593}, - {0.038588483925985, -0.11111111111111, 0.17023209847088, -0.20882058239687, - 0.22222222222222, -0.20882058239687, 0.17023209847088, -0.11111111111111, - 0.038588483925985}}; +PORTABLE_INLINE_FUNCTION void get_vmbox(Real *r) { + constexpr Real Vandermonde9[9][9] = { + {0.11111111111111, 0.11111111111111, 0.11111111111111, 0.11111111111111, + 0.11111111111111, 0.11111111111111, 0.11111111111111, 0.11111111111111, + 0.11111111111111}, + {-0.21884616733605, -0.19245008972988, -0.14284169104145, + -0.076004476294593, 0, 0.076004476294593, 0.14284169104145, + 0.19245008972988, 0.21884616733605}, + {0.20882058239687, 0.11111111111111, -0.038588483925985, + -0.17023209847088, -0.22222222222222, -0.17023209847088, + -0.038588483925985, 0.11111111111111, 0.20882058239687}, + {-0.19245008972988, 0, 0.19245008972988, 0.19245008972988, 0, + -0.19245008972988, -0.19245008972988, 0, 0.19245008972988}, + {0.17023209847088, -0.11111111111111, -0.20882058239687, + 0.038588483925985, 0.22222222222222, 0.038588483925985, + -0.20882058239687, -0.11111111111111, 0.17023209847088}, + {-0.14284169104145, 0.19245008972988, 0.076004476294593, + -0.21884616733605, 0, 0.21884616733605, -0.076004476294593, + -0.19245008972988, 0.14284169104145}, + {0.11111111111111, -0.22222222222222, 0.11111111111111, 0.11111111111111, + -0.22222222222222, 0.11111111111111, 0.11111111111111, -0.22222222222222, + 0.11111111111111}, + {-0.076004476294593, 0.19245008972988, -0.21884616733605, + 0.14284169104145, 0, -0.14284169104145, 0.21884616733605, + -0.19245008972988, 0.076004476294593}, + {0.038588483925985, -0.11111111111111, 0.17023209847088, + -0.20882058239687, 0.22222222222222, -0.20882058239687, 0.17023209847088, + -0.11111111111111, 0.038588483925985}}; - for(int i = 0; i < 9; ++i) - { - for (int j = 0; j < 9; ++j) - { - r[j + i * 9] = Vandermonde9[i][j]; - } - } + for (int i = 0; i < 9; ++i) { + for (int j = 0; j < 9; ++j) { + r[j + i * 9] = Vandermonde9[i][j]; + } } +} } // namespace chebyshev } // namespace singularity diff --git a/singularity-opac/constants/constants.hpp b/singularity-opac/constants/constants.hpp index 24fd739..8c4ddcd 100644 --- a/singularity-opac/constants/constants.hpp +++ b/singularity-opac/constants/constants.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_CONSTANTS_CONSTANTS_ #define SINGULARITY_OPAC_CONSTANTS_CONSTANTS_ diff --git a/singularity-opac/neutrinos/brt_neutrinos.hpp b/singularity-opac/neutrinos/brt_neutrinos.hpp index 6ac9cf1..e55f0da 100644 --- a/singularity-opac/neutrinos/brt_neutrinos.hpp +++ b/singularity-opac/neutrinos/brt_neutrinos.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_BRT_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_BRT_NEUTRINOS_ diff --git a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp index 8dae38d..5f7386d 100644 --- a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp +++ b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_GRAY_OPACITY_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_GRAY_OPACITY_NEUTRINOS_ diff --git a/singularity-opac/neutrinos/gray_s_opacity_neutrinos.hpp b/singularity-opac/neutrinos/gray_s_opacity_neutrinos.hpp index bd5071b..a633592 100644 --- a/singularity-opac/neutrinos/gray_s_opacity_neutrinos.hpp +++ b/singularity-opac/neutrinos/gray_s_opacity_neutrinos.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_GRAY_S_OPACITY_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_GRAY_S_OPACITY_NEUTRINOS_ diff --git a/singularity-opac/neutrinos/mean_neutrino_s_variant.hpp b/singularity-opac/neutrinos/mean_neutrino_s_variant.hpp index baeb39a..cbd0504 100644 --- a/singularity-opac/neutrinos/mean_neutrino_s_variant.hpp +++ b/singularity-opac/neutrinos/mean_neutrino_s_variant.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_MEAN_NEUTRINO_S_VARIANT_ #define SINGULARITY_OPAC_NEUTRINOS_MEAN_NEUTRINO_S_VARIANT_ diff --git a/singularity-opac/neutrinos/mean_neutrino_variant.hpp b/singularity-opac/neutrinos/mean_neutrino_variant.hpp index 4e9f5a0..a2d9a15 100644 --- a/singularity-opac/neutrinos/mean_neutrino_variant.hpp +++ b/singularity-opac/neutrinos/mean_neutrino_variant.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_MEAN_NEUTRINO_VARIANT_ #define SINGULARITY_OPAC_NEUTRINOS_MEAN_NEUTRINO_VARIANT_ diff --git a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp index bfb15b5..4ec3c68 100644 --- a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp +++ b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp @@ -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. @@ -12,9 +12,9 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_MEAN_OPACITY_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_MEAN_OPACITY_NEUTRINOS_ +// This file was made in part with generative AI. #include @@ -174,33 +174,32 @@ class MeanOpacity { lNuMin = toLog_(1.e-3 * PC::kb * fromLog_(lTMin) / PC::h); lNuMax = toLog_(1.e3 * PC::kb * fromLog_(lTMax) / PC::h); } - const Real dlnu = (lNuMax - lNuMin) / (NNu - 1); - // Integrate over frequency + const Real dlnu = (lNuMax - lNuMin) / NNu; + // Integrate over frequency using midpoint rule for (int inu = 0; inu < NNu; ++inu) { - const Real weight = - (inu == 0 || inu == NNu - 1) ? 0.5 : 1.; // Trapezoidal rule - const Real lnu = lNuMin + inu * dlnu; + const Real lnu = lNuMin + (inu + 0.5) * dlnu; const Real nu = fromLog_(lnu); const Real alpha = opac.AbsorptionCoefficient(rho, T, Ye, type, nu, lambda); const Real B = opac.ThermalDistributionOfTNu(T, type, nu); const Real dBdT = opac.DThermalDistributionOfTNuDT(T, type, nu); - kappaPlanckNum += weight * alpha / rho * B * nu * dlnu; - kappaPlanckDenom += weight * B * nu * dlnu; + kappaPlanckNum += alpha / rho * B * nu * dlnu; + kappaPlanckDenom += B * nu * dlnu; // Only contributions to integral from non-zero kappa if (alpha > singularity_opac::robust::SMALL()) { kappaRosselandNum += - weight * singularity_opac::robust::ratio(rho, alpha) * - dBdT * nu * dlnu; - kappaRosselandDenom += weight * dBdT * nu * dlnu; + singularity_opac::robust::ratio(rho, alpha) * dBdT * nu * + dlnu; + kappaRosselandDenom += dBdT * nu * dlnu; } } Real kappaPlanck = singularity_opac::robust::ratio( kappaPlanckNum, kappaPlanckDenom); Real kappaRosseland = - kappaPlanck > singularity_opac::robust::SMALL() + (kappaPlanck > singularity_opac::robust::SMALL() && + kappaRosselandNum > singularity_opac::robust::SMALL()) ? singularity_opac::robust::ratio(kappaRosselandDenom, kappaRosselandNum) : 0.; diff --git a/singularity-opac/neutrinos/mean_s_opacity_neutrinos.hpp b/singularity-opac/neutrinos/mean_s_opacity_neutrinos.hpp index 8665531..976be0f 100644 --- a/singularity-opac/neutrinos/mean_s_opacity_neutrinos.hpp +++ b/singularity-opac/neutrinos/mean_s_opacity_neutrinos.hpp @@ -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. @@ -12,9 +12,9 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_MEAN_S_OPACITY_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_MEAN_S_OPACITY_NEUTRINOS_ +// This file was made in part with generative AI. #include @@ -165,31 +165,30 @@ class MeanSOpacity { lNuMin = toLog_(1.e-3 * pc::kb * fromLog_(lTMin) / pc::h); lNuMax = toLog_(1.e3 * pc::kb * fromLog_(lTMax) / pc::h); } - const Real dlnu = (lNuMax - lNuMin) / (NNu - 1); - // Integrate over frequency + const Real dlnu = (lNuMax - lNuMin) / NNu; + // Integrate over frequency using midpoint rule for (int inu = 0; inu < NNu; ++inu) { - const Real weight = - (inu == 0 || inu == NNu - 1) ? 0.5 : 1.; // Trapezoidal rule - const Real lnu = lNuMin + inu * dlnu; + const Real lnu = lNuMin + (inu + 0.5) * dlnu; const Real nu = fromLog_(lnu); const Real alpha = s_opac.TotalScatteringCoefficient( rho, T, Ye, type, nu, lambda); const Real B = dist.ThermalDistributionOfTNu(T, type, nu); const Real dBdT = dist.DThermalDistributionOfTNuDT(T, type, nu); - kappaPlanckNum += weight * alpha / rho * B * nu * dlnu; - kappaPlanckDenom += weight * B * nu * dlnu; + kappaPlanckNum += alpha / rho * B * nu * dlnu; + kappaPlanckDenom += B * nu * dlnu; if (alpha > singularity_opac::robust::SMALL()) { kappaRosselandNum += - weight * singularity_opac::robust::ratio(rho, alpha) * - dBdT * nu * dlnu; - kappaRosselandDenom += weight * dBdT * nu * dlnu; + singularity_opac::robust::ratio(rho, alpha) * dBdT * nu * + dlnu; + kappaRosselandDenom += dBdT * nu * dlnu; } Real kappaPlanck = singularity_opac::robust::ratio( kappaPlanckNum, kappaPlanckDenom); Real kappaRosseland = - kappaPlanck > singularity_opac::robust::SMALL() + (kappaPlanck > singularity_opac::robust::SMALL() && + kappaRosselandNum > singularity_opac::robust::SMALL()) ? singularity_opac::robust::ratio(kappaRosselandDenom, kappaRosselandNum) : 0.; @@ -221,13 +220,10 @@ class MeanSOpacity { } // namespace impl -using MeanSOpacityScaleFree = - impl::MeanSOpacity, - PhysicalConstantsUnity>; -using MeanSOpacityCGS = +using MeanSOpacityBase = impl::MeanSOpacity, PhysicalConstantsCGS>; -using MeanSOpacity = impl::MeanSVariant>; +using MeanSOpacity = + impl::MeanSVariant>; } // namespace neutrinos } // namespace singularity diff --git a/singularity-opac/neutrinos/neutrino_s_variant.hpp b/singularity-opac/neutrinos/neutrino_s_variant.hpp index b3829de..cf22773 100644 --- a/singularity-opac/neutrinos/neutrino_s_variant.hpp +++ b/singularity-opac/neutrinos/neutrino_s_variant.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_NEUTRINO_S_VARIANT_ #define SINGULARITY_OPAC_NEUTRINOS_NEUTRINO_S_VARIANT_ diff --git a/singularity-opac/neutrinos/neutrino_variant.hpp b/singularity-opac/neutrinos/neutrino_variant.hpp index 2fc255a..6bdd445 100644 --- a/singularity-opac/neutrinos/neutrino_variant.hpp +++ b/singularity-opac/neutrinos/neutrino_variant.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_NEUTRINO_VARIANT_ #define SINGULARITY_OPAC_NEUTRINOS_NEUTRINO_VARIANT_ diff --git a/singularity-opac/neutrinos/non_cgs_neutrinos.hpp b/singularity-opac/neutrinos/non_cgs_neutrinos.hpp index 16a7b7f..7287741 100644 --- a/singularity-opac/neutrinos/non_cgs_neutrinos.hpp +++ b/singularity-opac/neutrinos/non_cgs_neutrinos.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_NON_CGS_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_NON_CGS_NEUTRINOS_ diff --git a/singularity-opac/neutrinos/non_cgs_s_neutrinos.hpp b/singularity-opac/neutrinos/non_cgs_s_neutrinos.hpp index ad29903..cd5084a 100644 --- a/singularity-opac/neutrinos/non_cgs_s_neutrinos.hpp +++ b/singularity-opac/neutrinos/non_cgs_s_neutrinos.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_NON_CGS_S_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_NON_CGS_S_NEUTRINOS_ diff --git a/singularity-opac/neutrinos/opac_neutrinos.hpp b/singularity-opac/neutrinos/opac_neutrinos.hpp index 9967168..f7faa52 100644 --- a/singularity-opac/neutrinos/opac_neutrinos.hpp +++ b/singularity-opac/neutrinos/opac_neutrinos.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_OPAC_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_OPAC_NEUTRINOS_ diff --git a/singularity-opac/neutrinos/s_opac_neutrinos.hpp b/singularity-opac/neutrinos/s_opac_neutrinos.hpp index 7e326fb..2baae71 100644 --- a/singularity-opac/neutrinos/s_opac_neutrinos.hpp +++ b/singularity-opac/neutrinos/s_opac_neutrinos.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_S_OPAC_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_S_OPAC_NEUTRINOS_ diff --git a/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp b/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp index ee8d1a0..05f67bc 100644 --- a/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp +++ b/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_SPINER_OPAC_NEUTRINOS_HPP_ #define SINGULARITY_OPAC_NEUTRINOS_SPINER_OPAC_NEUTRINOS_HPP_ diff --git a/singularity-opac/neutrinos/thermal_distributions_neutrinos.hpp b/singularity-opac/neutrinos/thermal_distributions_neutrinos.hpp index b66bc1e..4903076 100644 --- a/singularity-opac/neutrinos/thermal_distributions_neutrinos.hpp +++ b/singularity-opac/neutrinos/thermal_distributions_neutrinos.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_THERMAL_DISTRIBUTIONS_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_THERMAL_DISTRIBUTIONS_NEUTRINOS_ diff --git a/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp b/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp index 25df1e9..716da94 100644 --- a/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp +++ b/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_NEUTRINOS_TOPHAT_EMISSIVITY_NEUTRINOS_ #define SINGULARITY_OPAC_NEUTRINOS_TOPHAT_EMISSIVITY_NEUTRINOS_ @@ -56,7 +55,7 @@ class TophatEmissivity { const RadiationType type, const Real nu, Real *lambda = nullptr) const { return dist_.AbsorptionCoefficientFromKirchhoff(*this, rho, temp, Ye, type, - nu, lambda) / + nu, lambda) / (4. * M_PI); } diff --git a/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp b/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp index 04961cc..60dc84d 100644 --- a/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp +++ b/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_EPBREMSSTRAHLUNG_OPACITY_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_EPBREMSSTRAHLUNG_OPACITY_PHOTONS_ @@ -54,7 +53,7 @@ class EPBremsstrahlungOpacity { Real AbsorptionCoefficient(const Real rho, const Real temp, const Real nu, Real *lambda = nullptr) const { return dist_.AbsorptionCoefficientFromKirchhoff(*this, rho, temp, nu, - lambda); + lambda); } template diff --git a/singularity-opac/photons/example_ascii/kap_plaw.txt b/singularity-opac/photons/example_ascii/kap_plaw.txt deleted file mode 100644 index 59568c1..0000000 --- a/singularity-opac/photons/example_ascii/kap_plaw.txt +++ /dev/null @@ -1,19 +0,0 @@ - NRho 4 NT 4 - rho_min 1.00000000e-14 rho_max 7.94328235e-01 - T_min 1.00000000e+00 T_max 7.94328235e+06 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 -1.00000000e-03 1.00000000e-01 diff --git a/singularity-opac/photons/example_ascii/preproc_ascii_opac.py b/singularity-opac/photons/example_ascii/preproc_ascii_opac.py deleted file mode 100644 index 295c0ad..0000000 --- a/singularity-opac/photons/example_ascii/preproc_ascii_opac.py +++ /dev/null @@ -1,103 +0,0 @@ -# ====================================================================== -# © 2024. 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. -# ====================================================================== - -#--------------------------------------------------------------------------------------------------# -# Regrid (interpolate) Rosseland and Planck opacity data to a log-log rho-T grid. -#--------------------------------------------------------------------------------------------------# -import numpy as np -from scipy.interpolate import interp2d -import argparse - -#-- parse command line arguments -parser = argparse.ArgumentParser(description="Re-grid opacity data to be log-log in density-temperature.") -parser.add_argument("fname", type=str, help="Opacity file to re-grid.") -parser.add_argument("--is_loglog", action="store_true", help="Avoid regridding if already log-log.") -parser.add_argument("--fname_new", type=str, default="kappa_new.txt", help="File name for new file created.") -parser.add_argument("--plot", action="store_true", help="If not log-log, plot colormap of original and regridded data.") -args = parser.parse_args() - -NRho = -1 -NT = -1 -with open(args.fname, "r") as ff: - svec = ff.readline().split(" ") - NRho = int(svec[2]) - NT = int(svec[4]) - -print("NRho = ", NRho) -print("NT = ", NT) - -opac = np.loadtxt(args.fname, skiprows=1) - -print("opac.shape = ", opac.shape) -assert np.size(opac, 0) == NT * NRho, "np.size(opac, 0) != NT * NRho" - -#-- density, temperature grids -Rho = np.unique(opac[:,0]) -T = np.unique(opac[:,1]) - -assert np.size(Rho) == NRho, "np.size(Rho) != Rho" -assert np.size(T) == NT, "np.size(T) != NT" - -r = np.logspace(np.log10(Rho[0]), np.log10(Rho[NRho - 1]), NRho) -tt = np.logspace(np.log10(T[0]), np.log10(T[NT - 1]), NT) - -if (args.is_loglog): - assert np.max(abs(r - Rho) / Rho) < 1e-6, "np.max(abs(r - Rho) / Rho) >= 1e-6" - assert np.max(abs(tt - T) / T) < 1e-6, "np.max(abs(tt - T) / T) >= 1e-6" -else: - print() - print("Interpolating data to log-log rho-T grid...") - #-- reshape to rho-T grid - ross_old_opac = np.reshape(opac[:,2], (NRho, NT)) - plnk_old_opac = np.reshape(opac[:,3], (NRho, NT)) - #-- linearly interpolate in log-log rho-T - ross_f2d = interp2d(np.log10(T), np.log10(Rho), ross_old_opac, kind='linear') - plnk_f2d = interp2d(np.log10(T), np.log10(Rho), plnk_old_opac, kind='linear') - ross_new_opac = ross_f2d(np.log10(tt), np.log10(r)) - plnk_new_opac = plnk_f2d(np.log10(tt), np.log10(r)) - #-- plot 4-panel of data - if (args.plot): - import matplotlib.pyplot as plt - XOLD, YOLD = np.meshgrid(np.log10(T), np.log10(Rho)) - XNEW, YNEW = np.meshgrid(np.log10(tt), np.log10(r)) - fig, axes = plt.subplots(2, 2) - im1 = axes[0, 0].pcolormesh(XOLD, YOLD, ross_old_opac) - axes[0, 0].set_title('Old Rosseland') - fig.colorbar(im1, ax=axes[0, 0]) - im2 = axes[0, 1].pcolormesh(XOLD, YOLD, plnk_old_opac) - axes[0, 1].set_title('Old Planck') - fig.colorbar(im2, ax=axes[0, 1]) - im3 = axes[1, 0].pcolormesh(XNEW, YNEW, ross_new_opac) - axes[1, 0].set_title('New Rosseland') - fig.colorbar(im3, ax=axes[1, 0]) - im4 = axes[1, 1].pcolormesh(XNEW, YNEW, plnk_new_opac) - axes[1, 1].set_title('New Planck') - fig.colorbar(im4, ax=axes[1, 1]) - plt.tight_layout() - plt.show() - #-- reset opacity array - for i in range(NRho): - for j in range(NT): - k = j + NT * i - opac[k,0] = r[i] - opac[k,1] = tt[j] - opac[k,2] = ross_new_opac[i,j] - opac[k,3] = plnk_new_opac[i,j] - -#-- save with new header containing min/max rho, T -hdr = "NRho "+str(NRho)+" NT "+str(NT) + "\n" \ - "rho_min {:.8e}".format(r[0])+" rho_max {:.8e}".format(r[NRho-1]) + "\n" \ - "T_min {:.8e}".format(tt[0]) + " T_max {:.8e}".format(tt[NT-1]) -np.savetxt(args.fname_new, opac[:,2:], delimiter=" ", header=hdr, comments=" ", fmt="%.8e") diff --git a/singularity-opac/photons/gray_opacity_photons.hpp b/singularity-opac/photons/gray_opacity_photons.hpp index d0fa891..007eab4 100644 --- a/singularity-opac/photons/gray_opacity_photons.hpp +++ b/singularity-opac/photons/gray_opacity_photons.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_GRAY_OPACITY_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_GRAY_OPACITY_PHOTONS_ diff --git a/singularity-opac/photons/gray_s_opacity_photons.hpp b/singularity-opac/photons/gray_s_opacity_photons.hpp index d2f1a88..a4610ac 100644 --- a/singularity-opac/photons/gray_s_opacity_photons.hpp +++ b/singularity-opac/photons/gray_s_opacity_photons.hpp @@ -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. @@ -12,9 +12,9 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_GRAY_S_OPACITY_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_GRAY_S_OPACITY_PHOTONS_ +// This file was made in part with generative AI. #include #include @@ -22,6 +22,7 @@ #include #include +#include namespace singularity { namespace photons { diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp index 9e97557..4e9bad0 100644 --- a/singularity-opac/photons/mean_opacity_photons.hpp +++ b/singularity-opac/photons/mean_opacity_photons.hpp @@ -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. @@ -12,31 +12,33 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_MEAN_OPACITY_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_MEAN_OPACITY_PHOTONS_ +// This file was made in part with generative AI. +#include #include -#include -#include +#include +#include +#include #include -#include +#include +#include #include #include #include -#include - +#include +#include #include #include #include +#include namespace singularity { namespace photons { namespace impl { -#define EPS (10.0 * std::numeric_limits::min()) - // TODO(BRR) Note: It is assumed that lambda is constant for all densities and // temperatures @@ -44,84 +46,100 @@ template class MeanOpacity { public: using PC = pc; + using DataBox = Spiner::DataBox; MeanOpacity() = default; - template - MeanOpacity(const Opacity &opac, const Real lRhoMin, const Real lRhoMax, - const int NRho, const Real lTMin, const Real lTMax, const int NT, - Real *lambda = nullptr) { - MeanOpacityImpl_(opac, lRhoMin, lRhoMax, NRho, lTMin, lTMax, - NT, -1., -1., 100, lambda); - } - template + template MeanOpacity(const Opacity &opac, const Real lRhoMin, const Real lRhoMax, const int NRho, const Real lTMin, const Real lTMax, const int NT, - Real lNuMin, Real lNuMax, const int NNu, Real *lambda = nullptr) { - MeanOpacityImpl_(opac, lRhoMin, lRhoMax, NRho, lTMin, lTMax, - NT, lNuMin, lNuMax, NNu, lambda); + const GroupBoundsIndexer &group_bounds, const int ngroups, + const int NNuPerGroup = 64, Real *lambda = nullptr) { + MeanOpacityImpl_(opac, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, + group_bounds, ngroups, NNuPerGroup, lambda); } - // construct Planck/Rosseland DataBox from ascii file - MeanOpacity(const std::string &filename) : filename_(filename.c_str()) { - - // get number of density and temperature points - std::ifstream ff(filename.c_str()); - const bool fexists = ff.good(); - - if (fexists) { - - std::filesystem::path filePath(filename); - std::string extension = filePath.extension().string(); + template + MeanOpacity(const DataBox &kappaPlanck, const DataBox &kappaRosseland, + const GroupBoundsIndexer &group_bounds) { + // Table-backed multigroup opacities always carry explicit group bounds. + // To represent [nu_max, infinity), the final bound must literally be + // IEEE +infinity, not a large finite proxy value. + LoadOpacityTables_(kappaPlanck, kappaRosseland, group_bounds); + } - if (extension == ".txt") { - loadASCII(ff); #ifdef SPINER_USE_HDF - } else if (extension == ".hdf5" || extension == ".h5" || extension == ".sp5") { - herr_t status = H5_SUCCESS; - hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); - status += lkappa_.loadHDF(file, "Rosseland & Planck Mean Opacities"); - status += H5Fclose(file); - - if (status != H5_SUCCESS) { - OPAC_ERROR("photons::MeanOpacity: HDF5 error\n"); - } -#endif - } else { - OPAC_ERROR("photons::MeanOpacity: unrecognized file extension"); - } + MeanOpacity(const std::string &filename) { + // HDF-backed multigroup tables are expected to provide an ngroups + 1 + // "group bounds" dataset. If the last group is [nu_max, infinity), then + // the final stored bound must be IEEE +infinity. + DataBox kappaPlanck; + DataBox kappaRosseland; + DataBox groupBounds; + herr_t status = H5_SUCCESS; + hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + status += + kappaPlanck.loadHDF(file, SP5::MultigroupOpac::PlanckGroupOpacity); + status += kappaRosseland.loadHDF( + file, SP5::MultigroupOpac::RosselandGroupOpacity); + status += groupBounds.loadHDF(file, SP5::MultigroupOpac::GroupBounds); + status += H5Fclose(file); - } else { - OPAC_ERROR("photons::MeanOpacity: file does not exist"); + if (status != H5_SUCCESS) { + OPAC_ERROR("photons::MeanOpacity: HDF5 error\n"); } + + LoadOpacityTables_(kappaPlanck, kappaRosseland, groupBounds); + groupBounds.finalize(); + kappaPlanck.finalize(); + kappaRosseland.finalize(); } -#ifdef SPINER_USE_HDF void Save(const std::string &filename) const { + DataBox kappaPlanck; + DataBox kappaRosseland; + DataBox groupBounds; + ExportOpacityTables_(kappaPlanck, kappaRosseland); + ExportGroupBounds(groupBounds, groupBounds_, ngroups_); + herr_t status = H5_SUCCESS; hid_t file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - status += lkappa_.saveHDF(file, "Rosseland & Planck Mean Opacities"); + status += + kappaPlanck.saveHDF(file, SP5::MultigroupOpac::PlanckGroupOpacity); + status += kappaRosseland.saveHDF( + file, SP5::MultigroupOpac::RosselandGroupOpacity); + status += groupBounds.saveHDF(file, SP5::MultigroupOpac::GroupBounds); status += H5Fclose(file); + kappaPlanck.finalize(); + kappaRosseland.finalize(); + groupBounds.finalize(); + if (status != H5_SUCCESS) { OPAC_ERROR("photons::MeanOpacity: HDF5 error\n"); } } #endif - PORTABLE_INLINE_FUNCTION void PrintParams() const { - printf("Mean opacity\n"); + PORTABLE_INLINE_FUNCTION + void PrintParams() const { + printf("Photon multigroup opacity. ngroups = %d\n", ngroups_); } MeanOpacity GetOnDevice() { MeanOpacity other; - other.lkappa_ = Spiner::getOnDeviceDataBox(lkappa_); + other.lkappaPlanck_ = Spiner::getOnDeviceDataBox(lkappaPlanck_); + other.lkappaRosseland_ = Spiner::getOnDeviceDataBox(lkappaRosseland_); + other.groupBounds_ = Spiner::getOnDeviceDataBox(groupBounds_); + other.ngroups_ = ngroups_; return other; } void Finalize() { - lkappa_.finalize(); + lkappaPlanck_.finalize(); + lkappaRosseland_.finalize(); + groupBounds_.finalize(); } PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants @@ -129,205 +147,266 @@ class MeanOpacity { return RuntimePhysicalConstants(PC()); } + PORTABLE_INLINE_FUNCTION + int ngroups() const noexcept { return ngroups_; } + + PORTABLE_INLINE_FUNCTION + bool HasGroupBounds() const noexcept { return true; } + + std::vector GetGroupBounds() const { + std::vector bounds(ngroups_ + 1); + for (int group = 0; group <= ngroups_; ++group) { + bounds[group] = groupBounds_(group); + } + return bounds; + } + + // The group-index-less "mean" accessors only make sense when there is a + // single group. In that case the lone group spans the entire spectrum, so + // its group-integrated coefficient (stored at group 0) IS the traditional + // gray mean. We therefore require ngroups==1 and forward to group 0. This is + // not a distinguished "mean slot": for ngroups>1, group 0 is simply the + // lowest-frequency group and callers must use the group-index API. PORTABLE_INLINE_FUNCTION Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const { - Real lRho = toLog_(rho); - Real lT = toLog_(temp); - return rho * fromLog_(lkappa_.interpToReal(lRho, lT, Planck)); + PORTABLE_REQUIRE( + ngroups_ == 1, + "PlanckMeanAbsorptionCoefficient only valid for ngroups==1. " + "Use PlanckGroupAbsorptionCoefficient(rho, temp, group) for " + "multigroup."); + return PlanckGroupAbsorptionCoefficient(rho, temp, 0); } + // See PlanckMeanAbsorptionCoefficient: the gray mean is the single group 0. PORTABLE_INLINE_FUNCTION Real RosselandMeanAbsorptionCoefficient(const Real rho, const Real temp) const { - Real lRho = toLog_(rho); - Real lT = toLog_(temp); - return rho * fromLog_(lkappa_.interpToReal(lRho, lT, Rosseland)); + PORTABLE_REQUIRE( + ngroups_ == 1, + "RosselandMeanAbsorptionCoefficient only valid for ngroups==1. " + "Use RosselandGroupAbsorptionCoefficient(rho, temp, group) for " + "multigroup."); + return RosselandGroupAbsorptionCoefficient(rho, temp, 0); + } + + PORTABLE_INLINE_FUNCTION + Real PlanckGroupAbsorptionCoefficient(const Real rho, const Real temp, + const int group) const { + return GroupAbsorptionCoefficient_(lkappaPlanck_, rho, temp, group); + } + + PORTABLE_INLINE_FUNCTION + Real RosselandGroupAbsorptionCoefficient(const Real rho, const Real temp, + const int group) const { + return GroupAbsorptionCoefficient_(lkappaRosseland_, rho, temp, group); } - // 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)); + Real AbsorptionCoefficient(const Real rho, const Real temp, const int group, + const int gmode = Rosseland) const { + return (gmode == Planck) + ? PlanckGroupAbsorptionCoefficient(rho, temp, group) + : RosselandGroupAbsorptionCoefficient(rho, temp, group); } + // Like the mean accessors above, Emissivity has no group-index argument and + // so is only defined for ngroups==1, where group 0 is the whole-spectrum + // (gray) group. PORTABLE_INLINE_FUNCTION - Real Emissivity(const Real rho, const Real temp, - const int gmode = Rosseland, + 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; + if (ngroups_ != 1) { + OPAC_ERROR("photons::MeanOpacity: Emissivity only valid for ngroups==1"); + } + PlanckDistribution dist; + Real B = dist.ThermalDistributionOfT(temp, lambda); + return AbsorptionCoefficient(rho, temp, 0, gmode) * B; + } + + PORTABLE_INLINE_FUNCTION + int GroupOfNu(const Real nu) const { + if (!(nu >= GroupBoundAt(groupBounds_, 0) && + nu <= GroupBoundAt(groupBounds_, ngroups_))) { + OPAC_ERROR("photons::MeanOpacity: frequency is outside group bounds"); + } + return GroupOfNuImpl(groupBounds_, ngroups_, nu); + } + + PORTABLE_INLINE_FUNCTION + Real PlanckGroupAbsorptionCoefficientFromNu(const Real rho, const Real temp, + const Real nu) const { + return AbsorptionCoefficientFromNu(rho, temp, nu, Planck); + } + + PORTABLE_INLINE_FUNCTION + Real RosselandGroupAbsorptionCoefficientFromNu(const Real rho, + const Real temp, + const Real nu) const { + return AbsorptionCoefficientFromNu(rho, temp, nu, Rosseland); + } + + PORTABLE_INLINE_FUNCTION + Real AbsorptionCoefficientFromNu(const Real rho, const Real temp, + const Real nu, + const int gmode = Rosseland) const { + return AbsorptionCoefficient(rho, temp, GroupOfNu(nu), gmode); } private: - template + PORTABLE_INLINE_FUNCTION + Real GroupAbsorptionCoefficient_(const DataBox &lkappa, const Real rho, + const Real temp, const int group) const { + const Real lRho = ToLog(rho); + const Real lT = ToLog(temp); + return rho * FromLog(lkappa.interpToReal(lRho, lT, group)); + } + + void ValidateOpacityTables_(const DataBox &kappaPlanck, + const DataBox &kappaRosseland) const { + if (kappaPlanck.rank() != 3 || kappaRosseland.rank() != 3) { + OPAC_ERROR("photons::MeanOpacity: opacity tables must be rank 3"); + } + for (int dim = 1; dim <= 3; ++dim) { + if (kappaPlanck.dim(dim) != kappaRosseland.dim(dim)) { + OPAC_ERROR("photons::MeanOpacity: table dimensions do not match"); + } + } + if (kappaPlanck.dim(1) <= 0) { + OPAC_ERROR("photons::MeanOpacity: ngroups must be positive"); + } + if (kappaPlanck.range(1) != kappaRosseland.range(1) || + kappaPlanck.range(2) != kappaRosseland.range(2)) { + OPAC_ERROR("photons::MeanOpacity: table ranges do not match"); + } + } + + template + void LoadOpacityTables_(const DataBox &kappaPlanck, + const DataBox &kappaRosseland, + const GroupBoundsIndexer &group_bounds) { + ValidateOpacityTables_(kappaPlanck, kappaRosseland); + ngroups_ = kappaPlanck.dim(1); + ValidateGroupBounds(group_bounds, ngroups_); + SetGroupBounds(groupBounds_, group_bounds, ngroups_); + lkappaPlanck_.copyMetadata(kappaPlanck); + lkappaRosseland_.copyMetadata(kappaRosseland); + for (int i = 0; i < kappaPlanck.size(); ++i) { + lkappaPlanck_(i) = ToLog(kappaPlanck(i)); + lkappaRosseland_(i) = ToLog(kappaRosseland(i)); + } + } + + void ExportOpacityTables_(DataBox &kappaPlanck, + DataBox &kappaRosseland) const { + kappaPlanck.copyMetadata(lkappaPlanck_); + kappaRosseland.copyMetadata(lkappaRosseland_); + for (int i = 0; i < lkappaPlanck_.size(); ++i) { + kappaPlanck(i) = FromLog(lkappaPlanck_(i)); + kappaRosseland(i) = FromLog(lkappaRosseland_(i)); + } + } + + template void MeanOpacityImpl_(const Opacity &opac, const Real lRhoMin, const Real lRhoMax, const int NRho, const Real lTMin, - const Real lTMax, const int NT, Real lNuMin, - Real lNuMax, const int NNu, Real *lambda = nullptr) { + const Real lTMax, const int NT, + const GroupBoundsIndexer &group_bounds, + const int ngroups, const int NNuPerGroup, + Real *lambda = nullptr) { #ifndef NDEBUG auto RPC = RuntimePhysicalConstants(PC()); auto opc = opac.GetRuntimePhysicalConstants(); assert(RPC == opc && "Physical constants are the same"); #endif - lkappa_.resize(NRho, NT, 2); - lkappa_.setRange(1, lTMin, lTMax, NT); - lkappa_.setRange(2, lRhoMin, lRhoMax, NRho); - - // Fill tables - for (int iRho = 0; iRho < NRho; ++iRho) { - Real lRho = lkappa_.range(2).x(iRho); - Real rho = fromLog_(lRho); - for (int iT = 0; iT < NT; ++iT) { - Real lT = lkappa_.range(1).x(iT); - Real T = fromLog_(lT); - Real kappaPlanckNum = 0.; - Real kappaPlanckDenom = 0.; - Real kappaRosselandNum = 0.; - Real kappaRosselandDenom = 0.; - if (AUTOFREQ) { - lNuMin = toLog_(1.e-3 * PC::kb * fromLog_(lTMin) / PC::h); - lNuMax = toLog_(1.e3 * PC::kb * fromLog_(lTMax) / PC::h); - } - const Real dlnu = (lNuMax - lNuMin) / (NNu - 1); - // Integrate over frequency - for (int inu = 0; inu < NNu; ++inu) { - const Real weight = - (inu == 0 || inu == NNu - 1) ? 0.5 : 1.; // Trapezoidal rule - const Real lnu = lNuMin + inu * dlnu; - const Real nu = fromLog_(lnu); - const Real alpha = opac.AbsorptionCoefficient(rho, T, nu, lambda); - const Real B = opac.ThermalDistributionOfTNu(T, nu); - const Real dBdT = opac.DThermalDistributionOfTNuDT(T, nu); - kappaPlanckNum += weight * alpha / rho * B * nu * dlnu; - kappaPlanckDenom += weight * B * nu * dlnu; - - if (alpha > singularity_opac::robust::SMALL()) { - kappaRosselandNum += weight * - singularity_opac::robust::ratio(rho, alpha) * - dBdT * nu * dlnu; - kappaRosselandDenom += weight * dBdT * nu * dlnu; - } - } - - Real kappaPlanck = - singularity_opac::robust::ratio(kappaPlanckNum, kappaPlanckDenom); - Real kappaRosseland = kappaPlanck > singularity_opac::robust::SMALL() - ? singularity_opac::robust::ratio( - kappaRosselandDenom, kappaRosselandNum) - : 0.; - - lkappa_(iRho, iT, Rosseland) = toLog_(kappaRosseland); - lkappa_(iRho, iT, Planck) = toLog_(kappaPlanck); - if (std::isnan(lkappa_(iRho, iT, Rosseland)) || - std::isnan(lkappa_(iRho, iT, Planck))) { - OPAC_ERROR("photons::MeanOpacity: NAN in opacity evaluations"); - } - } + if (NNuPerGroup < 2) { + OPAC_ERROR("photons::MeanOpacity: NNuPerGroup must be at least 2"); } - } + ValidateGroupBounds(group_bounds, ngroups); + + ngroups_ = ngroups; + SetGroupBounds(groupBounds_, group_bounds, ngroups_); + lkappaPlanck_.resize(NRho, NT, ngroups_); + lkappaPlanck_.setRange(1, lTMin, lTMax, NT); + lkappaPlanck_.setRange(2, lRhoMin, lRhoMax, NRho); + lkappaRosseland_.copyMetadata(lkappaPlanck_); + + PlanckDistribution dist; + std::vector planckDenom(ngroups_, 0.); + std::vector rosselandDenom(ngroups_, 0.); + + for (int iT = 0; iT < NT; ++iT) { + const Real lT = lkappaPlanck_.range(1).x(iT); + const Real T = FromLog(lT); + + for (int group = 0; group < ngroups_; ++group) { + Real Baccum = 0.; + Real dBdTaccum = 0.; + const Real nuMin = GroupBoundAt(group_bounds, group); + const Real nuMax = GroupBoundAt(group_bounds, group + 1); + ForEachGroupFrequencySample( + T, nuMin, nuMax, NNuPerGroup, [&](const Real nu, const Real dnu) { + Real B = 0.; + Real dBdT = 0.; + ThermalWeightsAtNu(dist, T, nu, B, dBdT); + Baccum += B * dnu; + dBdTaccum += dBdT * dnu; + }); + + planckDenom[group] = Baccum; + rosselandDenom[group] = dBdTaccum; + } - // ASCII opacity file reader - void loadASCII(std::ifstream &ff) { - - int NRho = -1; - int NT = -1; - - // line read from file - std::string fline; - - // read 1st line of header to get sizes - std::getline(ff, fline); - - // tokenize fline - char* cfline = const_cast(fline.c_str()); - char* fl_tok = std::strtok(cfline, " "); - - // move to next token to get number of density points - fl_tok = std::strtok(nullptr, " "); - NRho = std::stoi(fl_tok); - - // move to next token to get number of temperature points - fl_tok = std::strtok(nullptr, " "); - fl_tok = std::strtok(nullptr, " "); - NT = std::stoi(fl_tok); - - // read 2nd line of header to get min/max density - std::getline(ff, fline); - // tokenize fline - cfline = const_cast(fline.c_str()); - fl_tok = std::strtok(cfline, " "); - fl_tok = std::strtok(nullptr, " "); - const Real RhoMin = std::stod(fl_tok); - fl_tok = std::strtok(nullptr, " "); - fl_tok = std::strtok(nullptr, " "); - const Real RhoMax = std::stod(fl_tok); - - // read 3nd line of header to get min/max temperature - std::getline(ff, fline); - // tokenize fline - cfline = const_cast(fline.c_str()); - fl_tok = std::strtok(cfline, " "); - fl_tok = std::strtok(nullptr, " "); - const Real TMin = std::stod(fl_tok); - fl_tok = std::strtok(nullptr, " "); - fl_tok = std::strtok(nullptr, " "); - const Real TMax = std::stod(fl_tok); - - // reseize the Planck and Rosseland databox - lkappa_.resize(NRho, NT, 2); - - // set rho-T ranges - const Real lTMin = toLog_(TMin); - const Real lTMax = toLog_(TMax); - const Real lRhoMin = toLog_(RhoMin); - const Real lRhoMax = toLog_(RhoMax); - lkappa_.setRange(1, lTMin, lTMax, NT); - lkappa_.setRange(2, lRhoMin, lRhoMax, NRho); - - // fill tables - for (int iRho = 0; iRho < NRho; ++iRho) { - const Real lRho_i = lkappa_.range(2).x(iRho); - for (int iT = 0; iT < NT; ++iT) { - - // get new file-line - std::getline(ff, fline); - cfline = const_cast(fline.c_str()); - fl_tok = std::strtok(cfline, " "); - - // populate Rosseland opacity [cm^2/g] - lkappa_(iRho, iT, Rosseland) = toLog_(std::stod(fl_tok)); - - // populate Planck opacity [cm^2/g] - fl_tok = std::strtok(nullptr, " "); - lkappa_(iRho, iT, Planck) = toLog_(std::stod(fl_tok)); - - if (std::isnan(lkappa_(iRho, iT, Rosseland)) || - std::isnan(lkappa_(iRho, iT, Planck))) { - OPAC_ERROR("photons::MeanOpacity: NAN in parsed ASCII opacity"); + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real lRho = lkappaPlanck_.range(2).x(iRho); + const Real rho = FromLog(lRho); + + for (int group = 0; group < ngroups_; ++group) { + Real kappaPlanckNum = 0.; + Real kappaRosselandNum = 0.; + const Real nuMin = GroupBoundAt(group_bounds, group); + const Real nuMax = GroupBoundAt(group_bounds, group + 1); + ForEachGroupFrequencySample( + T, nuMin, nuMax, NNuPerGroup, [&](const Real nu, const Real dnu) { + const Real alpha = + opac.AbsorptionCoefficient(rho, T, nu, lambda); + Real B = 0.; + Real dBdT = 0.; + ThermalWeightsAtNu(dist, T, nu, B, dBdT); + kappaPlanckNum += alpha / rho * B * dnu; + + if (alpha > singularity_opac::robust::SMALL()) { + kappaRosselandNum += + singularity_opac::robust::ratio(rho, alpha) * dBdT * dnu; + } + }); + + const Real kappaPlanck = singularity_opac::robust::ratio( + kappaPlanckNum, planckDenom[group]); + const Real kappaRosseland = + (rosselandDenom[group] > singularity_opac::robust::SMALL() && + kappaRosselandNum > singularity_opac::robust::SMALL()) + ? singularity_opac::robust::ratio(rosselandDenom[group], + kappaRosselandNum) + : 0.; + + lkappaPlanck_(iRho, iT, group) = ToLog(kappaPlanck); + lkappaRosseland_(iRho, iT, group) = ToLog(kappaRosseland); + if (std::isnan(lkappaPlanck_(iRho, iT, group)) || + std::isnan(lkappaRosseland_(iRho, iT, group))) { + OPAC_ERROR("photons::MeanOpacity: NAN in opacity evaluations"); + } } } } } - PORTABLE_INLINE_FUNCTION Real toLog_(const Real x) const { - return std::log10(std::abs(x) + EPS); - } - PORTABLE_INLINE_FUNCTION Real fromLog_(const Real lx) const { - return std::pow(10., lx); - } - Spiner::DataBox lkappa_; - const char *filename_; - PlanckDistribution dist_; + DataBox lkappaPlanck_; + DataBox lkappaRosseland_; + DataBox groupBounds_; + int ngroups_ = 0; }; -#undef EPS - } // namespace impl using MeanOpacityBase = impl::MeanOpacity; @@ -337,4 +416,4 @@ using MeanOpacity = } // namespace photons } // namespace singularity -#endif // SINGULARITY_OPAC_PHOTONS_MEAN_OPACITY_PHOTONS__ +#endif // SINGULARITY_OPAC_PHOTONS_MEAN_OPACITY_PHOTONS_ diff --git a/singularity-opac/photons/mean_photon_s_variant.hpp b/singularity-opac/photons/mean_photon_s_variant.hpp index 699e5e1..fafac87 100644 --- a/singularity-opac/photons/mean_photon_s_variant.hpp +++ b/singularity-opac/photons/mean_photon_s_variant.hpp @@ -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. @@ -12,9 +12,9 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_MEAN_PHOTON_S_VARIANT_ #define SINGULARITY_OPAC_PHOTONS_MEAN_PHOTON_S_VARIANT_ +// This file was made in part with generative AI. #include @@ -88,6 +88,70 @@ class MeanSVariant { s_opac_); } + PORTABLE_INLINE_FUNCTION + int ngroups() const noexcept { + return PortsOfCall::visit( + [](const auto &s_opac) { return s_opac.ngroups(); }, s_opac_); + } + + PORTABLE_INLINE_FUNCTION + bool HasGroupBounds() const noexcept { + return PortsOfCall::visit( + [](const auto &s_opac) { return s_opac.HasGroupBounds(); }, s_opac_); + } + + PORTABLE_INLINE_FUNCTION + Real PlanckGroupScatteringCoefficient(const Real rho, const Real temp, + const int group) const { + return ScatteringCoefficient(rho, temp, group, Planck); + } + + PORTABLE_INLINE_FUNCTION + Real RosselandGroupScatteringCoefficient(const Real rho, const Real temp, + const int group) const { + return ScatteringCoefficient(rho, temp, group, Rosseland); + } + + PORTABLE_INLINE_FUNCTION + Real ScatteringCoefficient(const Real rho, const Real temp, const int group, + const int gmode = Rosseland) const { + return PortsOfCall::visit( + [=](const auto &s_opac) { + return s_opac.ScatteringCoefficient(rho, temp, group, gmode); + }, + s_opac_); + } + + PORTABLE_INLINE_FUNCTION + int GroupOfNu(const Real nu) const { + return PortsOfCall::visit( + [=](const auto &s_opac) { return s_opac.GroupOfNu(nu); }, s_opac_); + } + + PORTABLE_INLINE_FUNCTION + Real PlanckGroupScatteringCoefficientFromNu(const Real rho, const Real temp, + const Real nu) const { + return ScatteringCoefficientFromNu(rho, temp, nu, Planck); + } + + PORTABLE_INLINE_FUNCTION + Real RosselandGroupScatteringCoefficientFromNu(const Real rho, + const Real temp, + const Real nu) const { + return ScatteringCoefficientFromNu(rho, temp, nu, Rosseland); + } + + PORTABLE_INLINE_FUNCTION + Real ScatteringCoefficientFromNu(const Real rho, const Real temp, + const Real nu, + const int gmode = Rosseland) const { + return PortsOfCall::visit( + [=](const auto &s_opac) { + return s_opac.ScatteringCoefficientFromNu(rho, temp, nu, gmode); + }, + s_opac_); + } + inline void Finalize() noexcept { return PortsOfCall::visit([](auto &s_opac) { return s_opac.Finalize(); }, s_opac_); diff --git a/singularity-opac/photons/mean_photon_types.hpp b/singularity-opac/photons/mean_photon_types.hpp index 763a3af..0606d41 100644 --- a/singularity-opac/photons/mean_photon_types.hpp +++ b/singularity-opac/photons/mean_photon_types.hpp @@ -1,5 +1,5 @@ // ====================================================================== -// © 2025. Triad National Security, LLC. All rights reserved. This +// © 2025-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. @@ -12,18 +12,24 @@ // 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_ +// This file was made in part with generative AI. + +#include namespace singularity { namespace photons { // mean-opacity mode -enum OpacityAveraging { - Rosseland = 0, - Planck = 1 -}; +enum OpacityAveraging { Rosseland = 0, Planck = 1 }; + +// Crossover into the Wien tail of the Planck function, in units of the +// dimensionless x = h * nu / (kb * temp). For x above this value the exact +// Planck form 1 / expm1(x) and the Wien approximation exp(-x) agree to well +// below machine precision (exp(-80) ~ 2e-35), so we switch to the cheaper, +// monotone Wien closed form. This is a regime switch, not an overflow guard. +constexpr Real wien_tail_x = 80.; } // namespace photons } // namespace singularity diff --git a/singularity-opac/photons/mean_photon_utils.hpp b/singularity-opac/photons/mean_photon_utils.hpp new file mode 100644 index 0000000..936cab6 --- /dev/null +++ b/singularity-opac/photons/mean_photon_utils.hpp @@ -0,0 +1,208 @@ +// ====================================================================== +// © 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 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_UTILS_ +#define SINGULARITY_OPAC_PHOTONS_MEAN_PHOTON_UTILS_ +// This file was made in part with generative AI. + +#include +#include + +#include +#include +#include +#include +#include + +namespace singularity { +namespace photons { +namespace impl { + +using MeanUtilsDataBox = Spiner::DataBox; + +// Log/anti-log transforms used to store and interpolate opacities. A small +// floor keeps toLog well-defined at zero. +PORTABLE_INLINE_FUNCTION +Real ToLog(const Real x) { + constexpr Real eps = 10.0 * std::numeric_limits::min(); + return std::log10(std::abs(x) + eps); +} + +PORTABLE_INLINE_FUNCTION +Real FromLog(const Real lx) { return std::pow(10., lx); } + +// Group-bound access, supporting both a generic indexer (used at construction +// time) and a stored DataBox (used after the bounds are cached). +template +PORTABLE_INLINE_FUNCTION Real +GroupBoundAt(const GroupBoundsIndexer &group_bounds, const int group) { + return group_bounds[group]; +} + +PORTABLE_INLINE_FUNCTION +Real GroupBoundAt(const MeanUtilsDataBox &group_bounds, const int group) { + return group_bounds(group); +} + +// Locate the group index containing nu via binary search over the cached +// bounds. Boundary cases are shortcut: nu at or below the first bound maps to +// group 0, and nu at or above the final bound maps to the last group. Callers +// are responsible for range-checking nu before invoking this helper. +PORTABLE_INLINE_FUNCTION +int GroupOfNuImpl(const MeanUtilsDataBox &groupBounds, const int ngroups, + const Real nu) { + // Shortcuts for boundary cases + if (nu <= GroupBoundAt(groupBounds, 0)) { + return 0; + } + if (nu >= GroupBoundAt(groupBounds, ngroups)) { + return ngroups - 1; + } + // Binary search to find group index containing nu + int lower = 0; + int upper = ngroups; + while (upper - lower > 1) { + const int middle = (lower + upper) / 2; + if (nu < GroupBoundAt(groupBounds, middle)) { + upper = middle; + } else { + lower = middle; + } + } + return lower; +} + +// Copy user-provided group bounds into the class-owned DataBox. +template +void SetGroupBounds(MeanUtilsDataBox &groupBounds, + const GroupBoundsIndexer &group_bounds, const int ngroups) { + groupBounds.resize(ngroups + 1); + for (int group = 0; group <= ngroups; ++group) { + groupBounds(group) = GroupBoundAt(group_bounds, group); + } +} + +// Copy the class-owned group bounds back out (e.g. for HDF5 export). +inline void ExportGroupBounds(MeanUtilsDataBox &groupBounds, + const MeanUtilsDataBox &storedBounds, + const int ngroups) { + groupBounds.resize(ngroups + 1); + for (int group = 0; group <= ngroups; ++group) { + groupBounds(group) = storedBounds(group); + } +} + +// Validate half-open group bounds [nu_g, nu_{g+1}): strictly increasing, +// nonnegative, with only the final bound permitted to be IEEE +infinity. +template +void ValidateGroupBounds(const GroupBoundsIndexer &group_bounds, + const int ngroups) { + if (ngroups <= 0) { + OPAC_ERROR("photons multigroup: ngroups must be positive"); + } + for (int group = 0; group <= ngroups; ++group) { + const Real bound = GroupBoundAt(group_bounds, group); + if (std::isnan(bound)) { + OPAC_ERROR("photons multigroup: group bounds must be finite " + "or IEEE +infinity"); + } + if (std::isinf(bound) && bound < 0.) { + OPAC_ERROR("photons multigroup: group bounds may not be -infinity"); + } + if (group == 0) { + if (!(bound >= 0.)) { + OPAC_ERROR("photons multigroup: first group bound must be " + "nonnegative"); + } + } else if (!(bound > GroupBoundAt(group_bounds, group - 1))) { + OPAC_ERROR("photons multigroup: group bounds must be strictly " + "increasing"); + } + if (!std::isfinite(bound) && group != ngroups) { + OPAC_ERROR("photons multigroup: only the final group bound may be " + "IEEE +infinity"); + } + } +} + +// Sample a group's frequency range on a logarithmic, midpoint-rule grid, +// invoking sample_op(nu, weight) for each of NNuPerGroup points. Extreme +// bounds (nu=0, +infinity, or far from the thermal peak) are clamped to a +// thermal-aware window so the integral stays well-conditioned. +template +void ForEachGroupFrequencySample(const Real temp, const Real nuMin, + const Real nuMax, const int NNuPerGroup, + SampleOp &&sample_op) { + // For [0, ∞) or very wide ranges, use thermal-aware sampling + // For reasonable finite ranges, integrate over the full group bounds + const Real nu_thermal_min = 1.e-3 * PC::kb * temp / PC::h; + const Real nu_thermal_max = 1.e3 * PC::kb * temp / PC::h; + + // Determine if we need special handling + const bool is_lower_extreme = (nuMin == 0.) || (nuMin < 0.1 * nu_thermal_min); + const bool is_upper_extreme = + !std::isfinite(nuMax) || (nuMax > 10. * nu_thermal_max); + + // Set integration bounds, but ensure they're valid + Real nu_sample_min = is_lower_extreme ? nu_thermal_min : nuMin; + Real nu_sample_max = is_upper_extreme ? nu_thermal_max : nuMax; + + // If thermal-aware bounds are invalid, use a small but valid range within + // group bounds + if (nu_sample_min >= nu_sample_max) { + if (std::isfinite(nuMax) && nuMax > 0.) { + // Group is [0 or small, nuMax]: sample near nuMax + nu_sample_min = 0.5 * nuMax; + nu_sample_max = nuMax; + } else { + // Group extends to infinity: sample around thermal peak + nu_sample_min = 0.1 * nu_thermal_max; + nu_sample_max = nu_thermal_max; + } + } + + // Use logarithmic spacing with midpoint rule + const Real lNuMin = ToLog(nu_sample_min); + const Real lNuMax = ToLog(nu_sample_max); + const Real dlnu = (lNuMax - lNuMin) / NNuPerGroup; + for (int inu = 0; inu < NNuPerGroup; ++inu) { + const Real lnu = lNuMin + (inu + 0.5) * dlnu; + const Real nu = FromLog(lnu); + sample_op(nu, nu * dlnu); + } +} + +// Evaluate the Planck function B_nu and its temperature derivative dB_nu/dT, +// switching to the Wien closed form deep in the tail (see wien_tail_x). +template +void ThermalWeightsAtNu(const PlanckDistribution &dist, const Real temp, + const Real nu, Real &B, Real &dBdT) { + const Real x = PC::h * nu / (PC::kb * temp); + if (x < wien_tail_x) { + B = dist.ThermalDistributionOfTNu(temp, nu); + dBdT = dist.DThermalDistributionOfTNuDT(temp, nu); + return; + } + + const Real expMinusX = std::exp(-x); + B = (2. * PC::h * nu * nu * nu / (PC::c * PC::c)) * expMinusX; + dBdT = 2. * PC::h * PC::h * nu * nu * nu * nu * expMinusX / + (temp * temp * PC::c * PC::c * PC::kb); +} + +} // namespace impl +} // namespace photons +} // namespace singularity + +#endif // SINGULARITY_OPAC_PHOTONS_MEAN_PHOTON_UTILS_ diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp index 7762a85..7740090 100644 --- a/singularity-opac/photons/mean_photon_variant.hpp +++ b/singularity-opac/photons/mean_photon_variant.hpp @@ -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. @@ -12,9 +12,9 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_MEAN_PHOTON_VARIANT_ #define SINGULARITY_OPAC_PHOTONS_MEAN_PHOTON_VARIANT_ +// This file was made in part with generative AI. #include @@ -96,19 +96,77 @@ class MeanVariant { PORTABLE_INLINE_FUNCTION Real AbsorptionCoefficient(const Real rho, const Real temp, const int gmode = Rosseland) const { + return AbsorptionCoefficient(rho, temp, 0, gmode); + } + + PORTABLE_INLINE_FUNCTION + Real Emissivity(const Real rho, const Real temp, const int gmode = Rosseland, + Real *lambda = nullptr) const { + return PortsOfCall::visit( + [=](const auto &opac) { return opac.Emissivity(rho, temp, gmode); }, + opac_); + } + + PORTABLE_INLINE_FUNCTION + int ngroups() const noexcept { + return PortsOfCall::visit([](const auto &opac) { return opac.ngroups(); }, + opac_); + } + + PORTABLE_INLINE_FUNCTION + bool HasGroupBounds() const noexcept { + return PortsOfCall::visit( + [](const auto &opac) { return opac.HasGroupBounds(); }, opac_); + } + + PORTABLE_INLINE_FUNCTION + Real PlanckGroupAbsorptionCoefficient(const Real rho, const Real temp, + const int group) const { + return AbsorptionCoefficient(rho, temp, group, Planck); + } + + PORTABLE_INLINE_FUNCTION + Real RosselandGroupAbsorptionCoefficient(const Real rho, const Real temp, + const int group) const { + return AbsorptionCoefficient(rho, temp, group, Rosseland); + } + + PORTABLE_INLINE_FUNCTION + Real AbsorptionCoefficient(const Real rho, const Real temp, const int group, + const int gmode) const { return PortsOfCall::visit( [=](const auto &opac) { - return opac.AbsorptionCoefficient(rho, temp, gmode); + return opac.AbsorptionCoefficient(rho, temp, group, gmode); }, opac_); } PORTABLE_INLINE_FUNCTION - Real Emissivity(const Real rho, const Real temp, const int gmode = Rosseland, - Real *lambda = nullptr) const { + int GroupOfNu(const Real nu) const { + return PortsOfCall::visit( + [=](const auto &opac) { return opac.GroupOfNu(nu); }, opac_); + } + + PORTABLE_INLINE_FUNCTION + Real PlanckGroupAbsorptionCoefficientFromNu(const Real rho, const Real temp, + const Real nu) const { + return AbsorptionCoefficientFromNu(rho, temp, nu, Planck); + } + + PORTABLE_INLINE_FUNCTION + Real RosselandGroupAbsorptionCoefficientFromNu(const Real rho, + const Real temp, + const Real nu) const { + return AbsorptionCoefficientFromNu(rho, temp, nu, Rosseland); + } + + PORTABLE_INLINE_FUNCTION + Real AbsorptionCoefficientFromNu(const Real rho, const Real temp, + const Real nu, + const int gmode = Rosseland) const { return PortsOfCall::visit( [=](const auto &opac) { - return opac.Emissivity(rho, temp, gmode, lambda); + return opac.AbsorptionCoefficientFromNu(rho, temp, nu, gmode); }, opac_); } diff --git a/singularity-opac/photons/mean_s_opacity_photons.hpp b/singularity-opac/photons/mean_s_opacity_photons.hpp index 7fb6ff3..40d13f3 100644 --- a/singularity-opac/photons/mean_s_opacity_photons.hpp +++ b/singularity-opac/photons/mean_s_opacity_photons.hpp @@ -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. @@ -12,207 +12,348 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_MEAN_S_OPACITY_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_MEAN_S_OPACITY_PHOTONS_ +// This file was made in part with generative AI. +#include #include +#include +#include +#include #include -#include +#include #include #include #include -#include - -#include +#include +#include #include #include +#include namespace singularity { namespace photons { namespace impl { -#define EPS (10.0 * std::numeric_limits::min()) - // TODO(BRR) Note: It is assumed that lambda is constant for all densities and // temperatures -template +template class MeanSOpacity { - public: + using PC = pc; using DataBox = Spiner::DataBox; MeanSOpacity() = default; - template + + template MeanSOpacity(const SOpacity &s_opac, const Real lRhoMin, const Real lRhoMax, const int NRho, const Real lTMin, const Real lTMax, const int NT, - Real *lambda = nullptr) { - MeanSOpacityImpl_(s_opac, lRhoMin, lRhoMax, NRho, lTMin, - lTMax, NT, -1., -1., 100, lambda); + const GroupBoundsIndexer &group_bounds, const int ngroups, + const int NNuPerGroup = 64, Real *lambda = nullptr) { + MeanSOpacityImpl_(s_opac, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, + group_bounds, ngroups, NNuPerGroup, lambda); } - template - MeanSOpacity(const SOpacity &s_opac, const Real lRhoMin, const Real lRhoMax, - const int NRho, const Real lTMin, const Real lTMax, const int NT, - Real lNuMin, Real lNuMax, const int NNu, - Real *lambda = nullptr) { - MeanSOpacityImpl_(s_opac, lRhoMin, lRhoMax, NRho, lTMin, - lTMax, NT, lNuMin, lNuMax, NNu, lambda); + template + MeanSOpacity(const DataBox &sigmaPlanck, const DataBox &sigmaRosseland, + const GroupBoundsIndexer &group_bounds) { + LoadScatteringTables_(sigmaPlanck, sigmaRosseland, group_bounds); } #ifdef SPINER_USE_HDF - MeanSOpacity(const std::string &filename) : filename_(filename.c_str()) { + MeanSOpacity(const std::string &filename) { + DataBox sigmaPlanck; + DataBox sigmaRosseland; + DataBox groupBounds; herr_t status = H5_SUCCESS; hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); - status += lkappaPlanck_.loadHDF(file, SP5::MeanSOpac::PlanckMeanSOpacity); status += - lkappaRosseland_.loadHDF(file, SP5::MeanSOpac::RosselandMeanSOpacity); + sigmaPlanck.loadHDF(file, SP5::MultigroupSOpac::PlanckGroupSOpacity); + status += sigmaRosseland.loadHDF( + file, SP5::MultigroupSOpac::RosselandGroupSOpacity); + status += groupBounds.loadHDF(file, SP5::MultigroupSOpac::GroupBounds); status += H5Fclose(file); if (status != H5_SUCCESS) { OPAC_ERROR("photons::MeanSOpacity: HDF5 error\n"); } + + LoadScatteringTables_(sigmaPlanck, sigmaRosseland, groupBounds); + groupBounds.finalize(); + sigmaPlanck.finalize(); + sigmaRosseland.finalize(); } void Save(const std::string &filename) const { + DataBox sigmaPlanck; + DataBox sigmaRosseland; + DataBox groupBounds; + ExportScatteringTables_(sigmaPlanck, sigmaRosseland); + ExportGroupBounds(groupBounds, groupBounds_, ngroups_); + herr_t status = H5_SUCCESS; hid_t file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - status += lkappaPlanck_.saveHDF(file, SP5::MeanSOpac::PlanckMeanSOpacity); status += - lkappaRosseland_.saveHDF(file, SP5::MeanSOpac::RosselandMeanSOpacity); + sigmaPlanck.saveHDF(file, SP5::MultigroupSOpac::PlanckGroupSOpacity); + status += sigmaRosseland.saveHDF( + file, SP5::MultigroupSOpac::RosselandGroupSOpacity); + status += groupBounds.saveHDF(file, SP5::MultigroupSOpac::GroupBounds); status += H5Fclose(file); + sigmaPlanck.finalize(); + sigmaRosseland.finalize(); + groupBounds.finalize(); + if (status != H5_SUCCESS) { OPAC_ERROR("photons::MeanSOpacity: HDF5 error\n"); } } #endif - PORTABLE_INLINE_FUNCTION void PrintParams() const { - printf("Mean scattering opacity\n"); + PORTABLE_INLINE_FUNCTION + void PrintParams() const { + printf("Photon multigroup scattering opacity. ngroups = %d\n", ngroups_); } MeanSOpacity GetOnDevice() { MeanSOpacity other; - other.lkappaPlanck_ = Spiner::getOnDeviceDataBox(lkappaPlanck_); - other.lkappaRosseland_ = Spiner::getOnDeviceDataBox(lkappaRosseland_); + other.lsigmaPlanck_ = Spiner::getOnDeviceDataBox(lsigmaPlanck_); + other.lsigmaRosseland_ = Spiner::getOnDeviceDataBox(lsigmaRosseland_); + other.groupBounds_ = Spiner::getOnDeviceDataBox(groupBounds_); + other.ngroups_ = ngroups_; return other; } void Finalize() { - lkappaPlanck_.finalize(); - lkappaRosseland_.finalize(); + lsigmaPlanck_.finalize(); + lsigmaRosseland_.finalize(); + groupBounds_.finalize(); + } + + PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants + GetRuntimePhysicalConstants() const { + return RuntimePhysicalConstants(PC()); + } + + PORTABLE_INLINE_FUNCTION + int ngroups() const noexcept { return ngroups_; } + + PORTABLE_INLINE_FUNCTION + bool HasGroupBounds() const noexcept { return true; } + + PORTABLE_INLINE_FUNCTION + Real PlanckGroupScatteringCoefficient(const Real rho, const Real temp, + const int group) const { + return GroupScatteringCoefficient_(lsigmaPlanck_, rho, temp, group); + } + + PORTABLE_INLINE_FUNCTION + Real RosselandGroupScatteringCoefficient(const Real rho, const Real temp, + const int group) const { + return GroupScatteringCoefficient_(lsigmaRosseland_, rho, temp, group); + } + + PORTABLE_INLINE_FUNCTION + Real ScatteringCoefficient(const Real rho, const Real temp, const int group, + const int gmode = Rosseland) const { + return (gmode == Planck) + ? PlanckGroupScatteringCoefficient(rho, temp, group) + : RosselandGroupScatteringCoefficient(rho, temp, group); + } + + PORTABLE_INLINE_FUNCTION + int GroupOfNu(const Real nu) const { + if (!(nu >= GroupBoundAt(groupBounds_, 0) && + nu <= GroupBoundAt(groupBounds_, ngroups_))) { + OPAC_ERROR("photons::MeanSOpacity: frequency is outside group bounds"); + } + return GroupOfNuImpl(groupBounds_, ngroups_, nu); } PORTABLE_INLINE_FUNCTION - Real PlanckMeanTotalScatteringCoefficient(const Real rho, - const Real temp) const { - Real lRho = toLog_(rho); - Real lT = toLog_(temp); - return rho * fromLog_(lkappaPlanck_.interpToReal(lRho, lT)); + Real PlanckGroupScatteringCoefficientFromNu(const Real rho, const Real temp, + const Real nu) const { + return ScatteringCoefficientFromNu(rho, temp, nu, Planck); } PORTABLE_INLINE_FUNCTION - Real RosselandMeanTotalScatteringCoefficient(const Real rho, - const Real temp) const { - Real lRho = toLog_(rho); - Real lT = toLog_(temp); - return rho * fromLog_(lkappaRosseland_.interpToReal(lRho, lT)); + Real RosselandGroupScatteringCoefficientFromNu(const Real rho, + const Real temp, + const Real nu) const { + return ScatteringCoefficientFromNu(rho, temp, nu, Rosseland); + } + + PORTABLE_INLINE_FUNCTION + Real ScatteringCoefficientFromNu(const Real rho, const Real temp, + const Real nu, + const int gmode = Rosseland) const { + return ScatteringCoefficient(rho, temp, GroupOfNu(nu), gmode); } private: - template + PORTABLE_INLINE_FUNCTION + Real GroupScatteringCoefficient_(const DataBox &lsigma, const Real rho, + const Real temp, const int group) const { + const Real lRho = ToLog(rho); + const Real lT = ToLog(temp); + return rho * FromLog(lsigma.interpToReal(lRho, lT, group)); + } + + void ValidateScatteringTables_(const DataBox &sigmaPlanck, + const DataBox &sigmaRosseland) const { + if (sigmaPlanck.rank() != 3 || sigmaRosseland.rank() != 3) { + OPAC_ERROR("photons::MeanSOpacity: scattering tables must be rank 3"); + } + for (int dim = 1; dim <= 3; ++dim) { + if (sigmaPlanck.dim(dim) != sigmaRosseland.dim(dim)) { + OPAC_ERROR("photons::MeanSOpacity: table dimensions do not match"); + } + } + if (sigmaPlanck.dim(1) <= 0) { + OPAC_ERROR("photons::MeanSOpacity: ngroups must be positive"); + } + if (sigmaPlanck.range(1) != sigmaRosseland.range(1) || + sigmaPlanck.range(2) != sigmaRosseland.range(2)) { + OPAC_ERROR("photons::MeanSOpacity: table ranges do not match"); + } + } + + template + void LoadScatteringTables_(const DataBox &sigmaPlanck, + const DataBox &sigmaRosseland, + const GroupBoundsIndexer &group_bounds) { + ValidateScatteringTables_(sigmaPlanck, sigmaRosseland); + ngroups_ = sigmaPlanck.dim(1); + ValidateGroupBounds(group_bounds, ngroups_); + SetGroupBounds(groupBounds_, group_bounds, ngroups_); + lsigmaPlanck_.copyMetadata(sigmaPlanck); + lsigmaRosseland_.copyMetadata(sigmaRosseland); + for (int i = 0; i < sigmaPlanck.size(); ++i) { + lsigmaPlanck_(i) = ToLog(sigmaPlanck(i)); + lsigmaRosseland_(i) = ToLog(sigmaRosseland(i)); + } + } + + void ExportScatteringTables_(DataBox &sigmaPlanck, + DataBox &sigmaRosseland) const { + sigmaPlanck.copyMetadata(lsigmaPlanck_); + sigmaRosseland.copyMetadata(lsigmaRosseland_); + for (int i = 0; i < lsigmaPlanck_.size(); ++i) { + sigmaPlanck(i) = FromLog(lsigmaPlanck_(i)); + sigmaRosseland(i) = FromLog(lsigmaRosseland_(i)); + } + } + + template void MeanSOpacityImpl_(const SOpacity &s_opac, const Real lRhoMin, const Real lRhoMax, const int NRho, const Real lTMin, - const Real lTMax, const int NT, Real lNuMin, - Real lNuMax, const int NNu, Real *lambda = nullptr) { - ThermalDistribution dist; - - lkappaPlanck_.resize(NRho, NT); - lkappaPlanck_.setRange(0, lTMin, lTMax, NT); - lkappaPlanck_.setRange(1, lRhoMin, lRhoMax, NRho); - lkappaRosseland_.copyMetadata(lkappaPlanck_); - - // Fill tables - for (int iRho = 0; iRho < NRho; ++iRho) { - Real lRho = lkappaPlanck_.range(1).x(iRho); - Real rho = fromLog_(lRho); - for (int iT = 0; iT < NT; ++iT) { - Real lT = lkappaPlanck_.range(0).x(iT); - Real T = fromLog_(lT); - Real kappaPlanckNum = 0.; - Real kappaPlanckDenom = 0.; - Real kappaRosselandNum = 0.; - Real kappaRosselandDenom = 0.; - if (AUTOFREQ) { - lNuMin = toLog_(1.e-3 * pc::kb * fromLog_(lTMin) / pc::h); - lNuMax = toLog_(1.e3 * pc::kb * fromLog_(lTMax) / pc::h); - } - const Real dlnu = (lNuMax - lNuMin) / (NNu - 1); - // Integrate over frequency - for (int inu = 0; inu < NNu; ++inu) { - const Real weight = - (inu == 0 || inu == NNu - 1) ? 0.5 : 1.; // Trapezoidal rule - const Real lnu = lNuMin + inu * dlnu; - const Real nu = fromLog_(lnu); - const Real alpha = - s_opac.TotalScatteringCoefficient(rho, T, nu, lambda); - const Real B = dist.ThermalDistributionOfTNu(T, nu); - const Real dBdT = dist.DThermalDistributionOfTNuDT(T, nu); - kappaPlanckNum += weight * alpha / rho * B * nu * dlnu; - kappaPlanckDenom += weight * B * nu * dlnu; - - if (alpha > singularity_opac::robust::SMALL()) { - kappaRosselandNum += weight * - singularity_opac::robust::ratio(rho, alpha) * - dBdT * nu * dlnu; - kappaRosselandDenom += weight * dBdT * nu * dlnu; - } - } + const Real lTMax, const int NT, + const GroupBoundsIndexer &group_bounds, + const int ngroups, const int NNuPerGroup, + Real *lambda = nullptr) { +#ifndef NDEBUG + auto RPC = RuntimePhysicalConstants(PC()); + auto opc = s_opac.GetRuntimePhysicalConstants(); + assert(RPC == opc && "Physical constants are the same"); +#endif + + if (NNuPerGroup < 2) { + OPAC_ERROR("photons::MeanSOpacity: NNuPerGroup must be at least 2"); + } + ValidateGroupBounds(group_bounds, ngroups); + + ngroups_ = ngroups; + SetGroupBounds(groupBounds_, group_bounds, ngroups_); + lsigmaPlanck_.resize(NRho, NT, ngroups_); + lsigmaPlanck_.setRange(1, lTMin, lTMax, NT); + lsigmaPlanck_.setRange(2, lRhoMin, lRhoMax, NRho); + lsigmaRosseland_.copyMetadata(lsigmaPlanck_); + + PlanckDistribution dist; + std::vector planckDenom(ngroups_, 0.); + std::vector rosselandDenom(ngroups_, 0.); + + for (int iT = 0; iT < NT; ++iT) { + const Real lT = lsigmaPlanck_.range(1).x(iT); + const Real T = FromLog(lT); + + for (int group = 0; group < ngroups_; ++group) { + Real Baccum = 0.; + Real dBdTaccum = 0.; + const Real nuMin = GroupBoundAt(group_bounds, group); + const Real nuMax = GroupBoundAt(group_bounds, group + 1); + ForEachGroupFrequencySample( + T, nuMin, nuMax, NNuPerGroup, [&](const Real nu, const Real dnu) { + Real B = 0.; + Real dBdT = 0.; + ThermalWeightsAtNu(dist, T, nu, B, dBdT); + Baccum += B * dnu; + dBdTaccum += dBdT * dnu; + }); + + planckDenom[group] = Baccum; + rosselandDenom[group] = dBdTaccum; + } - Real kappaPlanck = - singularity_opac::robust::ratio(kappaPlanckNum, kappaPlanckDenom); - Real kappaRosseland = kappaPlanck > singularity_opac::robust::SMALL() - ? singularity_opac::robust::ratio( - kappaRosselandDenom, kappaRosselandNum) - : 0.; - lkappaPlanck_(iRho, iT) = toLog_(kappaPlanck); - lkappaRosseland_(iRho, iT) = toLog_(kappaRosseland); - if (std::isnan(lkappaPlanck_(iRho, iT)) || - std::isnan(lkappaRosseland_(iRho, iT))) { - OPAC_ERROR("photons::MeanSOpacity: NAN in opacity evaluations"); + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real lRho = lsigmaPlanck_.range(2).x(iRho); + const Real rho = FromLog(lRho); + + for (int group = 0; group < ngroups_; ++group) { + Real sigmaPlanckNum = 0.; + Real sigmaRosselandNum = 0.; + const Real nuMin = GroupBoundAt(group_bounds, group); + const Real nuMax = GroupBoundAt(group_bounds, group + 1); + ForEachGroupFrequencySample( + T, nuMin, nuMax, NNuPerGroup, [&](const Real nu, const Real dnu) { + const Real sigma = + s_opac.TotalScatteringCoefficient(rho, T, nu, lambda); + Real B = 0.; + Real dBdT = 0.; + ThermalWeightsAtNu(dist, T, nu, B, dBdT); + sigmaPlanckNum += sigma / rho * B * dnu; + + if (sigma > singularity_opac::robust::SMALL()) { + sigmaRosselandNum += + singularity_opac::robust::ratio(rho, sigma) * dBdT * dnu; + } + }); + + const Real sigmaPlanck = singularity_opac::robust::ratio( + sigmaPlanckNum, planckDenom[group]); + const Real sigmaRosseland = + (rosselandDenom[group] > singularity_opac::robust::SMALL() && + sigmaRosselandNum > singularity_opac::robust::SMALL()) + ? singularity_opac::robust::ratio(rosselandDenom[group], + sigmaRosselandNum) + : 0.; + + lsigmaPlanck_(iRho, iT, group) = ToLog(sigmaPlanck); + lsigmaRosseland_(iRho, iT, group) = ToLog(sigmaRosseland); + if (std::isnan(lsigmaPlanck_(iRho, iT, group)) || + std::isnan(lsigmaRosseland_(iRho, iT, group))) { + OPAC_ERROR("photons::MeanSOpacity: NAN in opacity evaluations"); + } } } } } - PORTABLE_INLINE_FUNCTION Real toLog_(const Real x) const { - return std::log10(std::abs(x) + EPS); - } - PORTABLE_INLINE_FUNCTION Real fromLog_(const Real lx) const { - return std::pow(10., lx); - } - DataBox lkappaPlanck_; - DataBox lkappaRosseland_; - const char *filename_; -}; -#undef EPS + DataBox lsigmaPlanck_; + DataBox lsigmaRosseland_; + DataBox groupBounds_; + int ngroups_ = 0; +}; } // namespace impl -using MeanSOpacityScaleFree = - impl::MeanSOpacity, - PhysicalConstantsUnity>; -using MeanSOpacityCGS = - impl::MeanSOpacity, - PhysicalConstantsCGS>; -using MeanSOpacity = impl::MeanSVariant>; +using MeanSOpacityBase = impl::MeanSOpacity; } // namespace photons } // namespace singularity -#endif // SINGULARITY_OPAC_PHOTONS_MEAN_S_OPACITY_PHOTONS__ +#endif // SINGULARITY_OPAC_PHOTONS_MEAN_S_OPACITY_PHOTONS_ diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp index a7222b1..c871652 100644 --- a/singularity-opac/photons/non_cgs_photons.hpp +++ b/singularity-opac/photons/non_cgs_photons.hpp @@ -1,5 +1,5 @@ // ====================================================================== -// © 2021-2024. 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. @@ -12,13 +12,14 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_NON_CGS_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_NON_CGS_PHOTONS_ +// This file was made in part with generative AI. #include #include #include +#include #include #include @@ -242,79 +243,91 @@ class MeanNonCGSUnits { using PC = typename MeanOpac::PC; MeanNonCGSUnits() = default; - MeanNonCGSUnits(MeanOpac &&mean_opac, const Real time_unit, + MeanNonCGSUnits(MeanOpac &&multigroup_opac, const Real time_unit, const Real mass_unit, const Real length_unit, 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), + : multigroup_opac_(std::forward(multigroup_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_ * time_unit_ / - mass_unit_) {} + freq_unit_(1. / time_unit_) {} auto GetOnDevice() { - return MeanNonCGSUnits(mean_opac_.GetOnDevice(), time_unit_, + return MeanNonCGSUnits(multigroup_opac_.GetOnDevice(), time_unit_, mass_unit_, length_unit_, temp_unit_); } - inline void Finalize() noexcept { mean_opac_.Finalize(); } + inline void Finalize() noexcept { multigroup_opac_.Finalize(); } + + PORTABLE_INLINE_FUNCTION + int ngroups() const noexcept { return multigroup_opac_.ngroups(); } PORTABLE_INLINE_FUNCTION - int nlambda() const noexcept { return mean_opac_.nlambda(); } + bool HasGroupBounds() const noexcept { + return multigroup_opac_.HasGroupBounds(); + } + + PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants + GetRuntimePhysicalConstants() const { + return multigroup_opac_.GetRuntimePhysicalConstants(); + } #ifdef SPINER_USE_HDF void Save(const std::string &filename) const { - return mean_opac_.Save(filename); + return multigroup_opac_.Save(filename); } #endif PORTABLE_INLINE_FUNCTION - Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const { - const Real alpha = mean_opac_.PlanckMeanAbsorptionCoefficient( - rho_unit_ * rho, temp_unit_ * temp); - // alpha output in units of 1/cm. Want to convert out of CGS. - // multiplication by length_unit converts length to cm. - // division converts length from cm to unit system. - // thus multiplication converts (1/cm) to unit system. - return alpha * length_unit_; + Real PlanckGroupAbsorptionCoefficient(const Real rho, const Real temp, + const int group) const { + return AbsorptionCoefficient(rho, temp, group, Planck); } PORTABLE_INLINE_FUNCTION - Real RosselandMeanAbsorptionCoefficient(const Real rho, - const Real temp) const { - const Real alpha = mean_opac_.RosselandMeanAbsorptionCoefficient( - rho_unit_ * rho, temp_unit_ * temp); - // alpha output in units of 1/cm. Want to convert out of CGS. - // multiplication by length_unit converts length to cm. - // division converts length from cm to unit system. - // thus multiplication converts (1/cm) to unit system. - return alpha * length_unit_; + Real RosselandGroupAbsorptionCoefficient(const Real rho, const Real temp, + const int group) const { + return AbsorptionCoefficient(rho, temp, group, Rosseland); } PORTABLE_INLINE_FUNCTION - Real AbsorptionCoefficient(const Real rho, const Real temp, + Real AbsorptionCoefficient(const Real rho, const Real temp, const int group, const int gmode = Rosseland) const { - const Real alpha = mean_opac_.AbsorptionCoefficient( - rho_unit_ * rho, temp_unit_ * temp, gmode); + const Real alpha = multigroup_opac_.AbsorptionCoefficient( + rho_unit_ * rho, temp_unit_ * temp, group, 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_unit_ * rho, temp_unit_ * temp, gmode); - return J * inv_emiss_unit_; + int GroupOfNu(const Real nu) const { + return multigroup_opac_.GroupOfNu(nu * freq_unit_); } - PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants - GetRuntimePhysicalConstants() const { - return RuntimePhysicalConstants(PhysicalConstantsCGS(), time_unit_, - mass_unit_, length_unit_, temp_unit_); + PORTABLE_INLINE_FUNCTION + Real PlanckGroupAbsorptionCoefficientFromNu(const Real rho, const Real temp, + const Real nu) const { + return AbsorptionCoefficientFromNu(rho, temp, nu, Planck); + } + + PORTABLE_INLINE_FUNCTION + Real RosselandGroupAbsorptionCoefficientFromNu(const Real rho, + const Real temp, + const Real nu) const { + return AbsorptionCoefficientFromNu(rho, temp, nu, Rosseland); + } + + PORTABLE_INLINE_FUNCTION + Real AbsorptionCoefficientFromNu(const Real rho, const Real temp, + const Real nu, + const int gmode = Rosseland) const { + const Real alpha = multigroup_opac_.AbsorptionCoefficientFromNu( + rho_unit_ * rho, temp_unit_ * temp, nu * freq_unit_, gmode); + return alpha * length_unit_; } private: - MeanOpac mean_opac_; + MeanOpac multigroup_opac_; Real time_unit_, mass_unit_, length_unit_, temp_unit_; - Real rho_unit_, inv_emiss_unit_; + Real rho_unit_, freq_unit_; }; } // namespace photons diff --git a/singularity-opac/photons/non_cgs_s_photons.hpp b/singularity-opac/photons/non_cgs_s_photons.hpp index 44ae2fc..e9ab9b2 100644 --- a/singularity-opac/photons/non_cgs_s_photons.hpp +++ b/singularity-opac/photons/non_cgs_s_photons.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_NON_CGS_S_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_NON_CGS_S_PHOTONS_ diff --git a/singularity-opac/photons/opac_photons.hpp b/singularity-opac/photons/opac_photons.hpp index 50918a5..a9fd66b 100644 --- a/singularity-opac/photons/opac_photons.hpp +++ b/singularity-opac/photons/opac_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. @@ -12,19 +12,18 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_OPAC_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_OPAC_PHOTONS_ +// This file was made in part with generative AI. #include #include +#include #include #include #include #include -#include - namespace singularity { namespace photons { diff --git a/singularity-opac/photons/photon_s_variant.hpp b/singularity-opac/photons/photon_s_variant.hpp index a7abc2a..07a7e5d 100644 --- a/singularity-opac/photons/photon_s_variant.hpp +++ b/singularity-opac/photons/photon_s_variant.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_PHOTON_S_VARIANT_ #define SINGULARITY_OPAC_PHOTONS_PHOTON_S_VARIANT_ diff --git a/singularity-opac/photons/photon_variant.hpp b/singularity-opac/photons/photon_variant.hpp index 33df636..1151b2d 100644 --- a/singularity-opac/photons/photon_variant.hpp +++ b/singularity-opac/photons/photon_variant.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_PHOTON_VARIANT_ #define SINGULARITY_OPAC_PHOTONS_PHOTON_VARIANT_ @@ -267,7 +266,8 @@ class Variant { PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { - return PortsOfCall::visit([](const auto &opac) { return opac.nlambda(); }, opac_); + return PortsOfCall::visit([](const auto &opac) { return opac.nlambda(); }, + opac_); } template @@ -277,12 +277,13 @@ class Variant { PORTABLE_INLINE_FUNCTION void PrintParams() const noexcept { - return PortsOfCall::visit([](const auto &opac) { return opac.PrintParams(); }, - opac_); + return PortsOfCall::visit( + [](const auto &opac) { return opac.PrintParams(); }, opac_); } inline void Finalize() noexcept { - return PortsOfCall::visit([](auto &opac) { return opac.Finalize(); }, opac_); + return PortsOfCall::visit([](auto &opac) { return opac.Finalize(); }, + opac_); } }; diff --git a/singularity-opac/photons/powerlaw_opacity_photons.hpp b/singularity-opac/photons/powerlaw_opacity_photons.hpp index 69f5860..d156113 100644 --- a/singularity-opac/photons/powerlaw_opacity_photons.hpp +++ b/singularity-opac/photons/powerlaw_opacity_photons.hpp @@ -12,10 +12,8 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_ - // This file was made in part with generative AI. #include @@ -36,22 +34,23 @@ class PowerLawOpacity { PowerLawOpacity() = default; PowerLawOpacity(const Real kappa0, const Real rho_exp, const Real temp_exp, - const Real nu_exp = 0., const Real nu_ref = 1., const Real nu_off = 0., - const Real rho_ref = 1., const Real rho_off = 0., - const Real temp_ref = 1., const Real temp_off = 0., - const bool do_stim_emit = false) + const Real nu_exp = 0., const Real nu_ref = 1., + const Real nu_off = 0., const Real rho_ref = 1., + const Real rho_off = 0., const Real temp_ref = 1., + const Real temp_off = 0., const bool do_stim_emit = false) : PowerLawOpacity(PlanckDistribution{}, kappa0, rho_exp, temp_exp, nu_exp, nu_ref, nu_off, rho_ref, rho_off, temp_ref, temp_off, do_stim_emit) {} PowerLawOpacity(const PlanckDistribution &dist, const Real kappa0, const Real rho_exp, const Real temp_exp, - const Real nu_exp = 0., const Real nu_ref = 1., const Real nu_off = 0., - const Real rho_ref = 1., const Real rho_off = 0., - const Real temp_ref = 1., const Real temp_off = 0., - const bool do_stim_emit = false) + const Real nu_exp = 0., const Real nu_ref = 1., + const Real nu_off = 0., const Real rho_ref = 1., + const Real rho_off = 0., const Real temp_ref = 1., + const Real temp_off = 0., const bool do_stim_emit = false) : dist_(dist), kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp), - rho_ref_(rho_ref), rho_off_(rho_off), temp_ref_(temp_ref), temp_off_(temp_off), - nu_exp_(nu_exp), nu_ref_(nu_ref), nu_off_(nu_off), do_stim_emit_(do_stim_emit) { + rho_ref_(rho_ref), rho_off_(rho_off), temp_ref_(temp_ref), + temp_off_(temp_off), nu_exp_(nu_exp), nu_ref_(nu_ref), nu_off_(nu_off), + do_stim_emit_(do_stim_emit) { if (!(rho_ref_ > 0.)) { OPAC_ERROR("PowerLawOpacity: rho_ref must be positive"); } @@ -80,9 +79,8 @@ class PowerLawOpacity { printf("Power law opacity. kappa0 = %g rho_exp = %g temp_exp = %g " "nu_exp = %g nu_ref = %g nu_off = %g rho_ref = %g rho_off = %g" "temp_ref = %g temp_off = %g do_stim_emit = %d\n", - kappa0_, rho_exp_, temp_exp_, nu_exp_, nu_ref_, - nu_off_, rho_ref_, rho_off_, temp_ref_, temp_off_, - do_stim_emit_); + kappa0_, rho_exp_, temp_exp_, nu_exp_, nu_ref_, nu_off_, rho_ref_, + rho_off_, temp_ref_, temp_off_, do_stim_emit_); } inline void Finalize() noexcept {} @@ -235,7 +233,8 @@ class PowerLawOpacity { PORTABLE_INLINE_FUNCTION Real OpacityScale_(const Real rho, const Real temp, const Real nu) const { const Real freq_plaw = std::pow((nu + nu_off_) / nu_ref_, nu_exp_); - const Real stim_fact = do_stim_emit_ ? -std::expm1(-(pc::h * nu / (pc::kb * temp))) : 1.0; + const Real stim_fact = + do_stim_emit_ ? -std::expm1(-(pc::h * nu / (pc::kb * temp))) : 1.0; return OpacityPrefactor_(rho, temp) * freq_plaw * stim_fact; } @@ -246,17 +245,17 @@ class PowerLawOpacity { return kappa0_ * std::pow(rhom, rho_exp_) * std::pow(tempm, temp_exp_); } - Real kappa0_; // Opacity scale. Units depend on nu_exp and nu_ref. - Real rho_exp_; // Power law index of density - Real temp_exp_; // Power law index of temperature - Real rho_ref_; // Density normalization for rho_exp. Units of g/cm^3 - Real rho_off_; // Density offset (same units as rho_ref). Units of g/cm^3 - Real temp_ref_; // Temperature normalization for temp_exp. Units of K - Real temp_off_; // Temperature offset (same units as temp_ref). Units of K - Real nu_exp_; // Power law index of frequency - Real nu_ref_; // Frequency normalization for nu_exp. Units of 1/s - Real nu_off_; // Frequency offset (same units as nu_ref). Units of 1/s - bool do_stim_emit_; // indicator to use stimulated (LTE) emission factor + Real kappa0_; // Opacity scale. Units depend on nu_exp and nu_ref. + Real rho_exp_; // Power law index of density + Real temp_exp_; // Power law index of temperature + Real rho_ref_; // Density normalization for rho_exp. Units of g/cm^3 + Real rho_off_; // Density offset (same units as rho_ref). Units of g/cm^3 + Real temp_ref_; // Temperature normalization for temp_exp. Units of K + Real temp_off_; // Temperature offset (same units as temp_ref). Units of K + Real nu_exp_; // Power law index of frequency + Real nu_ref_; // Frequency normalization for nu_exp. Units of 1/s + Real nu_off_; // Frequency offset (same units as nu_ref). Units of 1/s + bool do_stim_emit_; // indicator to use stimulated (LTE) emission factor PlanckDistribution dist_; }; diff --git a/singularity-opac/photons/powerlaw_s_opacity_photons.hpp b/singularity-opac/photons/powerlaw_s_opacity_photons.hpp index 6185ea4..c801b0e 100644 --- a/singularity-opac/photons/powerlaw_s_opacity_photons.hpp +++ b/singularity-opac/photons/powerlaw_s_opacity_photons.hpp @@ -12,10 +12,8 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_POWERLAW_S_OPACITY_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_POWERLAW_S_OPACITY_PHOTONS_ - // This file was partly copied from a file made in part with generative AI. #include @@ -33,13 +31,14 @@ class PowerLawSOpacity { public: PowerLawSOpacity() = default; PowerLawSOpacity(const Real kappa0, const Real rho_exp, const Real temp_exp, - const Real nu_exp = 0., const Real nu_ref = 1., const Real nu_off = 0., - const Real rho_ref = 1., const Real rho_off = 0., - const Real temp_ref = 1., const Real temp_off = 0., - const Real avg_particle_mass = 1.) + const Real nu_exp = 0., const Real nu_ref = 1., + const Real nu_off = 0., const Real rho_ref = 1., + const Real rho_off = 0., const Real temp_ref = 1., + const Real temp_off = 0., const Real avg_particle_mass = 1.) : kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp), - rho_ref_(rho_ref), rho_off_(rho_off), temp_ref_(temp_ref), temp_off_(temp_off), - nu_exp_(nu_exp), nu_ref_(nu_ref), nu_off_(nu_off), apm_(avg_particle_mass) { + rho_ref_(rho_ref), rho_off_(rho_off), temp_ref_(temp_ref), + temp_off_(temp_off), nu_exp_(nu_exp), nu_ref_(nu_ref), nu_off_(nu_off), + apm_(avg_particle_mass) { if (!(nu_ref_ > 0.)) { OPAC_ERROR("PowerLawSOpacity: nu_ref must be positive"); } @@ -65,12 +64,12 @@ class PowerLawSOpacity { int nlambda() const noexcept { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const noexcept { - printf("Power law scattering opacity. kappa0 = %g rho_exp = %g temp_exp = %g " - "nu_exp = %g nu_ref = %g nu_off = %g rho_ref = %g rho_off = %g" - "temp_ref = %g temp_off = %g avg particle mass = %g\n", - kappa0_, rho_exp_, temp_exp_, nu_exp_, nu_ref_, - nu_off_, rho_ref_, rho_off_, temp_ref_, temp_off_, - apm_); + printf( + "Power law scattering opacity. kappa0 = %g rho_exp = %g temp_exp = %g " + "nu_exp = %g nu_ref = %g nu_off = %g rho_ref = %g rho_off = %g" + "temp_ref = %g temp_off = %g avg particle mass = %g\n", + kappa0_, rho_exp_, temp_exp_, nu_exp_, nu_ref_, nu_off_, rho_ref_, + rho_off_, temp_ref_, temp_off_, apm_); } inline void Finalize() noexcept {} @@ -96,7 +95,8 @@ class PowerLawSOpacity { private: PORTABLE_INLINE_FUNCTION Real OpacityScale_(const Real rho, const Real temp, const Real nu) const { - return OpacityPrefactor_(rho, temp) * std::pow((nu + nu_off_) / nu_ref_, nu_exp_); + return OpacityPrefactor_(rho, temp) * + std::pow((nu + nu_off_) / nu_ref_, nu_exp_); } PORTABLE_INLINE_FUNCTION @@ -106,17 +106,17 @@ class PowerLawSOpacity { return kappa0_ * std::pow(rhom, rho_exp_) * std::pow(tempm, temp_exp_); } - Real kappa0_; // Opacity scale. Units depend on nu_exp and nu_ref. - Real rho_exp_; // Power law index of density - Real temp_exp_; // Power law index of temperature - Real rho_ref_; // Density normalization for rho_exp. Units of g/cm^3 - Real rho_off_; // Density offset (same units as rho_ref). Units of g/cm^3 - Real temp_ref_; // Temperature normalization for temp_exp. Units of K - Real temp_off_; // Temperature offset (same units as temp_ref). Units of K - Real nu_exp_; // Power law index of frequency - Real nu_ref_; // Frequency normalization for nu_exp. Units of 1/s - Real nu_off_; // Frequency offset (same units as nu_ref). Units of 1/s - Real apm_; // Mean molecular weight. Units of g + Real kappa0_; // Opacity scale. Units depend on nu_exp and nu_ref. + Real rho_exp_; // Power law index of density + Real temp_exp_; // Power law index of temperature + Real rho_ref_; // Density normalization for rho_exp. Units of g/cm^3 + Real rho_off_; // Density offset (same units as rho_ref). Units of g/cm^3 + Real temp_ref_; // Temperature normalization for temp_exp. Units of K + Real temp_off_; // Temperature offset (same units as temp_ref). Units of K + Real nu_exp_; // Power law index of frequency + Real nu_ref_; // Frequency normalization for nu_exp. Units of 1/s + Real nu_off_; // Frequency offset (same units as nu_ref). Units of 1/s + Real apm_; // Mean molecular weight. Units of g }; } // namespace photons diff --git a/singularity-opac/photons/s_opac_photons.hpp b/singularity-opac/photons/s_opac_photons.hpp index 257837e..287dc04 100644 --- a/singularity-opac/photons/s_opac_photons.hpp +++ b/singularity-opac/photons/s_opac_photons.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_S_OPAC_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_S_OPAC_PHOTONS_ @@ -32,9 +31,9 @@ using GrayS = GraySOpacity<>; using ThomsonS = ThomsonSOpacity<>; using PowerLawS = PowerLawSOpacity<>; -using SOpacity = impl::S_Variant, NonCGSUnitsS, - NonCGSUnitsS>; +using SOpacity = + impl::S_Variant, + NonCGSUnitsS, NonCGSUnitsS>; } // namespace photons } // namespace singularity diff --git a/singularity-opac/photons/thermal_distributions_photons.hpp b/singularity-opac/photons/thermal_distributions_photons.hpp index 05d2565..592191e 100644 --- a/singularity-opac/photons/thermal_distributions_photons.hpp +++ b/singularity-opac/photons/thermal_distributions_photons.hpp @@ -12,10 +12,8 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_THERMAL_DISTRIBUTIONS_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_THERMAL_DISTRIBUTIONS_PHOTONS_ - // This file was made in part with generative AI. #include diff --git a/singularity-opac/photons/thomson_s_opacity_photons.hpp b/singularity-opac/photons/thomson_s_opacity_photons.hpp index e0a588c..6173583 100644 --- a/singularity-opac/photons/thomson_s_opacity_photons.hpp +++ b/singularity-opac/photons/thomson_s_opacity_photons.hpp @@ -12,7 +12,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. // ====================================================================== - #ifndef SINGULARITY_OPAC_PHOTONS_THOMSON_S_OPACITY_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_THOMSON_S_OPACITY_PHOTONS_ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f3152fb..21a754a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -56,15 +56,10 @@ PRIVATE test_thomson_s_opacities.cpp test_chebyshev.cpp test_spiner_opac_neutrinos.cpp - test_mean_opacities.cpp + test_multigroup_opacities.cpp test_variant.cpp ) -# Link data file directory for Zhu table test -file(CREATE_LINK - ${CMAKE_SOURCE_DIR}/singularity-opac/photons/example_ascii/kap_plaw.txt - ${CMAKE_CURRENT_BINARY_DIR}/kap_plaw.txt SYMBOLIC) - target_link_libraries(${PROJECT_NAME}_unit_tests PRIVATE catch2_define diff --git a/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp deleted file mode 100644 index 97f9415..0000000 --- a/test/test_mean_opacities.cpp +++ /dev/null @@ -1,872 +0,0 @@ -// ====================================================================== -// © 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. -// 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. -// ====================================================================== - -// This file was generated in part with the assistance of generative AI. - -#include -#include - -#include - -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace singularity; - -using pc = PhysicalConstantsCGS; - -#ifdef PORTABILITY_STRATEGY_KOKKOS -using atomic_view = Kokkos::MemoryTraits; -#endif - -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-100); -} -constexpr Real EPS_TEST = 1e-3; -template -PORTABLE_INLINE_FUNCTION bool IsWrong(const T &a, const T &b) { - // We only care if the data is different and significantly nonzero. - constexpr Real ZERO = - std::numeric_limits::epsilon() / (EPS_TEST * EPS_TEST); - return ((std::isnan(a) || std::isnan(b)) || - ((std::abs(a) > ZERO || std::abs(b) > ZERO) && - FractionalDifference(a, b) > EPS_TEST)); -} - -TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") { - const std::string grayname = "mean_gray.sp5"; - - WHEN("We initialize a mean neutrino opacity") { - constexpr Real MeV2K = 1e6 * pc::eV / pc::kb; - constexpr Real MeV2Hz = 1e6 * pc::eV / pc::h; - constexpr Real rho = 1e11; // g/cc - constexpr Real temp = 10 * MeV2K; // 10 MeV - constexpr Real Ye = 0.1; - constexpr RadiationType type = RadiationType::NU_ELECTRON; - constexpr Real nu = 1.25 * MeV2Hz; // 1 MeV - - constexpr int nT = 10; - const Real lRhoMin = std::log10(0.1 * rho); - const Real lRhoMax = std::log10(10. * rho); - constexpr int NRho = 2; - 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; - constexpr int NYe = 10; - - constexpr Real kappa = 1.e-20; - - neutrinos::Gray opac_host(kappa); - neutrinos::Opacity opac = opac_host.GetOnDevice(); - - neutrinos::MeanOpacity mean_opac_host = neutrinos::MeanOpacityBase( - 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") { - int n_wrong_h = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View n_wrong_d("wrong"); -#else - PortableMDArray n_wrong_d(&n_wrong_h, 1); -#endif - - portableFor( - "calc mean opacities", 0, 100, PORTABLE_LAMBDA(const int &i) { - Real alphaPlanck = - mean_opac.PlanckMeanAbsorptionCoefficient(rho, temp, Ye, type); - Real alphaRosseland = - mean_opac.PlanckMeanAbsorptionCoefficient(rho, temp, Ye, type); - if (FractionalDifference(kappa * rho, alphaPlanck) > EPS_TEST) { - n_wrong_d() += 1; - } - if (FractionalDifference(kappa * rho, alphaRosseland) > EPS_TEST) { - n_wrong_d() += 1; - } - }); - -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(n_wrong_h, n_wrong_d); -#endif - REQUIRE(n_wrong_h == 0); - } - -#ifdef SPINER_USE_HDF - THEN("We can save to disk and reload") { - mean_opac.Save(grayname); - neutrinos::MeanOpacity mean_opac_host_load(grayname); - AND_THEN("The reloaded table matches the gray opacities") { - - auto mean_opac_load = mean_opac_host_load.GetOnDevice(); - - int n_wrong = 0; - portableReduce( - "rebuilt table vs gray", 0, NRho, 0, NT, 0, NYe, 0, NEUTRINO_NTYPES, - PORTABLE_LAMBDA(const int iRho, const int iT, const int iYe, - const int itp, int &accumulate) { - const Real lRho = - lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho; - const Real rho = std::pow(10, lRho); - const Real lT = lTMin + (lTMax - lTMin) / (NT - 1) * iT; - const Real T = std::pow(10, lT); - const Real Ye = YeMin + (YeMax - YeMin) / (NYe - 1) * iYe; - const RadiationType type = Idx2RadType(itp); - - const Real kappaPgray = - mean_opac.PlanckMeanAbsorptionCoefficient(rho, T, Ye, type); - const Real kappaPload = - mean_opac_load.PlanckMeanAbsorptionCoefficient(rho, T, Ye, - type); - - const Real kappaRgray = - mean_opac.RosselandMeanAbsorptionCoefficient(rho, T, Ye, - type); - const Real kappaRload = - mean_opac_load.RosselandMeanAbsorptionCoefficient(rho, T, Ye, - type); - - if (IsWrong(kappaPgray, kappaPload)) { - accumulate += 1; - } - if (IsWrong(kappaRgray, kappaRload)) { - accumulate += 1; - } - }, - n_wrong); - REQUIRE(n_wrong == 0); - } - } -#endif - - THEN("We can create an opacity object with funny units") { - constexpr Real time_unit = 123; - constexpr Real mass_unit = 456; - constexpr Real length_unit = 789; - constexpr Real temp_unit = 276; - constexpr Real rho_unit = - mass_unit / (length_unit * length_unit * length_unit); - - auto mean_opac_host_base = - neutrinos::MeanOpacityBase(opac_host, lRhoMin, lRhoMax, NRho, lTMin, - lTMax, NT, YeMin, YeMax, NYe); - auto funny_units_host = - neutrinos::MeanNonCGSUnits( - std::forward(mean_opac_host_base), - time_unit, mass_unit, length_unit, temp_unit); - - auto funny_units = funny_units_host.GetOnDevice(); - - THEN("We can retrieve physical constants in code units") { - auto noncgs_rpc = funny_units.GetRuntimePhysicalConstants(); - REQUIRE(FractionalDifference(noncgs_rpc.length, length_unit) < - EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.time, time_unit) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.mass, mass_unit) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.temp, temp_unit) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.na, 6.022141e+23) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.alpha, 7.297353e-03) < - EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.h, 2.871060e-33) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.hbar, 4.569434e-34) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.kb, 2.030877e-18) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.r_gas, 4.431243e+03) < - EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.qe, 4.803205e-10) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.c, 4.673571e+09) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.g_newt, 4.508065e-15) < - EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.me, 1.997672e-30) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.mp, 3.668030e-27) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.mn, 3.673086e-27) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.amu, 3.641533e-27) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.sb, 1.342760e+09) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.ar, 1.149237e+00) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.eV, 8.538896e-17) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.Fc, 1.189435e-09) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.nu_sigma0, 2.829094e-74) < - EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.gA, -1.23) < EPS_TEST); - } - - int n_wrong_h = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View n_wrong_d("wrong"); -#else - PortableMDArray n_wrong_d(&n_wrong_h, 1); -#endif - - portableFor( - "compare different units", 0, 100, PORTABLE_LAMBDA(const int &i) { - Real alphaPlanck = - mean_opac.PlanckMeanAbsorptionCoefficient(rho, temp, Ye, type); - Real alphaRosseland = mean_opac.RosselandMeanAbsorptionCoefficient( - rho, temp, Ye, type); - Real alphaPlanckFunny = funny_units.PlanckMeanAbsorptionCoefficient( - rho / rho_unit, temp / temp_unit, Ye, type); - Real alphaRosselandFunny = - funny_units.RosselandMeanAbsorptionCoefficient( - rho / rho_unit, temp / temp_unit, Ye, type); - if (FractionalDifference(alphaPlanck, alphaPlanckFunny / - length_unit) > EPS_TEST) { - n_wrong_d() += 1; - } - if (FractionalDifference(alphaRosseland, - alphaRosselandFunny / length_unit) > - EPS_TEST) { - n_wrong_d() += 1; - } - }); - -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(n_wrong_h, n_wrong_d); -#endif - REQUIRE(n_wrong_h == 0); - } - - opac.Finalize(); - } -} - -TEST_CASE("Mean neutrino scattering opacities", "[MeanNeutrinosS]") { - const std::string grayname = "mean_gray_s.sp5"; - - WHEN("We initialize a mean neutrino scattering opacity") { - constexpr Real MeV2K = 1e6 * pc::eV / pc::kb; - constexpr Real MeV2Hz = 1e6 * pc::eV / pc::h; - constexpr Real rho = 1e11; // g/cc - constexpr Real temp = 10 * MeV2K; // 10 MeV - constexpr Real Ye = 0.1; - constexpr RadiationType type = RadiationType::NU_ELECTRON; - constexpr Real nu = 1.25 * MeV2Hz; // 1 MeV - - constexpr int nT = 10; - const Real lRhoMin = std::log10(0.1 * rho); - const Real lRhoMax = std::log10(10. * rho); - constexpr int NRho = 2; - 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; - constexpr int NYe = 10; - - constexpr Real sigma = 1.e-20; - - constexpr Real avg_particle_mass = pc::mp; - - neutrinos::GrayS opac_host(sigma, avg_particle_mass); - neutrinos::SOpacity opac = opac_host.GetOnDevice(); - - neutrinos::MeanSOpacityCGS mean_opac_host( - opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin, YeMax, NYe); - auto mean_opac = mean_opac_host.GetOnDevice(); - - 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"); -#else - PortableMDArray n_wrong_d(&n_wrong_h, 1); -#endif - - portableFor( - "calc mean opacities", 0, 100, PORTABLE_LAMBDA(const int &i) { - Real alphaPlanck = mean_opac.PlanckMeanTotalScatteringCoefficient( - rho, temp, Ye, type); - Real alphaRosseland = - mean_opac.PlanckMeanTotalScatteringCoefficient(rho, temp, Ye, - type); - if (FractionalDifference(sigma * rho / pc::mp, alphaPlanck) > - EPS_TEST) { - n_wrong_d() += 1; - } - if (FractionalDifference(sigma * rho / pc::mp, alphaRosseland) > - EPS_TEST) { - n_wrong_d() += 1; - } - }); - -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(n_wrong_h, n_wrong_d); -#endif - REQUIRE(n_wrong_h == 0); - } - -#ifdef SPINER_USE_HDF - THEN("We can save to disk and reload") { - mean_opac.Save(grayname); - neutrinos::MeanSOpacityCGS mean_opac_host_load(grayname); - AND_THEN("The reloaded table matches the gray opacities") { - - auto mean_opac_load = mean_opac_host_load.GetOnDevice(); - - int n_wrong = 0; - portableReduce( - "rebuilt table vs gray", 0, NRho, 0, NT, 0, NYe, 0, NEUTRINO_NTYPES, - PORTABLE_LAMBDA(const int iRho, const int iT, const int iYe, - const int itp, int &accumulate) { - const Real lRho = - lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho; - const Real rho = std::pow(10, lRho); - const Real lT = lTMin + (lTMax - lTMin) / (NT - 1) * iT; - const Real T = std::pow(10, lT); - const Real Ye = YeMin + (YeMax - YeMin) / (NYe - 1) * iYe; - const RadiationType type = Idx2RadType(itp); - - const Real kappaPgray = - mean_opac.PlanckMeanTotalScatteringCoefficient(rho, T, Ye, - type); - const Real kappaPload = - mean_opac_load.PlanckMeanTotalScatteringCoefficient(rho, T, - Ye, type); - - const Real kappaRgray = - mean_opac.RosselandMeanTotalScatteringCoefficient(rho, T, Ye, - type); - const Real kappaRload = - mean_opac_load.RosselandMeanTotalScatteringCoefficient( - rho, T, Ye, type); - - if (IsWrong(kappaPgray, kappaPload)) { - accumulate += 1; - } - if (IsWrong(kappaRgray, kappaRload)) { - accumulate += 1; - } - }, - n_wrong); - REQUIRE(n_wrong == 0); - } - } -#endif - - THEN("We can create an opacity object with funny units") { - constexpr Real time_unit = 123; - constexpr Real mass_unit = 456; - constexpr Real length_unit = 789; - constexpr Real temp_unit = 276; - constexpr Real rho_unit = - mass_unit / (length_unit * length_unit * length_unit); - - auto funny_units_host = - neutrinos::MeanNonCGSUnitsS( - std::forward(mean_opac_host), time_unit, - mass_unit, length_unit, temp_unit); - - auto funny_units = funny_units_host.GetOnDevice(); - - int n_wrong_h = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View n_wrong_d("wrong"); -#else - PortableMDArray n_wrong_d(&n_wrong_h, 1); -#endif - - portableFor( - "compare different units", 0, 100, PORTABLE_LAMBDA(const int &i) { - Real alphaPlanck = mean_opac.PlanckMeanTotalScatteringCoefficient( - rho, temp, Ye, type); - Real alphaRosseland = - mean_opac.RosselandMeanTotalScatteringCoefficient(rho, temp, Ye, - type); - Real alphaPlanckFunny = - funny_units.PlanckMeanTotalScatteringCoefficient( - rho / rho_unit, temp / temp_unit, Ye, type); - Real alphaRosselandFunny = - funny_units.RosselandMeanTotalScatteringCoefficient( - rho / rho_unit, temp / temp_unit, Ye, type); - if (FractionalDifference(alphaPlanck, alphaPlanckFunny / - length_unit) > EPS_TEST) { - n_wrong_d() += 1; - } - if (FractionalDifference(alphaRosseland, - alphaRosselandFunny / length_unit) > - EPS_TEST) { - n_wrong_d() += 1; - } - }); - -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(n_wrong_h, n_wrong_d); -#endif - REQUIRE(n_wrong_h == 0); - } - - opac.Finalize(); - } -} - -TEST_CASE("Mean photon opacities", "[MeanPhotons]") { - const std::string grayname = "mean_gray_photons.sp5"; - - WHEN("We initialize a mean photon opacity") { - constexpr Real MeV2K = 1e6 * pc::eV / pc::kb; - constexpr Real MeV2Hz = 1e6 * pc::eV / pc::h; - constexpr Real rho = 1e0; // g/cc - constexpr Real temp = 1.e5; // K - constexpr Real nu = 1.e15; // Hz - - constexpr int nT = 10; - const Real lRhoMin = std::log10(0.1 * rho); - const Real lRhoMax = std::log10(10. * rho); - constexpr int NRho = 2; - 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; - - // Load analytic Opacity Eos - photons::Gray opac_host(kappa); - // Put it on the GPU - photons::Opacity opac = opac_host.GetOnDevice(); - - // Load a table Opacity - photons::MeanOpacity mean_opac_host = photons::MeanOpacityBase( - opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT); - // Put it on the GPU - auto mean_opac = mean_opac_host.GetOnDevice(); - - // This is a part of the catch2 control flow for unit testing. - THEN("The emissivity per nu omega is consistent with the emissity per nu") { - int n_wrong_h = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - // Creates a device-side view with the atomic trait when using Kokkos! - // This allocates data on device. - Kokkos::View n_wrong_d("wrong"); -#else - // Create something that works like a View but is fine on CPU - PortableMDArray n_wrong_d(&n_wrong_h, 1); -#endif - - // This is a for loop that can be device side (default) or host side - portableFor( - // The for loop is called "calc mean opacities" - // It loops over an index from 0 to 100 - "calc mean opacities", 0, 100, - // Define a Lambda that takes an argument (we'll cal it i) - PORTABLE_LAMBDA(const int &i) { - // Inside of this lambda, we are calculating these opacities - Real alphaPlanck = - mean_opac.PlanckMeanAbsorptionCoefficient(rho, temp); - Real alphaRosseland = - mean_opac.RosselandMeanAbsorptionCoefficient(rho, temp); - // Check if they're wrong and if they are increase the wrong counter - if (FractionalDifference(kappa * rho, alphaPlanck) > 1.e-12) { - n_wrong_d() += 1; - } - if (FractionalDifference(kappa * rho, alphaRosseland) > 1.e-12) { - n_wrong_d() += 1; - } - }); - -#ifdef PORTABILITY_STRATEGY_KOKKOS - // Copies counter from device to host - Kokkos::deep_copy(n_wrong_h, n_wrong_d); -#endif - REQUIRE(n_wrong_h == 0); - } - -#ifdef SPINER_USE_HDF - THEN("We can save to disk and reload") { - mean_opac.Save(grayname); - photons::MeanOpacity mean_opac_host_load = - photons::MeanOpacityBase(grayname); - AND_THEN("The reloaded table matches the gray opacities") { - - auto mean_opac_load = mean_opac_host_load.GetOnDevice(); - - int n_wrong = 0; - portableReduce( - "rebuilt table vs gray", 0, NRho, 0, NT, 0, 0, - PORTABLE_LAMBDA(const int iRho, const int iT, const int igarbage, - int &accumulate) { - const Real lRho = - lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho; - const Real rho = std::pow(10, lRho); - const Real lT = lTMin + (lTMax - lTMin) / (NT - 1) * iT; - const Real T = std::pow(10, lT); - - const Real kappaPgray = - mean_opac.PlanckMeanAbsorptionCoefficient(rho, T); - const Real kappaPload = - mean_opac_load.PlanckMeanAbsorptionCoefficient(rho, T); - - const Real kappaRgray = - mean_opac.RosselandMeanAbsorptionCoefficient(rho, T); - const Real kappaRload = - mean_opac_load.RosselandMeanAbsorptionCoefficient(rho, T); - - if (IsWrong(kappaPgray, kappaPload)) { - accumulate += 1; - } - if (IsWrong(kappaRgray, kappaRload)) { - accumulate += 1; - } - }, - n_wrong); - REQUIRE(n_wrong == 0); - } - } -#endif - - THEN("We can create an opacity object with funny units") { - constexpr Real time_unit = 123; - constexpr Real mass_unit = 456; - constexpr Real length_unit = 789; - constexpr Real temp_unit = 276; - constexpr Real rho_unit = - mass_unit / (length_unit * length_unit * length_unit); - - auto mean_opac_host_base = photons::MeanOpacityBase( - opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT); - auto funny_units_host = - photons::MeanNonCGSUnits( - std::forward(mean_opac_host_base), - time_unit, mass_unit, length_unit, temp_unit); - - auto funny_units = funny_units_host.GetOnDevice(); - - THEN("We can retrieve physical constants in code units") { - auto noncgs_rpc = funny_units.GetRuntimePhysicalConstants(); - REQUIRE(FractionalDifference(noncgs_rpc.length, length_unit) < - EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.time, time_unit) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.mass, mass_unit) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.temp, temp_unit) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.na, 6.022141e+23) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.alpha, 7.297353e-03) < - EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.h, 2.871060e-33) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.hbar, 4.569434e-34) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.kb, 2.030877e-18) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.r_gas, 4.431243e+03) < - EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.qe, 4.803205e-10) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.c, 4.673571e+09) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.g_newt, 4.508065e-15) < - EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.me, 1.997672e-30) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.mp, 3.668030e-27) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.mn, 3.673086e-27) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.amu, 3.641533e-27) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.sb, 1.342760e+09) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.ar, 1.149237e+00) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.eV, 8.538896e-17) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.Fc, 1.189435e-09) < EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.nu_sigma0, 2.829094e-74) < - EPS_TEST); - REQUIRE(FractionalDifference(noncgs_rpc.gA, -1.23) < EPS_TEST); - } - - int n_wrong_h = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View n_wrong_d("wrong"); -#else - PortableMDArray n_wrong_d(&n_wrong_h, 1); -#endif - - portableFor( - "compare different units", 0, 100, PORTABLE_LAMBDA(const int &i) { - Real alphaPlanck = - mean_opac.PlanckMeanAbsorptionCoefficient(rho, temp); - Real alphaRosseland = - mean_opac.RosselandMeanAbsorptionCoefficient(rho, temp); - Real alphaPlanckFunny = funny_units.PlanckMeanAbsorptionCoefficient( - rho / rho_unit, temp / temp_unit); - Real alphaRosselandFunny = - funny_units.RosselandMeanAbsorptionCoefficient( - rho / rho_unit, temp / temp_unit); - if (FractionalDifference(alphaPlanck, alphaPlanckFunny / - length_unit) > EPS_TEST) { - n_wrong_d() += 1; - } - if (FractionalDifference(alphaRosseland, - alphaRosselandFunny / length_unit) > - EPS_TEST) { - n_wrong_d() += 1; - } - }); - -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(n_wrong_h, n_wrong_d); -#endif - REQUIRE(n_wrong_h == 0); - } - - opac.Finalize(); - } -} - -TEST_CASE("Mean photon scattering opacities", "[MeanPhotonS]") { - const std::string grayname = "mean_gray_s.sp5"; - - WHEN("We initialize a mean photon scattering opacity") { - constexpr Real MeV2K = 1e6 * pc::eV / pc::kb; - constexpr Real MeV2Hz = 1e6 * pc::eV / pc::h; - constexpr Real rho = 1e0; // g/cc - constexpr Real temp = 1.e5; // K - constexpr Real nu = 1.e15; // Hz - - constexpr int nT = 10; - const Real lRhoMin = std::log10(0.1 * rho); - const Real lRhoMax = std::log10(10. * rho); - constexpr int NRho = 2; - 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; - - constexpr Real avg_particle_mass = pc::mp / 2.; - - photons::GrayS opac_host(sigma, avg_particle_mass); - photons::SOpacity opac = opac_host.GetOnDevice(); - - photons::MeanSOpacityCGS mean_opac_host(opac_host, lRhoMin, lRhoMax, NRho, - lTMin, lTMax, NT); - auto mean_opac = mean_opac_host.GetOnDevice(); - - THEN("The emissivity per nu omega is consistent with the emissity per nu") { - int n_wrong_h = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View n_wrong_d("wrong"); -#else - PortableMDArray n_wrong_d(&n_wrong_h, 1); -#endif - - portableFor( - "calc mean opacities", 0, 100, PORTABLE_LAMBDA(const int &i) { - Real alphaPlanck = - mean_opac.PlanckMeanTotalScatteringCoefficient(rho, temp); - Real alphaRosseland = - mean_opac.PlanckMeanTotalScatteringCoefficient(rho, temp); - if (FractionalDifference(sigma * rho / avg_particle_mass, - alphaPlanck) > EPS_TEST) { - n_wrong_d() += 1; - } - if (FractionalDifference(sigma * rho / avg_particle_mass, - alphaRosseland) > EPS_TEST) { - n_wrong_d() += 1; - } - }); - -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(n_wrong_h, n_wrong_d); -#endif - REQUIRE(n_wrong_h == 0); - } - -#ifdef SPINER_USE_HDF - THEN("We can save to disk and reload") { - mean_opac.Save(grayname); - photons::MeanSOpacityCGS mean_opac_host_load(grayname); - AND_THEN("The reloaded table matches the gray opacities") { - - auto mean_opac_load = mean_opac_host_load.GetOnDevice(); - - int n_wrong = 0; - portableReduce( - "rebuilt table vs gray", 0, NRho, 0, NT, 0, 0, - PORTABLE_LAMBDA(const int iRho, const int iT, const int igarbage, - int &accumulate) { - const Real lRho = - lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho; - const Real rho = std::pow(10, lRho); - const Real lT = lTMin + (lTMax - lTMin) / (NT - 1) * iT; - const Real T = std::pow(10, lT); - - const Real kappaPgray = - mean_opac.PlanckMeanTotalScatteringCoefficient(rho, T); - const Real kappaPload = - mean_opac_load.PlanckMeanTotalScatteringCoefficient(rho, T); - - const Real kappaRgray = - mean_opac.RosselandMeanTotalScatteringCoefficient(rho, T); - const Real kappaRload = - mean_opac_load.RosselandMeanTotalScatteringCoefficient(rho, - T); - - if (IsWrong(kappaPgray, kappaPload)) { - accumulate += 1; - } - if (IsWrong(kappaRgray, kappaRload)) { - accumulate += 1; - } - }, - n_wrong); - REQUIRE(n_wrong == 0); - } - } -#endif - - THEN("We can create an opacity object with funny units") { - constexpr Real time_unit = 123; - constexpr Real mass_unit = 456; - constexpr Real length_unit = 789; - constexpr Real temp_unit = 276; - constexpr Real rho_unit = - mass_unit / (length_unit * length_unit * length_unit); - - auto funny_units_host = photons::MeanNonCGSUnitsS( - std::forward(mean_opac_host), time_unit, - mass_unit, length_unit, temp_unit); - - auto funny_units = funny_units_host.GetOnDevice(); - - int n_wrong_h = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View n_wrong_d("wrong"); -#else - PortableMDArray n_wrong_d(&n_wrong_h, 1); -#endif - - portableFor( - "compare different units", 0, 100, PORTABLE_LAMBDA(const int &i) { - Real alphaPlanck = - mean_opac.PlanckMeanTotalScatteringCoefficient(rho, temp); - Real alphaRosseland = - mean_opac.RosselandMeanTotalScatteringCoefficient(rho, temp); - Real alphaPlanckFunny = - funny_units.PlanckMeanTotalScatteringCoefficient( - rho / rho_unit, temp / temp_unit); - Real alphaRosselandFunny = - funny_units.RosselandMeanTotalScatteringCoefficient( - rho / rho_unit, temp / temp_unit); - if (FractionalDifference(alphaPlanck, alphaPlanckFunny / - length_unit) > EPS_TEST) { - n_wrong_d() += 1; - } - if (FractionalDifference(alphaRosseland, - alphaRosselandFunny / length_unit) > - EPS_TEST) { - n_wrong_d() += 1; - } - }); - -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(n_wrong_h, n_wrong_d); -#endif - REQUIRE(n_wrong_h == 0); - } - - opac.Finalize(); - } -} - -TEST_CASE("ASCII-parsed Mean photon opacities", "[MeanPhotons]") { - const std::string grayname = "kap_plaw.txt"; - - 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 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 - - photons::MeanOpacity mean_opac_host = photons::MeanOpacity(grayname); - auto mean_opac = mean_opac_host.GetOnDevice(); - - THEN("The emissivity per nu omega is consistent with the emissity per nu") { - int n_wrong_h = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View n_wrong_d("wrong"); -#else - PortableMDArray n_wrong_d(&n_wrong_h, 1); -#endif - - // get a test value at min rho-T point - Real mross = - mean_opac.RosselandMeanAbsorptionCoefficient(rho_min, temp_min); - Real mplnk = mean_opac.PlanckMeanAbsorptionCoefficient(rho_min, temp_min); - - // check min rho-T point - if (FractionalDifference(rho_min * ross_at_min, mross) > EPS_TEST) { - n_wrong_d() += 1; - } - if (FractionalDifference(rho_min * plnk_at_min, mplnk) > EPS_TEST) { - n_wrong_d() += 1; - } - - // get a test value at max rho-T point - mross = mean_opac.RosselandMeanAbsorptionCoefficient(rho_max, temp_max); - mplnk = mean_opac.PlanckMeanAbsorptionCoefficient(rho_max, temp_max); - - // check min rho-T point - if (FractionalDifference(rho_max * ross_at_max, mross) > EPS_TEST) { - n_wrong_d() += 1; - } - if (FractionalDifference(rho_max * plnk_at_max, mplnk) > EPS_TEST) { - 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, 1); - - // 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 - REQUIRE(n_wrong_h == 0); - } - - mean_opac.Finalize(); - } -} diff --git a/test/test_multigroup_opacities.cpp b/test/test_multigroup_opacities.cpp new file mode 100644 index 0000000..0b9efa9 --- /dev/null +++ b/test/test_multigroup_opacities.cpp @@ -0,0 +1,1128 @@ +// ====================================================================== +// © 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 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. +// ====================================================================== + +// This file was made in part with generative AI. + +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include + +using namespace singularity; + +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-100); +} + +constexpr Real EPS_TEST = 5e-3; + +TEST_CASE("Photon multigroup gray opacities are exact", "[MultigroupPhotons]") { + constexpr Real rho = 1.; + constexpr Real temp = 1.e5; + constexpr Real kappa = 3.e-2; + constexpr int NRho = 3; + constexpr int NT = 4; + constexpr int ngroups = 4; + constexpr int nnu_per_group = 96; + const Real lRhoMin = std::log10(0.1 * rho); + const Real lRhoMax = std::log10(10. * rho); + const Real lTMin = std::log10(0.25 * temp); + const Real lTMax = std::log10(4. * temp); + const std::array group_bounds = {1.e12, 3.e12, 1.e13, + 1.e14, 3.e15}; + + photons::Gray opac_host(kappa); + photons::MeanOpacityBase multigroup_host( + opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, group_bounds, + ngroups, nnu_per_group); + auto multigroup = multigroup_host.GetOnDevice(); + + REQUIRE(multigroup_host.ngroups() == ngroups); + + int n_wrong = 0; + portableReduce( + "check multigroup gray opacity", 0, ngroups, + PORTABLE_LAMBDA(const int group, int &accum) { + const Real alpha_planck = + multigroup.PlanckGroupAbsorptionCoefficient(rho, temp, group); + const Real alpha_rosseland = + multigroup.RosselandGroupAbsorptionCoefficient(rho, temp, group); + const Real alpha_default = + multigroup.AbsorptionCoefficient(rho, temp, group); + if (FractionalDifference(alpha_planck, rho * kappa) > EPS_TEST) { + accum += 1; + } + if (FractionalDifference(alpha_rosseland, rho * kappa) > EPS_TEST) { + accum += 1; + } + if (FractionalDifference(alpha_default, rho * kappa) > EPS_TEST) { + accum += 1; + } + }, + n_wrong); + REQUIRE(n_wrong == 0); + + multigroup.Finalize(); + multigroup_host.Finalize(); +} + +TEST_CASE("Photon multigroup can be constructed from pretabulated Spiner data", + "[MultigroupPhotons]") { + using DataBox = Spiner::DataBox; + + constexpr int NRho = 2; + constexpr int NT = 2; + constexpr int ngroups = 3; + constexpr Real kappaP0 = 1.5e-2; + constexpr Real kappaR0 = 4.0e-3; + constexpr Real rho_exp_p = 0.25; + constexpr Real temp_exp_p = -0.5; + constexpr Real rho_exp_r = -0.2; + constexpr Real temp_exp_r = 0.3; + const std::array group_bounds = {2.e11, 3.e11, 4.e11, + 5.e11}; + const Real lRhoMin = -4.; + const Real lRhoMax = 2.; + const Real lTMin = 2.; + const Real lTMax = 8.; + + DataBox kappa_planck(NRho, NT, ngroups); + kappa_planck.setRange(1, lTMin, lTMax, NT); + kappa_planck.setRange(2, lRhoMin, lRhoMax, NRho); + DataBox kappa_rosseland; + kappa_rosseland.copyMetadata(kappa_planck); + + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real rho = std::pow(10., kappa_planck.range(2).x(iRho)); + for (int iT = 0; iT < NT; ++iT) { + const Real temp = std::pow(10., kappa_planck.range(1).x(iT)); + for (int group = 0; group < ngroups; ++group) { + const Real group_factor = group + 1.; + kappa_planck(iRho, iT, group) = kappaP0 * std::pow(rho, rho_exp_p) * + std::pow(temp, temp_exp_p) * + group_factor; + kappa_rosseland(iRho, iT, group) = kappaR0 * std::pow(rho, rho_exp_r) * + std::pow(temp, temp_exp_r) / + group_factor; + } + } + } + + photons::MeanOpacityBase multigroup_host(kappa_planck, kappa_rosseland, + group_bounds); + + const Real rho_test = std::pow(10., 0.5 * (lRhoMin + lRhoMax)); + const Real temp_test = std::pow(10., 0.5 * (lTMin + lTMax)); + for (int group = 0; group < ngroups; ++group) { + const Real group_factor = group + 1.; + const Real kappa_planck_expected = kappaP0 * std::pow(rho_test, rho_exp_p) * + std::pow(temp_test, temp_exp_p) * + group_factor; + const Real kappa_rosseland_expected = + kappaR0 * std::pow(rho_test, rho_exp_r) * + std::pow(temp_test, temp_exp_r) / group_factor; + + REQUIRE( + FractionalDifference(multigroup_host.PlanckGroupAbsorptionCoefficient( + rho_test, temp_test, group), + rho_test * kappa_planck_expected) < EPS_TEST); + REQUIRE(FractionalDifference( + multigroup_host.RosselandGroupAbsorptionCoefficient( + rho_test, temp_test, group), + rho_test * kappa_rosseland_expected) < EPS_TEST); + } + + constexpr Real nu_min = 2.e11; + constexpr Real nu_max = 4.e11; + const std::array tail_group_bounds = { + 0., nu_min, nu_max, std::numeric_limits::infinity()}; + const std::array nu_probe = { + 0.5 * nu_min, std::sqrt(nu_min * nu_max), 2. * nu_max}; + photons::MeanOpacityBase with_tail_bounds(kappa_planck, kappa_rosseland, + tail_group_bounds); + + REQUIRE(with_tail_bounds.HasGroupBounds()); + REQUIRE(with_tail_bounds.GroupOfNu(0.) == 0); + REQUIRE(with_tail_bounds.GroupOfNu(0.5 * nu_min) == 0); + REQUIRE(with_tail_bounds.GroupOfNu(nu_min) == 1); + REQUIRE(with_tail_bounds.GroupOfNu(nu_max) == 2); + REQUIRE(with_tail_bounds.GroupOfNu(2. * nu_max) == 2); + + for (int group = 0; group < ngroups; ++group) { + REQUIRE(FractionalDifference( + with_tail_bounds.PlanckGroupAbsorptionCoefficientFromNu( + rho_test, temp_test, nu_probe[group]), + multigroup_host.PlanckGroupAbsorptionCoefficient( + rho_test, temp_test, group)) < EPS_TEST); + REQUIRE(FractionalDifference( + with_tail_bounds.RosselandGroupAbsorptionCoefficientFromNu( + rho_test, temp_test, nu_probe[group]), + multigroup_host.RosselandGroupAbsorptionCoefficient( + rho_test, temp_test, group)) < EPS_TEST); + } + + with_tail_bounds.Finalize(); + multigroup_host.Finalize(); + kappa_planck.finalize(); + kappa_rosseland.finalize(); +} + +#ifdef SPINER_USE_HDF +TEST_CASE("Photon multigroup tables can round-trip through SP5 HDF", + "[MultigroupPhotons]") { + constexpr int NRho = 2; + constexpr int NT = 2; + constexpr int ngroups = 3; + const Real lRhoMin = -4.; + const Real lRhoMax = 2.; + const Real lTMin = 2.; + const Real lTMax = 8.; + constexpr Real nu_min = 2.e11; + constexpr Real nu_max = 4.e11; + const std::array group_bounds = { + 0., nu_min, nu_max, std::numeric_limits::infinity()}; + const std::array nu_probe = { + 0.5 * nu_min, std::sqrt(nu_min * nu_max), 2. * nu_max}; + using DataBox = Spiner::DataBox; + + DataBox kappa_planck(NRho, NT, ngroups); + kappa_planck.setRange(1, lTMin, lTMax, NT); + kappa_planck.setRange(2, lRhoMin, lRhoMax, NRho); + DataBox kappa_rosseland; + kappa_rosseland.copyMetadata(kappa_planck); + + for (int iRho = 0; iRho < NRho; ++iRho) { + for (int iT = 0; iT < NT; ++iT) { + for (int group = 0; group < ngroups; ++group) { + kappa_planck(iRho, iT, group) = + 1.e-2 * (1. + iRho + 2. * iT + 3. * group); + kappa_rosseland(iRho, iT, group) = + 5.e-3 * (2. + 2. * iRho + iT + group); + } + } + } + + photons::MeanOpacityBase saved(kappa_planck, kappa_rosseland, + group_bounds); + const char *filename = "multigroup-photon-table.sp5"; + saved.Save(filename); + photons::MeanOpacityBase loaded(filename); + + REQUIRE(loaded.HasGroupBounds()); + REQUIRE(loaded.ngroups() == ngroups); + REQUIRE(loaded.GroupOfNu(0.) == 0); + REQUIRE(loaded.GroupOfNu(0.5 * nu_min) == 0); + REQUIRE(loaded.GroupOfNu(nu_min) == 1); + REQUIRE(loaded.GroupOfNu(nu_max) == ngroups - 1); + REQUIRE(loaded.GroupOfNu(2. * nu_max) == ngroups - 1); + + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real rho_test = + std::pow(10., lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho); + for (int iT = 0; iT < NT; ++iT) { + const Real temp_test = + std::pow(10., lTMin + (lTMax - lTMin) / (NT - 1) * iT); + for (int group = 0; group < ngroups; ++group) { + const Real alpha_planck_expected = + rho_test * kappa_planck(iRho, iT, group); + const Real alpha_rosseland_expected = + rho_test * kappa_rosseland(iRho, iT, group); + REQUIRE(FractionalDifference(loaded.PlanckGroupAbsorptionCoefficient( + rho_test, temp_test, group), + alpha_planck_expected) < EPS_TEST); + REQUIRE(FractionalDifference(loaded.RosselandGroupAbsorptionCoefficient( + rho_test, temp_test, group), + alpha_rosseland_expected) < EPS_TEST); + REQUIRE( + FractionalDifference(loaded.PlanckGroupAbsorptionCoefficientFromNu( + rho_test, temp_test, nu_probe[group]), + alpha_planck_expected) < EPS_TEST); + REQUIRE(FractionalDifference( + loaded.RosselandGroupAbsorptionCoefficientFromNu( + rho_test, temp_test, nu_probe[group]), + alpha_rosseland_expected) < EPS_TEST); + } + } + } + + loaded.Finalize(); + saved.Finalize(); + kappa_planck.finalize(); + kappa_rosseland.finalize(); +} +#endif + +TEST_CASE("Photon multigroup frequency lookup uses half-open group bounds", + "[MultigroupPhotons]") { + using DataBox = Spiner::DataBox; + + constexpr int NRho = 2; + constexpr int NT = 2; + constexpr int ngroups = 3; + constexpr Real rho = 3.; + constexpr Real temp = 5.e4; + const Real lRhoMin = std::log10(0.25 * rho); + const Real lRhoMax = std::log10(4. * rho); + const Real lTMin = std::log10(0.5 * temp); + const Real lTMax = std::log10(2. * temp); + const std::array group_bounds = {1.e12, 2.e12, 4.e12, + 8.e12}; + + DataBox kappa_planck(NRho, NT, ngroups); + kappa_planck.setRange(1, lTMin, lTMax, NT); + kappa_planck.setRange(2, lRhoMin, lRhoMax, NRho); + DataBox kappa_rosseland; + kappa_rosseland.copyMetadata(kappa_planck); + + for (int group = 0; group < ngroups; ++group) { + for (int iT = 0; iT < NT; ++iT) { + for (int iRho = 0; iRho < NRho; ++iRho) { + kappa_planck(iRho, iT, group) = 1.e-2 * (group + 1.); + kappa_rosseland(iRho, iT, group) = 2.e-2 * (group + 1.); + } + } + } + + photons::MeanOpacityBase multigroup_host(kappa_planck, kappa_rosseland, + group_bounds); + photons::MeanOpacity multigroup = multigroup_host; + + REQUIRE(multigroup.HasGroupBounds()); + REQUIRE(multigroup.GroupOfNu(group_bounds[0]) == 0); + REQUIRE(multigroup.GroupOfNu(1.5e12) == 0); + REQUIRE(multigroup.GroupOfNu(group_bounds[1]) == 1); + REQUIRE(multigroup.GroupOfNu(group_bounds[2]) == 2); + REQUIRE(multigroup.GroupOfNu(group_bounds[ngroups]) == ngroups - 1); + + const Real nu_mid = std::sqrt(group_bounds[1] * group_bounds[2]); + REQUIRE( + FractionalDifference( + multigroup.PlanckGroupAbsorptionCoefficientFromNu(rho, temp, nu_mid), + multigroup.PlanckGroupAbsorptionCoefficient(rho, temp, 1)) < + EPS_TEST); + REQUIRE(FractionalDifference( + multigroup.RosselandGroupAbsorptionCoefficientFromNu(rho, temp, + nu_mid), + multigroup.RosselandGroupAbsorptionCoefficient(rho, temp, 1)) < + EPS_TEST); + REQUIRE(FractionalDifference( + multigroup.AbsorptionCoefficientFromNu(rho, temp, nu_mid), + multigroup.AbsorptionCoefficient(rho, temp, 1, photons::Rosseland)) < EPS_TEST); + + multigroup.Finalize(); + kappa_planck.finalize(); + kappa_rosseland.finalize(); +} + +TEST_CASE("Photon multigroup non-CGS wrapper converts units", + "[MultigroupPhotons]") { + using DataBox = Spiner::DataBox; + + constexpr Real rho = 2.; + constexpr Real temp = 3.e4; + constexpr int NRho = 2; + constexpr int NT = 3; + constexpr int ngroups = 3; + constexpr Real time_unit = 13.; + constexpr Real mass_unit = 17.; + constexpr Real length_unit = 19.; + constexpr Real temp_unit = 23.; + constexpr Real rho_unit = + mass_unit / (length_unit * length_unit * length_unit); + const std::array group_bounds = {1.e11, 2.e11, 4.e11, + 8.e11}; + const Real lRhoMin = std::log10(0.25 * rho); + const Real lRhoMax = std::log10(4. * rho); + const Real lTMin = std::log10(0.5 * temp); + const Real lTMax = std::log10(2. * temp); + + DataBox kappa_planck(NRho, NT, ngroups); + kappa_planck.setRange(1, lTMin, lTMax, NT); + kappa_planck.setRange(2, lRhoMin, lRhoMax, NRho); + DataBox kappa_rosseland; + kappa_rosseland.copyMetadata(kappa_planck); + + for (int i = 0; i < kappa_planck.size(); ++i) { + kappa_planck(i) = 7.e-3; + kappa_rosseland(i) = 5.e-3; + } + + photons::MeanOpacityBase reference_host(kappa_planck, kappa_rosseland, + group_bounds); + photons::MeanOpacityBase multigroup_base(kappa_planck, kappa_rosseland, + group_bounds); + auto funny_host = + photons::MeanNonCGSUnits( + std::move(multigroup_base), time_unit, mass_unit, length_unit, + temp_unit); + + REQUIRE(funny_host.ngroups() == ngroups); + REQUIRE(funny_host.HasGroupBounds()); + + for (int group = 0; group < ngroups; ++group) { + const Real alpha_planck_cgs = funny_host.PlanckGroupAbsorptionCoefficient( + rho / rho_unit, temp / temp_unit, group) / + length_unit; + const Real alpha_rosseland_cgs = + funny_host.RosselandGroupAbsorptionCoefficient( + rho / rho_unit, temp / temp_unit, group) / + length_unit; + + REQUIRE( + FractionalDifference(alpha_planck_cgs, + reference_host.PlanckGroupAbsorptionCoefficient( + rho, temp, group)) < EPS_TEST); + REQUIRE( + FractionalDifference(alpha_rosseland_cgs, + reference_host.RosselandGroupAbsorptionCoefficient( + rho, temp, group)) < EPS_TEST); + REQUIRE(FractionalDifference(funny_host.AbsorptionCoefficient( + rho / rho_unit, temp / temp_unit, group) / + length_unit, + reference_host.AbsorptionCoefficient( + rho, temp, group)) < EPS_TEST); + + const Real nu_mid = + std::sqrt(group_bounds[group] * group_bounds[group + 1]); + REQUIRE(funny_host.GroupOfNu(nu_mid * time_unit) == group); + REQUIRE(FractionalDifference( + funny_host.PlanckGroupAbsorptionCoefficientFromNu( + rho / rho_unit, temp / temp_unit, nu_mid * time_unit) / + length_unit, + reference_host.PlanckGroupAbsorptionCoefficientFromNu( + rho, temp, nu_mid)) < EPS_TEST); + REQUIRE(FractionalDifference( + funny_host.RosselandGroupAbsorptionCoefficientFromNu( + rho / rho_unit, temp / temp_unit, nu_mid * time_unit) / + length_unit, + reference_host.RosselandGroupAbsorptionCoefficientFromNu( + rho, temp, nu_mid)) < EPS_TEST); + REQUIRE(FractionalDifference( + funny_host.AbsorptionCoefficientFromNu( + rho / rho_unit, temp / temp_unit, nu_mid * time_unit) / + length_unit, + reference_host.AbsorptionCoefficientFromNu(rho, temp, nu_mid)) < + EPS_TEST); + } + + funny_host.Finalize(); + reference_host.Finalize(); + kappa_planck.finalize(); + kappa_rosseland.finalize(); +} + +TEST_CASE("Photon multigroup non-CGS wrapper works for monochromatic-built " + "multigroup opacities", + "[MultigroupPhotons]") { + constexpr Real rho = 4.e-6; + constexpr Real temp = 7.e5; + constexpr Real kappa = 9.e-3; + constexpr int NRho = 3; + constexpr int NT = 4; + constexpr int ngroups = 4; + constexpr int nnu_per_group = 96; + constexpr Real time_unit = 11.; + constexpr Real mass_unit = 13.; + constexpr Real length_unit = 17.; + constexpr Real temp_unit = 19.; + constexpr Real rho_unit = + mass_unit / (length_unit * length_unit * length_unit); + const Real lRhoMin = std::log10(0.25 * rho); + const Real lRhoMax = std::log10(4. * rho); + const Real lTMin = std::log10(0.5 * temp); + const Real lTMax = std::log10(2. * temp); + const std::array group_bounds = {1.e12, 3.e12, 1.e13, + 1.e14, 3.e15}; + + photons::Gray opac_host(kappa); + photons::MeanOpacityBase reference_host( + opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, group_bounds, + ngroups, nnu_per_group); + photons::MeanOpacityBase multigroup_base( + opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, group_bounds, + ngroups, nnu_per_group); + auto funny_host = + photons::MeanNonCGSUnits( + std::move(multigroup_base), time_unit, mass_unit, length_unit, + temp_unit); + auto funny = funny_host.GetOnDevice(); + + REQUIRE(funny.ngroups() == ngroups); + REQUIRE(funny.HasGroupBounds()); + + int n_wrong_h = 0; +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View n_wrong_d("wrong"); +#else + PortableMDArray n_wrong_d(&n_wrong_h, 1); +#endif + + portableFor( + "compare multigroup non-cgs wrapper built from monochromatic opacity", 0, + ngroups, PORTABLE_LAMBDA(const int &group) { + const Real rho_funny = rho / rho_unit; + const Real temp_funny = temp / temp_unit; + const Real alpha_planck_funny = funny.PlanckGroupAbsorptionCoefficient( + rho_funny, temp_funny, group); + const Real alpha_rosseland_funny = + funny.RosselandGroupAbsorptionCoefficient(rho_funny, temp_funny, + group); + const Real alpha_default_funny = + funny.AbsorptionCoefficient(rho_funny, temp_funny, group); + const Real nu_mid = + std::sqrt(group_bounds[group] * group_bounds[group + 1]); + const Real alpha_planck_from_nu_funny = + funny.PlanckGroupAbsorptionCoefficientFromNu(rho_funny, temp_funny, + nu_mid * time_unit); + const Real alpha_rosseland_from_nu_funny = + funny.RosselandGroupAbsorptionCoefficientFromNu( + rho_funny, temp_funny, nu_mid * time_unit); + const Real alpha_default_from_nu_funny = + funny.AbsorptionCoefficientFromNu(rho_funny, temp_funny, + nu_mid * time_unit); + const Real alpha_planck = + reference_host.PlanckGroupAbsorptionCoefficient(rho, temp, group); + const Real alpha_rosseland = + reference_host.RosselandGroupAbsorptionCoefficient(rho, temp, + group); + + if (FractionalDifference(alpha_planck_funny / length_unit, + alpha_planck) > EPS_TEST) { + n_wrong_d() += 1; + } + if (FractionalDifference(alpha_rosseland_funny / length_unit, + alpha_rosseland) > EPS_TEST) { + n_wrong_d() += 1; + } + if (FractionalDifference(alpha_default_funny / length_unit, + alpha_rosseland) > EPS_TEST) { + n_wrong_d() += 1; + } + if (funny.GroupOfNu(nu_mid * time_unit) != group) { + n_wrong_d() += 1; + } + if (FractionalDifference(alpha_planck_from_nu_funny / length_unit, + alpha_planck) > EPS_TEST) { + n_wrong_d() += 1; + } + if (FractionalDifference(alpha_rosseland_from_nu_funny / length_unit, + alpha_rosseland) > EPS_TEST) { + n_wrong_d() += 1; + } + if (FractionalDifference(alpha_default_from_nu_funny / length_unit, + alpha_rosseland) > EPS_TEST) { + n_wrong_d() += 1; + } + }); + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(n_wrong_h, n_wrong_d); +#endif + REQUIRE(n_wrong_h == 0); + + funny_host.Finalize(); + reference_host.Finalize(); +} + +TEST_CASE("Photon multigroup with extreme bounds 0 to infinity recovers gray " + "solution", + "[MultigroupPhotons]") { + constexpr Real rho = 3.e-5; + constexpr Real temp = 2.e6; + constexpr Real kappa = 5.e-2; + constexpr int NRho = 4; + constexpr int NT = 5; + constexpr int ngroups = 1; + constexpr int nnu_per_group = 128; + const Real lRhoMin = std::log10(0.1 * rho); + const Real lRhoMax = std::log10(10. * rho); + const Real lTMin = std::log10(0.25 * temp); + const Real lTMax = std::log10(4. * temp); + const std::array full_spectrum_bounds = { + 0., std::numeric_limits::infinity()}; + + photons::Gray opac_host(kappa); + photons::MeanOpacityBase multigroup_host( + opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, full_spectrum_bounds, + ngroups, nnu_per_group); + + REQUIRE(multigroup_host.ngroups() == ngroups); + REQUIRE(multigroup_host.HasGroupBounds()); + + // Test at multiple densities and temperatures + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real test_rho = + std::pow(10., lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho); + for (int iT = 0; iT < NT; ++iT) { + const Real test_temp = + std::pow(10., lTMin + (lTMax - lTMin) / (NT - 1) * iT); + + const Real alpha_planck = + multigroup_host.PlanckGroupAbsorptionCoefficient(test_rho, test_temp, + 0); + const Real alpha_rosseland = + multigroup_host.RosselandGroupAbsorptionCoefficient(test_rho, + test_temp, 0); + const Real alpha_gray = test_rho * kappa; + + // Verify no NaNs or infinities + REQUIRE(std::isfinite(alpha_planck)); + REQUIRE(std::isfinite(alpha_rosseland)); + + // For a gray opacity, both Planck and Rosseland means should equal the + // gray opacity + REQUIRE(FractionalDifference(alpha_planck, alpha_gray) < EPS_TEST); + REQUIRE(FractionalDifference(alpha_rosseland, alpha_gray) < EPS_TEST); + } + } + + multigroup_host.Finalize(); +} + +TEST_CASE("Photon multigroup with tail groups 0 to nu_mid and nu_mid to " + "infinity are numerically stable", + "[MultigroupPhotons]") { + constexpr Real rho = 7.e-4; + constexpr Real temp = 5.e5; + constexpr Real kappa = 1.2e-2; + constexpr int NRho = 3; + constexpr int NT = 4; + constexpr int ngroups = 2; + constexpr int nnu_per_group = 96; + constexpr Real nu_mid = 1.e14; + const Real lRhoMin = std::log10(0.2 * rho); + const Real lRhoMax = std::log10(5. * rho); + const Real lTMin = std::log10(0.5 * temp); + const Real lTMax = std::log10(2. * temp); + const std::array tail_bounds = { + 0., nu_mid, std::numeric_limits::infinity()}; + + photons::Gray opac_host(kappa); + photons::MeanOpacityBase multigroup_host( + opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, tail_bounds, ngroups, + nnu_per_group); + + REQUIRE(multigroup_host.ngroups() == ngroups); + + // Test both groups at multiple state points + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real test_rho = + std::pow(10., lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho); + for (int iT = 0; iT < NT; ++iT) { + const Real test_temp = + std::pow(10., lTMin + (lTMax - lTMin) / (NT - 1) * iT); + + for (int group = 0; group < ngroups; ++group) { + const Real alpha_planck = + multigroup_host.PlanckGroupAbsorptionCoefficient(test_rho, + test_temp, group); + const Real alpha_rosseland = + multigroup_host.RosselandGroupAbsorptionCoefficient( + test_rho, test_temp, group); + const Real alpha_default = + multigroup_host.AbsorptionCoefficient(test_rho, test_temp, group); + + // Verify no NaNs or infinities + REQUIRE(std::isfinite(alpha_planck)); + REQUIRE(std::isfinite(alpha_rosseland)); + REQUIRE(std::isfinite(alpha_default)); + + // Verify positivity + REQUIRE(alpha_planck >= 0.); + REQUIRE(alpha_rosseland >= 0.); + REQUIRE(alpha_default >= 0.); + + // For gray opacity, both groups should give the gray value + const Real alpha_gray = test_rho * kappa; + REQUIRE(FractionalDifference(alpha_planck, alpha_gray) < EPS_TEST); + REQUIRE(FractionalDifference(alpha_rosseland, alpha_gray) < EPS_TEST); + } + } + } + + multigroup_host.Finalize(); +} + +TEST_CASE("Photon multigroup with asymmetric tail groups is numerically stable", + "[MultigroupPhotons]") { + constexpr Real rho = 2.e-3; + constexpr Real temp = 1.e6; + constexpr Real kappa = 3.5e-2; + constexpr int NRho = 3; + constexpr int NT = 3; + constexpr int ngroups = 4; + constexpr int nnu_per_group = 80; + constexpr Real nu_low = 5.e13; + constexpr Real nu_high = 2.e15; + const Real lRhoMin = std::log10(0.25 * rho); + const Real lRhoMax = std::log10(4. * rho); + const Real lTMin = std::log10(0.5 * temp); + const Real lTMax = std::log10(2. * temp); + // Groups: [0, nu_low], [nu_low, nu_mid], [nu_mid, nu_high], [nu_high, inf] + const std::array asymmetric_bounds = { + 0., nu_low, std::sqrt(nu_low * nu_high), nu_high, + std::numeric_limits::infinity()}; + + photons::Gray opac_host(kappa); + photons::MeanOpacityBase multigroup_host( + opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, asymmetric_bounds, + ngroups, nnu_per_group); + + REQUIRE(multigroup_host.ngroups() == ngroups); + + const Real test_rho = std::pow(10., 0.5 * (lRhoMin + lRhoMax)); + const Real test_temp = std::pow(10., 0.5 * (lTMin + lTMax)); + const Real alpha_gray = test_rho * kappa; + + for (int group = 0; group < ngroups; ++group) { + const Real alpha_planck = multigroup_host.PlanckGroupAbsorptionCoefficient( + test_rho, test_temp, group); + const Real alpha_rosseland = + multigroup_host.RosselandGroupAbsorptionCoefficient(test_rho, test_temp, + group); + + // Verify no NaNs or infinities + REQUIRE(std::isfinite(alpha_planck)); + REQUIRE(std::isfinite(alpha_rosseland)); + + // Verify positivity + REQUIRE(alpha_planck >= 0.); + REQUIRE(alpha_rosseland >= 0.); + + // For gray opacity, all groups should give the gray value + REQUIRE(FractionalDifference(alpha_planck, alpha_gray) < EPS_TEST); + REQUIRE(FractionalDifference(alpha_rosseland, alpha_gray) < EPS_TEST); + } + + multigroup_host.Finalize(); +} + +TEST_CASE("Photon multigroup with single group 0 to nuMax recovers gray " + "solution", + "[MultigroupPhotons]") { + constexpr Real rho = 4.e-4; + constexpr Real temp = 8.e5; + constexpr Real kappa = 7.e-3; + constexpr int NRho = 3; + constexpr int NT = 4; + constexpr int ngroups = 1; + constexpr int nnu_per_group = 96; + constexpr Real nu_max = 5.e15; + const Real lRhoMin = std::log10(0.2 * rho); + const Real lRhoMax = std::log10(5. * rho); + const Real lTMin = std::log10(0.25 * temp); + const Real lTMax = std::log10(4. * temp); + const std::array low_tail_bounds = {0., nu_max}; + + photons::Gray opac_host(kappa); + photons::MeanOpacityBase multigroup_host( + opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, low_tail_bounds, + ngroups, nnu_per_group); + + REQUIRE(multigroup_host.ngroups() == ngroups); + REQUIRE(multigroup_host.HasGroupBounds()); + + // Test at multiple state points + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real test_rho = + std::pow(10., lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho); + for (int iT = 0; iT < NT; ++iT) { + const Real test_temp = + std::pow(10., lTMin + (lTMax - lTMin) / (NT - 1) * iT); + + const Real alpha_planck = + multigroup_host.PlanckGroupAbsorptionCoefficient(test_rho, test_temp, + 0); + const Real alpha_rosseland = + multigroup_host.RosselandGroupAbsorptionCoefficient(test_rho, + test_temp, 0); + const Real alpha_gray = test_rho * kappa; + + // Verify no NaNs or infinities + REQUIRE(std::isfinite(alpha_planck)); + REQUIRE(std::isfinite(alpha_rosseland)); + + // For gray opacity integrated from [0, nuMax), should equal gray value + // (nuMax is large enough to capture essentially all Planck weight) + REQUIRE(FractionalDifference(alpha_planck, alpha_gray) < EPS_TEST); + REQUIRE(FractionalDifference(alpha_rosseland, alpha_gray) < EPS_TEST); + } + } + + // Verify GroupOfNu works correctly + REQUIRE(multigroup_host.GroupOfNu(0.) == 0); + REQUIRE(multigroup_host.GroupOfNu(0.5 * nu_max) == 0); + REQUIRE(multigroup_host.GroupOfNu(nu_max) == 0); + + multigroup_host.Finalize(); +} + +TEST_CASE("Photon multigroup with single group nuMin to infinity recovers " + "gray solution", + "[MultigroupPhotons]") { + constexpr Real rho = 6.e-3; + constexpr Real temp = 3.e5; + constexpr Real kappa = 2.5e-2; + constexpr int NRho = 3; + constexpr int NT = 4; + constexpr int ngroups = 1; + constexpr int nnu_per_group = 128; + constexpr Real nu_min = 1.e12; + const Real lRhoMin = std::log10(0.25 * rho); + const Real lRhoMax = std::log10(4. * rho); + const Real lTMin = std::log10(0.5 * temp); + const Real lTMax = std::log10(2. * temp); + const std::array high_tail_bounds = { + nu_min, std::numeric_limits::infinity()}; + + photons::Gray opac_host(kappa); + photons::MeanOpacityBase multigroup_host( + opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, high_tail_bounds, + ngroups, nnu_per_group); + + REQUIRE(multigroup_host.ngroups() == ngroups); + REQUIRE(multigroup_host.HasGroupBounds()); + + // Test at multiple state points + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real test_rho = + std::pow(10., lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho); + for (int iT = 0; iT < NT; ++iT) { + const Real test_temp = + std::pow(10., lTMin + (lTMax - lTMin) / (NT - 1) * iT); + + const Real alpha_planck = + multigroup_host.PlanckGroupAbsorptionCoefficient(test_rho, test_temp, + 0); + const Real alpha_rosseland = + multigroup_host.RosselandGroupAbsorptionCoefficient(test_rho, + test_temp, 0); + const Real alpha_gray = test_rho * kappa; + + // Verify no NaNs or infinities + REQUIRE(std::isfinite(alpha_planck)); + REQUIRE(std::isfinite(alpha_rosseland)); + + // For gray opacity integrated from [nuMin, infinity), should equal gray + // value (nuMin is small enough to capture essentially all Planck weight) + REQUIRE(FractionalDifference(alpha_planck, alpha_gray) < EPS_TEST); + REQUIRE(FractionalDifference(alpha_rosseland, alpha_gray) < EPS_TEST); + } + } + + // Verify GroupOfNu works correctly + REQUIRE(multigroup_host.GroupOfNu(nu_min) == 0); + REQUIRE(multigroup_host.GroupOfNu(2. * nu_min) == 0); + REQUIRE(multigroup_host.GroupOfNu(1.e20) == 0); + REQUIRE(multigroup_host.GroupOfNu(std::numeric_limits::infinity()) == + 0); + + multigroup_host.Finalize(); +} + +TEST_CASE("Photon multigroup GroupOfNu handles extreme bounds correctly", + "[MultigroupPhotons]") { + constexpr Real rho = 1.e-2; + constexpr Real temp = 3.e5; + constexpr Real kappa = 2.e-2; + constexpr int NRho = 2; + constexpr int NT = 2; + constexpr int ngroups = 3; + constexpr int nnu_per_group = 64; + constexpr Real nu_low = 1.e13; + constexpr Real nu_high = 1.e15; + const Real lRhoMin = std::log10(0.5 * rho); + const Real lRhoMax = std::log10(2. * rho); + const Real lTMin = std::log10(0.5 * temp); + const Real lTMax = std::log10(2. * temp); + const std::array extreme_bounds = { + 0., nu_low, nu_high, std::numeric_limits::infinity()}; + + photons::Gray opac_host(kappa); + photons::MeanOpacityBase multigroup_host( + opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, extreme_bounds, + ngroups, nnu_per_group); + + REQUIRE(multigroup_host.ngroups() == ngroups); + REQUIRE(multigroup_host.HasGroupBounds()); + + // Test GroupOfNu at various frequencies + REQUIRE(multigroup_host.GroupOfNu(0.) == 0); + REQUIRE(multigroup_host.GroupOfNu(1.e-10) == 0); + REQUIRE(multigroup_host.GroupOfNu(0.5 * nu_low) == 0); + REQUIRE(multigroup_host.GroupOfNu(nu_low) == 1); + REQUIRE(multigroup_host.GroupOfNu(std::sqrt(nu_low * nu_high)) == 1); + REQUIRE(multigroup_host.GroupOfNu(nu_high) == 2); + REQUIRE(multigroup_host.GroupOfNu(2. * nu_high) == 2); + REQUIRE(multigroup_host.GroupOfNu(1.e20) == 2); + REQUIRE(multigroup_host.GroupOfNu(std::numeric_limits::infinity()) == + ngroups - 1); + + // Verify that AbsorptionCoefficientFromNu works at extreme frequencies + const Real alpha_low = + multigroup_host.AbsorptionCoefficientFromNu(rho, temp, 0.); + const Real alpha_very_low = + multigroup_host.AbsorptionCoefficientFromNu(rho, temp, 1.e-10); + const Real alpha_high = + multigroup_host.AbsorptionCoefficientFromNu(rho, temp, 1.e20); + const Real alpha_inf = multigroup_host.AbsorptionCoefficientFromNu( + rho, temp, std::numeric_limits::infinity()); + + // All should be finite and positive + REQUIRE(std::isfinite(alpha_low)); + REQUIRE(std::isfinite(alpha_very_low)); + REQUIRE(std::isfinite(alpha_high)); + REQUIRE(std::isfinite(alpha_inf)); + REQUIRE(alpha_low > 0.); + REQUIRE(alpha_very_low > 0.); + REQUIRE(alpha_high > 0.); + REQUIRE(alpha_inf > 0.); + + // For gray opacity, all should equal the gray value + const Real alpha_gray = rho * kappa; + REQUIRE(FractionalDifference(alpha_low, alpha_gray) < EPS_TEST); + REQUIRE(FractionalDifference(alpha_very_low, alpha_gray) < EPS_TEST); + REQUIRE(FractionalDifference(alpha_high, alpha_gray) < EPS_TEST); + REQUIRE(FractionalDifference(alpha_inf, alpha_gray) < EPS_TEST); + + multigroup_host.Finalize(); +} + +TEST_CASE("Photon multigroup gray scattering opacities are exact", + "[MultigroupPhotons][MultigroupScattering]") { + constexpr Real rho = 1.e-3; + constexpr Real temp = 1.e6; + constexpr Real sigma = 0.665e-24; // Thomson cross section + constexpr Real apm = 2.0e-24; // avg particle mass ~ 2 protons + constexpr int NRho = 3; + constexpr int NT = 4; + constexpr int ngroups = 4; + constexpr int nnu_per_group = 96; + const Real lRhoMin = std::log10(0.1 * rho); + const Real lRhoMax = std::log10(10. * rho); + const Real lTMin = std::log10(0.25 * temp); + const Real lTMax = std::log10(4. * temp); + const std::array group_bounds = {1.e12, 3.e12, 1.e13, + 1.e14, 3.e15}; + + photons::GraySOpacity s_opac_host(sigma, apm); + photons::MeanSOpacityBase multigroup_host( + s_opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, group_bounds, + ngroups, nnu_per_group); + + REQUIRE(multigroup_host.ngroups() == ngroups); + + int n_wrong = 0; + for (int group = 0; group < ngroups; ++group) { + const Real sigma_planck = + multigroup_host.PlanckGroupScatteringCoefficient(rho, temp, group); + const Real sigma_rosseland = + multigroup_host.RosselandGroupScatteringCoefficient(rho, temp, group); + const Real sigma_default = + multigroup_host.ScatteringCoefficient(rho, temp, group); + const Real sigma_expected = (rho / apm) * sigma; + if (FractionalDifference(sigma_planck, sigma_expected) > EPS_TEST) { + n_wrong += 1; + } + if (FractionalDifference(sigma_rosseland, sigma_expected) > EPS_TEST) { + n_wrong += 1; + } + if (FractionalDifference(sigma_default, sigma_expected) > EPS_TEST) { + n_wrong += 1; + } + } + REQUIRE(n_wrong == 0); + + multigroup_host.Finalize(); +} + +TEST_CASE("Photon multigroup scattering with extreme bounds 0 to infinity " + "works correctly", + "[MultigroupPhotons][MultigroupScattering]") { + constexpr Real rho = 2.e-4; + constexpr Real temp = 5.e5; + constexpr Real sigma = 0.665e-24; + constexpr Real apm = 2.0e-24; + constexpr int NRho = 3; + constexpr int NT = 3; + constexpr int ngroups = 1; + constexpr int nnu_per_group = 128; + const Real lRhoMin = std::log10(0.2 * rho); + const Real lRhoMax = std::log10(5. * rho); + const Real lTMin = std::log10(0.5 * temp); + const Real lTMax = std::log10(2. * temp); + const std::array full_spectrum_bounds = { + 0., std::numeric_limits::infinity()}; + + photons::GraySOpacity s_opac_host(sigma, apm); + photons::MeanSOpacityBase multigroup_host( + s_opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, + full_spectrum_bounds, ngroups, nnu_per_group); + + REQUIRE(multigroup_host.ngroups() == ngroups); + REQUIRE(multigroup_host.HasGroupBounds()); + + const Real sigma_expected = (rho / apm) * sigma; + const Real sigma_planck = + multigroup_host.PlanckGroupScatteringCoefficient(rho, temp, 0); + const Real sigma_rosseland = + multigroup_host.RosselandGroupScatteringCoefficient(rho, temp, 0); + + REQUIRE(std::isfinite(sigma_planck)); + REQUIRE(std::isfinite(sigma_rosseland)); + REQUIRE(FractionalDifference(sigma_planck, sigma_expected) < EPS_TEST); + REQUIRE(FractionalDifference(sigma_rosseland, sigma_expected) < EPS_TEST); + + multigroup_host.Finalize(); +} + +TEST_CASE("Photon multigroup scattering GroupOfNu handles extreme bounds", + "[MultigroupPhotons][MultigroupScattering]") { + constexpr Real rho = 5.e-3; + constexpr Real temp = 2.e6; + constexpr Real sigma = 0.665e-24; + constexpr Real apm = 2.0e-24; + constexpr int NRho = 2; + constexpr int NT = 2; + constexpr int ngroups = 3; + constexpr int nnu_per_group = 64; + constexpr Real nu_low = 1.e13; + constexpr Real nu_high = 1.e15; + const Real lRhoMin = std::log10(0.5 * rho); + const Real lRhoMax = std::log10(2. * rho); + const Real lTMin = std::log10(0.5 * temp); + const Real lTMax = std::log10(2. * temp); + const std::array extreme_bounds = { + 0., nu_low, nu_high, std::numeric_limits::infinity()}; + + photons::GraySOpacity s_opac_host(sigma, apm); + photons::MeanSOpacityBase multigroup_host( + s_opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, extreme_bounds, + ngroups, nnu_per_group); + + REQUIRE(multigroup_host.ngroups() == ngroups); + REQUIRE(multigroup_host.HasGroupBounds()); + + // Test GroupOfNu at various frequencies + REQUIRE(multigroup_host.GroupOfNu(0.) == 0); + REQUIRE(multigroup_host.GroupOfNu(0.5 * nu_low) == 0); + REQUIRE(multigroup_host.GroupOfNu(nu_low) == 1); + REQUIRE(multigroup_host.GroupOfNu(std::sqrt(nu_low * nu_high)) == 1); + REQUIRE(multigroup_host.GroupOfNu(nu_high) == 2); + REQUIRE(multigroup_host.GroupOfNu(2. * nu_high) == 2); + REQUIRE(multigroup_host.GroupOfNu(std::numeric_limits::infinity()) == + ngroups - 1); + + // Verify ScatteringCoefficientFromNu works at extremes + const Real sigma_expected = (rho / apm) * sigma; + const Real sigma_low = + multigroup_host.ScatteringCoefficientFromNu(rho, temp, 0.); + const Real sigma_high = + multigroup_host.ScatteringCoefficientFromNu(rho, temp, 1.e20); + + REQUIRE(std::isfinite(sigma_low)); + REQUIRE(std::isfinite(sigma_high)); + REQUIRE(FractionalDifference(sigma_low, sigma_expected) < EPS_TEST); + REQUIRE(FractionalDifference(sigma_high, sigma_expected) < EPS_TEST); + + multigroup_host.Finalize(); +} + +#ifdef SPINER_USE_HDF +TEST_CASE("Photon multigroup scattering tables can round-trip through SP5 HDF", + "[MultigroupPhotons][MultigroupScattering]") { + using DataBox = Spiner::DataBox; + + constexpr int NRho = 2; + constexpr int NT = 2; + constexpr int ngroups = 3; + constexpr Real sigma0 = 1.5e-24; + const Real lRhoMin = -4.; + const Real lRhoMax = 2.; + const Real lTMin = 2.; + const Real lTMax = 8.; + constexpr Real nu_min = 2.e11; + constexpr Real nu_max = 4.e11; + const std::array group_bounds = { + 0., nu_min, nu_max, std::numeric_limits::infinity()}; + + DataBox sigma_planck(NRho, NT, ngroups); + sigma_planck.setRange(1, lTMin, lTMax, NT); + sigma_planck.setRange(2, lRhoMin, lRhoMax, NRho); + DataBox sigma_rosseland; + sigma_rosseland.copyMetadata(sigma_planck); + + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real rho = std::pow(10., sigma_planck.range(2).x(iRho)); + for (int iT = 0; iT < NT; ++iT) { + const Real temp = std::pow(10., sigma_planck.range(1).x(iT)); + for (int group = 0; group < ngroups; ++group) { + const Real group_factor = group + 1.; + sigma_planck(iRho, iT, group) = sigma0 * rho * group_factor; + sigma_rosseland(iRho, iT, group) = sigma0 * rho / group_factor; + } + } + } + + photons::MeanSOpacityBase saved(sigma_planck, sigma_rosseland, + group_bounds); + const char *filename = "multigroup-photon-scattering-table.sp5"; + saved.Save(filename); + photons::MeanSOpacityBase loaded(filename); + + REQUIRE(loaded.HasGroupBounds()); + REQUIRE(loaded.ngroups() == ngroups); + REQUIRE(loaded.GroupOfNu(0.) == 0); + REQUIRE(loaded.GroupOfNu(nu_min) == 1); + REQUIRE(loaded.GroupOfNu(nu_max) == ngroups - 1); + + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real rho_test = + std::pow(10., lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho); + for (int iT = 0; iT < NT; ++iT) { + const Real temp_test = + std::pow(10., lTMin + (lTMax - lTMin) / (NT - 1) * iT); + for (int group = 0; group < ngroups; ++group) { + const Real sigma_planck_expected = + rho_test * sigma_planck(iRho, iT, group) / rho_test; + const Real sigma_rosseland_expected = + rho_test * sigma_rosseland(iRho, iT, group) / rho_test; + REQUIRE(FractionalDifference(loaded.PlanckGroupScatteringCoefficient( + rho_test, temp_test, group), + sigma_planck_expected * rho_test) < + EPS_TEST); + REQUIRE(FractionalDifference(loaded.RosselandGroupScatteringCoefficient( + rho_test, temp_test, group), + sigma_rosseland_expected * rho_test) < + EPS_TEST); + } + } + } + + loaded.Finalize(); + saved.Finalize(); + sigma_planck.finalize(); + sigma_rosseland.finalize(); +} +#endif