diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 0d8189caa1d..f3cace3125b 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -994,6 +994,15 @@ class LibMesh : public UnstructuredMesh { LibMesh(const std::string& filename, double length_multiplier = 1.0); LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0); + //! Create a mesh from an externally constructed libMesh mesh, transferring + //! ownership of the mesh to OpenMC + // + //! \param[in] input_mesh Externally built mesh (must be replicated) + //! \param[in] length_multiplier Multiplier applied to mesh coordinates + //! \param[in] filename Name of the file the mesh was read from, if any + LibMesh(unique_ptr input_mesh, + double length_multiplier = 1.0, const std::string& filename = ""); + static const std::string mesh_lib_type; // Overridden Methods diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index a5d404133ce..6e04fc690f3 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -244,6 +244,11 @@ void apply_weight_window(Particle& p, WeightWindow weight_window); //! Free memory associated with weight windows void free_memory_weight_windows(); +//! Build a WeightWindows object from multigroup adjoint flux stored as +//! elemental data in an Exodus II file (requires libMesh support) +//! \param[in] node XML node for in settings.xml +void read_weight_windows_exodus(pugi::xml_node node); + //! Search weight window that apply to a particle //! \param[in] p Particle to search weight window for std::pair search_weight_window(const Particle& p); diff --git a/openmc/settings.py b/openmc/settings.py index 8120eb073e6..3b6e3e54951 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -17,7 +17,8 @@ from .source import SourceBase, MeshSource, IndependentSource from .utility_funcs import input_path from .volume import VolumeCalculation -from .weight_windows import WeightWindows, WeightWindowGenerator, WeightWindowsList +from .weight_windows import (WeightWindows, WeightWindowGenerator, + WeightWindowsList, WeightWindowsExodus) class RunMode(Enum): @@ -395,6 +396,12 @@ class Settings: Path to a weight window file to load during simulation initialization .. versionadded::0.14.0 + weight_windows_exodus : openmc.WeightWindowsExodus + Specification for building weight windows from multigroup adjoint flux + stored as elemental data in an Exodus file. Requires OpenMC to be + built with libMesh support. + + .. versionadded:: 0.15.4 write_initial_source : bool Indicate whether to write the initial source distribution to file """ @@ -493,6 +500,7 @@ def __init__(self, **kwargs): self._weight_windows_on = None self._shared_secondary_bank = None self._weight_windows_file = None + self._weight_windows_exodus = None self._weight_window_checkpoints = {} self._max_history_splits = None self._max_tracks = None @@ -1383,6 +1391,16 @@ def weight_windows_file(self, value: PathLike | None): cv.check_type('weight windows file', value, PathLike) self._weight_windows_file = input_path(value) + @property + def weight_windows_exodus(self) -> WeightWindowsExodus | None: + return self._weight_windows_exodus + + @weight_windows_exodus.setter + def weight_windows_exodus(self, value: WeightWindowsExodus | None): + if value is not None: + cv.check_type('weight windows exodus', value, WeightWindowsExodus) + self._weight_windows_exodus = value + @property def weight_window_generators(self) -> list[WeightWindowGenerator]: return self._weight_window_generators @@ -1973,6 +1991,10 @@ def _create_weight_windows_file_element(self, root): element.text = str(self.weight_windows_file) root.append(element) + def _create_weight_windows_exodus_subelement(self, root): + if self._weight_windows_exodus is not None: + root.append(self._weight_windows_exodus.to_xml_element()) + def _create_weight_window_checkpoints_subelement(self, root): if not self._weight_window_checkpoints: return @@ -2464,6 +2486,11 @@ def _weight_windows_file_from_xml_element(self, root): if text is not None: self.weight_windows_file = text + def _weight_windows_exodus_from_xml_element(self, root): + elem = root.find('weight_windows_exodus') + if elem is not None: + self.weight_windows_exodus = WeightWindowsExodus.from_xml_element(elem) + def _weight_window_checkpoints_from_xml_element(self, root): elem = root.find('weight_window_checkpoints') if elem is None: @@ -2629,6 +2656,7 @@ def to_xml_element(self, mesh_memo=None): self._create_shared_secondary_bank_subelement(element) self._create_weight_window_generators_subelement(element, mesh_memo) self._create_weight_windows_file_element(element) + self._create_weight_windows_exodus_subelement(element) self._create_weight_window_checkpoints_subelement(element) self._create_max_history_splits_subelement(element) self._create_max_tracks_subelement(element) @@ -2746,6 +2774,7 @@ def from_xml_element(cls, elem, meshes=None): settings._weight_windows_on_from_xml_element(elem) settings._shared_secondary_bank_from_xml_element(elem) settings._weight_windows_file_from_xml_element(elem) + settings._weight_windows_exodus_from_xml_element(elem) settings._weight_window_generators_from_xml_element(elem, meshes) settings._weight_window_checkpoints_from_xml_element(elem) settings._max_history_splits_from_xml_element(elem) diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 63af2596efc..7f9f1ca4e97 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -17,6 +17,7 @@ from ._xml import get_elem_list, get_text, clean_indentation from .mixin import IDManagerMixin from .particle_type import ParticleType +from .utility_funcs import input_path class WeightWindows(IDManagerMixin): @@ -800,6 +801,314 @@ def from_xml_element(cls, elem: ET.Element, meshes: dict) -> Self: return wwg + +class WeightWindowsExodus: + """Specification for building weight windows from an Exodus file. + + The Exodus file is expected to contain multigroup adjoint flux stored as + CONSTANT MONOMIAL elemental variables (one variable per energy group). + At simulation initialization, OpenMC reads the mesh and flux, registers the + mesh as an unstructured (libMesh) mesh, applies FW-CADIS-style normalization + and creates the corresponding weight windows. An instance of this class can + be assigned to the :attr:`openmc.Settings.weight_windows_exodus` attribute. + + Requires OpenMC to be built with libMesh support. + + .. versionadded:: 0.15.4 + + Parameters + ---------- + file : path-like + Path to the Exodus file containing the mesh and adjoint flux + adjoint_flux_variables : iterable of str + Names of the elemental variables containing the adjoint flux, one per + energy group, ordered consistently with `energy_bounds` (ascending + energy). Solvers that write group 0 as the fastest group (e.g. + Griffin) require the variables to be listed thermal-first. + energy_bounds : iterable of float or openmc.mgxs.EnergyGroups + Monotonically increasing energy group boundaries in [eV]. The number + of boundaries must be one more than the number of flux variables. An + :class:`openmc.mgxs.EnergyGroups` instance may be passed directly. + timestep : int, optional + Zero-based index of the Exodus time step to read the flux from. If + not given, the last time step in the file is used. + particle_type : str or int or openmc.ParticleType + Particle type the weight windows apply to (default: 'neutron') + survival_ratio : float, optional + Ratio of the survival weight to the lower weight window bound for + rouletting. If not given, the default of the transport code (3.0) + applies. + upper_bound_ratio : float, optional + Ratio of the upper to lower weight window bounds. If not given, the + default of the transport code (5.0) applies. + max_split : int, optional + Maximum allowable number of particles when splitting. If not given, + the default of the transport code (10) applies. + + Attributes + ---------- + file : pathlib.Path + Path to the Exodus file containing the mesh and adjoint flux + adjoint_flux_variables : list of str + Names of the elemental variables containing the adjoint flux + energy_bounds : numpy.ndarray of float + Monotonically increasing energy group boundaries in [eV] + timestep : int or None + Zero-based index of the Exodus time step to read the flux from + particle_type : openmc.ParticleType + Particle type the weight windows apply to + survival_ratio : float or None + Ratio of the survival weight to the lower weight window bound + upper_bound_ratio : float or None + Ratio of the upper to lower weight window bounds + max_split : int or None + Maximum allowable number of particles when splitting + + See Also + -------- + openmc.Settings.weight_windows_exodus + + """ + + def __init__( + self, + file: PathLike, + adjoint_flux_variables: Iterable[str], + energy_bounds, + timestep: int | None = None, + particle_type: str | int | openmc.ParticleType = 'neutron', + survival_ratio: float | None = None, + upper_bound_ratio: float | None = None, + max_split: int | None = None + ): + self.file = file + self.adjoint_flux_variables = adjoint_flux_variables + self.energy_bounds = energy_bounds + self.timestep = timestep + self.particle_type = particle_type + self.survival_ratio = survival_ratio + self.upper_bound_ratio = upper_bound_ratio + self.max_split = max_split + self._check_consistency() + + def _check_consistency(self): + """Cross-attribute checks mirroring those performed by the C++ layer""" + n_groups = len(self.adjoint_flux_variables) + if self.energy_bounds.size != n_groups + 1: + raise ValueError( + f'Number of energy bounds ({self.energy_bounds.size}) must be ' + f'one more than the number of adjoint flux variables ' + f'({n_groups}).') + # compare using the transport code defaults when a value is unset + survival = 3.0 if self.survival_ratio is None else self.survival_ratio + upper = 5.0 if self.upper_bound_ratio is None else self.upper_bound_ratio + if upper <= survival: + raise ValueError( + f'Upper bound ratio ({upper}) must be larger than the ' + f'survival ratio ({survival}).') + + def __repr__(self) -> str: + string = type(self).__name__ + '\n' + string += f'\t{"File":<20}=\t{self.file}\n' + string += f'\t{"Flux variables":<20}=\t{self.adjoint_flux_variables}\n' + string += f'\t{"Energy bounds":<20}=\t{self.energy_bounds}\n' + string += f'\t{"Timestep":<20}=\t{self.timestep}\n' + string += f'\t{"Particle":<20}=\t{str(self.particle_type)}\n' + string += f'\t{"Survival ratio":<20}=\t{self.survival_ratio}\n' + string += f'\t{"Upper bound ratio":<20}=\t{self.upper_bound_ratio}\n' + string += f'\t{"Max split":<20}=\t{self.max_split}\n' + return string + + def __eq__(self, other) -> bool: + if not isinstance(other, WeightWindowsExodus): + return False + attrs = ('file', 'adjoint_flux_variables', 'timestep', + 'particle_type', 'survival_ratio', 'upper_bound_ratio', + 'max_split') + for attr in attrs: + if getattr(self, attr) != getattr(other, attr): + return False + return np.array_equal(self.energy_bounds, other.energy_bounds) + + @property + def file(self) -> Path: + return self._file + + @file.setter + def file(self, value: PathLike): + cv.check_type('Exodus weight windows file', value, PathLike) + self._file = input_path(value) + + @property + def adjoint_flux_variables(self) -> list[str]: + return self._adjoint_flux_variables + + @adjoint_flux_variables.setter + def adjoint_flux_variables(self, variables: Iterable[str]): + cv.check_type('adjoint flux variables', variables, Iterable, str) + variables = list(variables) + cv.check_greater_than( + 'number of adjoint flux variables', len(variables), 0) + self._adjoint_flux_variables = variables + + @property + def energy_bounds(self) -> np.ndarray: + return self._energy_bounds + + @energy_bounds.setter + def energy_bounds(self, bounds): + # accept an openmc.mgxs.EnergyGroups directly; local import avoids a + # circular import between openmc.weight_windows and openmc.mgxs + from openmc.mgxs import EnergyGroups + if isinstance(bounds, EnergyGroups): + bounds = bounds.group_edges + cv.check_type('energy bounds', bounds, Iterable, Real) + bounds = np.asarray(bounds, dtype=float) + if bounds.ndim != 1 or bounds.size < 2: + raise ValueError('At least two energy bounds must be provided.') + if np.any(np.diff(bounds) <= 0.0) + raise ValueError('Energy bounds must be strictly increasing.') + self._energy_bounds = bounds + + @property + def timestep(self) -> int | None: + return self._timestep + + @timestep.setter + def timestep(self, value: int | None): + if value is not None: + cv.check_type('timestep', value, Integral) + cv.check_greater_than('timestep', value, 0, equality=True) + self._timestep = value + + @property + def particle_type(self) -> ParticleType: + return self._particle_type + + @particle_type.setter + def particle_type(self, pt): + ptype = ParticleType(pt) + if ptype not in {ParticleType.NEUTRON, ParticleType.PHOTON}: + raise ValueError( + 'Weight windows can only be applied for neutrons or photons') + self._particle_type = ptype + + @property + def survival_ratio(self) -> float | None: + return self._survival_ratio + + @survival_ratio.setter + def survival_ratio(self, value: float | None): + if value is not None: + cv.check_type('survival ratio', value, Real) + cv.check_greater_than('survival ratio', value, 1.0) + self._survival_ratio = value + + @property + def upper_bound_ratio(self) -> float | None: + return self._upper_bound_ratio + + @upper_bound_ratio.setter + def upper_bound_ratio(self, value: float | None): + if value is not None: + cv.check_type('upper bound ratio', value, Real) + cv.check_greater_than('upper bound ratio', value, 1.0) + self._upper_bound_ratio = value + + @property + def max_split(self) -> int | None: + return self._max_split + + @max_split.setter + def max_split(self, value: int | None): + if value is not None: + cv.check_type('max split', value, Integral) + cv.check_greater_than('max split', value, 1) + self._max_split = value + + def to_xml_element(self) -> ET.Element: + """Create a 'weight_windows_exodus' element to be written to an XML file. + """ + self._check_consistency() + + element = ET.Element('weight_windows_exodus') + + subelement = ET.SubElement(element, 'file') + subelement.text = str(self.file) + + subelement = ET.SubElement(element, 'adjoint_flux_variables') + subelement.text = ' '.join(self.adjoint_flux_variables) + + subelement = ET.SubElement(element, 'energy_bounds') + subelement.text = ' '.join(str(e) for e in self.energy_bounds) + + if self.timestep is not None: + subelement = ET.SubElement(element, 'timestep') + subelement.text = str(self.timestep) + + subelement = ET.SubElement(element, 'particle_type') + subelement.text = str(self.particle_type) + + # optional values are omitted so that the transport code defaults apply + if self.survival_ratio is not None: + subelement = ET.SubElement(element, 'survival_ratio') + subelement.text = str(self.survival_ratio) + + if self.upper_bound_ratio is not None: + subelement = ET.SubElement(element, 'upper_bound_ratio') + subelement.text = str(self.upper_bound_ratio) + + if self.max_split is not None: + subelement = ET.SubElement(element, 'max_split') + subelement.text = str(self.max_split) + + clean_indentation(element) + + return element + + @classmethod + def from_xml_element(cls, elem: ET.Element) -> Self: + """Create a WeightWindowsExodus object from an XML element + + Parameters + ---------- + elem : lxml.etree._Element + XML element + + Returns + ------- + openmc.WeightWindowsExodus + """ + file = get_text(elem, 'file') + variables = get_elem_list(elem, 'adjoint_flux_variables', str) + energy_bounds = get_elem_list(elem, 'energy_bounds', float) + + wwe = cls(file, variables, energy_bounds) + + timestep = get_text(elem, 'timestep') + if timestep is not None: + wwe.timestep = int(timestep) + + particle_type = get_text(elem, 'particle_type') + if particle_type is not None: + wwe.particle_type = particle_type + + survival_ratio = get_text(elem, 'survival_ratio') + if survival_ratio is not None: + wwe.survival_ratio = float(survival_ratio) + + upper_bound_ratio = get_text(elem, 'upper_bound_ratio') + if upper_bound_ratio is not None: + wwe.upper_bound_ratio = float(upper_bound_ratio) + + max_split = get_text(elem, 'max_split') + if max_split is not None: + wwe.max_split = int(max_split) + + wwe._check_consistency() + return wwe + + def hdf5_to_wws(path='weight_windows.h5') -> WeightWindowsList: """Create a WeightWindowsList from a weight windows HDF5 file diff --git a/src/mesh.cpp b/src/mesh.cpp index 181af846694..4d8907583f6 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -3587,6 +3587,23 @@ LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier) initialize(); } +// create the mesh from an externally constructed libMesh mesh, transferring +// ownership to OpenMC +LibMesh::LibMesh(unique_ptr input_mesh, + double length_multiplier, const std::string& filename) +{ + if (!input_mesh->is_replicated()) { + fatal_error("At present LibMesh tallies require a replicated mesh. Please " + "ensure 'input_mesh' is a libMesh::ReplicatedMesh."); + } + + unique_m_ = std::move(input_mesh); + m_ = unique_m_.get(); + filename_ = filename; + set_length_multiplier(length_multiplier); + initialize(); +} + // create the mesh from an input file LibMesh::LibMesh(const std::string& filename, double length_multiplier) { diff --git a/src/settings.cpp b/src/settings.cpp index 8ae252ae1ab..e2d9b294cd0 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -1240,6 +1240,11 @@ void read_settings_xml(pugi::xml_node root) std::make_unique(node_ww)); } + // Weight windows built from adjoint flux in an Exodus file (libMesh) + if (check_for_node(root, "weight_windows_exodus")) { + read_weight_windows_exodus(root.child("weight_windows_exodus")); + } + // Enable weight windows by default if one or more are present if (variance_reduction::weight_windows.size() > 0) settings::weight_windows_on = true; diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 0614110cd32..7474dcf8b26 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -30,6 +30,17 @@ #include +#ifdef OPENMC_LIBMESH_ENABLED +#include "libmesh/dof_map.h" +#include "libmesh/elem.h" +#include "libmesh/equation_systems.h" +#include "libmesh/exodusII_io.h" +#include "libmesh/explicit_system.h" +#include "libmesh/mesh_communication.h" +#include "libmesh/numeric_vector.h" +#include "libmesh/replicated_mesh.h" +#endif + namespace openmc { //============================================================================== @@ -929,6 +940,280 @@ void WeightWindowsGenerator::update() const // Non-member functions //============================================================================== +//! Compute FW-CADIS weight window bounds from multigroup adjoint flux +// +//! Mirrors the FW_CADIS branch of WeightWindows::update_weights(): positive +//! flux values are inverted and normalized by twice the global maximum of the +//! inverted values. Elements with non-positive flux keep the sentinel -1.0. +//! \param[in] flux flux[g][e]: adjoint flux for group g, element e +//! \param[in] upper_bound_ratio ratio of upper to lower ww bounds +//! \param[out] flat_lower lower bounds, flat layout [g * n_elem + e] +//! \param[out] flat_upper upper bounds, same layout +//! \return false if no positive flux value exists anywhere +static bool fw_cadis_bounds(const vector>& flux, + double upper_bound_ratio, vector& flat_lower, + vector& flat_upper) +{ + const size_t n_groups = flux.size(); + const size_t n_elem = n_groups ? flux[0].size() : 0; + + flat_lower.assign(n_groups * n_elem, -1.0); + flat_upper.assign(n_groups * n_elem, -1.0); + + // Invert positive flux values and track the global maximum + double inv_max = 0.0; + for (size_t g = 0; g < n_groups; ++g) { + for (size_t e = 0; e < n_elem; ++e) { + if (flux[g][e] > 0.0) { + double inv = 1.0 / flux[g][e]; + flat_lower[g * n_elem + e] = inv; + inv_max = std::max(inv_max, inv); + } + } + } + + if (inv_max <= 0.0) + return false; + + const double norm_factor = 1.0 / (2.0 * inv_max); + for (size_t i = 0; i < n_groups * n_elem; ++i) { + if (flat_lower[i] >= 0.0) { + flat_lower[i] *= norm_factor; + flat_upper[i] = flat_lower[i] * upper_bound_ratio; + } + } + return true; +} + +void read_weight_windows_exodus(pugi::xml_node node) +{ +#ifndef OPENMC_LIBMESH_ENABLED + (void)node; + fatal_error(" requires OpenMC to be compiled " + "with libMesh support (-DOPENMC_USE_LIBMESH=on)."); +#else + // Make sure required elements are present + const vector required_elems { + "file", "adjoint_flux_variables", "energy_bounds"}; + for (const auto& elem : required_elems) { + if (!check_for_node(node, elem.c_str())) { + fatal_error(fmt::format( + "Must specify <{}> for .", elem)); + } + } + + const std::string file = get_node_value(node, "file", true); + if (!file_exists(file)) + fatal_error(fmt::format( + ": mesh file '{}' does not exist.", file)); + + // One elemental variable per energy group, ordered by ascending energy + // consistently with + const vector flux_vars = + get_node_array(node, "adjoint_flux_variables"); + if (flux_vars.empty()) + fatal_error(": must " + "list at least one variable."); + const int n_groups = static_cast(flux_vars.size()); + + const vector e_bounds = + get_node_array(node, "energy_bounds"); + if (static_cast(e_bounds.size()) != n_groups + 1) + fatal_error(fmt::format( + ": must have exactly {} values " + "for {} group(s), but {} were provided.", + n_groups + 1, n_groups, e_bounds.size())); + for (int g = 0; g < n_groups; ++g) { + if (e_bounds[g] >= e_bounds[g + 1]) + fatal_error(fmt::format( + ": must be strictly " + "increasing; bounds[{}] = {} >= bounds[{}] = {}.", + g, e_bounds[g], g + 1, e_bounds[g + 1])); + } + + // is 0-based; default -1 selects the last step in the file + const int ts_user = check_for_node(node, "timestep") + ? std::stoi(get_node_value(node, "timestep", true)) + : -1; + + const std::string p_type_str = check_for_node(node, "particle_type") + ? get_node_value(node, "particle_type", true) + : "neutron"; + + const double survival_ratio = + check_for_node(node, "survival_ratio") + ? std::stod(get_node_value(node, "survival_ratio", true)) + : 3.0; + if (survival_ratio <= 1) + fatal_error("Survival to lower weight window ratio must bigger than 1 " + "and less than the upper to lower weight window ratio."); + + const double upper_bound_ratio = + check_for_node(node, "upper_bound_ratio") + ? std::stod(get_node_value(node, "upper_bound_ratio", true)) + : 5.0; + if (upper_bound_ratio <= survival_ratio) + fatal_error(fmt::format( + ": ({}) must be larger " + "than ({}).", + upper_bound_ratio, survival_ratio)); + + const int max_split = + check_for_node(node, "max_split") + ? std::stoi(get_node_value(node, "max_split", true)) + : 10; + if (max_split <= 1) + fatal_error("max split must be larger than 1"); + + // Read the mesh and all group flux variables in a single pass. Note that + // copy_elemental_solution() must be called on the same ExodusII_IO object + // that performed read(), and allow_renumbering(false) must be set before + // read() so that element IDs match the Exodus element block entries. + if (!settings::libmesh_comm) + fatal_error(": no libMesh communicator is " + "initialized."); + + auto mesh = make_unique(*settings::libmesh_comm, 3); + mesh->allow_renumbering(false); + + libMesh::ExodusII_IO exo(*mesh); + exo.read(file); + + // The reader only populates rank 0, so replicate the mesh to the other MPI + // ranks before use (no-op in serial) + libMesh::MeshCommunication().broadcast(*mesh); + mesh->prepare_for_use(); + + const int n_elem = static_cast(mesh->n_active_elem()); + if (n_elem == 0) + fatal_error(fmt::format( + ": mesh file '{}' has no elements.", file)); + + // Resolve the requested time step (Exodus steps are 1-based internally). + // The file is only open on rank 0, so query metadata there and broadcast. + const auto& comm = mesh->comm(); + int n_steps = 0; + if (comm.rank() == 0) + n_steps = static_cast(exo.get_time_steps().size()); + comm.broadcast(n_steps); + + const int ts_1based = (ts_user < 0) ? n_steps : (ts_user + 1); + if (ts_1based < 1 || ts_1based > n_steps) + fatal_error(fmt::format( + ": requested timestep {} is out of range " + "[0, {}) for file '{}'.", + (ts_user < 0 ? n_steps - 1 : ts_user), n_steps, file)); + + // Verify every requested variable exists before reading any of them + if (comm.rank() == 0) { + const auto& exo_elem_vars = exo.get_elem_var_names(); + for (const auto& vname : flux_vars) { + if (std::find(exo_elem_vars.begin(), exo_elem_vars.end(), vname) == + exo_elem_vars.end()) { + std::string available; + for (size_t vi = 0; vi < exo_elem_vars.size(); ++vi) { + if (vi) + available += ", "; + available += exo_elem_vars[vi]; + } + fatal_error(fmt::format( + ": variable '{}' not found in '{}'.\n" + " Available element variables: [{}]", + vname, file, available)); + } + } + } + + // Index flux arrays by elem->id() - first_id so that the flux index matches + // the mesh bin computed by LibMesh::get_bin_from_element() + const auto first_id = (*mesh->elements_begin())->id(); + + // flux[g][e] matches the (n_energy_bins, n_mesh_bins) layout of lower_ww_ + vector> flux(n_groups); + + for (int g = 0; g < n_groups; ++g) { + // Use a fresh EquationSystems per group to avoid DOF conflicts from + // multiple active variables + libMesh::EquationSystems eq_sys(*mesh); + auto& sys = eq_sys.add_system("adjoint_ww"); + sys.add_variable(flux_vars[g], libMesh::CONSTANT, libMesh::MONOMIAL); + eq_sys.init(); + + exo.copy_elemental_solution(sys, flux_vars[g], flux_vars[g], ts_1based); + + // Under MPI the solution vector is distributed; gather the full vector + // onto every rank (collective, no-op in serial) + std::vector soln_local; + sys.solution->localize(soln_local); + + const libMesh::DofMap& dof_map = sys.get_dof_map(); + flux[g].assign(n_elem, 0.0); + for (const auto* elem : mesh->active_element_ptr_range()) { + std::vector dofs; + dof_map.dof_indices(elem, dofs); + if (dofs.size() != 1) + fatal_error(fmt::format( + ": expected one DOF per element but found " + "{} for element {}.", + dofs.size(), elem->id())); + const auto bin = elem->id() - first_id; + if (bin >= static_cast(n_elem)) + fatal_error(fmt::format( + ": element IDs in '{}' are not contiguous " + "(element {} with first ID {}).", + file, elem->id(), first_id)); + flux[g][bin] = soln_local[dofs[0]]; + } + } + + // Register the mesh with OpenMC, transferring ownership + int32_t mesh_id = 1; + for (const auto& m : model::meshes) + mesh_id = std::max(mesh_id, m->id_ + 1); + + model::meshes.push_back( + make_unique(std::move(mesh), 1.0, file)); + model::meshes.back()->set_id(mesh_id); + + // Normalize (FW-CADIS) and build the WeightWindows object + vector flat_lower; + vector flat_upper; + if (!fw_cadis_bounds(flux, upper_bound_ratio, flat_lower, flat_upper)) + fatal_error(fmt::format( + ": all adjoint flux values across all {} " + "group(s) in '{}' are zero or negative -- cannot compute FW-CADIS " + "weight windows.", + n_groups, file)); + + // set_mesh() and set_energy_bounds() must precede set_bounds() since both + // trigger allocate_ww_bounds() + WeightWindows* wws = WeightWindows::create(); + wws->set_mesh(model::mesh_map.at(mesh_id)); + wws->set_particle_type(ParticleType {p_type_str}); + wws->set_energy_bounds( + span(e_bounds.data(), e_bounds.size())); + wws->survival_ratio() = survival_ratio; + wws->max_split() = max_split; + wws->set_bounds( + span(flat_lower.data(), flat_lower.size()), + span(flat_upper.data(), flat_upper.size())); + + std::string varlist; + for (int g = 0; g < n_groups; ++g) { + if (g) + varlist += ", "; + varlist += flux_vars[g]; + } + write_message( + fmt::format("Loaded {}-group adjoint weight windows from '{}':\n" + " {} elements, variables [{}], timestep {}, " + "upper_bound_ratio {:g}.", + n_groups, file, n_elem, varlist, + (ts_user < 0 ? n_steps - 1 : ts_user), upper_bound_ratio), + 5); +#endif // OPENMC_LIBMESH_ENABLED +} + std::pair search_weight_window(const Particle& p) { // TODO: this is a linear search - should do something more clever @@ -1399,4 +1684,4 @@ extern "C" int openmc_weight_windows_import(const char* filename) return 0; } -} // namespace openmc +} // namespace openmc \ No newline at end of file