Towards Multigroup API

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<Real>::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
 

plan_histories/69.md

@@ -1,3 +1,5 @@
+<!-- This file was made in part with generative AI. -->
+
 ## Goal
 
 Correct the photon and neutrino `dB_\nu/dT` implementations used for Rosseland

plan_histories/70.md

@@ -1,3 +1,5 @@
+<!-- This file was made in part with generative AI. -->
+
 ## Goal
 
 Extend photon `PowerLaw` opacities to support an optional frequency power law

plan_histories/71.md

@@ -0,0 +1,38 @@
+<!-- This file was made in part with generative AI. -->
+
+## 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).

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 <fast-math/logs.hpp>
-#include <spiner/databox.hpp>
 #include <ports-of-call/variant.hpp>
+#include <spiner/databox.hpp>
 
 #include <singularity-opac/chebyshev/chebyshev.hpp>
 
@@ -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);
   }
 

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_
 

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_
 

singularity-opac/base/robust_utils.hpp

@@ -11,47 +11,46 @@
 // 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_
 
 #include <limits>
 #include <ports-of-call/portability.hpp>
 
 namespace singularity_opac {
-  namespace robust {
-
-  template <typename T = Real>
-  PORTABLE_FORCEINLINE_FUNCTION constexpr auto SMALL() {
-    return 10 * std::numeric_limits<T>::min();
-  }
-
-  template <typename T = Real>
-  PORTABLE_FORCEINLINE_FUNCTION constexpr auto EPS() {
-    return 10 * std::numeric_limits<T>::epsilon();
-  }
-
-  template <typename T>
-  PORTABLE_FORCEINLINE_FUNCTION auto make_positive(const T val) {
-    return std::max(val, EPS<T>());
-  }
-
-  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 <typename T>
-  PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val) {
-    return (T(0) <= val) - (val < T(0));
-  }
-
-  template <typename A, typename B>
-  PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) {
-    return a / (b + sgn(b) * SMALL<B>());
-  }
-
-  } // namespace robust
+namespace robust {
+
+template <typename T = Real>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto SMALL() {
+  return 10 * std::numeric_limits<T>::min();
+}
+
+template <typename T = Real>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto EPS() {
+  return 10 * std::numeric_limits<T>::epsilon();
+}
+
+template <typename T>
+PORTABLE_FORCEINLINE_FUNCTION auto make_positive(const T val) {
+  return std::max(val, EPS<T>());
+}
+
+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 <typename T>
+PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val) {
+  return (T(0) <= val) - (val < T(0));
+}
+
+template <typename A, typename B>
+PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) {
+  return a / (b + sgn(b) * SMALL<B>());
+}
+
+} // namespace robust
 } // namespace singularity_opac
 
 #endif // SINGULARITY_OPAC_BASE_ROBUST_UTILS_HPP_

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_

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_
 

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