Add support for OpenMC MGXS delayed neutron data.

doc/source/userguide/discrete_ordinates_problems.rst

@@ -208,10 +208,18 @@ The most important options for everyday use are:
 * ``ags_convergence_check``
 * ``field_function_prefix_option``
 * ``field_function_prefix``
+* ``use_precursors``
 
 There are also restart, precursor, and message-size options for more
 specialized workflows.
 
+``use_precursors`` controls precursor storage and transient delayed-neutron
+evolution. Steady-state and k-eigenvalue solves include delayed neutron
+production whenever the active cross-section data contains it, independent of
+this option. In transient mode, ``use_precursors=False`` selects a prompt-only
+transient, while ``use_precursors=True`` evolves precursor concentrations and
+includes delayed production through those precursor equations.
+
 Restart-related options include:
 
 * ``restart_writes_enabled``: enable restart dump writes.

doc/source/userguide/materials_xs.rst

@@ -267,11 +267,88 @@ fission production. This is useful for advanced data preparation, but most user
 inputs are clearer when written in terms of ``SIGMA_F``, ``NU`` or
 ``NU_PROMPT``/``NU_DELAYED``, and the corresponding spectra.
 
+Delayed-Neutron Data and Solver Behavior
+----------------------------------------
+
+Delayed-neutron data has two related but distinct roles in OpenSn:
+
+* it contributes delayed fission production to steady-state and k-eigenvalue
+  fission sources, and
+* it can define precursor concentrations for time-dependent delayed-neutron
+  evolution.
+
+These roles are intentionally not controlled by the same switch. If a
+steady-state or k-eigenvalue problem uses cross sections with delayed-neutron
+data, OpenSn includes delayed neutron production in the steady fission source
+and in fission-production normalization. The steady eigenvalue therefore
+corresponds to total fission production, prompt plus delayed. This is true even
+when the problem option ``use_precursors`` is ``False``.
+
+For transient problems, ``use_precursors`` controls whether delayed-neutron
+precursor state is stored and evolved:
+
+* ``use_precursors=True`` uses prompt production and precursor-mediated delayed
+  production.
+* ``use_precursors=False`` is a prompt-only transient, even if the cross-section
+  file contains delayed-neutron data.
+
+This distinction is useful when a steady k-eigenvalue solve is used only to
+prepare an initial condition. The steady solve still computes the physical
+total-production eigenvalue, while a later transient can deliberately choose
+whether to evolve delayed precursors.
+
+.. important::
+
+   If delayed-neutron production is active and any fissionable material in the
+   active cross-section map contains precursor data, all fissionable materials
+   in that map must contain compatible precursor data. Non-fissionable materials
+   may have zero precursors. This avoids solving a mixed prompt-only and
+   prompt-plus-delayed fission model by accident.
+
+Native OpenSn delayed data is interpreted as follows:
+
+* ``NU_PROMPT`` and ``CHI_PROMPT`` define the prompt fission production matrix.
+* ``NU_DELAYED`` defines the delayed fission production by source group.
+* ``PRECURSOR_FRACTIONAL_YIELDS`` partitions delayed production among precursor
+  families. The yields are normalized when read.
+* ``CHI_DELAYED`` defines each precursor family's delayed emission spectrum by
+  destination group.
+* ``PRECURSOR_DECAY_CONSTANTS`` are used for precursor concentrations and
+  transient evolution; they do not change the steady-state k-eigenvalue.
+
+For steady-state source assembly and k-eigenvalue normalization, OpenSn adds the
+delayed term
+
+.. math::
+
+   \sum_j \chi_{d,j,g}\,\gamma_j\,\nu_{d,g'}\Sigma_{f,g'}\phi_{g'}
+
+to the prompt production matrix contribution. Here ``j`` is the precursor family
+index, ``g`` is the destination group, and ``g'`` is the source group.
+
+When the native ``PRODUCTION_MATRIX`` form is used, delayed-neutron precursor
+data is currently not allowed in the same material. This avoids ambiguity over
+whether the matrix already contains delayed production.
+
+Practical guidance:
+
+* For production steady-state or k-eigenvalue calculations, leave delayed data in
+  the XS when it exists. Do not disable ``use_precursors`` to try to remove
+  delayed neutron production from a steady eigenvalue solve.
+* For production transients with delayed-neutron physics, use
+  ``options={"use_precursors": True, "save_angular_flux": True}``.
+* Use ``use_precursors=False`` with delayed XS data only when a prompt-only
+  transient is intentional.
+* If a k-eigenvalue changes when delayed data is removed from the XS, that is
+  expected. If it changes only because ``use_precursors`` was toggled in a
+  steady-state calculation, that is not expected.
+
 Loading OpenMC MGXS Files
 =========================
 
 Use :py:meth:`pyopensn.xs.MultiGroupXS.LoadFromOpenMC` to load cross sections
-from an OpenMC multi-group HDF5 library.
+from an OpenMC multi-group HDF5 library. Delayed-neutron data is imported when
+it is present in the OpenMC file.
 
 .. code-block:: python
 
@@ -309,6 +386,13 @@ Example:
    absorption dataset directly. This makes the imported data less sensitive to
    small inconsistencies or statistical noise.
 
+   When delayed-neutron data is present, OpenSn also imports precursor-group
+   information, prompt and delayed fission production, precursor decay
+   constants, and delayed emission spectra. If ``chi-prompt`` is not present,
+   OpenSn falls back to the steady fission spectrum ``chi`` as the prompt
+   spectrum. Once imported, the same steady-state and transient delayed-neutron
+   behavior described above applies.
+
 What OpenMC Import Expects
 --------------------------
 
@@ -326,6 +410,14 @@ The imported object may include:
 * fission cross section
 * fission production data
 * fission spectrum
+* delayed-neutron data, including:
+
+  - ``prompt-nu-fission``
+  - ``delayed-nu-fission``
+  - ``decay-rate``
+  - ``chi-delayed``
+  - optional ``chi-prompt``
+
 * any requested named custom one-dimensional XS in ``extra_xs_names``
 
 If the requested dataset or temperature is not present, the load fails with an

doc/source/userguide/solvers.rst

@@ -178,10 +178,18 @@ Problem options available in Python include:
 * ``field_function_prefix_option``
 * ``field_function_prefix``
 
-``use_precursors`` controls whether delayed-neutron precursor treatment is kept
-active for the problem. The default is ``True``. This should usually stay
-enabled for transient and k-eigen workflows unless you explicitly want a
-prompt-only model.
+``use_precursors`` controls delayed-neutron precursor storage and transient
+precursor evolution. It does not remove delayed neutron production from
+steady-state or k-eigenvalue fission sources. If the active cross sections
+contain delayed-neutron data, steady-state and k-eigenvalue solves use total
+fission production, prompt plus delayed, regardless of ``use_precursors``.
+
+For transient solves, ``use_precursors=True`` evolves precursor concentrations
+and includes delayed production through the transient precursor equations.
+``use_precursors=False`` is a prompt-only transient, even if the cross sections
+contain delayed-neutron data. Keep the default ``True`` for production
+delayed-neutron transients; set it to ``False`` only when a prompt-only transient
+is intentional.
 
 ``read_restart_path`` reads a full restart for continuing a compatible solve.
 ``read_initial_condition_path`` reads restart data as an initial condition; this
@@ -205,9 +213,11 @@ new families start at zero and removed families are discarded. If a cell passes
 through a material with zero precursors, its precursor history is dropped and
 any later reintroduced precursor families restart from zero.
 
-If any fissionable material in the active map contains precursor data and
-``use_precursors=True``, then all fissionable materials in that map must
-contain precursor data. Non-fissionable materials may have zero precursors.
+If delayed-neutron production is active and any fissionable material in the
+active map contains precursor data, then all fissionable materials in that map
+must contain precursor data. Non-fissionable materials may have zero precursors.
+In steady-state mode, delayed-neutron production is active whenever delayed data
+exists. In transient mode, it is active only when ``use_precursors=True``.
 
 .. note::
 
@@ -516,6 +526,15 @@ Meaning of ``initial_state``:
    ``initial_state="existing"`` is the usual choice when a transient begins
    from a previously computed steady-state or eigenvalue solution.
 
+.. note::
+
+   A steady-state or k-eigenvalue initial condition computed with delayed
+   cross-section data uses total fission production. A following transient uses
+   delayed precursors only if the problem option ``use_precursors`` is ``True``.
+   If ``use_precursors`` is ``False``, the transient is prompt-only even though
+   the steady initial condition may have included delayed production in its
+   eigenvalue or source balance.
+
 .. note::
 
    ``initial_state="zero"`` is mostly useful for controlled studies where the
@@ -730,6 +749,12 @@ At a high level, the solver:
 This solver uses the groupset inner solver and AGS settings configured on the
 problem.
 
+If the cross sections contain delayed-neutron data, the fission source and
+``k_eff`` update use total fission production, prompt plus delayed. This
+steady-state behavior is independent of ``use_precursors``. The
+``use_precursors`` option only controls whether precursor concentrations are
+stored and computed after the solve.
+
 .. note::
 
    Power iteration is usually the first k-eigen method to try because it uses

framework/materials/multi_group_xs/openmc_import.cc

@@ -8,6 +8,7 @@
 #include "framework/utils/utils.h"
 #include <numeric>
 #include <algorithm>
+#include <cstdint>
 
 namespace opensn
 {
@@ -71,7 +72,11 @@ MultiGroupXS::LoadFromOpenMC(const std::string& file_name,
     throw std::runtime_error(file_name + " has an unsupported scatter shape");
 
   // Number of precursors
-  mgxs.num_precursors_ = 0;
+  unsigned int delayed_groups = 0;
+  if (not H5ReadAttribute<unsigned int>(file.Id(), "delayed_groups", delayed_groups))
+    delayed_groups = 0;
+  OpenSnLogicalErrorIf(delayed_groups < 0, "OpenMC number of delayed_groups must be non-negative.");
+  mgxs.num_precursors_ = delayed_groups;
 
   // Inverse velocity
   H5ReadDataset1D<double>(file.Id(), path + "inverse-velocity", mgxs.inv_velocity_);
@@ -162,6 +167,96 @@ MultiGroupXS::LoadFromOpenMC(const std::string& file_name,
       if (mgxs.sigma_f_[g] > 0.0)
         nu[g] /= mgxs.sigma_f_[g];
 
+    if (mgxs.num_precursors_ > 0)
+    {
+      std::vector<double> nu_prompt_sigma_f;
+      if (H5Has(file.Id(), path + "prompt-nu-fission"))
+        H5ReadDataset1D<double>(file.Id(), path + "prompt-nu-fission", nu_prompt_sigma_f);
+
+      std::vector<std::vector<double>> delayed_nu_sigma_f_by_prec;
+      if (not H5ReadDataset2D<double>(
+            file.Id(), path + "delayed-nu-fission", delayed_nu_sigma_f_by_prec))
+      {
+        std::vector<double> delayed_nu_sigma_f_1d;
+        OpenSnLogicalErrorIf(not H5ReadDataset1D<double>(
+                               file.Id(), path + "delayed-nu-fission", delayed_nu_sigma_f_1d),
+                             "Failed to read dataset \"" + path + "delayed-nu-fission\".");
+        OpenSnLogicalErrorIf(
+          mgxs.num_precursors_ != 1 or delayed_nu_sigma_f_1d.size() != mgxs.num_groups_,
+          "Dataset \"" + path + "delayed-nu-fission\" has incorrect dimensions.");
+        delayed_nu_sigma_f_by_prec.push_back(std::move(delayed_nu_sigma_f_1d));
+      }
+      std::vector<double> nu_delayed_sigma_f(mgxs.num_groups_, 0.0);
+      for (const auto& delayed_family : delayed_nu_sigma_f_by_prec)
+        for (unsigned int g = 0; g < mgxs.num_groups_; ++g)
+          nu_delayed_sigma_f[g] += delayed_family[g];
+
+      if (nu_prompt_sigma_f.empty())
+      {
+        nu_prompt_sigma_f = mgxs.nu_sigma_f_;
+        for (unsigned int g = 0; g < mgxs.num_groups_; ++g)
+          nu_prompt_sigma_f[g] -= nu_delayed_sigma_f[g];
+      }
+
+      OpenSnLogicalErrorIf(not IsNonNegative(nu_prompt_sigma_f),
+                           "OpenMC delayed fission data imply a negative prompt nu-sigma-f.");
+      OpenSnLogicalErrorIf(not IsNonNegative(nu_delayed_sigma_f),
+                           "OpenMC delayed nu-sigma-f must be non-negative.");
+
+      std::vector<double> precursor_decay_constants;
+      H5ReadDataset1D<double>(file.Id(), path + "decay-rate", precursor_decay_constants);
+      OpenSnLogicalErrorIf(precursor_decay_constants.size() != mgxs.num_precursors_,
+                           "OpenMC decay-rate data has incorrect size.");
+
+      std::vector<std::vector<double>> chi_delayed_by_prec;
+      if (not H5ReadDataset2D<double>(file.Id(), path + "chi-delayed", chi_delayed_by_prec))
+      {
+        std::vector<double> chi_delayed_1d;
+        OpenSnLogicalErrorIf(
+          not H5ReadDataset1D<double>(file.Id(), path + "chi-delayed", chi_delayed_1d),
+          "Failed to read dataset \"" + path + "chi-delayed\".");
+        OpenSnLogicalErrorIf(mgxs.num_precursors_ != 1 or chi_delayed_1d.size() != mgxs.num_groups_,
+                             "Dataset \"" + path + "chi-delayed\" has incorrect dimensions.");
+        chi_delayed_by_prec.push_back(std::move(chi_delayed_1d));
+      }
+
+      std::vector<double> chi_prompt = mgxs.chi_;
+      if (H5Has(file.Id(), path + "chi-prompt"))
+      {
+        H5ReadDataset1D<double>(file.Id(), path + "chi-prompt", chi_prompt);
+        OpenSnLogicalErrorIf(chi_prompt.size() != mgxs.num_groups_,
+                             "OpenMC chi-prompt data has incorrect size.");
+      }
+      else
+      {
+        log.Log0Warning() << "OpenMC delayed-neutron data found in \"" << file_name
+                          << "\" without chi-prompt. Using chi as the prompt spectrum.";
+      }
+
+      mgxs.nu_prompt_sigma_f_ = std::move(nu_prompt_sigma_f);
+      mgxs.nu_delayed_sigma_f_ = std::move(nu_delayed_sigma_f);
+      mgxs.precursors_.resize(mgxs.num_precursors_);
+
+      double delayed_total = 0.0;
+      for (const auto& delayed_family : delayed_nu_sigma_f_by_prec)
+        delayed_total += std::accumulate(delayed_family.begin(), delayed_family.end(), 0.0);
+
+      for (unsigned int j = 0; j < mgxs.num_precursors_; ++j)
+      {
+        auto& precursor = mgxs.precursors_[j];
+        precursor.decay_constant = precursor_decay_constants[j];
+        precursor.emission_spectrum = chi_delayed_by_prec[j];
+        const double family_total = std::accumulate(
+          delayed_nu_sigma_f_by_prec[j].begin(), delayed_nu_sigma_f_by_prec[j].end(), 0.0);
+        precursor.fractional_yield = delayed_total > 0.0 ? family_total / delayed_total : 0.0;
+      }
+
+      mgxs.production_matrix_.assign(mgxs.num_groups_, std::vector<double>(mgxs.num_groups_, 0.0));
+      for (unsigned int gp = 0; gp < mgxs.num_groups_; ++gp)
+        for (unsigned int g = 0; g < mgxs.num_groups_; ++g)
+          mgxs.production_matrix_[g][gp] = chi_prompt[g] * mgxs.nu_prompt_sigma_f_[gp];
+    }
+
     // Production matrix (computed)
     // TODO: This path uses chi * nu_sigma_f (total production). If delayed-neutron data is
     // introduced here in the future, ensure LBS source/fission-production routines remain
@@ -189,6 +284,7 @@ MultiGroupXS::LoadFromOpenMC(const std::string& file_name,
   else
   {
     // Clear fission data if not fissionable
+    mgxs.num_precursors_ = 0;
     mgxs.sigma_f_.clear();
     mgxs.nu_sigma_f_.clear();
     mgxs.nu_prompt_sigma_f_.clear();

framework/utils/hdf_utils.h

@@ -312,6 +312,52 @@ H5ReadDataset1D(hid_t id, const std::string& name, std::vector<T>& data)
   return success;
 }
 
+template <typename T>
+bool
+H5ReadDataset2D(hid_t id, const std::string& name, std::vector<std::vector<T>>& data)
+{
+  static_assert(
+    !std::is_same_v<T, bool>,
+    "H5ReadDataset2D does not support std::vector<std::vector<bool>>. Use unsigned char.");
+
+  bool success = false;
+  data.clear();
+
+  auto dataset = H5Dopen2(id, name.c_str(), H5P_DEFAULT);
+  if (dataset == H5I_INVALID_HID)
+    return false;
+
+  auto dataspace = H5Dget_space(dataset);
+  if (dataspace != H5I_INVALID_HID)
+  {
+    const int rank = H5Sget_simple_extent_ndims(dataspace);
+    if (rank == 2)
+    {
+      hsize_t dims[2] = {0, 0};
+      if (H5Sget_simple_extent_dims(dataspace, dims, nullptr) >= 0)
+      {
+        auto n0 = static_cast<size_t>(dims[0]);
+        auto n1 = static_cast<size_t>(dims[1]);
+        std::vector<T> flat_data(n0 * n1);
+        if (flat_data.empty() or
+            H5Dread(dataset, get_datatype<T>(), H5S_ALL, H5S_ALL, H5P_DEFAULT, flat_data.data()) >=
+              0)
+        {
+          data.assign(n0, std::vector<T>(n1));
+          for (size_t i = 0; i < n0; ++i)
+            for (size_t j = 0; j < n1; ++j)
+              data[i][j] = flat_data[i * n1 + j];
+          success = true;
+        }
+      }
+    }
+    H5Sclose(dataspace);
+  }
+
+  H5Dclose(dataset);
+  return success;
+}
+
 template <typename T>
 bool
 H5ReadAttribute(hid_t id, const std::string& name, T& value)