Merge 670fee26cb into f1fb6721f0

include/openmc/nuclide.h

@@ -84,6 +84,10 @@ public:
   double collapse_rate(int MT, double temperature, span<const double> energy,
     span<const double> flux) const;
 
+  void set_xs(int MT, int T_index, const std::vector<double>& values);
+  void rebuild_derived_xs();
+  std::vector<double> get_xs(int MT, int T_index) const;
+  std::vector<double> get_energy_grid(int T_index) const;
   //! Return a ParticleType object representing this nuclide
   ParticleType particle_type() const { return {Z_, A_, metastable_}; }
 

openmc/lib/core.py

@@ -11,6 +11,7 @@ import traceback as tb
 
 import numpy as np
 from numpy.ctypeslib import as_array
+from ctypes import byref
 
 from . import _dll
 from .error import _error_handler
@@ -115,6 +116,36 @@ _dll.openmc_global_bounding_box.errcheck = _error_handler
 _dll.openmc_sample_external_source.argtypes = [c_size_t, POINTER(c_uint64), POINTER(_SourceSite)]
 _dll.openmc_sample_external_source.restype = c_int
 _dll.openmc_sample_external_source.errcheck = _error_handler
+#----------------------------------------------------------------------------
+# --- Nuclide functions -----------------------------------------------------
+#----------------------------------------------------------------------------
+_dll.openmc_nuclide_rebuild_derived_xs.argtypes = [c_int]
+_dll.openmc_nuclide_rebuild_derived_xs.restype  = c_int
+_dll.openmc_nuclide_get_energy_grid.argtypes = [c_int, c_int, POINTER(c_double), c_int]
+_dll.openmc_nuclide_get_energy_grid.restype = c_int
+_dll.openmc_nuclide_get_energy_grid_size.argtypes = [c_int, c_int, POINTER(c_int)]
+_dll.openmc_nuclide_get_energy_grid_size.restype = c_int
+_dll.openmc_nuclide_get_xs_size.argtypes = [c_int, c_int, c_int, POINTER(c_int)]
+_dll.openmc_nuclide_get_xs_size.restype = c_int
+_dll.openmc_nuclide_get_xs.argtypes = [c_int, c_int, c_int, POINTER(c_double), c_int]
+_dll.openmc_nuclide_get_xs.restype = c_int
+_dll.openmc_nuclide_get_reaction_threshold_energy.argtypes = [c_int, c_int, c_int, POINTER(c_double)]
+_dll.openmc_nuclide_get_reaction_threshold_energy.restype = c_int
+_dll.openmc_nuclide_get_mt_numbers.argtypes = [c_int, POINTER(c_int), POINTER(c_int)]
+_dll.openmc_nuclide_get_mt_numbers.restype = c_int
+_dll.openmc_nuclide_set_reaction_xs_with_threshold.argtypes = [
+    c_int,                     # index
+    c_int,                     # mt
+    c_int,                     # T_index
+    POINTER(c_double),         # energy
+    POINTER(c_double),         # values
+    c_int,                     # n
+    c_double                   # threshold_energy
+]
+_dll.openmc_nuclide_set_reaction_xs_with_threshold.restype = c_int
+#----------------------------------------------------------------------------
+#----------------------------------------------------------------------------
+
 
 def global_bounding_box():
     """Calculate a global bounding box for the model"""
@@ -827,3 +858,201 @@ def quiet_dll(output=True):
             sys.stdout = initial_stdout
             sys.stdout.flush()
             os.dup2(initial_stdout_fno, 1)
+
+# --- Nuclide functions -----------------------------------------------------
+
+def nuclide_rebuild_derived_xs(index):
+    # C function: int openmc_nuclide_rebuild_derived_xs(int index)
+    """Rebuild derived cross sections for nuclide at given index."""
+    err = _dll.openmc_nuclide_rebuild_derived_xs(index)
+    if err != 0:
+        raise RuntimeError(f"OpenMC core error in nuclide_rebuild_derived_xs (index={index}: {err}")
+
+def nuclide_energy_grid(index, temperature_index=0):
+    """Return the energy grid for a nuclide at a given temperature index."""
+
+    n = c_int()
+
+    err = _dll.openmc_nuclide_get_energy_grid_size(index, temperature_index, byref(n))
+    if err != 0:
+        raise RuntimeError(f"OpenMC error getting energy grid size (index={index}): {err}")
+
+    E = np.empty(n.value, dtype=np.float64)
+
+    err = _dll.openmc_nuclide_get_energy_grid(index, temperature_index, E.ctypes.data_as(POINTER(c_double)),n.value)
+    if err != 0:
+        raise RuntimeError(f"OpenMC core error in nuclide_energy_grid (index={index}: {err}")
+
+    return E
+
+def get_nuclide_xs(index, mt, temperature=0):
+    # Get XS size
+    n = c_int()
+
+    err = _dll.openmc_nuclide_get_xs_size(index, mt, temperature, byref(n))
+    if err != 0:
+        raise RuntimeError(f"OpenMC core error in get_nuclide_xs ref (index={index}, mt={mt}): {err}")
+
+    xs = np.zeros(n.value, dtype=np.float64)
+
+    err = _dll.openmc_nuclide_get_xs(index, mt, temperature, xs.ctypes.data_as(POINTER(c_double)), n.value)
+    if err != 0:
+        raise RuntimeError(f"OpenMC core error in get_nuclide_xs (index={index}, mt={mt}): {err}")
+    
+    return xs
+
+def get_reaction_threshold_energy(index, mt, temperature=0):
+    """Return reaction threshold energy in eV."""
+    threshold = c_double()
+
+    err = _dll.openmc_nuclide_get_reaction_threshold_energy(index, mt, temperature, byref(threshold))
+
+    if err != 0:
+        raise RuntimeError(f"OpenMC core error in get_reaction_threshold_energy (index={index}, mt={mt}): {err}")
+    return threshold.value
+
+def get_nuclide_mt_numbers(index):
+    # First, find how many MTs there are
+
+    n_mts = c_int()
+
+    # We don’t know the array size yet, but C++ sets *n_mts
+    # So call once with a null pointer
+    err = _dll.openmc_nuclide_get_mt_numbers(index, None, byref(n_mts))
+    if err != 0:
+        raise RuntimeError(f"OpenMC error getting MT count: {err}")
+    # Allocate array
+    mts = np.zeros(n_mts.value, dtype=np.int32)
+
+    # Call again to fill it
+    err = _dll.openmc_nuclide_get_mt_numbers(
+        index,
+        mts.ctypes.data_as(POINTER(c_int)),
+        byref(n_mts)
+    )
+    if err != 0:
+        raise RuntimeError(f"OpenMC error getting MT list: {err}")
+
+    return mts
+    
+
+def new_nuclide_xs_with_threshold(
+    index,
+    mt,
+    energy,
+    xs_perturbed,
+    threshold_energy,
+    temperature=0
+):
+    """
+    Set a reaction cross section with a modified threshold.
+
+    Parameters
+    ----------
+    index : int
+        Nuclide index in OpenMC.
+    mt : int
+        MT number of the reaction.
+    energy : numpy.ndarray
+        Full nuclide energy grid (must match OpenMC's grid).
+    xs_perturbed : numpy.ndarray
+        Cross section values defined on the full grid.
+        Must be zero below threshold_energy.
+    threshold_energy : float
+        Threshold energy in eV.
+    temperature : int, optional
+        Temperature index (default is 0).
+    """
+
+    import numpy as np
+
+    # --- Basic validation on Python side ---
+    if not isinstance(energy, np.ndarray) or not isinstance(xs_perturbed, np.ndarray):
+        raise TypeError("energy and xs_perturbed must be numpy arrays")
+
+    if energy.dtype != np.float64:
+        energy = np.asarray(energy, dtype=np.float64)
+
+    if xs_perturbed.dtype != np.float64:
+        xs_perturbed = np.asarray(xs_perturbed, dtype=np.float64)
+
+    if energy.size != xs_perturbed.size:
+        raise ValueError("energy and xs_perturbed must have the same length")
+
+    # --- Call C API ---
+    err = _dll.openmc_nuclide_set_reaction_xs_with_threshold(
+        index,
+        mt,
+        temperature,
+        energy.ctypes.data_as(POINTER(c_double)),
+        xs_perturbed.ctypes.data_as(POINTER(c_double)),
+        energy.size,
+        threshold_energy
+    )
+
+    if err != 0:
+        raise RuntimeError(
+            f"OpenMC core error in new_nuclide_xs_with_threshold "
+            f"(index={index}, mt={mt}, T_index={temperature}): {err}"
+        )
+
+def new_nuclide_xs_with_threshold_old(
+    index,
+    mt,
+    energy,
+    xs_perturbed,
+    threshold_energy,
+    temperature=0
+):
+    """
+    Set a reaction cross section with a modified threshold.
+
+    Parameters
+    ----------
+    index : int
+        Nuclide index in OpenMC.
+    mt : int
+        MT number of the reaction.
+    energy : numpy.ndarray
+        Full nuclide energy grid (must match OpenMC's grid).
+    xs_perturbed : numpy.ndarray
+        Cross section values defined on the full grid.
+        Must be zero below threshold_energy.
+    threshold_energy : float
+        Threshold energy in eV.
+    temperature : int, optional
+        Temperature index (default is 0).
+    """
+
+
+
+    # --- Basic validation on Python side ---
+    if not isinstance(energy, np.ndarray) or not isinstance(xs_perturbed, np.ndarray):
+        raise TypeError("energy and xs_perturbed must be numpy arrays")
+
+    if energy.dtype != np.float64:
+        energy = np.asarray(energy, dtype=np.float64)
+
+    if xs_perturbed.dtype != np.float64:
+        xs_perturbed = np.asarray(xs_perturbed, dtype=np.float64)
+    '''
+    if energy.size != xs_perturbed.size:
+        raise ValueError("energy and xs_perturbed must have the same length")
+    '''
+
+    # --- Call C API ---
+    err = _dll.openmc_nuclide_set_reaction_xs_with_threshold(
+        index,
+        mt,
+        temperature,
+        energy.ctypes.data_as(POINTER(c_double)),
+        xs_perturbed.ctypes.data_as(POINTER(c_double)),
+        energy.size,
+        threshold_energy
+    )
+
+    if err != 0:
+        raise RuntimeError(
+            f"OpenMC core error in new_nuclide_xs_with_threshold "
+            f"(index={index}, mt={mt}, T_index={temperature}): {err}"
+        )

src/nuclide.cpp

@@ -493,6 +493,111 @@ void Nuclide::create_derived(
   }
 }
 
+//======================================================
+//===================== get current xs ========
+
+std::vector<double> Nuclide::get_xs(int MT, int T_index) const
+{
+  // Find reaction index
+  size_t i_rx = reaction_index_.at(MT);
+  if (i_rx == C_NONE) {
+    throw std::runtime_error("No reaction with MT = " + std::to_string(MT));
+  }
+
+  const auto& rx = *reactions_[i_rx];
+  if (T_index >= rx.xs_.size()) {
+    throw std::runtime_error("Temperature index out of range");
+  }
+
+  return rx.xs_[T_index].value;
+}
+
+//===================== end of code get current xs =====
+//======================================================
+
+//======================================================
+//== Rebuild derived xs_ from reactions_ data =
+
+void Nuclide::rebuild_derived_xs()
+{
+  // Loop over temperature indices
+  for (int t = 0; t < static_cast<int>(kTs_.size()); ++t) {
+
+    auto& xs_t = xs_[t];
+
+    // Reset derived cross sections
+    xs_t.fill(0.0);
+
+    // Loop over all reactions
+    for (const auto& rx_ptr : reactions_) {
+      const auto& rx = *rx_ptr;
+
+      // Skip redundant reactions
+      if (rx.redundant_) {
+        continue;
+      }
+
+      const int j = rx.xs_[t].threshold;
+      const int n = static_cast<int>(rx.xs_[t].value.size());
+
+      const auto& vals = rx.xs_[t].value;
+
+      // -----------------------------
+      // TOTAL CROSS SECTION
+      // -----------------------------
+      for (int k = 0; k < n; ++k) {
+        xs_t(j + k, XS_TOTAL) += vals[k];
+      }
+
+      // -----------------------------
+      // ABSORPTION (disappearance)
+      // -----------------------------
+      if (is_disappearance(rx.mt_)) {
+        for (int k = 0; k < n; ++k) {
+          xs_t(j + k, XS_ABSORPTION) += vals[k];
+        }
+      }
+
+      // -----------------------------
+      // FISSION
+      // -----------------------------
+      if (is_fission(rx.mt_)) {
+        for (int k = 0; k < n; ++k) {
+          xs_t(j + k, XS_FISSION) += vals[k];
+          xs_t(j + k, XS_ABSORPTION) += vals[k];
+        }
+      }
+
+      // -----------------------------
+      // PHOTON
+      // -----------------------------
+      for (const auto& p : rx.products_) {
+        if (p.particle_.is_photon()) {
+          for (int k = 0; k < n; ++k) {
+            double E = grid_[t].energy[j + k];
+            xs_t(j + k, XS_PHOTON_PROD) += vals[k] * (*p.yield_)(E);
+          }
+        }
+      }
+    }
+
+    // -----------------------------
+    // NU-FISSION RECONSTRUCTION
+    // -----------------------------
+    if (fissionable_) {
+      const int ngrid = static_cast<int>(grid_[t].energy.size());
+      for (int i = 0; i < ngrid; ++i) {
+        double E = grid_[t].energy[i];
+        xs_t(i, XS_NU_FISSION) =
+          nu(E, EmissionMode::total) * xs_t(i, XS_FISSION);
+      }
+    }
+  }
+}
+
+//=end of code Rebuild derived xs_ from reactions_ data=
+//======================================================
+
 void Nuclide::init_grid()
 {
   int neutron = ParticleType::neutron().transport_index();
@@ -1226,4 +1331,274 @@ bool multipole_in_range(const Nuclide& nuc, double E)
   return E >= nuc.multipole_->E_min_ && E <= nuc.multipole_->E_max_;
 }
 
+extern "C" int openmc_nuclide_get_xs(
+  int i_nuclide, int MT, int T_index, double* xs_out, int n_values)
+{
+  using namespace openmc;
+
+  if (i_nuclide < 0 || i_nuclide >= static_cast<int>(data::nuclides.size()))
+    return OPENMC_E_INVALID_ID;
+
+  try {
+    const Nuclide* nuc = data::nuclides[i_nuclide].get();
+    auto xs = nuc->get_xs(MT, T_index);
+    if (xs.size() != static_cast<size_t>(n_values)) {
+      set_errmsg("Requested length does not match stored XS length");
+      return OPENMC_E_DATA;
+    }
+
+    // std::memcpy(xs_out, xs.data(), n_values * sizeof(double));
+    for (int i = 0; i < n_values; ++i) {
+      xs_out[i] = xs[i];
+    }
+
+    return 0;
+  } catch (const std::exception& e) {
+    set_errmsg(e.what());
+    return OPENMC_E_DATA;
+  }
+}
+
+//======================================================
+//===================== new code get current xs size ===
+
+extern "C" int openmc_nuclide_get_xs_size(
+  int index, int MT, int T_index, int* n)
+{
+  if (index < 0 || index >= data::nuclides.size()) {
+    set_errmsg("Index in nuclides vector is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  try {
+    auto& nuc = *data::nuclides[index];
+    size_t i_rx = nuc.reaction_index_[MT];
+    if (i_rx == C_NONE) {
+      set_errmsg("No reaction with MT = " + std::to_string(MT));
+      return OPENMC_E_DATA;
+    }
+
+    const auto& rx = *nuc.reactions_[i_rx];
+    if (T_index >= rx.xs_.size()) {
+      set_errmsg("Temperature index out of range");
+      return OPENMC_E_DATA;
+    }
+
+    *n = static_cast<int>(rx.xs_[T_index].value.size());
+  } catch (const std::exception& e) {
+    set_errmsg(e.what());
+    return OPENMC_E_DATA;
+  }
+
+  return 0;
+}
+//===================== end of code get current sizexs ===
+//========================================================
+extern "C" int openmc_nuclide_get_energy_grid(
+  int index, int T_index, double* E, int n)
+{
+  using namespace openmc;
+
+  if (index < 0 || index >= data::nuclides.size())
+    return OPENMC_E_INVALID_ARGUMENT;
+
+  try {
+    const auto& nuc = *data::nuclides[index];
+    const auto& grid = nuc.grid_.at(T_index).energy;
+
+    if (static_cast<int>(grid.size()) != n)
+      return OPENMC_E_INVALID_ARGUMENT;
+
+    std::memcpy(E, grid.data(), n * sizeof(double));
+  } catch (const std::exception& e) {
+    set_errmsg(e.what());
+    return OPENMC_E_DATA;
+  }
+
+  return 0;
+}
+
+extern "C" int openmc_nuclide_rebuild_derived_xs(int index)
+{
+  if (index < 0 || index >= data::nuclides.size()) {
+    set_errmsg("Index in nuclides vector out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  try {
+    data::nuclides[index]->rebuild_derived_xs();
+  } catch (const std::exception& e) {
+    set_errmsg(e.what());
+    return OPENMC_E_DATA;
+  }
+  return 0;
+}
+
+extern "C" int openmc_nuclide_get_energy_grid_size(
+  int index, int T_index, int* n)
+{
+  using namespace openmc;
+
+  if (index < 0 || index >= data::nuclides.size())
+    return OPENMC_E_INVALID_ARGUMENT;
+
+  try {
+    const auto& nuc = *data::nuclides[index];
+    const auto& grid = nuc.grid_.at(T_index).energy;
+
+    *n = static_cast<int>(grid.size());
+  } catch (const std::exception& e) {
+    set_errmsg(e.what());
+    return OPENMC_E_DATA;
+  }
+
+  return 0;
+}
+
+extern "C" int openmc_nuclide_get_reaction_threshold_energy(
+  int index, int mt, int T_index, double* threshold)
+{
+  using namespace openmc;
+
+  if (index < 0 || index >= data::nuclides.size())
+    return OPENMC_E_INVALID_ARGUMENT;
+  if (!threshold)
+    return OPENMC_E_INVALID_ARGUMENT;
+
+  try {
+    const auto& nuc = *data::nuclides[index];
+
+    // Find reaction
+    const Reaction* r = nullptr;
+    for (const auto& rx : nuc.reactions_) {
+      if (rx && rx->mt_ == mt) {
+        r = rx.get();
+        break;
+      }
+    }
+    if (!r)
+      return OPENMC_E_INVALID_ARGUMENT;
+
+    // Access threshold index in NuclideMicroXS
+    if (T_index < 0 || T_index >= static_cast<int>(r->xs_.size()))
+      return OPENMC_E_INVALID_ARGUMENT;
+
+    int idx = r->xs_[T_index].threshold;
+
+    // Energy grid for this temperature
+    const auto& grid = nuc.grid_.at(T_index).energy;
+
+    if (idx < 0 || idx >= static_cast<int>(grid.size()))
+      return OPENMC_E_INVALID_ARGUMENT;
+
+    *threshold = grid[idx]; // <-- actual energy in eV
+  } catch (const std::exception& e) {
+    set_errmsg(e.what());
+    return OPENMC_E_DATA;
+  }
+
+  return 0;
+}
+
+extern "C" int openmc_nuclide_get_mt_numbers(int index, int* mts, int* n_mts)
+{
+  using namespace openmc;
+
+  if (index < 0 || index >= data::nuclides.size())
+    return OPENMC_E_INVALID_ARGUMENT;
+
+  try {
+    const auto& nuc = *data::nuclides[index];
+    int n = static_cast<int>(nuc.reactions_.size());
+    *n_mts = n;
+
+    if (mts == nullptr)
+      return 0;
+
+    for (int i = 0; i < n; ++i) {
+      mts[i] = nuc.reactions_[i]->mt_;
+    }
+  } catch (const std::exception& e) {
+    set_errmsg(e.what());
+    return OPENMC_E_DATA;
+  }
+
+  return 0;
+}
+
+extern "C" int openmc_nuclide_set_reaction_xs_with_threshold(int index, int mt,
+  int T_index, const double* energy, const double* values, int n,
+  double threshold_energy)
+{
+  using namespace openmc;
+
+  if (index < 0 || index >= static_cast<int>(data::nuclides.size()))
+    return OPENMC_E_INVALID_ARGUMENT;
+  if (!energy || !values || n <= 0)
+    return OPENMC_E_INVALID_ARGUMENT;
+
+  try {
+    auto& nuc = *data::nuclides[index];
+
+    // --- Find reaction ---
+    Reaction* rx = nullptr;
+    for (auto& rx_ptr : nuc.reactions_) {
+      if (rx_ptr && rx_ptr->mt_ == mt) {
+        rx = rx_ptr.get();
+        break;
+      }
+    }
+    if (!rx)
+      return OPENMC_E_INVALID_ARGUMENT;
+
+    if (T_index < 0 || T_index >= static_cast<int>(rx->xs_.size()))
+      return OPENMC_E_INVALID_ARGUMENT;
+
+    auto& xs = rx->xs_[T_index];
+
+    // --- Validate energy grid ---
+    const auto& nuc_grid = nuc.grid_.at(T_index).energy;
+
+    if (static_cast<int>(nuc_grid.size()) != n)
+      return OPENMC_E_INVALID_ARGUMENT;
+
+    for (int i = 0; i < n; ++i) {
+      if (energy[i] != nuc_grid[i])
+        return OPENMC_E_INVALID_ARGUMENT;
+    }
+
+    // --- Compute threshold index ---
+    int thr_idx = -1;
+    for (int i = 0; i < n; ++i) {
+      if (energy[i] >= threshold_energy) {
+        thr_idx = i;
+        break;
+      }
+    }
+    if (thr_idx < 0)
+      return OPENMC_E_INVALID_ARGUMENT;
+
+    // --- Enforce zero below threshold ---
+    for (int i = 0; i < thr_idx; ++i) {
+      if (values[i] != 0.0)
+        return OPENMC_E_INVALID_ARGUMENT;
+    }
+
+    // --- Trim XS ---
+    std::vector<double> new_values;
+    new_values.reserve(n - thr_idx);
+
+    for (int i = thr_idx; i < n; ++i)
+      new_values.push_back(values[i]);
+
+    // --- Update reaction --
+    xs.value = std::move(new_values);
+    xs.threshold = thr_idx;
+  } catch (const std::exception& e) {
+    set_errmsg(e.what());
+    return OPENMC_E_DATA;
+  }
+
+  return 0;
+}
+
 } // namespace openmc