From c2e9f89aad7c9a618658440aacd97f445af67d2a Mon Sep 17 00:00:00 2001 From: Foraejee Date: Tue, 30 Jun 2026 11:44:44 -0600 Subject: [PATCH 1/4] read 1 group adjoint from exodus using libmesh -no MPI --- src/settings.cpp | 358 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 358 insertions(+) diff --git a/src/settings.cpp b/src/settings.cpp index 8ae252ae1ab..d1883b4f1fa 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -36,6 +36,28 @@ #include "openmc/weight_windows.h" #include "openmc/xml_interface.h" +// ───────────────────────────────────────────────────────────────────────────── +// BEGIN INSERT A ─ new #include directives +// ───────────────────────────────────────────────────────────────────────────── + +#ifdef OPENMC_LIBMESH_ENABLED +// libMesh headers needed only by read_weight_windows_from_exodus(). +// All libMesh symbols are confined to this guard so non-libMesh builds +// compile without change. +#include "libmesh/replicated_mesh.h" // ReplicatedMesh – reads .exo files +#include "libmesh/exodusII_io.h" // ExodusII_IO::read + read_elemental_var +#include "libmesh/equation_systems.h" // required by ExodusII_IO for var reading +#include "libmesh/explicit_system.h" // ExplicitSystem – lightweight container +#include "libmesh/numeric_vector.h" // NumericVector used by ExplicitSystem +#include "libmesh/dof_map.h" // DoFMap for per-element DOF lookup +#endif // OPENMC_LIBMESH_ENABLED + +// ───────────────────────────────────────────────────────────────────────────── +// END INSERT A +// ───────────────────────────────────────────────────────────────────────────── + + + namespace openmc { //============================================================================== @@ -381,6 +403,335 @@ void get_run_parameters(pugi::xml_node node_base) } } + +// ----------------------------------------------------------------------------- +// read_weight_windows_from_exodus +// +// Called by read_settings_xml() when a node is +// present in settings.xml. The function: +// 1. Parses XML fields. +// 2. Loads the Exodus mesh via openmc::LibMesh; registers it in model::meshes. +// 3. Reads the named elemental flux variable from the same Exodus file. +// 4. Normalises the flux → lower_ww_bounds; scales → upper_ww_bounds. +// 5. Synthesises an in-memory pugi node and delegates +// construction to the existing WeightWindows::from_xml() path so that +// all ID management and registration logic is reused. +// 6. Sets settings::weight_windows_on = true. +// ----------------------------------------------------------------------------- +static void read_weight_windows_from_exodus(pugi::xml_node node) +{ +#ifndef OPENMC_LIBMESH_ENABLED + // Hard error: the stanza was present but libMesh was not compiled in. + // Failing loudly prevents silent no-op behaviour. + fatal_error(" requires OpenMC to be compiled " + "with libMesh support (-DOPENMC_USE_LIBMESH=ON)."); +#else + + // ── Step 1: Parse XML fields ────────────────────────────────────────────── + + // – path to the Exodus (.exo) file produced by the adjoint + // solver. Required; fatal_error if absent or empty. + const std::string mesh_file = + get_node_value(node, "mesh_file", /*strip_whitespace=*/true); + if (mesh_file.empty()) + fatal_error(": is required."); + + // – the element-variable name in the Exodus file + // that holds the adjoint scalar flux. Case-sensitive. + const std::string flux_var = + get_node_value(node, "adjoint_flux_variable", /*strip=*/true); + if (flux_var.empty()) + fatal_error(": " + " is required."); + + // – at least two space-separated energy values in eV. + // For a single energy group: two values (E_lo E_hi). + const std::vector energy_bounds = + get_node_array(node, "energy_bounds"); + if (energy_bounds.size() < 2) + fatal_error(": " + " must contain at least two values."); + + // – 0-based Exodus time step index. Sentinel -1 means "last". + // An absent element defaults to the last available time step. + const int ts_user = check_for_node(node, "timestep") + ? std::stoi(get_node_value(node, "timestep", true)) + : -1; // -1 → resolve to last step after opening the file + + // – weight-window roulette ratio (default 3.0). + const double survival_ratio = check_for_node(node, "survival_ratio") + ? std::stod(get_node_value(node, "survival_ratio", true)) + : 3.0; + + // – upper_ww = lower_ww * this factor (default 5.0). + const double upper_bound_ratio = check_for_node(node, "upper_bound_ratio") + ? std::stod(get_node_value(node, "upper_bound_ratio", true)) + : 5.0; + + // ── Step 2: Load the Exodus mesh; register as openmc::LibMesh ──────────── + // + // openmc::LibMesh (src/mesh.cpp) wraps a libMesh::ReplicatedMesh, builds a + // PointLocator, and is the canonical unstructured-mesh type for OpenMC. We + // construct it exactly as read_meshes() does for + // nodes, so the mesh is a proper first-class citizen: it appears in the + // statepoint, can be used as a tally filter, and its point-locator is ready + // for particle tracking. + + // Verify the file exists early; the libMesh error for a missing file is + // cryptic, so we give a cleaner message. + if (!file_exists(mesh_file)) + fatal_error(fmt::format( + ": mesh file '{}' does not exist.", + mesh_file)); + + // ── Step 2a: Read flux variable using a standalone ReplicatedMesh ───────── + // + // openmc::LibMesh does NOT expose its internal libMesh::MeshBase (there is no + // libmesh_mesh() accessor on that class). The clean solution is to open the + // Exodus file once in a fully independent libMesh::ReplicatedMesh that we own, + // extract the adjoint flux values from it, then discard it. After that we + // pass the same file path to openmc::LibMesh, which re-reads it for transport. + // The two ReplicatedMesh instances are independent objects that both represent + // the same Exodus geometry; the element traversal order is identical between + // them, which preserves the bin-index correspondence (see ordering guarantee + // below). + // + // Ordering guarantee + // ────────────────── + // ExodusII_IO::copy_elemental_solution() populates DOFs by iterating + // active_element_ptr_range() in ascending element-ID order on a + // ReplicatedMesh. openmc::LibMesh::get_bin() identifies bins by the same + // traversal (it stores the first active element ID and maps + // bin = elem->id() - first_elem_id). + // Because both meshes are loaded from the same Exodus file without any + // renumbering, their element IDs are identical and the traversal order matches. + // Therefore DOF index k from our standalone mesh == weight-window bin k in + // OpenMC, with no re-ordering step. + + std::vector flux; // filled below; size = n_active_elements + int n_elem = 0; + int n_steps_saved = 1; // set inside the scope block below + { + // Standalone ReplicatedMesh — scoped so it is destroyed before the + // openmc::LibMesh object is created, freeing memory early. + // ReplicatedMesh requires a libMesh::Parallel::Communicator reference. + // settings::libmesh_comm is a const libMesh::Parallel::Communicator* set + // during initialize.cpp and used by all mesh construction in openmc. + // This is identical to how openmc::LibMesh constructs its own mesh. + + libMesh::ReplicatedMesh standalone_mesh(*settings::libmesh_comm); + + // allow_renumbering(false) MUST be called before read(). + // copy_elemental_solution maps Exodus element-block entries to DOFs by + // element ID. If the mesh is renumbered those IDs change and the mapping + // is wrong (segfault or silent wrong values). Setting this flag before + // read() prevents renumbering during both read() and prepare_for_use(). + standalone_mesh.allow_renumbering(false); + + // Use a single ExodusII_IO object for both the mesh read and the later + // copy_elemental_solution call. Constructing a second ExodusII_IO and + // calling read() again on an already-populated mesh causes a segfault + // because the internal Exodus file handle and element maps are rebuilt + // inconsistently. One object, one read, one copy — that is the correct + // libMesh pattern. + libMesh::ExodusII_IO exo_reader(standalone_mesh); + exo_reader.read(mesh_file); + standalone_mesh.prepare_for_use(); + + // Count active elements and verify the mesh is non-empty. + n_elem = static_cast(standalone_mesh.n_active_elem()); + if (n_elem == 0) + fatal_error(fmt::format( + ": mesh file '{}' has no elements.", + mesh_file)); + + // Attach a throw-away EquationSystems / ExplicitSystem so that + // ExodusII_IO::copy_elemental_solution() has somewhere to write the data. + libMesh::EquationSystems eq_sys(standalone_mesh); + auto& sys = eq_sys.add_system("adjoint_ww"); + + // CONSTANT MONOMIAL: one scalar DOF per active element — matches how + // MOOSE/Griffin writes element-averaged scalar fluxes. + sys.add_variable(flux_var, libMesh::CONSTANT, libMesh::MONOMIAL); + eq_sys.init(); // allocate the DOF vectors + + // Verify the variable is present in the file. + const auto& exo_elem_vars = exo_reader.get_elem_var_names(); + if (std::find(exo_elem_vars.begin(), exo_elem_vars.end(), flux_var) + == exo_elem_vars.end()) { + // Build comma-separated list of available variable names for the error msg. + std::string available_vars; + for (std::size_t vi = 0; vi < exo_elem_vars.size(); ++vi) { + if (vi) available_vars += ", "; + available_vars += exo_elem_vars[vi]; + } + fatal_error(fmt::format( + ": variable '{}' not found in '{}'.\n" + " Available element variables: [{}]", + flux_var, mesh_file, available_vars)); + } + + // Resolve time step. Exodus uses 1-based step indices internally. + const int n_steps = static_cast(exo_reader.get_time_steps().size()); + // ts_user is 0-based (-1 = last). Convert to 1-based for libMesh. + 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, mesh_file)); + + // Populate the ExplicitSystem solution with the adjoint flux values. + exo_reader.copy_elemental_solution(sys, flux_var, flux_var, ts_1based); + + // Extract per-element values in ascending element-ID order. + // CONSTANT MONOMIAL → exactly one DOF per element → dof_indices[0]. + const libMesh::DofMap& dof_map = sys.get_dof_map(); + flux.resize(n_elem, 0.0); + int bin = 0; + for (const auto* elem : standalone_mesh.active_element_ptr_range()) { + std::vector dofs; + dof_map.dof_indices(elem, dofs); + flux[bin++] = sys.solution->el(dofs[0]); + } + // Save n_steps for later use in write_message (outside this scope). + n_steps_saved = n_steps; + // standalone_mesh and eq_sys destruct here; memory is freed. + } + + // ── Step 2b: Register the Exodus mesh with OpenMC ───────────────────────── + // + // Now construct the openmc::LibMesh wrapper (which re-reads the same Exodus + // file and builds the PointLocator for transport). We assign it a fresh ID + // and register it in model::meshes exactly as read_meshes() does in mesh.cpp. + + // Choose a fresh mesh ID: one beyond the current maximum. + int mesh_id = 1; + for (const auto& m : model::meshes) + mesh_id = std::max(mesh_id, m->id_ + 1); + + // openmc::LibMesh constructors (from include/openmc/mesh.h line 994): + // LibMesh(const std::string& filename, double length_multiplier = 1.0) + // LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0) + // There is no 3-argument constructor; ID is assigned separately via set_id(). + auto lm_ptr = std::make_unique( + mesh_file, // Exodus file path + 1.0 // length_multiplier (cm → cm, no conversion needed) + ); + // Verify element count matches what we extracted in Step 2a. + if (static_cast(lm_ptr->n_bins()) != n_elem) + fatal_error(fmt::format( + ": element count mismatch between " + "standalone read ({}) and openmc::LibMesh ({}).", + n_elem, lm_ptr->n_bins())); + + // Register in model::meshes first (set_id searches model::meshes for 'this'), + // then assign the ID via set_id() — which also writes model::mesh_map. + // This exactly mirrors the Mesh::create() pattern in mesh.cpp. + model::meshes.push_back(std::move(lm_ptr)); + model::meshes.back()->set_id(mesh_id); + + // ── Step 4: FW-CADIS normalisation → lower_ww_bounds; scale → upper_ww ──── + // + // This exactly mirrors src/weight_windows.cpp WeightWindows::update_weights() + // for the FW_CADIS branch (lines ~850-890 of that file). + // + // FW-CADIS: weight windows are INVERSELY proportional to the adjoint flux. + // Step A invert: importance[e] = 1 / phi_adj[e] + // Step B normalize: lower_ww[e] = importance[e] / (2 * max(importance)) + // = min(phi_adj) / (2 * phi_adj[e]) + // upper_ww[e] = lower_ww[e] * upper_bound_ratio + // + // The factor of 2 in the denominator is OpenMC's convention: it centres the + // particle's nominal weight inside the window when upper_bound_ratio = 5 + // (geometric mean ≈ 2.24 × lower_ww), matching what update_weights() does. + // + // Elements where phi_adj <= 0 have no importance for the detector; OpenMC + // stores -1.0 as the sentinel meaning "no window here" (the same value that + // update_weights() writes for bins where sum <= 0 or rel_err > threshold). + // The transport kernel skips windows whose lower bound is < 0, so those + // elements are simply left as analog Monte Carlo. + + const double phi_max = *std::max_element(flux.begin(), flux.end()); + if (phi_max <= 0.0) + fatal_error(fmt::format( + ": max value of variable '{}' in '{}' " + "is {:g} <= 0. Check variable name and time step index.", + flux_var, mesh_file, phi_max)); + + // Step A: invert adjoint flux; record max of inverted values. + // Elements with phi_adj <= 0 are flagged with inv = -1 (sentinel). + std::vector inv(n_elem); + double inv_max = 0.0; + for (int i = 0; i < n_elem; ++i) { + if (flux[i] > 0.0) { + inv[i] = 1.0 / flux[i]; // importance ∝ 1/phi_adj + if (inv[i] > inv_max) + inv_max = inv[i]; + } else { + inv[i] = -1.0; // sentinel: no window in this element + } + } + + if (inv_max <= 0.0) + fatal_error(fmt::format( + ": all values of variable '{}' in '{}' " + "are zero or negative — cannot compute FW-CADIS weight windows.", + flux_var, mesh_file)); + + // Step B: normalise by (2 * inv_max), matching OpenMC's update_weights(). + const double norm_factor = 1.0 / (2.0 * inv_max); + + std::vector lower_ww(n_elem), upper_ww(n_elem); + for (int i = 0; i < n_elem; ++i) { + if (inv[i] < 0.0) { + // Zero-adjoint-flux element: no window (OpenMC sentinel -1). + lower_ww[i] = -1.0; + upper_ww[i] = -1.0; + } else { + lower_ww[i] = inv[i] * norm_factor; // = 1/phi_adj / (2*max(1/phi_adj)) + upper_ww[i] = lower_ww[i] * upper_bound_ratio; // upper bound + } + } + + // ── Step 5: Build the WeightWindows object ──────────────────────────────── + // + // From include/openmc/weight_windows.h (verified at commit fd1bc26a): + // static WeightWindows* create(int32_t id = -1) + // → pushes to variance_reduction::weight_windows, sets index_, registers + // in ww_map, returns raw ptr + // double& survival_ratio() → non-const ref accessor, writable + // void set_mesh(int32_t mesh_idx) + // void set_particle_type(ParticleType) + // void set_energy_bounds(span) + // void set_bounds(span lower, span upper) + { + WeightWindows* wws = WeightWindows::create(); + wws->set_mesh(model::mesh_map.at(mesh_id)); + wws->set_particle_type(ParticleType {"neutron"}); + wws->set_energy_bounds( + span(energy_bounds.data(), energy_bounds.size())); + // survival_ratio() returns a non-const reference (weight_windows.h line ~160) + wws->survival_ratio() = survival_ratio; + wws->set_bounds( + span(lower_ww.data(), lower_ww.size()), + span(upper_ww.data(), upper_ww.size())); + } + + // ── Step 6: Enable weight windows globally ──────────────────────────────── + // read_settings_xml() sets this flag inside the normal loop; + // our path bypasses that loop, so we set it explicitly here. + settings::weight_windows_on = true; + + write_message(fmt::format( + "Loaded adjoint weight windows from Exodus file '{}':\n" + " {} elements, variable '{}', timestep {}, upper_bound_ratio {:g}.", + mesh_file, n_elem, flux_var, + (ts_user < 0 ? n_steps_saved - 1 : ts_user), upper_bound_ratio), 5); + +#endif // OPENMC_LIBMESH_ENABLED +} + void read_settings_xml() { using namespace settings; @@ -1240,6 +1591,13 @@ void read_settings_xml(pugi::xml_node root) std::make_unique(node_ww)); } + // ──────────────────────────────────────────────────────── + // If the user provided a stanza, load the + // adjoint flux from the Exodus file and build WeightWindows in memory. + if (auto wwe = root.child("weight_windows_from_exodus")) + read_weight_windows_from_exodus(wwe); + // ──────────────────────────────────────────────────────────────────────── + // Enable weight windows by default if one or more are present if (variance_reduction::weight_windows.size() > 0) settings::weight_windows_on = true; From a5c4874c170348a23bf39496fcba74243c2a2e88 Mon Sep 17 00:00:00 2001 From: Foraejee Date: Mon, 13 Jul 2026 09:22:02 -0600 Subject: [PATCH 2/4] multigroup adjoint yet to be tested --- src/settings.cpp | 401 +++++++++++++++++++++-------------------------- 1 file changed, 181 insertions(+), 220 deletions(-) diff --git a/src/settings.cpp b/src/settings.cpp index d1883b4f1fa..49df9324e68 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -404,334 +404,295 @@ void get_run_parameters(pugi::xml_node node_base) } -// ----------------------------------------------------------------------------- + // ----------------------------------------------------------------------------- // read_weight_windows_from_exodus // -// Called by read_settings_xml() when a node is -// present in settings.xml. The function: -// 1. Parses XML fields. -// 2. Loads the Exodus mesh via openmc::LibMesh; registers it in model::meshes. -// 3. Reads the named elemental flux variable from the same Exodus file. -// 4. Normalises the flux → lower_ww_bounds; scales → upper_ww_bounds. -// 5. Synthesises an in-memory pugi node and delegates -// construction to the existing WeightWindows::from_xml() path so that -// all ID management and registration logic is reused. -// 6. Sets settings::weight_windows_on = true. +// Builds a WeightWindows object from adjoint flux variables stored as elemental +// data in an Exodus II file. Single-group and multi-group are both supported. +// +// Tensor layout (verified from weight_windows.cpp bounds_size() and set_bounds): +// lower_ww_ and upper_ww_ are shaped (n_energy_bins, n_mesh_bins). +// Energy is the OUTER (slow) index; mesh element is the INNER (fast) index. +// set_bounds(span, span) expects a flat array in that order: +// [ group0_elem0, group0_elem1, ..., group1_elem0, group1_elem1, ... ] +// +// FW-CADIS normalisation (verified from update_weights() FW_CADIS branch): +// 1. Invert all values across all groups: importance = 1/phi_adj +// 2. Find the GLOBAL maximum of the inverted values (across all groups) +// 3. Normalise everything by 1/(2*global_max) +// This single normalisation across groups is what the source code does — +// not per-group normalisation (that is the MAGIC method). // ----------------------------------------------------------------------------- static void read_weight_windows_from_exodus(pugi::xml_node node) { #ifndef OPENMC_LIBMESH_ENABLED - // Hard error: the stanza was present but libMesh was not compiled in. - // Failing loudly prevents silent no-op behaviour. fatal_error(" requires OpenMC to be compiled " "with libMesh support (-DOPENMC_USE_LIBMESH=ON)."); #else // ── Step 1: Parse XML fields ────────────────────────────────────────────── - // – path to the Exodus (.exo) file produced by the adjoint - // solver. Required; fatal_error if absent or empty. const std::string mesh_file = - get_node_value(node, "mesh_file", /*strip_whitespace=*/true); + get_node_value(node, "mesh_file", /*strip=*/true); if (mesh_file.empty()) fatal_error(": is required."); - // – the element-variable name in the Exodus file - // that holds the adjoint scalar flux. Case-sensitive. - const std::string flux_var = - get_node_value(node, "adjoint_flux_variable", /*strip=*/true); - if (flux_var.empty()) + if (!file_exists(mesh_file)) + fatal_error(fmt::format( + ": mesh file '{}' does not exist.", + mesh_file)); + + // – space-separated list of Exodus element-variable + // names, one per energy group, ordered from group 0 (highest energy) to + // group N-1 (lowest energy), matching the order of . + const std::vector flux_vars = + get_node_array(node, "adjoint_flux_variables"); + if (flux_vars.empty()) fatal_error(": " - " is required."); + " must list at least one variable."); + + const int n_groups = static_cast(flux_vars.size()); - // – at least two space-separated energy values in eV. - // For a single energy group: two values (E_lo E_hi). + // – must have exactly n_groups + 1 values. const std::vector energy_bounds = get_node_array(node, "energy_bounds"); - if (energy_bounds.size() < 2) - fatal_error(": " - " must contain at least two values."); + if (static_cast(energy_bounds.size()) != n_groups + 1) + fatal_error(fmt::format( + ": must have exactly " + "{} values for {} group(s), but {} were provided.", + n_groups + 1, n_groups, energy_bounds.size())); - // – 0-based Exodus time step index. Sentinel -1 means "last". - // An absent element defaults to the last available time step. + // – 0-based; -1 means last step. const int ts_user = check_for_node(node, "timestep") ? std::stoi(get_node_value(node, "timestep", true)) - : -1; // -1 → resolve to last step after opening the file + : -1; - // – weight-window roulette ratio (default 3.0). const double survival_ratio = check_for_node(node, "survival_ratio") ? std::stod(get_node_value(node, "survival_ratio", true)) : 3.0; - // – upper_ww = lower_ww * this factor (default 5.0). const double upper_bound_ratio = check_for_node(node, "upper_bound_ratio") ? std::stod(get_node_value(node, "upper_bound_ratio", true)) : 5.0; - // ── Step 2: Load the Exodus mesh; register as openmc::LibMesh ──────────── + const double max_split = check_for_node(node, "max_split") + ? std::stod(get_node_value(node, "max_split", true)) + : 10; + + // ── Step 2a: Read all group flux variables from the Exodus file ─────────── // - // openmc::LibMesh (src/mesh.cpp) wraps a libMesh::ReplicatedMesh, builds a - // PointLocator, and is the canonical unstructured-mesh type for OpenMC. We - // construct it exactly as read_meshes() does for - // nodes, so the mesh is a proper first-class citizen: it appears in the - // statepoint, can be used as a tally filter, and its point-locator is ready - // for particle tracking. - - // Verify the file exists early; the libMesh error for a missing file is - // cryptic, so we give a cleaner message. - if (!file_exists(mesh_file)) - fatal_error(fmt::format( - ": mesh file '{}' does not exist.", - mesh_file)); - - // ── Step 2a: Read flux variable using a standalone ReplicatedMesh ───────── + // Use a single ExodusII_IO object for both read() and copy_elemental_solution(). + // Constructing a second ExodusII_IO on an already-populated mesh and calling + // read() again causes a segfault (internal element maps are rebuilt + // inconsistently). // - // openmc::LibMesh does NOT expose its internal libMesh::MeshBase (there is no - // libmesh_mesh() accessor on that class). The clean solution is to open the - // Exodus file once in a fully independent libMesh::ReplicatedMesh that we own, - // extract the adjoint flux values from it, then discard it. After that we - // pass the same file path to openmc::LibMesh, which re-reads it for transport. - // The two ReplicatedMesh instances are independent objects that both represent - // the same Exodus geometry; the element traversal order is identical between - // them, which preserves the bin-index correspondence (see ordering guarantee - // below). + // allow_renumbering(false) MUST be called before read(). copy_elemental_solution + // maps Exodus element-block entries to DOFs by element ID; renumbering changes + // those IDs and produces wrong values or a segfault. // - // Ordering guarantee - // ────────────────── - // ExodusII_IO::copy_elemental_solution() populates DOFs by iterating - // active_element_ptr_range() in ascending element-ID order on a - // ReplicatedMesh. openmc::LibMesh::get_bin() identifies bins by the same - // traversal (it stores the first active element ID and maps - // bin = elem->id() - first_elem_id). - // Because both meshes are loaded from the same Exodus file without any - // renumbering, their element IDs are identical and the traversal order matches. - // Therefore DOF index k from our standalone mesh == weight-window bin k in - // OpenMC, with no re-ordering step. + // flux[g][e] = adjoint flux for group g, element e. + // Outer index = group (energy), inner index = element (mesh bin). + // This matches the (n_energy_bins, n_mesh_bins) layout of lower_ww_. - std::vector flux; // filled below; size = n_active_elements int n_elem = 0; - int n_steps_saved = 1; // set inside the scope block below + int n_steps_saved = 1; + // flux[group][element] + std::vector> flux(n_groups); + { - // Standalone ReplicatedMesh — scoped so it is destroyed before the - // openmc::LibMesh object is created, freeing memory early. - // ReplicatedMesh requires a libMesh::Parallel::Communicator reference. - // settings::libmesh_comm is a const libMesh::Parallel::Communicator* set - // during initialize.cpp and used by all mesh construction in openmc. - // This is identical to how openmc::LibMesh constructs its own mesh. - libMesh::ReplicatedMesh standalone_mesh(*settings::libmesh_comm); + standalone_mesh.allow_renumbering(false); // must be before read() - // allow_renumbering(false) MUST be called before read(). - // copy_elemental_solution maps Exodus element-block entries to DOFs by - // element ID. If the mesh is renumbered those IDs change and the mapping - // is wrong (segfault or silent wrong values). Setting this flag before - // read() prevents renumbering during both read() and prepare_for_use(). - standalone_mesh.allow_renumbering(false); - - // Use a single ExodusII_IO object for both the mesh read and the later - // copy_elemental_solution call. Constructing a second ExodusII_IO and - // calling read() again on an already-populated mesh causes a segfault - // because the internal Exodus file handle and element maps are rebuilt - // inconsistently. One object, one read, one copy — that is the correct - // libMesh pattern. libMesh::ExodusII_IO exo_reader(standalone_mesh); exo_reader.read(mesh_file); standalone_mesh.prepare_for_use(); - // Count active elements and verify the mesh is non-empty. n_elem = static_cast(standalone_mesh.n_active_elem()); if (n_elem == 0) fatal_error(fmt::format( ": mesh file '{}' has no elements.", mesh_file)); - // Attach a throw-away EquationSystems / ExplicitSystem so that - // ExodusII_IO::copy_elemental_solution() has somewhere to write the data. - libMesh::EquationSystems eq_sys(standalone_mesh); - auto& sys = eq_sys.add_system("adjoint_ww"); - - // CONSTANT MONOMIAL: one scalar DOF per active element — matches how - // MOOSE/Griffin writes element-averaged scalar fluxes. - sys.add_variable(flux_var, libMesh::CONSTANT, libMesh::MONOMIAL); - eq_sys.init(); // allocate the DOF vectors - - // Verify the variable is present in the file. - const auto& exo_elem_vars = exo_reader.get_elem_var_names(); - if (std::find(exo_elem_vars.begin(), exo_elem_vars.end(), flux_var) - == exo_elem_vars.end()) { - // Build comma-separated list of available variable names for the error msg. - std::string available_vars; - for (std::size_t vi = 0; vi < exo_elem_vars.size(); ++vi) { - if (vi) available_vars += ", "; - available_vars += exo_elem_vars[vi]; - } - fatal_error(fmt::format( - ": variable '{}' not found in '{}'.\n" - " Available element variables: [{}]", - flux_var, mesh_file, available_vars)); - } - - // Resolve time step. Exodus uses 1-based step indices internally. + // Resolve time step (Exodus is 1-based internally). const int n_steps = static_cast(exo_reader.get_time_steps().size()); - // ts_user is 0-based (-1 = last). Convert to 1-based for libMesh. + n_steps_saved = 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 '{}'.", + ": requested timestep {} is out of " + "range [0, {}) for file '{}'.", (ts_user < 0 ? n_steps - 1 : ts_user), n_steps, mesh_file)); - // Populate the ExplicitSystem solution with the adjoint flux values. - exo_reader.copy_elemental_solution(sys, flux_var, flux_var, ts_1based); + // Verify all requested variable names exist in the file before reading any. + const auto& exo_elem_vars = exo_reader.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 (std::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, mesh_file, available)); + } + } - // Extract per-element values in ascending element-ID order. - // CONSTANT MONOMIAL → exactly one DOF per element → dof_indices[0]. - const libMesh::DofMap& dof_map = sys.get_dof_map(); - flux.resize(n_elem, 0.0); - int bin = 0; - for (const auto* elem : standalone_mesh.active_element_ptr_range()) { - std::vector dofs; - dof_map.dof_indices(elem, dofs); - flux[bin++] = sys.solution->el(dofs[0]); - } - // Save n_steps for later use in write_message (outside this scope). - n_steps_saved = n_steps; - // standalone_mesh and eq_sys destruct here; memory is freed. + // Read each group variable into its own EquationSystems instance. + // We reuse the same ExodusII_IO object (exo_reader) for all groups — + // copy_elemental_solution only needs the file handle that read() opened, + // not a fresh system. We reinitialise eq_sys between groups to avoid + // DOF conflicts from having multiple variables active simultaneously. + for (int g = 0; g < n_groups; ++g) { + libMesh::EquationSystems eq_sys(standalone_mesh); + auto& sys = eq_sys.add_system("adjoint_ww"); + // CONSTANT MONOMIAL: one scalar DOF per active element. + sys.add_variable(flux_vars[g], libMesh::CONSTANT, libMesh::MONOMIAL); + eq_sys.init(); + + exo_reader.copy_elemental_solution( + sys, flux_vars[g], flux_vars[g], ts_1based); + + const libMesh::DofMap& dof_map = sys.get_dof_map(); + flux[g].resize(n_elem, 0.0); + int bin = 0; + for (const auto* elem : standalone_mesh.active_element_ptr_range()) { + std::vector dofs; + dof_map.dof_indices(elem, dofs); + // CONSTANT MONOMIAL → exactly one DOF per element. + flux[g][bin++] = sys.solution->el(dofs[0]); + } + // eq_sys destructs here; frees the DOF vectors for this group. + } + // standalone_mesh destructs here. } // ── Step 2b: Register the Exodus mesh with OpenMC ───────────────────────── - // - // Now construct the openmc::LibMesh wrapper (which re-reads the same Exodus - // file and builds the PointLocator for transport). We assign it a fresh ID - // and register it in model::meshes exactly as read_meshes() does in mesh.cpp. - // Choose a fresh mesh ID: one beyond the current maximum. int mesh_id = 1; for (const auto& m : model::meshes) mesh_id = std::max(mesh_id, m->id_ + 1); - // openmc::LibMesh constructors (from include/openmc/mesh.h line 994): + // openmc::LibMesh constructors (include/openmc/mesh.h): // LibMesh(const std::string& filename, double length_multiplier = 1.0) - // LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0) - // There is no 3-argument constructor; ID is assigned separately via set_id(). - auto lm_ptr = std::make_unique( - mesh_file, // Exodus file path - 1.0 // length_multiplier (cm → cm, no conversion needed) - ); - // Verify element count matches what we extracted in Step 2a. + // No 3-arg constructor; ID is set via set_id() after push_back. + auto lm_ptr = std::make_unique(mesh_file, 1.0); + if (static_cast(lm_ptr->n_bins()) != n_elem) fatal_error(fmt::format( ": element count mismatch between " "standalone read ({}) and openmc::LibMesh ({}).", n_elem, lm_ptr->n_bins())); - // Register in model::meshes first (set_id searches model::meshes for 'this'), - // then assign the ID via set_id() — which also writes model::mesh_map. - // This exactly mirrors the Mesh::create() pattern in mesh.cpp. model::meshes.push_back(std::move(lm_ptr)); model::meshes.back()->set_id(mesh_id); - // ── Step 4: FW-CADIS normalisation → lower_ww_bounds; scale → upper_ww ──── + // ── Step 3: FW-CADIS normalisation ──────────────────────────────────────── + // + // From weight_windows.cpp update_weights() FW_CADIS branch (verified): + // 1. Invert all values across ALL groups simultaneously. + // 2. Find the GLOBAL maximum of the inverted values (single value for all + // groups — NOT per-group; per-group normalisation is the MAGIC method). + // 3. Normalise by 1/(2*global_max). + // 4. Elements with phi_adj <= 0 → sentinel -1.0 (no window). // - // This exactly mirrors src/weight_windows.cpp WeightWindows::update_weights() - // for the FW_CADIS branch (lines ~850-890 of that file). + // ── Energy group ordering ──────────────────────────────────────────────── // - // FW-CADIS: weight windows are INVERSELY proportional to the adjoint flux. - // Step A invert: importance[e] = 1 / phi_adj[e] - // Step B normalize: lower_ww[e] = importance[e] / (2 * max(importance)) - // = min(phi_adj) / (2 * phi_adj[e]) - // upper_ww[e] = lower_ww[e] * upper_bound_ratio + // The user lists in the SAME order as the energy + // intervals implied by : + // energy_bounds[0..1] → flux_vars[0] (lowest-energy group) + // energy_bounds[1..2] → flux_vars[1] + // ... + // energy_bounds[N-1..N] → flux_vars[N-1] (highest-energy group) // - // The factor of 2 in the denominator is OpenMC's convention: it centres the - // particle's nominal weight inside the window when upper_bound_ratio = 5 - // (geometric mean ≈ 2.24 × lower_ww), matching what update_weights() does. + // This matches OpenMC's internal storage directly: + // lower_ww_(energy_bin, mesh_bin) where energy_bin=0 = lowest energy + // set_bounds(span) copies flat[g * n_elem + e] → lower_ww_(g, e) // - // Elements where phi_adj <= 0 have no importance for the detector; OpenMC - // stores -1.0 as the sentinel meaning "no window here" (the same value that - // update_weights() writes for bins where sum <= 0 or rel_err > threshold). - // The transport kernel skips windows whose lower bound is < 0, so those - // elements are simply left as analog Monte Carlo. + // If your adjoint solver uses a different ordering (e.g. Griffin writes + // g0=fast), list the variable names in ascending-energy order in the XML, + // i.e. thermal group variable first. - const double phi_max = *std::max_element(flux.begin(), flux.end()); - if (phi_max <= 0.0) - fatal_error(fmt::format( - ": max value of variable '{}' in '{}' " - "is {:g} <= 0. Check variable name and time step index.", - flux_var, mesh_file, phi_max)); + // Step A: invert; find global max of inverted values across all groups. + std::vector flat_lower(n_groups * n_elem, -1.0); + std::vector flat_upper(n_groups * n_elem, -1.0); - // Step A: invert adjoint flux; record max of inverted values. - // Elements with phi_adj <= 0 are flagged with inv = -1 (sentinel). - std::vector inv(n_elem); double inv_max = 0.0; - for (int i = 0; i < n_elem; ++i) { - if (flux[i] > 0.0) { - inv[i] = 1.0 / flux[i]; // importance ∝ 1/phi_adj - if (inv[i] > inv_max) - inv_max = inv[i]; - } else { - inv[i] = -1.0; // sentinel: no window in this element + for (int g = 0; g < n_groups; ++g) { + for (int 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; // temporary; normalised below + if (inv > inv_max) inv_max = inv; + } + // else: remains -1.0 (sentinel) } } if (inv_max <= 0.0) fatal_error(fmt::format( - ": all values of variable '{}' in '{}' " - "are zero or negative — cannot compute FW-CADIS weight windows.", - flux_var, mesh_file)); + ": all adjoint flux values across all " + "{} group(s) in '{}' are zero or negative — cannot compute " + "FW-CADIS weight windows.", n_groups, mesh_file)); - // Step B: normalise by (2 * inv_max), matching OpenMC's update_weights(). + // Step B: normalise by global 1/(2*inv_max) and set upper bounds. const double norm_factor = 1.0 / (2.0 * inv_max); - - std::vector lower_ww(n_elem), upper_ww(n_elem); - for (int i = 0; i < n_elem; ++i) { - if (inv[i] < 0.0) { - // Zero-adjoint-flux element: no window (OpenMC sentinel -1). - lower_ww[i] = -1.0; - upper_ww[i] = -1.0; - } else { - lower_ww[i] = inv[i] * norm_factor; // = 1/phi_adj / (2*max(1/phi_adj)) - upper_ww[i] = lower_ww[i] * upper_bound_ratio; // upper bound + for (int 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; } + // else: both remain -1.0 (sentinel — transport kernel skips these) } - // ── Step 5: Build the WeightWindows object ──────────────────────────────── + // ── Step 4: Build the WeightWindows object ──────────────────────────────── + // + // Verified API from include/openmc/weight_windows.h at fd1bc26a: + // WeightWindows::create() → allocates, registers, returns ptr + // set_mesh(int32_t mesh_idx) → takes vector index + // set_particle_type(ParticleType) + // set_energy_bounds(span)→ also calls allocate_ww_bounds() + // set_bounds(span, span) + // double& survival_ratio() → non-const ref, writable // - // From include/openmc/weight_windows.h (verified at commit fd1bc26a): - // static WeightWindows* create(int32_t id = -1) - // → pushes to variance_reduction::weight_windows, sets index_, registers - // in ww_map, returns raw ptr - // double& survival_ratio() → non-const ref accessor, writable - // void set_mesh(int32_t mesh_idx) - // void set_particle_type(ParticleType) - // void set_energy_bounds(span) - // void set_bounds(span lower, span upper) + // Order matters: set_mesh() and set_energy_bounds() must be called before + // set_bounds() because both call allocate_ww_bounds() which sizes the tensors. + // set_bounds() then checks that span size == n_energy_bins * n_mesh_bins and + // copies the data in. { WeightWindows* wws = WeightWindows::create(); wws->set_mesh(model::mesh_map.at(mesh_id)); wws->set_particle_type(ParticleType {"neutron"}); wws->set_energy_bounds( span(energy_bounds.data(), energy_bounds.size())); - // survival_ratio() returns a non-const reference (weight_windows.h line ~160) wws->survival_ratio() = survival_ratio; + wws->max_split() = max_split; wws->set_bounds( - span(lower_ww.data(), lower_ww.size()), - span(upper_ww.data(), upper_ww.size())); + span(flat_lower.data(), flat_lower.size()), + span(flat_upper.data(), flat_upper.size())); } - // ── Step 6: Enable weight windows globally ──────────────────────────────── - // read_settings_xml() sets this flag inside the normal loop; - // our path bypasses that loop, so we set it explicitly here. + // ── Step 5: Enable weight windows globally ──────────────────────────────── settings::weight_windows_on = true; - write_message(fmt::format( - "Loaded adjoint weight windows from Exodus file '{}':\n" - " {} elements, variable '{}', timestep {}, upper_bound_ratio {:g}.", - mesh_file, n_elem, flux_var, - (ts_user < 0 ? n_steps_saved - 1 : ts_user), upper_bound_ratio), 5); + { + 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, mesh_file, n_elem, varlist, + (ts_user < 0 ? n_steps_saved - 1 : ts_user), upper_bound_ratio), 5); + } #endif // OPENMC_LIBMESH_ENABLED } - + void read_settings_xml() { using namespace settings; From b125f5f6ab782e09ce4aa4e9a3dc545878227118 Mon Sep 17 00:00:00 2001 From: Foraejee Date: Tue, 11 Aug 2026 13:08:02 -0600 Subject: [PATCH 3/4] read adjoint solution from exodus file --- include/openmc/mesh.h | 9 + include/openmc/weight_windows.h | 5 + src/mesh.cpp | 17 ++ src/settings.cpp | 322 +------------------------------- src/weight_windows.cpp | 287 +++++++++++++++++++++++++++- 5 files changed, 321 insertions(+), 319 deletions(-) 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/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 49df9324e68..e2d9b294cd0 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -36,28 +36,6 @@ #include "openmc/weight_windows.h" #include "openmc/xml_interface.h" -// ───────────────────────────────────────────────────────────────────────────── -// BEGIN INSERT A ─ new #include directives -// ───────────────────────────────────────────────────────────────────────────── - -#ifdef OPENMC_LIBMESH_ENABLED -// libMesh headers needed only by read_weight_windows_from_exodus(). -// All libMesh symbols are confined to this guard so non-libMesh builds -// compile without change. -#include "libmesh/replicated_mesh.h" // ReplicatedMesh – reads .exo files -#include "libmesh/exodusII_io.h" // ExodusII_IO::read + read_elemental_var -#include "libmesh/equation_systems.h" // required by ExodusII_IO for var reading -#include "libmesh/explicit_system.h" // ExplicitSystem – lightweight container -#include "libmesh/numeric_vector.h" // NumericVector used by ExplicitSystem -#include "libmesh/dof_map.h" // DoFMap for per-element DOF lookup -#endif // OPENMC_LIBMESH_ENABLED - -// ───────────────────────────────────────────────────────────────────────────── -// END INSERT A -// ───────────────────────────────────────────────────────────────────────────── - - - namespace openmc { //============================================================================== @@ -403,296 +381,6 @@ void get_run_parameters(pugi::xml_node node_base) } } - - // ----------------------------------------------------------------------------- -// read_weight_windows_from_exodus -// -// Builds a WeightWindows object from adjoint flux variables stored as elemental -// data in an Exodus II file. Single-group and multi-group are both supported. -// -// Tensor layout (verified from weight_windows.cpp bounds_size() and set_bounds): -// lower_ww_ and upper_ww_ are shaped (n_energy_bins, n_mesh_bins). -// Energy is the OUTER (slow) index; mesh element is the INNER (fast) index. -// set_bounds(span, span) expects a flat array in that order: -// [ group0_elem0, group0_elem1, ..., group1_elem0, group1_elem1, ... ] -// -// FW-CADIS normalisation (verified from update_weights() FW_CADIS branch): -// 1. Invert all values across all groups: importance = 1/phi_adj -// 2. Find the GLOBAL maximum of the inverted values (across all groups) -// 3. Normalise everything by 1/(2*global_max) -// This single normalisation across groups is what the source code does — -// not per-group normalisation (that is the MAGIC method). -// ----------------------------------------------------------------------------- -static void read_weight_windows_from_exodus(pugi::xml_node node) -{ -#ifndef OPENMC_LIBMESH_ENABLED - fatal_error(" requires OpenMC to be compiled " - "with libMesh support (-DOPENMC_USE_LIBMESH=ON)."); -#else - - // ── Step 1: Parse XML fields ────────────────────────────────────────────── - - const std::string mesh_file = - get_node_value(node, "mesh_file", /*strip=*/true); - if (mesh_file.empty()) - fatal_error(": is required."); - - if (!file_exists(mesh_file)) - fatal_error(fmt::format( - ": mesh file '{}' does not exist.", - mesh_file)); - - // – space-separated list of Exodus element-variable - // names, one per energy group, ordered from group 0 (highest energy) to - // group N-1 (lowest energy), matching the order of . - const std::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()); - - // – must have exactly n_groups + 1 values. - const std::vector energy_bounds = - get_node_array(node, "energy_bounds"); - if (static_cast(energy_bounds.size()) != n_groups + 1) - fatal_error(fmt::format( - ": must have exactly " - "{} values for {} group(s), but {} were provided.", - n_groups + 1, n_groups, energy_bounds.size())); - - // – 0-based; -1 means last step. - const int ts_user = check_for_node(node, "timestep") - ? std::stoi(get_node_value(node, "timestep", true)) - : -1; - - const double survival_ratio = check_for_node(node, "survival_ratio") - ? std::stod(get_node_value(node, "survival_ratio", true)) - : 3.0; - - const double upper_bound_ratio = check_for_node(node, "upper_bound_ratio") - ? std::stod(get_node_value(node, "upper_bound_ratio", true)) - : 5.0; - - const double max_split = check_for_node(node, "max_split") - ? std::stod(get_node_value(node, "max_split", true)) - : 10; - - // ── Step 2a: Read all group flux variables from the Exodus file ─────────── - // - // Use a single ExodusII_IO object for both read() and copy_elemental_solution(). - // Constructing a second ExodusII_IO on an already-populated mesh and calling - // read() again causes a segfault (internal element maps are rebuilt - // inconsistently). - // - // allow_renumbering(false) MUST be called before read(). copy_elemental_solution - // maps Exodus element-block entries to DOFs by element ID; renumbering changes - // those IDs and produces wrong values or a segfault. - // - // flux[g][e] = adjoint flux for group g, element e. - // Outer index = group (energy), inner index = element (mesh bin). - // This matches the (n_energy_bins, n_mesh_bins) layout of lower_ww_. - - int n_elem = 0; - int n_steps_saved = 1; - // flux[group][element] - std::vector> flux(n_groups); - - { - libMesh::ReplicatedMesh standalone_mesh(*settings::libmesh_comm); - standalone_mesh.allow_renumbering(false); // must be before read() - - libMesh::ExodusII_IO exo_reader(standalone_mesh); - exo_reader.read(mesh_file); - standalone_mesh.prepare_for_use(); - - n_elem = static_cast(standalone_mesh.n_active_elem()); - if (n_elem == 0) - fatal_error(fmt::format( - ": mesh file '{}' has no elements.", - mesh_file)); - - // Resolve time step (Exodus is 1-based internally). - const int n_steps = static_cast(exo_reader.get_time_steps().size()); - n_steps_saved = 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, mesh_file)); - - // Verify all requested variable names exist in the file before reading any. - const auto& exo_elem_vars = exo_reader.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 (std::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, mesh_file, available)); - } - } - - // Read each group variable into its own EquationSystems instance. - // We reuse the same ExodusII_IO object (exo_reader) for all groups — - // copy_elemental_solution only needs the file handle that read() opened, - // not a fresh system. We reinitialise eq_sys between groups to avoid - // DOF conflicts from having multiple variables active simultaneously. - for (int g = 0; g < n_groups; ++g) { - libMesh::EquationSystems eq_sys(standalone_mesh); - auto& sys = eq_sys.add_system("adjoint_ww"); - // CONSTANT MONOMIAL: one scalar DOF per active element. - sys.add_variable(flux_vars[g], libMesh::CONSTANT, libMesh::MONOMIAL); - eq_sys.init(); - - exo_reader.copy_elemental_solution( - sys, flux_vars[g], flux_vars[g], ts_1based); - - const libMesh::DofMap& dof_map = sys.get_dof_map(); - flux[g].resize(n_elem, 0.0); - int bin = 0; - for (const auto* elem : standalone_mesh.active_element_ptr_range()) { - std::vector dofs; - dof_map.dof_indices(elem, dofs); - // CONSTANT MONOMIAL → exactly one DOF per element. - flux[g][bin++] = sys.solution->el(dofs[0]); - } - // eq_sys destructs here; frees the DOF vectors for this group. - } - // standalone_mesh destructs here. - } - - // ── Step 2b: Register the Exodus mesh with OpenMC ───────────────────────── - - int mesh_id = 1; - for (const auto& m : model::meshes) - mesh_id = std::max(mesh_id, m->id_ + 1); - - // openmc::LibMesh constructors (include/openmc/mesh.h): - // LibMesh(const std::string& filename, double length_multiplier = 1.0) - // No 3-arg constructor; ID is set via set_id() after push_back. - auto lm_ptr = std::make_unique(mesh_file, 1.0); - - if (static_cast(lm_ptr->n_bins()) != n_elem) - fatal_error(fmt::format( - ": element count mismatch between " - "standalone read ({}) and openmc::LibMesh ({}).", - n_elem, lm_ptr->n_bins())); - - model::meshes.push_back(std::move(lm_ptr)); - model::meshes.back()->set_id(mesh_id); - - // ── Step 3: FW-CADIS normalisation ──────────────────────────────────────── - // - // From weight_windows.cpp update_weights() FW_CADIS branch (verified): - // 1. Invert all values across ALL groups simultaneously. - // 2. Find the GLOBAL maximum of the inverted values (single value for all - // groups — NOT per-group; per-group normalisation is the MAGIC method). - // 3. Normalise by 1/(2*global_max). - // 4. Elements with phi_adj <= 0 → sentinel -1.0 (no window). - // - // ── Energy group ordering ──────────────────────────────────────────────── - // - // The user lists in the SAME order as the energy - // intervals implied by : - // energy_bounds[0..1] → flux_vars[0] (lowest-energy group) - // energy_bounds[1..2] → flux_vars[1] - // ... - // energy_bounds[N-1..N] → flux_vars[N-1] (highest-energy group) - // - // This matches OpenMC's internal storage directly: - // lower_ww_(energy_bin, mesh_bin) where energy_bin=0 = lowest energy - // set_bounds(span) copies flat[g * n_elem + e] → lower_ww_(g, e) - // - // If your adjoint solver uses a different ordering (e.g. Griffin writes - // g0=fast), list the variable names in ascending-energy order in the XML, - // i.e. thermal group variable first. - - // Step A: invert; find global max of inverted values across all groups. - std::vector flat_lower(n_groups * n_elem, -1.0); - std::vector flat_upper(n_groups * n_elem, -1.0); - - double inv_max = 0.0; - for (int g = 0; g < n_groups; ++g) { - for (int 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; // temporary; normalised below - if (inv > inv_max) inv_max = inv; - } - // else: remains -1.0 (sentinel) - } - } - - if (inv_max <= 0.0) - 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, mesh_file)); - - // Step B: normalise by global 1/(2*inv_max) and set upper bounds. - const double norm_factor = 1.0 / (2.0 * inv_max); - for (int 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; - } - // else: both remain -1.0 (sentinel — transport kernel skips these) - } - - // ── Step 4: Build the WeightWindows object ──────────────────────────────── - // - // Verified API from include/openmc/weight_windows.h at fd1bc26a: - // WeightWindows::create() → allocates, registers, returns ptr - // set_mesh(int32_t mesh_idx) → takes vector index - // set_particle_type(ParticleType) - // set_energy_bounds(span)→ also calls allocate_ww_bounds() - // set_bounds(span, span) - // double& survival_ratio() → non-const ref, writable - // - // Order matters: set_mesh() and set_energy_bounds() must be called before - // set_bounds() because both call allocate_ww_bounds() which sizes the tensors. - // set_bounds() then checks that span size == n_energy_bins * n_mesh_bins and - // copies the data in. - { - WeightWindows* wws = WeightWindows::create(); - wws->set_mesh(model::mesh_map.at(mesh_id)); - wws->set_particle_type(ParticleType {"neutron"}); - wws->set_energy_bounds( - span(energy_bounds.data(), energy_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())); - } - - // ── Step 5: Enable weight windows globally ──────────────────────────────── - settings::weight_windows_on = true; - - { - 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, mesh_file, n_elem, varlist, - (ts_user < 0 ? n_steps_saved - 1 : ts_user), upper_bound_ratio), 5); - } - -#endif // OPENMC_LIBMESH_ENABLED -} - void read_settings_xml() { using namespace settings; @@ -1552,12 +1240,10 @@ void read_settings_xml(pugi::xml_node root) std::make_unique(node_ww)); } - // ──────────────────────────────────────────────────────── - // If the user provided a stanza, load the - // adjoint flux from the Exodus file and build WeightWindows in memory. - if (auto wwe = root.child("weight_windows_from_exodus")) - read_weight_windows_from_exodus(wwe); - // ──────────────────────────────────────────────────────────────────────── + // 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) 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 From 488f4fda31f4edc47a1cdc942f50fd14b90ff55e Mon Sep 17 00:00:00 2001 From: Foraejee Date: Thu, 13 Aug 2026 12:19:20 -0600 Subject: [PATCH 4/4] added weight window exodus to python api --- openmc/settings.py | 31 +++- openmc/weight_windows.py | 309 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 339 insertions(+), 1 deletion(-) 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