diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index ae39a53ddf5..a58c2d5e62c 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -84,6 +84,10 @@ class Nuclide { double collapse_rate(int MT, double temperature, span energy, span flux) const; + void set_xs(int MT, int T_index, const std::vector& values); + void rebuild_derived_xs(); + std::vector get_xs(int MT, int T_index) const; + std::vector get_energy_grid(int T_index) const; //! Return a ParticleType object representing this nuclide ParticleType particle_type() const { return {Z_, A_, metastable_}; } diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 0ee8c3ff06b..a9659b20c13 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -11,6 +11,7 @@ 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 @@ class _SourceSite(Structure): _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}" + ) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 17d6e952c39..51ffbc0d824 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -493,6 +493,111 @@ void Nuclide::create_derived( } } +//====================================================== +//===================== get current xs ======== + +std::vector 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(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(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(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(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(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(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(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(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(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(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(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(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(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(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 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