From 92ce1e78a60f7c6fbb379e87421edce294c2dca4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 3 Aug 2026 14:08:41 -0500 Subject: [PATCH 1/5] Initial implementation of C++ based WWList export --- include/openmc/capi.h | 9 ++ include/openmc/mesh.h | 18 ++- include/openmc/weight_windows.h | 3 + openmc/checkvalue.py | 12 ++ openmc/lib/mesh.py | 95 ++++++++++++++- openmc/lib/weight_windows.py | 15 ++- openmc/weight_windows.py | 72 +++++++++-- src/mesh.cpp | 115 ++++++++++++++++-- src/weight_windows.cpp | 24 +++- .../unit_tests/weightwindows/test_ww_list.py | 78 +++++++++++- 10 files changed, 410 insertions(+), 31 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 9f6987d74ed..d577f49c78b 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -101,6 +101,9 @@ int openmc_get_n_batches(int* n_batches, bool get_max_batches); int openmc_get_nuclide_index(const char name[], int* index); int openmc_add_unstructured_mesh( const char filename[], const char library[], int* id); +int openmc_add_unstructured_mesh_with_properties(const char filename[], + const char library[], double length_multiplier, const char options[], + int32_t id, int32_t* index); int64_t openmc_get_seed(); uint64_t openmc_get_stride(); int openmc_get_tally_index(int32_t id, int32_t* index); @@ -140,12 +143,18 @@ int openmc_mesh_filter_get_translation(int32_t index, double translation[3]); int openmc_mesh_filter_set_translation(int32_t index, double translation[3]); int openmc_mesh_get_id(int32_t index, int32_t* id); int openmc_mesh_set_id(int32_t index, int32_t id); +int openmc_mesh_get_name(int32_t index, const char** name); +int openmc_mesh_set_name(int32_t index, const char* name); int openmc_mesh_get_n_elements(int32_t index, size_t* n); int openmc_mesh_get_volumes(int32_t index, double* volumes); int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz, int max_mats, int32_t* materials, double* volumes, double* bboxes); int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh); int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh); +int openmc_cylindrical_mesh_get_origin(int32_t index, double origin[3]); +int openmc_cylindrical_mesh_set_origin(int32_t index, const double origin[3]); +int openmc_spherical_mesh_get_origin(int32_t index, double origin[3]); +int openmc_spherical_mesh_set_origin(int32_t index, const double origin[3]); int openmc_new_filter(const char* type, int32_t* index); int openmc_next_batch(int* status); int openmc_nuclide_name(int index, const char** name); diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index be97b09d4ec..94354109466 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -209,6 +209,8 @@ class Mesh { const std::string& name() const { return name_; } + void set_name(const std::string& name) { name_ = name; } + //! Set the mesh ID void set_id(int32_t id = -1); @@ -477,6 +479,16 @@ class PeriodicStructuredMesh : public StructuredMesh { return r - origin_; }; + const Position& origin() const { return origin_; } + + virtual int set_grid() = 0; + + int set_origin(Position origin) + { + origin_ = origin; + return set_grid(); + } + // Data members Position origin_ {0.0, 0.0, 0.0}; //!< Origin of the mesh }; @@ -834,7 +846,8 @@ class MOABMesh : public UnstructuredMesh { MOABMesh() = default; MOABMesh(pugi::xml_node); MOABMesh(hid_t group); - MOABMesh(const std::string& filename, double length_multiplier = 1.0); + MOABMesh(const std::string& filename, double length_multiplier = 1.0, + const std::string& options = {}); MOABMesh(std::shared_ptr external_mbi); static const std::string mesh_lib_type; @@ -1004,7 +1017,8 @@ class LibMesh : public UnstructuredMesh { // Constructors LibMesh(pugi::xml_node node); LibMesh(hid_t group); - LibMesh(const std::string& filename, double length_multiplier = 1.0); + LibMesh(const std::string& filename, double length_multiplier = 1.0, + const std::string& options = {}); LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0); static const std::string mesh_lib_type; diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index d0b385d169d..8d6692d01c3 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -109,6 +109,9 @@ class WeightWindows { //! Ready the weight window class for use void set_defaults(); + //! Replace the energy grid with defaults for the selected particle type + void reset_energy_bounds(); + //! Ensure the weight window lower bounds are properly allocated void allocate_ww_bounds(); diff --git a/openmc/checkvalue.py b/openmc/checkvalue.py index 5ff2cf9ac5a..cc747230bcb 100644 --- a/openmc/checkvalue.py +++ b/openmc/checkvalue.py @@ -1,6 +1,7 @@ import copy import os from collections.abc import Iterable +from numbers import Real import numpy as np @@ -80,7 +81,18 @@ def check_iterable_type(name, value, expected_type, min_depth=1, max_depth=1): max_depth : int The maximum number of layers of nested iterables there should be before reaching the ultimately contained items + + Notes + ----- + For NumPy floating-point arrays with an allowed number of dimensions, the + dtype guarantees the element type and the per-element scan is skipped when + *expected_type* is :class:`numbers.Real` or :class:`float`. """ + if (isinstance(value, np.ndarray) and value.dtype.kind == 'f' + and min_depth <= value.ndim <= max_depth + and expected_type in (Real, float)): + return + # Initialize the tree at the very first item. tree = [value] index = [0] diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 19e6f74d7ad..4485df71478 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -36,6 +36,12 @@ _dll.openmc_mesh_set_id.argtypes = [c_int32, c_int32] _dll.openmc_mesh_set_id.restype = c_int _dll.openmc_mesh_set_id.errcheck = _error_handler +_dll.openmc_mesh_get_name.argtypes = [c_int32, POINTER(c_char_p)] +_dll.openmc_mesh_get_name.restype = c_int +_dll.openmc_mesh_get_name.errcheck = _error_handler +_dll.openmc_mesh_set_name.argtypes = [c_int32, c_char_p] +_dll.openmc_mesh_set_name.restype = c_int +_dll.openmc_mesh_set_name.errcheck = _error_handler _dll.openmc_mesh_get_n_elements.argtypes = [c_int32, POINTER(c_size_t)] _dll.openmc_mesh_get_n_elements.restype = c_int _dll.openmc_mesh_get_n_elements.errcheck = _error_handler @@ -97,6 +103,12 @@ c_int, POINTER(c_double), c_int, POINTER(c_double), c_int] _dll.openmc_cylindrical_mesh_set_grid.restype = c_int _dll.openmc_cylindrical_mesh_set_grid.errcheck = _error_handler +_dll.openmc_cylindrical_mesh_get_origin.argtypes = [c_int32, POINTER(c_double)] +_dll.openmc_cylindrical_mesh_get_origin.restype = c_int +_dll.openmc_cylindrical_mesh_get_origin.errcheck = _error_handler +_dll.openmc_cylindrical_mesh_set_origin.argtypes = [c_int32, POINTER(c_double)] +_dll.openmc_cylindrical_mesh_set_origin.restype = c_int +_dll.openmc_cylindrical_mesh_set_origin.errcheck = _error_handler _dll.openmc_spherical_mesh_get_grid.argtypes = [c_int32, POINTER(POINTER(c_double)), POINTER(c_int), POINTER(POINTER(c_double)), @@ -107,6 +119,17 @@ c_int, POINTER(c_double), c_int, POINTER(c_double), c_int] _dll.openmc_spherical_mesh_set_grid.restype = c_int _dll.openmc_spherical_mesh_set_grid.errcheck = _error_handler +_dll.openmc_spherical_mesh_get_origin.argtypes = [c_int32, POINTER(c_double)] +_dll.openmc_spherical_mesh_get_origin.restype = c_int +_dll.openmc_spherical_mesh_get_origin.errcheck = _error_handler +_dll.openmc_spherical_mesh_set_origin.argtypes = [c_int32, POINTER(c_double)] +_dll.openmc_spherical_mesh_set_origin.restype = c_int +_dll.openmc_spherical_mesh_set_origin.errcheck = _error_handler + +_dll.openmc_add_unstructured_mesh_with_properties.argtypes = [ + c_char_p, c_char_p, c_double, c_char_p, c_int32, POINTER(c_int32)] +_dll.openmc_add_unstructured_mesh_with_properties.restype = c_int +_dll.openmc_add_unstructured_mesh_with_properties.errcheck = _error_handler class Mesh(_FortranObjectWithID): @@ -155,6 +178,16 @@ def id(self): def id(self, mesh_id): _dll.openmc_mesh_set_id(self._index, mesh_id) + @property + def name(self): + name = c_char_p() + _dll.openmc_mesh_get_name(self._index, name) + return name.value.decode() + + @name.setter + def name(self, name): + _dll.openmc_mesh_set_name(self._index, name.encode()) + @property def n_elements(self) -> int: n = c_size_t() @@ -621,6 +654,21 @@ def set_grid(self, r_grid, phi_grid, z_grid): _dll.openmc_cylindrical_mesh_set_grid(self._index, r_grid, nr, phi_grid, nphi, z_grid, nz) + @property + def origin(self): + origin = np.empty(3) + _dll.openmc_cylindrical_mesh_get_origin( + self._index, origin.ctypes.data_as(POINTER(c_double))) + return origin + + @origin.setter + def origin(self, origin): + origin = np.ascontiguousarray(origin, dtype=np.float64) + if origin.shape != (3,): + raise ValueError('Mesh origin must have three coordinates') + _dll.openmc_cylindrical_mesh_set_origin( + self._index, origin.ctypes.data_as(POINTER(c_double))) + class SphericalMesh(Mesh): """SphericalMesh stored internally. @@ -726,9 +774,54 @@ def set_grid(self, r_grid, theta_grid, phi_grid): _dll.openmc_spherical_mesh_set_grid(self._index, r_grid, nr, theta_grid, ntheta, phi_grid, nphi) + @property + def origin(self): + origin = np.empty(3) + _dll.openmc_spherical_mesh_get_origin( + self._index, origin.ctypes.data_as(POINTER(c_double))) + return origin + + @origin.setter + def origin(self, origin): + origin = np.ascontiguousarray(origin, dtype=np.float64) + if origin.shape != (3,): + raise ValueError('Mesh origin must have three coordinates') + _dll.openmc_spherical_mesh_set_origin( + self._index, origin.ctypes.data_as(POINTER(c_double))) + class UnstructuredMesh(Mesh): - pass + @classmethod + def from_file(cls, filename, library, uid=None, length_multiplier=1.0, + options=None): + """Create an unstructured mesh from a file. + + Parameters + ---------- + filename : path-like + Path to the unstructured mesh file. + library : {'libmesh', 'moab'} + Library used to load the mesh. + uid : int, optional + Unique ID for the mesh. If omitted, an ID is assigned. + length_multiplier : float, optional + Multiplicative factor applied to mesh coordinates. + options : str, optional + Options used to construct spatial search data structures. + + Returns + ------- + openmc.lib.UnstructuredMesh + The newly allocated mesh. + + """ + index = c_int32() + mesh_id = -1 if uid is None else uid + options = None if options is None else options.encode() + _dll.openmc_add_unstructured_mesh_with_properties( + str(filename).encode(), library.encode(), length_multiplier, + options, mesh_id, index) + return cls(index=index.value) _MESH_TYPE_MAP = { diff --git a/openmc/lib/weight_windows.py b/openmc/lib/weight_windows.py index 2b26d3b55f5..559be4d9719 100644 --- a/openmc/lib/weight_windows.py +++ b/openmc/lib/weight_windows.py @@ -194,10 +194,15 @@ def energy_bounds(self): @energy_bounds.setter def energy_bounds(self, e_bounds): - e_bounds_arr = np.asarray(e_bounds, dtype=float) - e_bounds_ptr = e_bounds_arr.ctypes.data_as(POINTER(c_double)) + if e_bounds is None: + e_bounds_ptr = None + size = 0 + else: + e_bounds_arr = np.ascontiguousarray(e_bounds, dtype=np.float64) + e_bounds_ptr = e_bounds_arr.ctypes.data_as(POINTER(c_double)) + size = e_bounds_arr.size _dll.openmc_weight_windows_set_energy_bounds( - self._index, e_bounds_ptr, e_bounds_arr.size) + self._index, e_bounds_ptr, size) @property def particle(self): @@ -222,8 +227,8 @@ def bounds(self): @bounds.setter def bounds(self, bounds): - lower = np.asarray(bounds[0]) - upper = np.asarray(bounds[1]) + lower = np.ascontiguousarray(bounds[0], dtype=np.float64) + upper = np.ascontiguousarray(bounds[1], dtype=np.float64) lower_p = lower.ctypes.data_as(POINTER(c_double)) upper_p = upper.ctypes.data_as(POINTER(c_double)) diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 63af2596efc..6ad160c8b0c 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -140,9 +140,7 @@ def __init__( "upper_bound_ratio must be present.") if upper_bound_ratio: - self.upper_ww_bounds = [ - lb * upper_bound_ratio for lb in self.lower_ww_bounds - ] + self.upper_ww_bounds = self.lower_ww_bounds * upper_bound_ratio if upper_ww_bounds is not None: self.upper_ww_bounds = upper_ww_bounds @@ -1074,18 +1072,78 @@ def export_to_hdf5(self, path: PathLike = 'weight_windows.h5', **init_kwargs): import openmc.lib cv.check_type('path', path, PathLike) - # Create a temporary model with the weight windows + # Create a minimal model without serializing the weight window bounds to + # XML. The meshes and weight windows are added through openmc.lib after + # the shared library has been initialized. model = openmc.Model() sph = openmc.Sphere(boundary_type='vacuum') cell = openmc.Cell(region=-sph) model.geometry = openmc.Geometry([cell]) - model.settings.weight_windows = self - model.settings.particles = 100 + model.settings.particles = 1 model.settings.batches = 1 + model.settings.output = {'summary': False} # Get absolute path before moving to temporary directory path = Path(path).resolve() + original_dir = Path.cwd() - # Load the model with openmc.lib and then export it to an HDF5 file + # Populate the C++ model directly and use its existing HDF5 writer. with openmc.lib.TemporarySession(model, **init_kwargs): + lib_meshes = {} + for ww in self: + mesh = ww.mesh + if mesh.id not in lib_meshes: + lib_meshes[mesh.id] = _create_lib_mesh(mesh, original_dir) + + lib_ww = openmc.lib.WeightWindows(ww.id) + lib_ww.particle = ww.particle_type + lib_ww.mesh = lib_meshes[mesh.id] + lib_ww.energy_bounds = ww.energy_bounds + + lower = np.ascontiguousarray( + ww.lower_ww_bounds.ravel(order='F'), dtype=np.float64) + upper = np.ascontiguousarray( + ww.upper_ww_bounds.ravel(order='F'), dtype=np.float64) + lib_ww.bounds = lower, upper + + lib_ww.survival_ratio = ww.survival_ratio + if ww.max_lower_bound_ratio is not None: + lib_ww.max_lower_bound_ratio = ww.max_lower_bound_ratio + lib_ww.max_split = ww.max_split + lib_ww.weight_cutoff = ww.weight_cutoff + openmc.lib.export_weight_windows(path) + + +def _create_lib_mesh(mesh: MeshBase, original_dir: Path): + """Create an openmc.lib mesh corresponding to a Python API mesh.""" + import openmc.lib + + if isinstance(mesh, openmc.RegularMesh): + lib_mesh = openmc.lib.RegularMesh(uid=mesh.id) + lib_mesh.dimension = mesh.dimension + lib_mesh.set_parameters( + lower_left=mesh.lower_left, upper_right=mesh.upper_right) + elif isinstance(mesh, openmc.RectilinearMesh): + lib_mesh = openmc.lib.RectilinearMesh(uid=mesh.id) + lib_mesh.set_grid(mesh.x_grid, mesh.y_grid, mesh.z_grid) + elif isinstance(mesh, openmc.CylindricalMesh): + lib_mesh = openmc.lib.CylindricalMesh(uid=mesh.id) + lib_mesh.set_grid(mesh.r_grid, mesh.phi_grid, mesh.z_grid) + lib_mesh.origin = mesh.origin + elif isinstance(mesh, openmc.SphericalMesh): + lib_mesh = openmc.lib.SphericalMesh(uid=mesh.id) + lib_mesh.set_grid(mesh.r_grid, mesh.theta_grid, mesh.phi_grid) + lib_mesh.origin = mesh.origin + elif isinstance(mesh, openmc.UnstructuredMesh): + filename = Path(mesh.filename) + if not filename.is_absolute(): + filename = original_dir / filename + lib_mesh = openmc.lib.UnstructuredMesh.from_file( + filename.resolve(), mesh.library, uid=mesh.id, + length_multiplier=mesh.length_multiplier, options=mesh.options) + else: + raise TypeError(f'Unsupported weight window mesh type: {type(mesh)}') + + lib_mesh.name = mesh.name + return lib_mesh diff --git a/src/mesh.cpp b/src/mesh.cpp index a0e497613c0..01f0df1c6e1 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2505,24 +2505,28 @@ extern "C" int openmc_extend_meshes( return 0; } -//! Adds a new unstructured mesh to OpenMC -extern "C" int openmc_add_unstructured_mesh( - const char filename[], const char library[], int* id) +//! Adds a new unstructured mesh to OpenMC with all supported properties +extern "C" int openmc_add_unstructured_mesh_with_properties( + const char filename[], const char library[], double length_multiplier, + const char options[], int32_t id, int32_t* index) { std::string lib_name(library); std::string mesh_file(filename); + std::string mesh_options(options ? options : ""); bool valid_lib = false; #ifdef OPENMC_DAGMC_ENABLED if (lib_name == MOABMesh::mesh_lib_type) { - model::meshes.push_back(std::move(make_unique(mesh_file))); + model::meshes.push_back( + make_unique(mesh_file, length_multiplier, mesh_options)); valid_lib = true; } #endif #ifdef OPENMC_LIBMESH_ENABLED if (lib_name == LibMesh::mesh_lib_type) { - model::meshes.push_back(std::move(make_unique(mesh_file))); + model::meshes.push_back( + make_unique(mesh_file, length_multiplier, mesh_options)); valid_lib = true; } #endif @@ -2534,13 +2538,26 @@ extern "C" int openmc_add_unstructured_mesh( return OPENMC_E_INVALID_ARGUMENT; } - // auto-assign new ID - model::meshes.back()->set_id(-1); - *id = model::meshes.back()->id_; + model::meshes.back()->set_id(id); + *index = model::meshes.size() - 1; return 0; } +//! Adds a new unstructured mesh to OpenMC +extern "C" int openmc_add_unstructured_mesh( + const char filename[], const char library[], int* id) +{ + int32_t index; + int err = openmc_add_unstructured_mesh_with_properties( + filename, library, 1.0, nullptr, C_NONE, &index); + if (err) + return err; + + *id = model::meshes[index]->id_; + return 0; +} + //! Return the index in the meshes array of a mesh with a given ID extern "C" int openmc_get_mesh_index(int32_t id, int32_t* index) { @@ -2567,8 +2584,25 @@ extern "C" int openmc_mesh_set_id(int32_t index, int32_t id) { if (int err = check_mesh(index)) return err; - model::meshes[index]->id_ = id; - model::mesh_map[id] = index; + model::meshes[index]->set_id(id); + return 0; +} + +//! Return the name of a mesh +extern "C" int openmc_mesh_get_name(int32_t index, const char** name) +{ + if (int err = check_mesh(index)) + return err; + *name = model::meshes[index]->name().c_str(); + return 0; +} + +//! Set the name of a mesh +extern "C" int openmc_mesh_set_name(int32_t index, const char* name) +{ + if (int err = check_mesh(index)) + return err; + model::meshes[index]->set_name(name); return 0; } @@ -2888,6 +2922,59 @@ extern "C" int openmc_spherical_mesh_set_grid(int32_t index, index, grid_x, nx, grid_y, ny, grid_z, nz); } +template +int openmc_periodic_mesh_get_origin_impl(int32_t index, double origin[3]) +{ + if (int err = check_mesh(index)) + return err; + T* mesh = dynamic_cast(model::meshes[index].get()); + if (!mesh) { + set_errmsg("This mesh is not of the expected type."); + return OPENMC_E_INVALID_TYPE; + } + const auto& mesh_origin = mesh->origin(); + origin[0] = mesh_origin.x; + origin[1] = mesh_origin.y; + origin[2] = mesh_origin.z; + return 0; +} + +template +int openmc_periodic_mesh_set_origin_impl(int32_t index, const double origin[3]) +{ + if (int err = check_mesh(index)) + return err; + T* mesh = dynamic_cast(model::meshes[index].get()); + if (!mesh) { + set_errmsg("This mesh is not of the expected type."); + return OPENMC_E_INVALID_TYPE; + } + return mesh->set_origin({origin[0], origin[1], origin[2]}); +} + +extern "C" int openmc_cylindrical_mesh_get_origin( + int32_t index, double origin[3]) +{ + return openmc_periodic_mesh_get_origin_impl(index, origin); +} + +extern "C" int openmc_cylindrical_mesh_set_origin( + int32_t index, const double origin[3]) +{ + return openmc_periodic_mesh_set_origin_impl(index, origin); +} + +extern "C" int openmc_spherical_mesh_get_origin(int32_t index, double origin[3]) +{ + return openmc_periodic_mesh_get_origin_impl(index, origin); +} + +extern "C" int openmc_spherical_mesh_set_origin( + int32_t index, const double origin[3]) +{ + return openmc_periodic_mesh_set_origin_impl(index, origin); +} + #ifdef OPENMC_DAGMC_ENABLED const std::string MOABMesh::mesh_lib_type = "moab"; @@ -2902,11 +2989,13 @@ MOABMesh::MOABMesh(hid_t group) : UnstructuredMesh(group) initialize(); } -MOABMesh::MOABMesh(const std::string& filename, double length_multiplier) +MOABMesh::MOABMesh(const std::string& filename, double length_multiplier, + const std::string& options) : UnstructuredMesh() { n_dimension_ = 3; filename_ = filename; + options_ = options; set_length_multiplier(length_multiplier); initialize(); } @@ -3633,9 +3722,11 @@ LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier) } // create the mesh from an input file -LibMesh::LibMesh(const std::string& filename, double length_multiplier) +LibMesh::LibMesh(const std::string& filename, double length_multiplier, + const std::string& options) { n_dimension_ = 3; + options_ = options; set_mesh_pointer_from_filename(filename); set_length_multiplier(length_multiplier); initialize(); diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index d3565eaaf63..081c49a7b8d 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -185,11 +185,25 @@ void WeightWindows::set_defaults() if (p_type == C_NONE) { fatal_error("Weight windows particle is not supported for transport."); } - energy_bounds_.push_back(data::energy_min[p_type]); - energy_bounds_.push_back(data::energy_max[p_type]); + double energy_min = data::energy_min[p_type]; + double energy_max = data::energy_max[p_type]; + if (energy_min >= energy_max) { + energy_min = 0.0; + energy_max = INFTY; + } + energy_bounds_.push_back(energy_min); + energy_bounds_.push_back(energy_max); } } +void WeightWindows::reset_energy_bounds() +{ + energy_bounds_.clear(); + set_defaults(); + if (mesh_idx_ != C_NONE) + allocate_ww_bounds(); +} + void WeightWindows::allocate_ww_bounds() { auto shape = bounds_size(); @@ -1156,7 +1170,11 @@ extern "C" int openmc_weight_windows_set_energy_bounds( if (int err = verify_ww_index(ww_idx)) return err; const auto& wws = variance_reduction::weight_windows.at(ww_idx); - wws->set_energy_bounds({e_bounds, e_bounds_size}); + if (e_bounds_size == 0) { + wws->reset_energy_bounds(); + } else { + wws->set_energy_bounds({e_bounds, e_bounds_size}); + } return 0; } diff --git a/tests/unit_tests/weightwindows/test_ww_list.py b/tests/unit_tests/weightwindows/test_ww_list.py index d148f382a53..ce94c973472 100644 --- a/tests/unit_tests/weightwindows/test_ww_list.py +++ b/tests/unit_tests/weightwindows/test_ww_list.py @@ -1,11 +1,21 @@ +import h5py +import numpy as np +import pytest + import openmc +import openmc.lib -def test_ww_roundtrip(request, run_in_tmpdir): +def test_ww_roundtrip(request, run_in_tmpdir, monkeypatch): # Load weight windows from a wwinp file wwinp_file = request.path.with_name('wwinp_n') wws = openmc.WeightWindowsList.from_wwinp(wwinp_file) + def fail_xml_export(*args, **kwargs): + pytest.fail('Weight windows should not be serialized to XML') + + monkeypatch.setattr(openmc.WeightWindows, 'to_xml_element', fail_xml_export) + # Roundtrip them, writing to HDF5 and reading back in wws.export_to_hdf5('ww.h5') wws_new = openmc.WeightWindowsList.from_hdf5('ww.h5') @@ -21,3 +31,69 @@ def test_ww_roundtrip(request, run_in_tmpdir): assert ww.max_split == ww_new.max_split assert ww.weight_cutoff == ww_new.weight_cutoff assert ww.mesh.id == ww_new.mesh.id + + +def test_export_hdf5_format(request, run_in_tmpdir): + # openmc_weight_windows_import expects this on-disk layout. + wws = openmc.WeightWindowsList.from_wwinp(request.path.with_name('wwinp_n')) + wws.export_to_hdf5('ww.h5') + + with h5py.File('ww.h5') as f: + assert f.attrs['filetype'] == b'weight_windows' + assert list(f.attrs['version']) == [1, 0] + wws_group = f['weight_windows'] + assert int(wws_group.attrs['n_weight_windows']) == len(wws) + for ww in wws: + group = wws_group[f'weight_windows_{ww.id}'] + assert group['lower_ww_bounds'].ndim == 2 + assert group['lower_ww_bounds'].shape[0] == ww.num_energy_bins + assert 'max_lower_bound_ratio' in group + + +def test_export_periodic_mesh_metadata(run_in_tmpdir): + cylindrical = openmc.CylindricalMesh( + r_grid=[0.0, 1.0, 2.0], phi_grid=[0.0, np.pi], + z_grid=[-1.0, 1.0], origin=(1.0, 2.0, 3.0), mesh_id=10, + name='cylindrical') + spherical = openmc.SphericalMesh( + r_grid=[0.0, 2.0], theta_grid=[0.0, np.pi], + phi_grid=[0.0, 2.0 * np.pi], origin=(-1.0, -2.0, -3.0), mesh_id=11, + name='spherical') + windows = openmc.WeightWindowsList([ + openmc.WeightWindows(cylindrical, [1.0, 2.0], + upper_bound_ratio=5.0, id=10), + openmc.WeightWindows(spherical, [3.0], + upper_bound_ratio=5.0, id=11), + ]) + + windows.export_to_hdf5('ww.h5') + roundtrip = openmc.WeightWindowsList.from_hdf5('ww.h5') + meshes = {ww.mesh.id: ww.mesh for ww in roundtrip} + + assert meshes[10].name == cylindrical.name + assert meshes[11].name == spherical.name + assert np.allclose(meshes[10].origin, cylindrical.origin) + assert np.allclose(meshes[11].origin, spherical.origin) + for ww in roundtrip: + assert ww.energy_bounds[0] == 0.0 + assert ww.energy_bounds[-1] == np.finfo(float).max + + +@pytest.mark.parametrize('library', ('libmesh', 'moab')) +def test_export_hdf5_unstructured_mesh(request, run_in_tmpdir, library): + if library == 'libmesh' and not openmc.lib._libmesh_enabled(): + pytest.skip('LibMesh not enabled in this build.') + if library == 'moab' and not openmc.lib._dagmc_enabled(): + pytest.skip('DAGMC (and MOAB) not enabled in this build.') + + mesh = openmc.UnstructuredMesh( + request.path.with_name('test_mesh_tets.exo'), library, mesh_id=20, + name='unstructured', length_multiplier=2.0) + ww = openmc.WeightWindows(mesh, np.ones(12_000), upper_bound_ratio=5.0) + openmc.WeightWindowsList([ww]).export_to_hdf5('ww.h5') + + with h5py.File('ww.h5') as f: + mesh_group = f['meshes'][f'mesh {mesh.id}'] + assert mesh_group['type'][()] == b'unstructured' + assert mesh_group['name'][()] == b'unstructured' + assert mesh_group['length_multiplier'][()] == 2.0 From 4aeabad5d71f9e530f0725ccec54081dba72b99b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 3 Aug 2026 14:57:23 -0500 Subject: [PATCH 2/5] Expose MeshBase.to_lib_object --- openmc/mesh.py | 69 +++++++++++++++++++ openmc/weight_windows.py | 37 +--------- .../unit_tests/weightwindows/test_ww_list.py | 25 +++++++ 3 files changed, 96 insertions(+), 35 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index bcb3ad2538e..4f12afd2c01 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -265,6 +265,75 @@ def axis_labels(self): """tuple of str : Names of the mesh axes, one per dimension.""" pass + def to_lib_object(self, uid: int | None = None, + base_dir: PathLike | None = None): + """Create a corresponding :mod:`openmc.lib` mesh. + + The OpenMC shared library must be initialized before calling this + method. The returned object is tied to the active library session and + becomes invalid when that session is finalized. + + Parameters + ---------- + uid : int, optional + ID to assign to the library mesh. If omitted, the ID of this mesh + is used. + base_dir : path-like, optional + Directory used to resolve relative filenames for unstructured + meshes. If omitted, the current working directory is used. + + Returns + ------- + openmc.lib.Mesh + The corresponding library mesh. The concrete type depends on the + type of this mesh. + + Raises + ------ + RuntimeError + If the OpenMC shared library has not been initialized. + + """ + import openmc.lib + + if not openmc.lib.is_initialized: + raise RuntimeError( + 'The OpenMC shared library must be initialized before ' + 'creating a library mesh.') + + if uid is None: + uid = self.id + base_dir = Path.cwd() if base_dir is None else Path(base_dir) + + if isinstance(self, RegularMesh): + lib_mesh = openmc.lib.RegularMesh(uid=uid) + lib_mesh.dimension = self.dimension + lib_mesh.set_parameters( + lower_left=self.lower_left, upper_right=self.upper_right) + elif isinstance(self, RectilinearMesh): + lib_mesh = openmc.lib.RectilinearMesh(uid=uid) + lib_mesh.set_grid(self.x_grid, self.y_grid, self.z_grid) + elif isinstance(self, CylindricalMesh): + lib_mesh = openmc.lib.CylindricalMesh(uid=uid) + lib_mesh.set_grid(self.r_grid, self.phi_grid, self.z_grid) + lib_mesh.origin = self.origin + elif isinstance(self, SphericalMesh): + lib_mesh = openmc.lib.SphericalMesh(uid=uid) + lib_mesh.set_grid(self.r_grid, self.theta_grid, self.phi_grid) + lib_mesh.origin = self.origin + elif isinstance(self, UnstructuredMesh): + filename = Path(self.filename) + if not filename.is_absolute(): + filename = base_dir / filename + lib_mesh = openmc.lib.UnstructuredMesh.from_file( + filename.resolve(), self.library, uid=uid, + length_multiplier=self.length_multiplier, options=self.options) + else: + raise TypeError(f'Unsupported mesh type: {type(self)}') + + lib_mesh.name = self.name + return lib_mesh + def __repr__(self): string = type(self).__name__ + '\n' string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 6ad160c8b0c..238dcf9e61c 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -1093,7 +1093,8 @@ def export_to_hdf5(self, path: PathLike = 'weight_windows.h5', **init_kwargs): for ww in self: mesh = ww.mesh if mesh.id not in lib_meshes: - lib_meshes[mesh.id] = _create_lib_mesh(mesh, original_dir) + lib_meshes[mesh.id] = mesh.to_lib_object( + base_dir=original_dir) lib_ww = openmc.lib.WeightWindows(ww.id) lib_ww.particle = ww.particle_type @@ -1113,37 +1114,3 @@ def export_to_hdf5(self, path: PathLike = 'weight_windows.h5', **init_kwargs): lib_ww.weight_cutoff = ww.weight_cutoff openmc.lib.export_weight_windows(path) - - -def _create_lib_mesh(mesh: MeshBase, original_dir: Path): - """Create an openmc.lib mesh corresponding to a Python API mesh.""" - import openmc.lib - - if isinstance(mesh, openmc.RegularMesh): - lib_mesh = openmc.lib.RegularMesh(uid=mesh.id) - lib_mesh.dimension = mesh.dimension - lib_mesh.set_parameters( - lower_left=mesh.lower_left, upper_right=mesh.upper_right) - elif isinstance(mesh, openmc.RectilinearMesh): - lib_mesh = openmc.lib.RectilinearMesh(uid=mesh.id) - lib_mesh.set_grid(mesh.x_grid, mesh.y_grid, mesh.z_grid) - elif isinstance(mesh, openmc.CylindricalMesh): - lib_mesh = openmc.lib.CylindricalMesh(uid=mesh.id) - lib_mesh.set_grid(mesh.r_grid, mesh.phi_grid, mesh.z_grid) - lib_mesh.origin = mesh.origin - elif isinstance(mesh, openmc.SphericalMesh): - lib_mesh = openmc.lib.SphericalMesh(uid=mesh.id) - lib_mesh.set_grid(mesh.r_grid, mesh.theta_grid, mesh.phi_grid) - lib_mesh.origin = mesh.origin - elif isinstance(mesh, openmc.UnstructuredMesh): - filename = Path(mesh.filename) - if not filename.is_absolute(): - filename = original_dir / filename - lib_mesh = openmc.lib.UnstructuredMesh.from_file( - filename.resolve(), mesh.library, uid=mesh.id, - length_multiplier=mesh.length_multiplier, options=mesh.options) - else: - raise TypeError(f'Unsupported weight window mesh type: {type(mesh)}') - - lib_mesh.name = mesh.name - return lib_mesh diff --git a/tests/unit_tests/weightwindows/test_ww_list.py b/tests/unit_tests/weightwindows/test_ww_list.py index ce94c973472..274d2a136d0 100644 --- a/tests/unit_tests/weightwindows/test_ww_list.py +++ b/tests/unit_tests/weightwindows/test_ww_list.py @@ -79,6 +79,31 @@ def test_export_periodic_mesh_metadata(run_in_tmpdir): assert ww.energy_bounds[-1] == np.finfo(float).max +def test_mesh_to_lib_object(run_in_tmpdir): + mesh = openmc.RegularMesh(mesh_id=17, name='runtime mesh') + mesh.dimension = (2, 3) + mesh.lower_left = (0.0, 1.0) + mesh.upper_right = (2.0, 4.0) + + with pytest.raises(RuntimeError, match='must be initialized'): + mesh.to_lib_object() + + model = openmc.Model() + sphere = openmc.Sphere(boundary_type='vacuum') + model.geometry = openmc.Geometry([openmc.Cell(region=-sphere)]) + model.settings.particles = 1 + model.settings.batches = 1 + + with openmc.lib.TemporarySession(model): + lib_mesh = mesh.to_lib_object() + assert isinstance(lib_mesh, openmc.lib.RegularMesh) + assert lib_mesh.id == mesh.id + assert lib_mesh.name == mesh.name + assert tuple(lib_mesh.dimension) == mesh.dimension + assert np.allclose(lib_mesh.lower_left, mesh.lower_left) + assert np.allclose(lib_mesh.upper_right, mesh.upper_right) + + @pytest.mark.parametrize('library', ('libmesh', 'moab')) def test_export_hdf5_unstructured_mesh(request, run_in_tmpdir, library): if library == 'libmesh' and not openmc.lib._libmesh_enabled(): From 9708ec0b12fbb89178a6a65765ff420168060d1a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 3 Aug 2026 15:13:57 -0500 Subject: [PATCH 3/5] Simplify Model use --- openmc/weight_windows.py | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 238dcf9e61c..57cd2eb0fe7 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -1072,23 +1072,12 @@ def export_to_hdf5(self, path: PathLike = 'weight_windows.h5', **init_kwargs): import openmc.lib cv.check_type('path', path, PathLike) - # Create a minimal model without serializing the weight window bounds to - # XML. The meshes and weight windows are added through openmc.lib after - # the shared library has been initialized. - model = openmc.Model() - sph = openmc.Sphere(boundary_type='vacuum') - cell = openmc.Cell(region=-sph) - model.geometry = openmc.Geometry([cell]) - model.settings.particles = 1 - model.settings.batches = 1 - model.settings.output = {'summary': False} - # Get absolute path before moving to temporary directory path = Path(path).resolve() original_dir = Path.cwd() # Populate the C++ model directly and use its existing HDF5 writer. - with openmc.lib.TemporarySession(model, **init_kwargs): + with openmc.lib.TemporarySession(**init_kwargs): lib_meshes = {} for ww in self: mesh = ww.mesh From ed9b48c738df75c5805a51c6d82e35316adeccf4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 3 Aug 2026 15:40:48 -0500 Subject: [PATCH 4/5] Use existing openmc_add_unstructured_mesh function --- include/openmc/capi.h | 7 ++----- openmc/lib/mesh.py | 8 ++++---- src/mesh.cpp | 20 +++----------------- 3 files changed, 9 insertions(+), 26 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index d577f49c78b..8332a9160ca 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -99,11 +99,8 @@ int openmc_get_material_index(int32_t id, int32_t* index); int openmc_get_mesh_index(int32_t id, int32_t* index); int openmc_get_n_batches(int* n_batches, bool get_max_batches); int openmc_get_nuclide_index(const char name[], int* index); -int openmc_add_unstructured_mesh( - const char filename[], const char library[], int* id); -int openmc_add_unstructured_mesh_with_properties(const char filename[], - const char library[], double length_multiplier, const char options[], - int32_t id, int32_t* index); +int openmc_add_unstructured_mesh(const char filename[], const char library[], + double length_multiplier, const char options[], int32_t id, int32_t* index); int64_t openmc_get_seed(); uint64_t openmc_get_stride(); int openmc_get_tally_index(int32_t id, int32_t* index); diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 4485df71478..10f377f42a5 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -126,10 +126,10 @@ _dll.openmc_spherical_mesh_set_origin.restype = c_int _dll.openmc_spherical_mesh_set_origin.errcheck = _error_handler -_dll.openmc_add_unstructured_mesh_with_properties.argtypes = [ +_dll.openmc_add_unstructured_mesh.argtypes = [ c_char_p, c_char_p, c_double, c_char_p, c_int32, POINTER(c_int32)] -_dll.openmc_add_unstructured_mesh_with_properties.restype = c_int -_dll.openmc_add_unstructured_mesh_with_properties.errcheck = _error_handler +_dll.openmc_add_unstructured_mesh.restype = c_int +_dll.openmc_add_unstructured_mesh.errcheck = _error_handler class Mesh(_FortranObjectWithID): @@ -818,7 +818,7 @@ def from_file(cls, filename, library, uid=None, length_multiplier=1.0, index = c_int32() mesh_id = -1 if uid is None else uid options = None if options is None else options.encode() - _dll.openmc_add_unstructured_mesh_with_properties( + _dll.openmc_add_unstructured_mesh( str(filename).encode(), library.encode(), length_multiplier, options, mesh_id, index) return cls(index=index.value) diff --git a/src/mesh.cpp b/src/mesh.cpp index 01f0df1c6e1..577fbc54224 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2506,9 +2506,9 @@ extern "C" int openmc_extend_meshes( } //! Adds a new unstructured mesh to OpenMC with all supported properties -extern "C" int openmc_add_unstructured_mesh_with_properties( - const char filename[], const char library[], double length_multiplier, - const char options[], int32_t id, int32_t* index) +extern "C" int openmc_add_unstructured_mesh(const char filename[], + const char library[], double length_multiplier, const char options[], + int32_t id, int32_t* index) { std::string lib_name(library); std::string mesh_file(filename); @@ -2544,20 +2544,6 @@ extern "C" int openmc_add_unstructured_mesh_with_properties( return 0; } -//! Adds a new unstructured mesh to OpenMC -extern "C" int openmc_add_unstructured_mesh( - const char filename[], const char library[], int* id) -{ - int32_t index; - int err = openmc_add_unstructured_mesh_with_properties( - filename, library, 1.0, nullptr, C_NONE, &index); - if (err) - return err; - - *id = model::meshes[index]->id_; - return 0; -} - //! Return the index in the meshes array of a mesh with a given ID extern "C" int openmc_get_mesh_index(int32_t id, int32_t* index) { From b4633d04edd8e8400a63191d1d44d80a37f929ce Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 7 Aug 2026 14:42:16 -0500 Subject: [PATCH 5/5] Move/add tests --- tests/unit_tests/test_checkvalue.py | 26 +++++++++++++++++++ tests/unit_tests/test_mesh.py | 25 ++++++++++++++++++ .../unit_tests/weightwindows/test_ww_list.py | 25 ------------------ 3 files changed, 51 insertions(+), 25 deletions(-) create mode 100644 tests/unit_tests/test_checkvalue.py diff --git a/tests/unit_tests/test_checkvalue.py b/tests/unit_tests/test_checkvalue.py new file mode 100644 index 00000000000..5e8081c71f0 --- /dev/null +++ b/tests/unit_tests/test_checkvalue.py @@ -0,0 +1,26 @@ +from numbers import Real + +import numpy as np +import pytest + +from openmc.checkvalue import check_iterable_type + + +@pytest.mark.parametrize('dtype', (np.float16, np.float32, np.float64)) +@pytest.mark.parametrize('expected_type', (Real, float)) +def test_check_iterable_type_float_array(dtype, expected_type): + values = np.ones((2, 3, 4), dtype=dtype) + check_iterable_type( + 'values', values, expected_type, min_depth=1, max_depth=3) + + +def test_check_iterable_type_float_array_depth(): + values = np.ones((2, 3)) + with pytest.raises(TypeError, match='maximum depth'): + check_iterable_type('values', values, Real, max_depth=1) + + +def test_check_iterable_type_nonfloat_array(): + values = np.ones(3, dtype=np.complex128) + with pytest.raises(TypeError, match='Items must be of type'): + check_iterable_type('values', values, Real) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 9b1469fc590..7c08d70dc02 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -117,6 +117,31 @@ def test_spherical_mesh_bounding_box(): np.testing.assert_array_equal(bb.upper_right, (4, 6, 8)) +def test_mesh_to_lib_object(run_in_tmpdir): + mesh = openmc.RegularMesh(mesh_id=17, name='runtime mesh') + mesh.dimension = (2, 3) + mesh.lower_left = (0.0, 1.0) + mesh.upper_right = (2.0, 4.0) + + with pytest.raises(RuntimeError, match='must be initialized'): + mesh.to_lib_object() + + model = openmc.Model() + sphere = openmc.Sphere(boundary_type='vacuum') + model.geometry = openmc.Geometry([openmc.Cell(region=-sphere)]) + model.settings.particles = 1 + model.settings.batches = 1 + + with openmc.lib.TemporarySession(model): + lib_mesh = mesh.to_lib_object() + assert isinstance(lib_mesh, openmc.lib.RegularMesh) + assert lib_mesh.id == mesh.id + assert lib_mesh.name == mesh.name + assert tuple(lib_mesh.dimension) == mesh.dimension + assert np.allclose(lib_mesh.lower_left, mesh.lower_left) + assert np.allclose(lib_mesh.upper_right, mesh.upper_right) + + def test_SphericalMesh_initiation(): # test defaults mesh = openmc.SphericalMesh(r_grid=(0, 10)) diff --git a/tests/unit_tests/weightwindows/test_ww_list.py b/tests/unit_tests/weightwindows/test_ww_list.py index 274d2a136d0..ce94c973472 100644 --- a/tests/unit_tests/weightwindows/test_ww_list.py +++ b/tests/unit_tests/weightwindows/test_ww_list.py @@ -79,31 +79,6 @@ def test_export_periodic_mesh_metadata(run_in_tmpdir): assert ww.energy_bounds[-1] == np.finfo(float).max -def test_mesh_to_lib_object(run_in_tmpdir): - mesh = openmc.RegularMesh(mesh_id=17, name='runtime mesh') - mesh.dimension = (2, 3) - mesh.lower_left = (0.0, 1.0) - mesh.upper_right = (2.0, 4.0) - - with pytest.raises(RuntimeError, match='must be initialized'): - mesh.to_lib_object() - - model = openmc.Model() - sphere = openmc.Sphere(boundary_type='vacuum') - model.geometry = openmc.Geometry([openmc.Cell(region=-sphere)]) - model.settings.particles = 1 - model.settings.batches = 1 - - with openmc.lib.TemporarySession(model): - lib_mesh = mesh.to_lib_object() - assert isinstance(lib_mesh, openmc.lib.RegularMesh) - assert lib_mesh.id == mesh.id - assert lib_mesh.name == mesh.name - assert tuple(lib_mesh.dimension) == mesh.dimension - assert np.allclose(lib_mesh.lower_left, mesh.lower_left) - assert np.allclose(lib_mesh.upper_right, mesh.upper_right) - - @pytest.mark.parametrize('library', ('libmesh', 'moab')) def test_export_hdf5_unstructured_mesh(request, run_in_tmpdir, library): if library == 'libmesh' and not openmc.lib._libmesh_enabled():