diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b8c14279b92..c3da79c0482 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -177,7 +177,7 @@ jobs: - name: Setup tmate debug session continue-on-error: true - if: ${{ contains(env.COMMIT_MESSAGE, '[gha-debug]') }} + if: ${{ failure() && contains(env.COMMIT_MESSAGE, '[gha-debug]') }} uses: mxschmitt/action-tmate@v3 timeout-minutes: 10 diff --git a/docs/source/io_formats/collision_track.rst b/docs/source/io_formats/collision_track.rst index 8e1e00ffb8e..4123fdac6b0 100644 --- a/docs/source/io_formats/collision_track.rst +++ b/docs/source/io_formats/collision_track.rst @@ -10,7 +10,7 @@ may also be written after each batch when multiple files are requested (``collision_track.N.h5``) or when the run is performed in parallel. The file contains the information needed to reconstruct each recorded collision. -The current revision of the collision track file format is 1.1. +The current revision of the collision track file format is 1.2. **/** @@ -33,7 +33,7 @@ The current revision of the collision track file format is 1.1. - ``event_mt`` (*int*) -- ENDF MT number identifying the reaction. - ``delayed_group`` (*int*) -- Delayed neutron group index (non-zero for delayed events). - ``cell_id`` (*int*) -- ID of the cell in which the collision occurred. - - ``nuclide_id`` (*int*) -- ZA identifier of the nuclide (ZZZAAAM format). + - ``nuclide_id`` (*int*) -- PDG number of the nuclide (100ZZZAAAM). - ``material_id`` (*int*) -- ID of the material containing the collision site. - ``universe_id`` (*int*) -- ID of the universe containing the collision site. - ``n_collision`` (*int*) -- Collision counter for the particle history. diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 1a436a75ede..fb02159169e 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -98,6 +98,11 @@ sub-elements: A list of strings representing the nuclide, to define specific define specific target nuclide collisions to be banked. + .. note:: + Electron and positron collision-track events are not associated with + a specific nuclide. If a ``nuclides`` entry is specified, these events + are omitted. + *Default*: None :reactions: @@ -606,30 +611,30 @@ found in the :ref:`random ray user guide `. *Default*: None :adjoint_source: - Specifies an adjoint fixed source for adjoint transport simulations, and - follows the format for :ref:`source_element`. The distributions which make - up the adjoint source are subject to the same restrictions as forward + Specifies an adjoint fixed source for adjoint transport simulations, and + follows the format for :ref:`source_element`. The distributions which make + up the adjoint source are subject to the same restrictions as forward fixed sources in Random Ray mode. *Default*: None - + :adjoint: - Specifies whether to perform adjoint transport. The default is 'False', + Specifies whether to perform adjoint transport. The default is 'False', corresponding to forward transport. *Default*: None - + :volume_estimator: - Specifies choice of volume estimator for the random ray solver. Options + Specifies choice of volume estimator for the random ray solver. Options are 'naive', 'simulation_averaged', or 'hybrid'. The default is 'hybrid'. *Default*: None :volume_normalized_flux_tallies: - Specifies whether to normalize flux tallies by volume (bool). The - default is 'False'. When enabled, flux tallies will be reported in units - of cm/cm^3. When disabled, flux tallies will be reported in units of cm - (i.e., total distance traveled by neutrons in the spatial tally + Specifies whether to normalize flux tallies by volume (bool). The + default is 'False'. When enabled, flux tallies will be reported in units + of cm/cm^3. When disabled, flux tallies will be reported in units of cm + (i.e., total distance traveled by neutrons in the spatial tally region). *Default*: None @@ -1757,11 +1762,11 @@ mesh-based weight windows. The ratio of the lower to upper weight window bounds. *Default*: 5.0 - + For FW-CADIS: :targets: - A sequence of IDs corresponding to the tallies which cover phase + A sequence of IDs corresponding to the tallies which cover phase space regions of interest for local variance reduction. *Default*: None diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 8ac07f3c892..e8514b8561a 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -792,6 +792,11 @@ collision_track.h5 file at the end of the simulation. The file contains 300 recorded collisions that occurred in materials with IDs 1 or 2, involving fission or (n,2n) reactions on the nuclides U-238 or O-16, within cells with IDs 5 and 12. + +.. note:: + Electron and positron collision-track events are not associated with a + specific nuclide. If a ``nuclides`` entry is specified, these events are omitted. + The file can be read using :func:`openmc.read_collision_track_file`. The example below shows how to extract the data from the collision_track feature and displays the fields stored in the file: diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 911654d318f..6b78145a4c2 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -123,8 +123,13 @@ 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); int openmc_plot_geometry(); +// Deprecated; use openmc_slice_data. int openmc_id_map(const void* slice, int32_t* data_out); +// Deprecated; use openmc_slice_data. int openmc_property_map(const void* slice, double* data_out); +int openmc_slice_data(const double origin[3], const double u_span[3], + const double v_span[3], const size_t pixels[2], bool show_overlaps, int level, + int32_t filter_index, int32_t* geom_data, double* property_data); int openmc_get_plot_index(int32_t id, int32_t* index); int openmc_plot_get_id(int32_t index, int32_t* id); int openmc_plot_set_id(int32_t index, int32_t id); diff --git a/include/openmc/constants.h b/include/openmc/constants.h index c49700f85fe..ccdc7ff4319 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -35,7 +35,7 @@ constexpr array VERSION_VOXEL {2, 0}; constexpr array VERSION_MGXS_LIBRARY {1, 0}; constexpr array VERSION_PROPERTIES {1, 1}; constexpr array VERSION_WEIGHT_WINDOWS {1, 0}; -constexpr array VERSION_COLLISION_TRACK {1, 1}; +constexpr array VERSION_COLLISION_TRACK {1, 2}; // ============================================================================ // ADJUSTABLE PARAMETERS diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 7a8b2acadd9..ae39a53ddf5 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -84,6 +84,9 @@ class Nuclide { double collapse_rate(int MT, double temperature, span energy, span flux) const; + //! Return a ParticleType object representing this nuclide + ParticleType particle_type() const { return {Z_, A_, metastable_}; } + //============================================================================ // Data members std::string name_; //!< Name of nuclide, e.g. "U235" diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 6c00fc2f6ba..f97d313847f 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -18,6 +18,8 @@ #include "openmc/position.h" #include "openmc/random_lcg.h" #include "openmc/ray.h" +#include "openmc/tallies/filter.h" +#include "openmc/tallies/filter_match.h" #include "openmc/xml_interface.h" namespace openmc { @@ -148,10 +150,11 @@ class PlottableInterface { struct IdData { // Constructor - IdData(size_t h_res, size_t v_res); + IdData(size_t h_res, size_t v_res, bool include_filter = false); // Methods - void set_value(size_t y, size_t x, const GeometryState& p, int level); + void set_value(size_t y, size_t x, const Particle& p, int level, + Filter* filter = nullptr, FilterMatch* match = nullptr); void set_overlap(size_t y, size_t x); // Members @@ -160,16 +163,34 @@ struct IdData { struct PropertyData { // Constructor - PropertyData(size_t h_res, size_t v_res); + PropertyData(size_t h_res, size_t v_res, bool include_filter = false); // Methods - void set_value(size_t y, size_t x, const GeometryState& p, int level); + void set_value(size_t y, size_t x, const Particle& p, int level, + Filter* filter = nullptr, FilterMatch* match = nullptr); void set_overlap(size_t y, size_t x); // Members tensor::Tensor data_; //!< 2D array of temperature & density data }; +struct RasterData { + // Constructor + RasterData(size_t h_res, size_t v_res, bool include_filter = false); + + // Methods + void set_value(size_t y, size_t x, const Particle& p, int level, + Filter* filter = nullptr, FilterMatch* match = nullptr); + void set_overlap(size_t y, size_t x); + + // Members + tensor::Tensor + id_data_; //!< [v_res, h_res, 3 or 4]: cell, instance, mat, [filter_bin] + tensor::Tensor + property_data_; //!< [v_res, h_res, 2]: temperature, density + bool include_filter_; //!< Whether filter bin index is included +}; + //=============================================================================== // Plot class //=============================================================================== @@ -177,7 +198,7 @@ struct PropertyData { class SlicePlotBase { public: template - T get_map() const; + T get_map(int32_t filter_index = -1) const; enum class PlotBasis { xy = 1, xz = 2, yz = 3 }; @@ -188,70 +209,65 @@ class SlicePlotBase { // Members public: - Position origin_; //!< Plot origin in geometry - Position width_; //!< Plot width in geometry - PlotBasis basis_; //!< Plot basis (XY/XZ/YZ) - array pixels_; //!< Plot size in pixels - bool slice_color_overlaps_; //!< Show overlapping cells? - int slice_level_ {-1}; //!< Plot universe level + Position origin_; //!< Plot origin in geometry + Direction u_span_; //!< Full-width span vector in geometry + Direction v_span_; //!< Full-height span vector in geometry + array pixels_; //!< Plot size in pixels + bool show_overlaps_; //!< Show overlapping cells? + int slice_level_ {-1}; //!< Plot universe level private: }; template -T SlicePlotBase::get_map() const +T SlicePlotBase::get_map(int32_t filter_index) const { size_t width = pixels_[0]; size_t height = pixels_[1]; - // get pixel size - double in_pixel = (width_[0]) / static_cast(width); - double out_pixel = (width_[1]) / static_cast(height); + // Determine if filter is being used + bool include_filter = (filter_index >= 0); + Filter* filter = nullptr; + if (include_filter) { + filter = model::tally_filters[filter_index].get(); + } // size data array - T data(width, height); - - // setup basis indices and initial position centered on pixel - int in_i, out_i; - Position xyz = origin_; - switch (basis_) { - case PlotBasis::xy: - in_i = 0; - out_i = 1; - break; - case PlotBasis::xz: - in_i = 0; - out_i = 2; - break; - case PlotBasis::yz: - in_i = 1; - out_i = 2; - break; - default: - UNREACHABLE(); - } + T data(width, height, include_filter); - // set initial position - xyz[in_i] = origin_[in_i] - width_[0] / 2. + in_pixel / 2.; - xyz[out_i] = origin_[out_i] + width_[1] / 2. - out_pixel / 2.; + // compute pixel steps and top-left pixel center + Direction u_step = u_span_ / static_cast(width); + Direction v_step = v_span_ / static_cast(height); + + Position start = + origin_ - 0.5 * u_span_ + 0.5 * v_span_ + 0.5 * u_step - 0.5 * v_step; + + // Validate that span vectors define a valid plane + Position cross = u_span_.cross(v_span_); + if (cross.norm() == 0.0) { + fatal_error("Slice span vectors are invalid (zero area)."); + } - // arbitrary direction - Direction dir = {1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.0}; + // Use an arbitrary direction that is not aligned with any coordinate axis. + // The direction has no physical meaning for plotting but is used by + // Surface::sense() to break ties when a pixel is coincident with a surface. + Direction dir = {1.0 / std::sqrt(2.0), 1.0 / std::sqrt(2.0), 0.0}; #pragma omp parallel { - GeometryState p; - p.r() = xyz; + Particle p; + p.r() = start; p.u() = dir; p.coord(0).universe() = model::root_universe; int level = slice_level_; int j {}; + FilterMatch match; #pragma omp for for (int y = 0; y < height; y++) { - p.r()[out_i] = xyz[out_i] - out_pixel * y; + Position row = start - v_step * static_cast(y); for (int x = 0; x < width; x++) { - p.r()[in_i] = xyz[in_i] + in_pixel * x; + p.r() = row + u_step * static_cast(x); p.n_coord() = 1; // local variables bool found_cell = exhaustive_find_cell(p); @@ -260,9 +276,9 @@ T SlicePlotBase::get_map() const j = level; } if (found_cell) { - data.set_value(y, x, p, j); + data.set_value(y, x, p, j, filter, &match); } - if (slice_color_overlaps_ && check_cell_overlap(p, false)) { + if (show_overlaps_ && check_cell_overlap(p, false)) { data.set_overlap(y, x); } } // inner for @@ -297,6 +313,8 @@ class Plot : public PlottableInterface, public SlicePlotBase { void print_info() const override; PlotType type_; //!< Plot type (Slice/Voxel) + Position width_; //!< Axis-aligned width from plot.xml + PlotBasis basis_; //!< Basis from plot.xml for slice plots int meshlines_width_; //!< Width of lines added to the plot int index_meshlines_mesh_ {-1}; //!< Index of the mesh to draw on the plot RGBColor meshlines_color_; //!< Color of meshlines on the plot diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index c12863956b9..bdfc32763be 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -339,8 +339,12 @@ def get_material_rates(self, mat_index, nuc_index, react_index): for i_nuc in nuc_index: nuc = self.nuc_ind_map[i_nuc] + if nuc not in xs._index_nuc: + continue for i_rx in react_index: rx = self.rx_ind_map[i_rx] + if rx not in xs._index_rx: + continue # Determine reaction rate by multiplying xs in [b] by flux # in [n-cm/src] to give [(reactions/src)*b-cm/atom] diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 687cf646f29..42bb958caf7 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -84,7 +84,12 @@ def get_microxs_and_flux( reactions listed in the depletion chain file are used. energies : iterable of float or str Energy group boundaries in [eV] or the name of the group structure. - If left as None energies will default to [0.0, 100e6] + If left as None, no energy filter is applied to the flux tally. When + `reaction_rate_mode` is "direct", these boundaries define the output + flux and microscopic cross section energy group structure. When + `reaction_rate_mode` is "flux", these boundaries define the multigroup + flux tally used to collapse continuous-energy cross sections; returned + fluxes and microscopic cross sections are one-group. reaction_rate_mode : {"direct", "flux"}, optional The "direct" method tallies reaction rates directly (per energy group). The "flux" method tallies a multigroup flux spectrum and then @@ -110,7 +115,9 @@ def get_microxs_and_flux( reaction_rate_opts : dict, optional When `reaction_rate_mode="flux"`, allows selecting a subset of nuclide/reaction pairs to be computed via direct reaction-rate tallies - (per energy group). Supported keys: "nuclides", "reactions". + over one energy bin spanning the full `energies` range. Supported keys: + "nuclides", "reactions". If "reactions" are specified without + "nuclides", all selected nuclides are used. Returns ------- @@ -139,10 +146,14 @@ def get_microxs_and_flux( nuclides = [nuc.name for nuc in chain.nuclides if nuc.name in nuclides_with_data] - # Set up the reaction rate and flux tallies + # Set up the reaction rate and flux tallies. When energies are omitted, no + # energy filter is needed for the transport calculation. A one-group energy + # range is still needed later if flux collapse is requested. + collapse_energies = energies if energies is None: - energies = [0.0, 100.0e6] - if isinstance(energies, str): + energy_filter = None + collapse_energies = [0.0, 100.0e6] + elif isinstance(energies, str): energy_filter = openmc.EnergyFilter.from_group_structure(energies) else: energy_filter = openmc.EnergyFilter(energies) @@ -172,8 +183,11 @@ def get_microxs_and_flux( rr_reactions = list(reactions) elif reaction_rate_mode == 'flux' and reaction_rate_opts: opts = reaction_rate_opts or {} - rr_nuclides = list(opts.get('nuclides', [])) rr_reactions = list(opts.get('reactions', [])) + if rr_reactions: + rr_nuclides = list(opts.get('nuclides', nuclides)) + else: + rr_nuclides = list(opts.get('nuclides', [])) # Keep only requested pairs within overall sets if rr_nuclides: rr_nuclides = [n for n in rr_nuclides if n in set(nuclides)] @@ -182,7 +196,7 @@ def get_microxs_and_flux( # Use 1-group energy filter for RR in flux mode has_rr = bool(rr_nuclides and rr_reactions) - if has_rr and reaction_rate_mode == 'flux': + if has_rr and reaction_rate_mode == 'flux' and energy_filter is not None: rr_energy_filter = openmc.EnergyFilter( [energy_filter.values[0], energy_filter.values[-1]]) else: @@ -194,14 +208,18 @@ def get_microxs_and_flux( model.tallies = [] for i, domain_filter in enumerate(domain_filters): flux_tally = openmc.Tally(name=f'MicroXS flux {i}') - flux_tally.filters = [domain_filter, energy_filter] + flux_tally.filters = [domain_filter] + if energy_filter is not None: + flux_tally.filters.append(energy_filter) flux_tally.scores = ['flux'] model.tallies.append(flux_tally) flux_tallies.append(flux_tally) if has_rr: rr_tally = openmc.Tally(name=f'MicroXS RR {i}') - rr_tally.filters = [domain_filter, rr_energy_filter] + rr_tally.filters = [domain_filter] + if rr_energy_filter is not None: + rr_tally.filters.append(rr_energy_filter) rr_tally.nuclides = rr_nuclides rr_tally.multiply_density = False rr_tally.scores = rr_reactions @@ -255,8 +273,12 @@ def get_microxs_and_flux( all_flux_arrays = [] for flux_tally in flux_tallies: # Get flux values and make energy groups last dimension - flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1) - flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups) + flux = flux_tally.get_reshaped_data() + if energy_filter is None: + flux = flux[..., np.newaxis] # (domains, 1, 1, groups) + else: + # (domains, groups, 1, 1) -> (domains, 1, 1, groups) + flux = np.moveaxis(flux, 1, -1) all_flux_arrays.append(flux) fluxes.extend(flux.squeeze((1, 2))) @@ -266,8 +288,15 @@ def get_microxs_and_flux( for flux_arr, rr_tally in zip(all_flux_arrays, rr_tallies): flux = flux_arr # Get reaction rates and make energy groups last dimension - reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions) - reaction_rates = np.moveaxis(reaction_rates, 1, -1) # (domains, nuclides, reactions, groups) + reaction_rates = rr_tally.get_reshaped_data() + if rr_energy_filter is None: + # (domains, nuclides, reactions) -> + # (domains, nuclides, reactions, groups) + reaction_rates = reaction_rates[..., np.newaxis] + else: + # (domains, groups, nuclides, reactions) -> + # (domains, nuclides, reactions, groups) + reaction_rates = np.moveaxis(reaction_rates, 1, -1) # If RR is 1-group, sum flux over groups if reaction_rate_mode == "flux": @@ -279,16 +308,20 @@ def get_microxs_and_flux( direct_micros.extend( MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs) - # If using flux mode, compute flux-collapsed microscopic XS if reaction_rate_mode == 'flux': + # Compute flux-collapsed microscopic XS flux_micros = [MicroXS.from_multigroup_flux( - energies=energies, + energies=collapse_energies, multigroup_flux=flux_i, chain_file=chain_file, nuclides=nuclides, reactions=reactions ) for flux_i in fluxes] + # We need to return one-group fluxes to match the microscopic cross + # sections, which are always one-group by virtue of the collapse + fluxes = [flux.sum(keepdims=True) for flux in fluxes] + # Decide which micros to use and merge if needed if reaction_rate_mode == 'flux' and rr_tallies: micros = [m1.merge(m2) for m1, m2 in zip(flux_micros, direct_micros)] diff --git a/openmc/examples.py b/openmc/examples.py index 94f1668211d..80ab98e2cb8 100644 --- a/openmc/examples.py +++ b/openmc/examples.py @@ -150,7 +150,7 @@ def pwr_core() -> openmc.Model: rpv_steel.add_nuclide('Ni60', 0.0026776, 'wo') rpv_steel.add_nuclide('Mn55', 0.01, 'wo') rpv_steel.add_nuclide('Cr52', 0.002092475, 'wo') - rpv_steel.add_nuclide('C0', 0.0025, 'wo') + rpv_steel.add_element('C', 0.0025, 'wo') rpv_steel.add_nuclide('Cu63', 0.0013696, 'wo') lower_rad_ref = openmc.Material(6, name='Lower radial reflector') diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py index 90af80d5b76..44d6ac273d3 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -9,6 +9,7 @@ from .error import _error_handler import numpy as np +import warnings class _Position(Structure): @@ -51,218 +52,209 @@ def __repr__(self): return f"({self.x}, {self.y}, {self.z})" -class _PlotBase(Structure): - """A structure defining a 2-D geometry slice with underlying c-types +def _extract_slice_data_args(plot): + """Convert a legacy plot-like object into slice_data keyword arguments.""" + try: + kwargs = { + 'origin': tuple(plot.origin), + 'width': (plot.width, plot.height), + 'basis': plot.basis, + 'pixels': (plot.h_res, plot.v_res), + 'show_overlaps': getattr(plot, 'color_overlaps', False), + 'level': getattr(plot, 'level', -1), + } + except AttributeError as exc: + raise TypeError( + "plot must be a legacy plot-like object with origin, width, " + "height, basis, h_res, and v_res attributes." + ) from exc + return kwargs + + +_dll.openmc_slice_data.argtypes = [ + POINTER(c_double * 3), # origin + POINTER(c_double * 3), # u_span + POINTER(c_double * 3), # v_span + POINTER(c_size_t * 2), # pixels + c_bool, # show_overlaps + c_int, # level + c_int32, # filter_index + POINTER(c_int32), # geom_data + POINTER(c_double), # property_data (can be None) +] +_dll.openmc_slice_data.restype = c_int +_dll.openmc_slice_data.errcheck = _error_handler + + +def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None, + pixels=None, show_overlaps=False, level=-1, filter=None, + include_properties=True): + """Generate a 2D raster of geometry and property data for plotting. - C-Type Attributes - ----------------- - origin_ : openmc.lib.plot._Position - A position defining the origin of the plot. - width_ : openmc.lib.plot._Position - The width of the plot along the x, y, and z axes, respectively - basis_ : c_int - The axes basis of the plot view. - pixels_ : c_size_t[3] - The resolution of the plot in the horizontal and vertical dimensions - color_overlaps_ : c_bool - Whether to assign unique IDs (-3) to overlapping regions. - level_ : c_int - The universe level for the plot view - - Attributes + Parameters ---------- - origin : tuple or list of ndarray - Origin (center) of the plot - width : float - The horizontal dimension of the plot in geometry units (cm) - height : float - The vertical dimension of the plot in geometry units (cm) - basis : string - One of {'xy', 'xz', 'yz'} indicating the horizontal and vertical - axes of the plot. - h_res : int - The horizontal resolution of the plot in pixels - v_res : int - The vertical resolution of the plot in pixels - level : int - The universe level for the plot (default: -1 -> all universes shown) - """ - _fields_ = [('origin_', _Position), - ('width_', _Position), - ('basis_', c_int), - ('pixels_', 3*c_size_t), - ('color_overlaps_', c_bool), - ('level_', c_int)] - - def __init__(self): - self.level_ = -1 - self.basis_ = 1 - self.color_overlaps_ = False - - @property - def origin(self): - return self.origin_ - - @origin.setter - def origin(self, origin): - self.origin_.x = origin[0] - self.origin_.y = origin[1] - self.origin_.z = origin[2] - - @property - def width(self): - return self.width_.x - - @width.setter - def width(self, width): - self.width_.x = width - - @property - def height(self): - return self.width_.y - - @height.setter - def height(self, height): - self.width_.y = height + origin : sequence of float + Center position of the plot [x, y, z] + width : sequence of float + Width of the plot [horizontal, vertical]. Mutually exclusive with + u_span/v_span. + basis : {'xy', 'xz', 'yz'} or int + Plot basis. Ignored if u_span/v_span are provided. + u_span : sequence of float, optional + Full-width span vector for the horizontal axis (3 values). Mutually + exclusive with width. + v_span : sequence of float, optional + Full-height span vector for the vertical axis (3 values). Mutually + exclusive with width. + pixels : sequence of int + Number of pixels [horizontal, vertical] + show_overlaps : bool, optional + Whether to detect overlapping cells + level : int, optional + Universe level (-1 for deepest) + filter : openmc.lib.Filter, optional + Filter for bin index lookup + include_properties : bool, optional + Whether to compute temperature/density - @property - def basis(self): - if self.basis_ == 1: - return 'xy' - elif self.basis_ == 2: - return 'xz' - elif self.basis_ == 3: - return 'yz' - - raise ValueError(f"Plot basis {self.basis_} is invalid") - - @basis.setter - def basis(self, basis): + Returns + ------- + geom_data : numpy.ndarray + Array of shape (v_res, h_res, 3) or (v_res, h_res, 4) with int32 dtype. + Contains [cell_id, cell_instance, material_id] when no filter is provided, + or [cell_id, cell_instance, material_id, filter_bin] when a filter is provided. + property_data : numpy.ndarray or None + Array of shape (v_res, h_res, 2) with float64 dtype containing + [temperature, density], or None if include_properties=False + """ + if pixels is None: + raise ValueError("pixels must be specified.") + if len(pixels) != 2: + raise ValueError("pixels must be a length-2 sequence.") + + if width is not None and (u_span is not None or v_span is not None): + raise ValueError("width is mutually exclusive with u_span/v_span.") + + if u_span is not None or v_span is not None: + if u_span is None or v_span is None: + raise ValueError("Both u_span and v_span must be provided.") + u_span = np.asarray(u_span, dtype=float) + v_span = np.asarray(v_span, dtype=float) + if u_span.shape != (3,) or v_span.shape != (3,): + raise ValueError("u_span and v_span must be length-3 sequences.") + u_norm = np.linalg.norm(u_span) + v_norm = np.linalg.norm(v_span) + if u_norm == 0.0 or v_norm == 0.0: + raise ValueError("u_span and v_span must be non-zero vectors.") + dot = float(np.dot(u_span, v_span)) + ortho_tol = 1.0e-10 * u_norm * v_norm + if abs(dot) > ortho_tol: + raise ValueError("u_span and v_span must be orthogonal.") + else: + if width is None: + raise ValueError("width must be provided when u_span/v_span are not set.") + if len(width) != 2: + raise ValueError("width must be a length-2 sequence.") + basis_map = {'xy': 1, 'xz': 2, 'yz': 3} if isinstance(basis, str): - valid_bases = ('xy', 'xz', 'yz') basis = basis.lower() - if basis not in valid_bases: + if basis not in basis_map: raise ValueError(f"{basis} is not a valid plot basis.") - - if basis == 'xy': - self.basis_ = 1 - elif basis == 'xz': - self.basis_ = 2 - elif basis == 'yz': - self.basis_ = 3 - return - - if isinstance(basis, int): - valid_bases = (1, 2, 3) - if basis not in valid_bases: + basis = basis_map[basis] + elif isinstance(basis, int): + if basis not in basis_map.values(): raise ValueError(f"{basis} is not a valid plot basis.") - self.basis_ = basis - return - - raise ValueError(f"{basis} of type {type(basis)} is an invalid plot basis") - - @property - def h_res(self): - return self.pixels_[0] - - @h_res.setter - def h_res(self, h_res): - self.pixels_[0] = h_res - - @property - def v_res(self): - return self.pixels_[1] - - @v_res.setter - def v_res(self, v_res): - self.pixels_[1] = v_res - - @property - def level(self): - return int(self.level_) - - @level.setter - def level(self, level): - self.level_ = level - - @property - def color_overlaps(self): - return self.color_overlaps_ - - @color_overlaps.setter - def color_overlaps(self, color_overlaps): - self.color_overlaps_ = color_overlaps - - def __repr__(self): - out_str = ["-----", - "Plot:", - "-----", - f"Origin: {self.origin}", - f"Width: {self.width}", - f"Height: {self.height}", - f"Basis: {self.basis}", - f"HRes: {self.h_res}", - f"VRes: {self.v_res}", - f"Color Overlaps: {self.color_overlaps}", - f"Level: {self.level}"] - return '\n'.join(out_str) - - -_dll.openmc_id_map.argtypes = [POINTER(_PlotBase), POINTER(c_int32)] -_dll.openmc_id_map.restype = c_int -_dll.openmc_id_map.errcheck = _error_handler + else: + raise ValueError(f"{basis} is not a valid plot basis.") + + if basis == 1: + u_span = np.array([width[0], 0.0, 0.0], dtype=float) + v_span = np.array([0.0, width[1], 0.0], dtype=float) + elif basis == 2: + u_span = np.array([width[0], 0.0, 0.0], dtype=float) + v_span = np.array([0.0, 0.0, width[1]], dtype=float) + else: + u_span = np.array([0.0, width[0], 0.0], dtype=float) + v_span = np.array([0.0, 0.0, width[1]], dtype=float) + + origin = np.asarray(origin, dtype=float) + if origin.shape != (3,): + raise ValueError("origin must be a length-3 sequence.") + + # Prepare ctypes arrays + origin_arr = (c_double * 3)(*origin) + u_span_arr = (c_double * 3)(*u_span) + v_span_arr = (c_double * 3)(*v_span) + pixels_arr = (c_size_t * 2)(*pixels) + + # Get internal filter index from filter ID if filter is provided + if filter is not None: + filter_index = c_int32() + _dll.openmc_get_filter_index(filter.id, filter_index) + filter_index = filter_index.value + else: + filter_index = -1 + + # Allocate output arrays with dynamic size based on filter + n_geom_fields = 4 if filter is not None else 3 + geom_data = np.zeros((pixels[1], pixels[0], n_geom_fields), dtype=np.int32) + if include_properties: + property_data = np.zeros((pixels[1], pixels[0], 2), dtype=np.float64) + prop_ptr = property_data.ctypes.data_as(POINTER(c_double)) + else: + property_data = None + prop_ptr = None + + _dll.openmc_slice_data( + origin_arr, + u_span_arr, + v_span_arr, + pixels_arr, + show_overlaps, + level, + filter_index, + geom_data.ctypes.data_as(POINTER(c_int32)), + prop_ptr + ) + + return geom_data, property_data def id_map(plot): - """ - Generate a 2-D map of cell and material IDs. Used for in-memory image - generation. - - Parameters - ---------- - plot : openmc.lib.plot._PlotBase - Object describing the slice of the model to be generated - - Returns - ------- - id_map : numpy.ndarray - A NumPy array with shape (vertical pixels, horizontal pixels, 3) of - OpenMC property ids with dtype int32. The last dimension of the array - contains, in order, cell IDs, cell instances, and material IDs. + """Deprecated compatibility wrapper for geometry ID maps. + This function is kept for compatibility and will be removed in a future + release. Use `slice_data(..., include_properties=False)` instead. """ - img_data = np.zeros((plot.v_res, plot.h_res, 3), - dtype=np.dtype('int32')) - _dll.openmc_id_map(plot, img_data.ctypes.data_as(POINTER(c_int32))) - return img_data + warnings.warn( + "openmc.lib.id_map is deprecated and will be removed in a future " + "release; use openmc.lib.slice_data(..., include_properties=False).", + FutureWarning, + ) - -_dll.openmc_property_map.argtypes = [POINTER(_PlotBase), POINTER(c_double)] -_dll.openmc_property_map.restype = c_int -_dll.openmc_property_map.errcheck = _error_handler + kwargs = _extract_slice_data_args(plot) + geom_data, _ = slice_data(include_properties=False, **kwargs) + return geom_data[:, :, :3] def property_map(plot): - """ - Generate a 2-D map of cell temperatures and material densities. Used for - in-memory image generation. - - Parameters - ---------- - plot : openmc.lib.plot._PlotBase - Object describing the slice of the model to be generated - - Returns - ------- - property_map : numpy.ndarray - A NumPy array with shape (vertical pixels, horizontal pixels, 2) of - OpenMC property ids with dtype float + """Deprecated compatibility wrapper for temperature/density maps. + This function is kept for compatibility and will be removed in a future + release. Use `slice_data(..., include_properties=True)` instead. """ - prop_data = np.zeros((plot.v_res, plot.h_res, 2)) - _dll.openmc_property_map(plot, prop_data.ctypes.data_as(POINTER(c_double))) + warnings.warn( + "openmc.lib.property_map is deprecated and will be removed in a " + "future release; use openmc.lib.slice_data(..., " + "include_properties=True).", + FutureWarning, + ) + + kwargs = _extract_slice_data_args(plot) + _, prop_data = slice_data(include_properties=True, **kwargs) return prop_data + _dll.openmc_get_plot_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_plot_index.restype = c_int _dll.openmc_get_plot_index.errcheck = _error_handler diff --git a/openmc/model/model.py b/openmc/model/model.py index da578159e30..82b0ae4dad8 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -292,7 +292,7 @@ def _assign_fw_cadis_tally_IDs(self): id_next = reference_tal.id break - if id_next == None: + if id_next is None: raise RuntimeError( f'Local FW-CADIS target tally {tal.id} not found on model.tallies!') else: @@ -1127,28 +1127,148 @@ def id_map( array contains cell IDs, cell instances, and material IDs (in that order). """ + ids, _ = self.slice_data( + origin=origin, + width=width, + pixels=pixels, + basis=basis, + show_overlaps=color_overlaps, + level=-1, + include_properties=False, + **init_kwargs, + ) + return ids + + def slice_data( + self, + origin: Sequence[float] | None = None, + width: Sequence[float] | None = None, + pixels: int | Sequence[int] = 40000, + basis: str = 'xy', + u_span: Sequence[float] | None = None, + v_span: Sequence[float] | None = None, + show_overlaps: bool = False, + level: int = -1, + filter: openmc.Filter | None = None, + include_properties: bool = True, + **init_kwargs + ) -> tuple[np.ndarray, np.ndarray | None]: + """Generate geometry and property data for a 2D plot slice. + + This method combines the functionality of :meth:`id_map` and property + mapping into a single call, avoiding duplicate geometry lookups. It also + supports filter bin index lookup for tally visualization. + + .. versionadded:: 0.16.0 + + Parameters + ---------- + origin : Sequence[float], optional + Origin of the plot. If unspecified, this argument defaults to the + center of the bounding box if the bounding box does not contain inf + values for the provided basis, otherwise (0.0, 0.0, 0.0). + width : Sequence[float], optional + Width of the plot. If unspecified, this argument defaults to the + width of the bounding box if the bounding box does not contain inf + values for the provided basis, otherwise (10.0, 10.0). + pixels : int | Sequence[int], optional + If an iterable of ints is provided then this directly sets the + number of pixels to use in each basis direction. If a single int is + provided then this sets the total number of pixels in the plot and + the number of pixels in each basis direction is calculated from this + total and the image aspect ratio based on the width argument. + basis : {'xy', 'yz', 'xz'}, optional + Basis of the plot. + u_span : Sequence[float], optional + Full-width span vector for an oriented slice (3 values). Mutually + exclusive with width. + v_span : Sequence[float], optional + Full-height span vector for an oriented slice (3 values). Mutually + exclusive with width. + show_overlaps : bool, optional + Whether to identify and assign unique IDs (-3) to overlapping + regions. If False, overlapping regions will be assigned the ID of + the lowest-numbered cell that occupies that region. Defaults to + False. + level : int, optional + Universe level to plot (-1 for deepest). Defaults to -1. + filter : openmc.Filter, optional + If provided, the information for each pixel also includes an index + in the filter corresponding to the pixel position. + include_properties : bool, optional + Whether to include temperature/density data. Defaults to True. + **init_kwargs + Keyword arguments passed to :meth:`Model.init_lib`. + + Returns + ------- + geom_data : numpy.ndarray + Shape (v_res, h_res, 3) or (v_res, h_res, 4) int32 array. Contains + [cell_id, cell_instance, material_id] when no filter, or [cell_id, + cell_instance, material_id, filter_bin] with filter. + property_data : numpy.ndarray or None + Shape (v_res, h_res, 2) float64 array with [temperature, density], + or None if include_properties=False. + """ import openmc.lib - origin, width, pixels = self._set_plot_defaults( - origin, width, pixels, basis) + if width is not None and (u_span is not None or v_span is not None): + raise ValueError("width is mutually exclusive with u_span/v_span.") - # initialize the openmc.lib.plot._PlotBase object - plot_obj = openmc.lib.plot._PlotBase() - plot_obj.origin = origin - plot_obj.width = width[0] - plot_obj.height = width[1] - plot_obj.h_res = pixels[0] - plot_obj.v_res = pixels[1] - plot_obj.basis = basis - plot_obj.color_overlaps = color_overlaps + if u_span is not None or v_span is not None: + if u_span is None or v_span is None: + raise ValueError("Both u_span and v_span must be provided.") + if origin is None: + origin = (0.0, 0.0, 0.0) + if isinstance(pixels, int): + u_norm = np.linalg.norm(u_span) + v_norm = np.linalg.norm(v_span) + aspect_ratio = u_norm / v_norm + pixels_y = math.sqrt(pixels / aspect_ratio) + pixels = (int(pixels / pixels_y), int(pixels_y)) + else: + origin, width, pixels = self._set_plot_defaults( + origin, width, pixels, basis) # Silence output by default. Also set arguments to start in volume # calculation mode to avoid loading cross sections init_kwargs.setdefault('output', False) init_kwargs.setdefault('args', ['-c']) + # If filter does not already appear in the model, temporarily add a + # tally with the filter + original_length = len(self.tallies) + if filter is not None: + filter_ids = {f.id for t in self.tallies for f in t.filters} + if filter.id not in filter_ids: + # Create temporary tally while preserving ID assignment + next_id = openmc.Tally.next_id + temp_tally = openmc.Tally() + temp_tally.filters = [filter] + temp_tally.scores = ['flux'] + self.tallies.append(temp_tally) + openmc.Tally.used_ids.remove(temp_tally.id) + openmc.Tally.next_id = next_id + with openmc.lib.TemporarySession(self, **init_kwargs): - return openmc.lib.id_map(plot_obj) + geom_data, property_data = openmc.lib.slice_data( + origin=origin, + width=width, + basis=basis, + u_span=u_span, + v_span=v_span, + pixels=pixels, + show_overlaps=show_overlaps, + level=level, + filter=filter, + include_properties=include_properties, + ) + + # If filter was temporarily added, remove it + if len(self.tallies) > original_length: + self.tallies.pop() + + return geom_data, property_data @add_plot_params def plot( @@ -1216,13 +1336,14 @@ def plot( "openmc.config before plotting.") break - # Get ID map from the C API - id_map = self.id_map( + # Get plot IDs from the C API + id_map, _ = self.slice_data( origin=origin, width=width, pixels=pixels, basis=basis, - color_overlaps=show_overlaps + show_overlaps=show_overlaps, + include_properties=False, ) # Generate colors if not provided @@ -1748,7 +1869,7 @@ def differentiate_mats(self, diff_volume_method: str = None, depletable_only: bo @staticmethod def _auto_generate_mgxs_lib( - model: openmc.model.model, + model: openmc.model.Model, groups: openmc.mgxs.EnergyGroups, correction: str | None, directory: PathLike, @@ -2712,12 +2833,15 @@ def convert_to_multigroup( self.settings.run_mode = original_run_mode break - # Make sure all materials have a name, and that the name is a valid HDF5 - # dataset name + # Temporarily replace each material's name with a unique, valid HDF5 + # dataset name (its name plus ID) for use as its MGXS library entry + # and macroscopic. The ID keeps the name unique even when materials + # share a name; the original names are restored at the end. + original_names = [material.name for material in self.materials] for material in self.materials: - if not material.name or not material.name.strip(): - material.name = f"material {material.id}" - material.name = re.sub(r'[^a-zA-Z0-9]', '_', material.name) + base = material.name if material.name and material.name.strip() \ + else "material" + material.name = re.sub(r'[^a-zA-Z0-9]', '_', base) + f"_{material.id}" # If needed, generate the needed MGXS data library file if not Path(mgxs_path).is_file() or overwrite_mgxs_library: @@ -2749,6 +2873,10 @@ def convert_to_multigroup( self.settings.energy_mode = 'multi-group' + # Restore the user's original material names. + for material, name in zip(self.materials, original_names): + material.name = name + def convert_to_random_ray(self): """Convert a multigroup model to use random ray. diff --git a/src/collision_track.cpp b/src/collision_track.cpp index 03cbc32b7b2..2e749007ff5 100644 --- a/src/collision_track.cpp +++ b/src/collision_track.cpp @@ -200,8 +200,13 @@ void collision_track_record(Particle& particle) return; int cell_id = model::cells[cell_index]->id_; - const auto* nuclide_ptr = data::nuclides[particle.event_nuclide()].get(); - std::string nuclide = nuclide_ptr->name_; + std::string nuclide {}; + int nuclide_id = 0; + if (particle.event_nuclide() != NUCLIDE_NONE) { + const auto* nuclide_ptr = data::nuclides[particle.event_nuclide()].get(); + nuclide = nuclide_ptr->name_; + nuclide_id = nuclide_ptr->particle_type().pdg_number(); + } int universe_id = model::universes[particle.lowest_coord().universe()]->id_; double delta_E = particle.E_last() - particle.E(); int material_index = particle.material(); @@ -224,8 +229,7 @@ void collision_track_record(Particle& particle) site.event_mt = particle.event_mt(); site.delayed_group = particle.delayed_group(); site.cell_id = cell_id; - site.nuclide_id = - 10000 * nuclide_ptr->Z_ + 10 * nuclide_ptr->A_ + nuclide_ptr->metastable_; + site.nuclide_id = nuclide_id; site.material_id = material_id; site.universe_id = universe_id; site.n_collision = particle.n_collision(); diff --git a/src/mesh.cpp b/src/mesh.cpp index 5ab7ac3988b..181af846694 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -3670,7 +3670,7 @@ Position LibMesh::sample_element(int32_t bin, uint64_t* seed) const // Get tet vertex coordinates from LibMesh std::array tet_verts; for (int i = 0; i < elem.n_nodes(); i++) { - auto node_ref = elem.node_ref(i); + const auto& node_ref = elem.node_ref(i); tet_verts[i] = {node_ref(0), node_ref(1), node_ref(2)}; } // Samples position within tet using Barycentric coordinates @@ -3700,7 +3700,7 @@ int LibMesh::n_vertices() const Position LibMesh::vertex(int vertex_id) const { - const auto node_ref = m_->node_ref(vertex_id); + const auto& node_ref = m_->node_ref(vertex_id); if (length_multiplier_ > 0.0) { return length_multiplier_ * Position(node_ref(0), node_ref(1), node_ref(2)); } else { diff --git a/src/plot.cpp b/src/plot.cpp index 707d53dc2cb..b9a1136afdb 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -33,6 +33,7 @@ #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/string_utils.h" +#include "openmc/tallies/filter.h" namespace openmc { @@ -44,10 +45,12 @@ constexpr int PLOT_LEVEL_LOWEST {-1}; //!< lower bound on plot universe level constexpr int32_t NOT_FOUND {-2}; constexpr int32_t OVERLAP {-3}; -IdData::IdData(size_t h_res, size_t v_res) : data_({v_res, h_res, 3}, NOT_FOUND) +IdData::IdData(size_t h_res, size_t v_res, bool /*include_filter*/) + : data_({v_res, h_res, 3}, NOT_FOUND) {} -void IdData::set_value(size_t y, size_t x, const GeometryState& p, int level) +void IdData::set_value(size_t y, size_t x, const Particle& p, int level, + Filter* /*filter*/, FilterMatch* /*match*/) { // set cell data if (p.n_coord() <= level) { @@ -64,7 +67,6 @@ void IdData::set_value(size_t y, size_t x, const GeometryState& p, int level) Cell* c = model::cells.at(p.lowest_coord().cell()).get(); if (p.material() == MATERIAL_VOID) { data_(y, x, 2) = MATERIAL_VOID; - return; } else if (c->type_ == Fill::MATERIAL) { Material* m = model::materials.at(p.material()).get(); data_(y, x, 2) = m->id_; @@ -77,12 +79,12 @@ void IdData::set_overlap(size_t y, size_t x) data_(y, x, k) = OVERLAP; } -PropertyData::PropertyData(size_t h_res, size_t v_res) +PropertyData::PropertyData(size_t h_res, size_t v_res, bool /*include_filter*/) : data_({v_res, h_res, 2}, NOT_FOUND) {} -void PropertyData::set_value( - size_t y, size_t x, const GeometryState& p, int level) +void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level, + Filter* /*filter*/, FilterMatch* /*match*/) { Cell* c = model::cells.at(p.lowest_coord().cell()).get(); data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN; @@ -97,6 +99,74 @@ void PropertyData::set_overlap(size_t y, size_t x) data_(y, x) = OVERLAP; } +//============================================================================== +// RasterData implementation +//============================================================================== + +RasterData::RasterData(size_t h_res, size_t v_res, bool include_filter) + : id_data_({v_res, h_res, include_filter ? 4u : 3u}, NOT_FOUND), + property_data_({v_res, h_res, 2}, static_cast(NOT_FOUND)), + include_filter_(include_filter) +{} + +void RasterData::set_value(size_t y, size_t x, const Particle& p, int level, + Filter* filter, FilterMatch* match) +{ + // set cell data + if (p.n_coord() <= level) { + id_data_(y, x, 0) = NOT_FOUND; + id_data_(y, x, 1) = NOT_FOUND; + } else { + id_data_(y, x, 0) = model::cells.at(p.coord(level).cell())->id_; + id_data_(y, x, 1) = level == p.n_coord() - 1 + ? p.cell_instance() + : cell_instance_at_level(p, level); + } + + // set material data + Cell* c = model::cells.at(p.lowest_coord().cell()).get(); + if (p.material() == MATERIAL_VOID) { + id_data_(y, x, 2) = MATERIAL_VOID; + } else if (c->type_ == Fill::MATERIAL) { + Material* m = model::materials.at(p.material()).get(); + id_data_(y, x, 2) = m->id_; + } + + // set filter index (only if filter is being used) + if (include_filter_ && filter) { + filter->get_all_bins(p, TallyEstimator::COLLISION, *match); + if (match->bins_.empty()) { + id_data_(y, x, 3) = -1; + } else { + id_data_(y, x, 3) = match->bins_[0]; + } + match->bins_.clear(); + match->weights_.clear(); + } + + // set temperature (in K) + property_data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN; + + // set density (g/cm³) + if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) { + Material* m = model::materials.at(p.material()).get(); + property_data_(y, x, 1) = m->density_gpcc_; + } +} + +void RasterData::set_overlap(size_t y, size_t x) +{ + // Set cell, instance, and material to OVERLAP, but preserve filter bin + id_data_(y, x, 0) = OVERLAP; + id_data_(y, x, 1) = OVERLAP; + id_data_(y, x, 2) = OVERLAP; + // Note: id_data_(y, x, 3) is NOT overwritten - preserves filter bin for tally + // plotting + + property_data_(y, x, 0) = OVERLAP; + property_data_(y, x, 1) = OVERLAP; +} + //============================================================================== // Global variables //============================================================================== @@ -450,6 +520,22 @@ void Plot::set_width(pugi::xml_node plot_node) if (pl_width.size() == 2) { width_.x = pl_width[0]; width_.y = pl_width[1]; + switch (basis_) { + case PlotBasis::xy: + u_span_ = {width_.x, 0.0, 0.0}; + v_span_ = {0.0, width_.y, 0.0}; + break; + case PlotBasis::xz: + u_span_ = {width_.x, 0.0, 0.0}; + v_span_ = {0.0, 0.0, width_.y}; + break; + case PlotBasis::yz: + u_span_ = {0.0, width_.x, 0.0}; + v_span_ = {0.0, 0.0, width_.y}; + break; + default: + UNREACHABLE(); + } } else { fatal_error( fmt::format(" must be length 2 in slice plot {}", id())); @@ -765,7 +851,7 @@ Plot::Plot(pugi::xml_node plot_node, PlotType type) set_width(plot_node); set_meshlines(plot_node); slice_level_ = level_; // Copy level employed in SlicePlotBase::get_map - slice_color_overlaps_ = color_overlaps_; + show_overlaps_ = color_overlaps_; } //============================================================================== @@ -862,23 +948,39 @@ void Plot::draw_mesh_lines(ImageData& data) const rgb = meshlines_color_; int ax1, ax2; + Position expected_u {}; + Position expected_v {}; switch (basis_) { case PlotBasis::xy: ax1 = 0; ax2 = 1; + expected_u = {width_[0], 0.0, 0.0}; + expected_v = {0.0, width_[1], 0.0}; break; case PlotBasis::xz: ax1 = 0; ax2 = 2; + expected_u = {width_[0], 0.0, 0.0}; + expected_v = {0.0, 0.0, width_[1]}; break; case PlotBasis::yz: ax1 = 1; ax2 = 2; + expected_u = {0.0, width_[0], 0.0}; + expected_v = {0.0, 0.0, width_[1]}; break; default: UNREACHABLE(); } + // Meshlines rely on axis-aligned indexing in global coordinates. + constexpr double rel_tol {1e-12}; + double span_tol = rel_tol * (1.0 + u_span_.norm() + v_span_.norm()); + if ((u_span_ - expected_u).norm() > span_tol || + (v_span_ - expected_v).norm() > span_tol) { + fatal_error("Meshlines are only supported for axis-aligned slice plots."); + } + Position ll_plot {origin_}; Position ur_plot {origin_}; @@ -1008,11 +1110,11 @@ void Plot::create_voxel() const voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace); SlicePlotBase pltbase; - pltbase.width_ = width_; pltbase.origin_ = origin_; - pltbase.basis_ = PlotBasis::xy; + pltbase.u_span_ = {width_.x, 0.0, 0.0}; + pltbase.v_span_ = {0.0, width_.y, 0.0}; pltbase.pixels() = pixels(); - pltbase.slice_color_overlaps_ = color_overlaps_; + pltbase.show_overlaps_ = color_overlaps_; ProgressBar pb; for (int z = 0; z < pixels()[2]; z++) { @@ -1794,6 +1896,12 @@ void PhongRay::on_intersection() extern "C" int openmc_id_map(const void* plot, int32_t* data_out) { + static bool warned {false}; + if (!warned) { + warning("openmc_id_map is deprecated and will be removed in a future " + "release. Use openmc_slice_data."); + warned = true; + } auto plt = reinterpret_cast(plot); if (!plt) { @@ -1801,7 +1909,7 @@ extern "C" int openmc_id_map(const void* plot, int32_t* data_out) return OPENMC_E_INVALID_ARGUMENT; } - if (plt->slice_color_overlaps_ && model::overlap_check_count.size() == 0) { + if (plt->show_overlaps_ && model::overlap_check_count.size() == 0) { model::overlap_check_count.resize(model::cells.size()); } @@ -1815,14 +1923,20 @@ extern "C" int openmc_id_map(const void* plot, int32_t* data_out) extern "C" int openmc_property_map(const void* plot, double* data_out) { + static bool warned {false}; + if (!warned) { + warning("openmc_property_map is deprecated and will be removed in a future " + "release. Use openmc_slice_data."); + warned = true; + } auto plt = reinterpret_cast(plot); if (!plt) { - set_errmsg("Invalid slice pointer passed to openmc_id_map"); + set_errmsg("Invalid slice pointer passed to openmc_property_map"); return OPENMC_E_INVALID_ARGUMENT; } - if (plt->slice_color_overlaps_ && model::overlap_check_count.size() == 0) { + if (plt->show_overlaps_ && model::overlap_check_count.size() == 0) { model::overlap_check_count.resize(model::cells.size()); } @@ -1834,6 +1948,68 @@ extern "C" int openmc_property_map(const void* plot, double* data_out) return 0; } +extern "C" int openmc_slice_data(const double origin[3], const double u_span[3], + const double v_span[3], const size_t pixels[2], bool color_overlaps, + int level, int32_t filter_index, int32_t* geom_data, double* property_data) +{ + // Validate span vectors + Direction u_span_pos {u_span[0], u_span[1], u_span[2]}; + Direction v_span_pos {v_span[0], v_span[1], v_span[2]}; + double u_norm = u_span_pos.norm(); + double v_norm = v_span_pos.norm(); + if (u_norm == 0.0 || v_norm == 0.0) { + set_errmsg("Slice span vectors must be non-zero."); + return OPENMC_E_INVALID_ARGUMENT; + } + + constexpr double ORTHO_REL_TOL = 1e-10; + double dot = u_span_pos.dot(v_span_pos); + if (std::abs(dot) > ORTHO_REL_TOL * u_norm * v_norm) { + set_errmsg("Slice span vectors must be orthogonal."); + return OPENMC_E_INVALID_ARGUMENT; + } + + // Validate filter index if provided + if (filter_index >= 0) { + if (int err = verify_filter(filter_index)) + return err; + } + + // Initialize overlap check vector if needed + if (color_overlaps && model::overlap_check_count.size() == 0) { + model::overlap_check_count.resize(model::cells.size()); + } + + try { + // Create a temporary SlicePlotBase object to reuse get_map logic + SlicePlotBase plot_params; + plot_params.origin_ = Position {origin[0], origin[1], origin[2]}; + plot_params.u_span_ = u_span_pos; + plot_params.v_span_ = v_span_pos; + plot_params.pixels_[0] = pixels[0]; + plot_params.pixels_[1] = pixels[1]; + plot_params.show_overlaps_ = color_overlaps; + plot_params.slice_level_ = level; + + // Use get_map to generate data + auto data = plot_params.get_map(filter_index); + + // Copy geometry data + std::copy(data.id_data_.begin(), data.id_data_.end(), geom_data); + + // Copy property data if requested + if (property_data != nullptr) { + std::copy( + data.property_data_.begin(), data.property_data_.end(), property_data); + } + } catch (const std::exception& e) { + set_errmsg(e.what()); + return OPENMC_E_UNASSIGNED; + } + + return 0; +} + extern "C" int openmc_get_plot_index(int32_t id, int32_t* index) { auto it = model::plot_map.find(id); diff --git a/tests/regression_tests/collision_track/case_1_Reactions/collision_track_true.h5 b/tests/regression_tests/collision_track/case_1_Reactions/collision_track_true.h5 new file mode 100644 index 00000000000..c3154621342 Binary files /dev/null and b/tests/regression_tests/collision_track/case_1_Reactions/collision_track_true.h5 differ diff --git a/tests/regression_tests/collision_track/case_1_Reactions/inputs_true.dat b/tests/regression_tests/collision_track/case_1_Reactions/inputs_true.dat index 7533616c059..005a5602058 100644 --- a/tests/regression_tests/collision_track/case_1_Reactions/inputs_true.dat +++ b/tests/regression_tests/collision_track/case_1_Reactions/inputs_true.dat @@ -38,9 +38,9 @@ eigenvalue - 100 + 80 5 - 1 + 4 -2.0 -2.0 -2.0 2.0 2.0 2.0 @@ -51,7 +51,7 @@ (n,fission) 101 - 300 + 100 1 diff --git a/tests/regression_tests/collision_track/case_1_Reactions/results_true.dat b/tests/regression_tests/collision_track/case_1_Reactions/results_true.dat deleted file mode 100644 index d4d1d1e5ad4..00000000000 --- a/tests/regression_tests/collision_track/case_1_Reactions/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_2_Cell_ID/collision_track_true.h5 b/tests/regression_tests/collision_track/case_2_Cell_ID/collision_track_true.h5 new file mode 100644 index 00000000000..850bba8415e Binary files /dev/null and b/tests/regression_tests/collision_track/case_2_Cell_ID/collision_track_true.h5 differ diff --git a/tests/regression_tests/collision_track/case_2_Cell_ID/inputs_true.dat b/tests/regression_tests/collision_track/case_2_Cell_ID/inputs_true.dat index 55fb835de03..8d46c3d3b97 100644 --- a/tests/regression_tests/collision_track/case_2_Cell_ID/inputs_true.dat +++ b/tests/regression_tests/collision_track/case_2_Cell_ID/inputs_true.dat @@ -38,9 +38,9 @@ eigenvalue - 100 + 80 5 - 1 + 4 -2.0 -2.0 -2.0 2.0 2.0 2.0 @@ -51,7 +51,7 @@ 22 - 300 + 100 1 diff --git a/tests/regression_tests/collision_track/case_2_Cell_ID/results_true.dat b/tests/regression_tests/collision_track/case_2_Cell_ID/results_true.dat deleted file mode 100644 index d4d1d1e5ad4..00000000000 --- a/tests/regression_tests/collision_track/case_2_Cell_ID/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_3_Material_ID/collision_track_true.h5 b/tests/regression_tests/collision_track/case_3_Material_ID/collision_track_true.h5 new file mode 100644 index 00000000000..777d91f0888 Binary files /dev/null and b/tests/regression_tests/collision_track/case_3_Material_ID/collision_track_true.h5 differ diff --git a/tests/regression_tests/collision_track/case_3_Material_ID/inputs_true.dat b/tests/regression_tests/collision_track/case_3_Material_ID/inputs_true.dat index 61890414bab..322e74f4938 100644 --- a/tests/regression_tests/collision_track/case_3_Material_ID/inputs_true.dat +++ b/tests/regression_tests/collision_track/case_3_Material_ID/inputs_true.dat @@ -38,9 +38,9 @@ eigenvalue - 100 + 80 5 - 1 + 4 -2.0 -2.0 -2.0 2.0 2.0 2.0 @@ -51,7 +51,7 @@ 1 - 300 + 100 1 diff --git a/tests/regression_tests/collision_track/case_3_Material_ID/results_true.dat b/tests/regression_tests/collision_track/case_3_Material_ID/results_true.dat deleted file mode 100644 index d4d1d1e5ad4..00000000000 --- a/tests/regression_tests/collision_track/case_3_Material_ID/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_4_Nuclide_ID/collision_track_true.h5 b/tests/regression_tests/collision_track/case_4_Nuclide_ID/collision_track_true.h5 new file mode 100644 index 00000000000..527705297fe Binary files /dev/null and b/tests/regression_tests/collision_track/case_4_Nuclide_ID/collision_track_true.h5 differ diff --git a/tests/regression_tests/collision_track/case_4_Nuclide_ID/inputs_true.dat b/tests/regression_tests/collision_track/case_4_Nuclide_ID/inputs_true.dat index 8960dde5cbd..6202bcacc75 100644 --- a/tests/regression_tests/collision_track/case_4_Nuclide_ID/inputs_true.dat +++ b/tests/regression_tests/collision_track/case_4_Nuclide_ID/inputs_true.dat @@ -38,9 +38,9 @@ eigenvalue - 100 + 80 5 - 1 + 4 -2.0 -2.0 -2.0 2.0 2.0 2.0 @@ -51,7 +51,7 @@ O16 U235 - 300 + 100 1 diff --git a/tests/regression_tests/collision_track/case_4_Nuclide_ID/results_true.dat b/tests/regression_tests/collision_track/case_4_Nuclide_ID/results_true.dat deleted file mode 100644 index d4d1d1e5ad4..00000000000 --- a/tests/regression_tests/collision_track/case_4_Nuclide_ID/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_5_Universe_ID/collision_track_true.h5 b/tests/regression_tests/collision_track/case_5_Universe_ID/collision_track_true.h5 new file mode 100644 index 00000000000..65419d5ead9 Binary files /dev/null and b/tests/regression_tests/collision_track/case_5_Universe_ID/collision_track_true.h5 differ diff --git a/tests/regression_tests/collision_track/case_5_Universe_ID/inputs_true.dat b/tests/regression_tests/collision_track/case_5_Universe_ID/inputs_true.dat index 8c0d7aa8ee6..c06061210c3 100644 --- a/tests/regression_tests/collision_track/case_5_Universe_ID/inputs_true.dat +++ b/tests/regression_tests/collision_track/case_5_Universe_ID/inputs_true.dat @@ -38,9 +38,9 @@ eigenvalue - 100 + 80 5 - 1 + 4 -2.0 -2.0 -2.0 2.0 2.0 2.0 @@ -52,7 +52,7 @@ 22 77 - 300 + 100 1 diff --git a/tests/regression_tests/collision_track/case_5_Universe_ID/results_true.dat b/tests/regression_tests/collision_track/case_5_Universe_ID/results_true.dat deleted file mode 100644 index d4d1d1e5ad4..00000000000 --- a/tests/regression_tests/collision_track/case_5_Universe_ID/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/collision_track_true.h5 b/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/collision_track_true.h5 new file mode 100644 index 00000000000..507665a60fd Binary files /dev/null and b/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/collision_track_true.h5 differ diff --git a/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/inputs_true.dat b/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/inputs_true.dat index 5173dc35cf1..434f6b434bc 100644 --- a/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/inputs_true.dat +++ b/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/inputs_true.dat @@ -38,9 +38,9 @@ eigenvalue - 100 + 80 5 - 1 + 4 -2.0 -2.0 -2.0 2.0 2.0 2.0 @@ -51,7 +51,7 @@ 550000.0 - 300 + 100 1 diff --git a/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/results_true.dat b/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/results_true.dat deleted file mode 100644 index d4d1d1e5ad4..00000000000 --- a/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_7_all_parameters_used_together/collision_track_true.h5 b/tests/regression_tests/collision_track/case_7_all_parameters_used_together/collision_track_true.h5 new file mode 100644 index 00000000000..b7ff67c48e1 Binary files /dev/null and b/tests/regression_tests/collision_track/case_7_all_parameters_used_together/collision_track_true.h5 differ diff --git a/tests/regression_tests/collision_track/case_7_all_parameters_used_together/inputs_true.dat b/tests/regression_tests/collision_track/case_7_all_parameters_used_together/inputs_true.dat index 005d9feb271..1449c3596fb 100644 --- a/tests/regression_tests/collision_track/case_7_all_parameters_used_together/inputs_true.dat +++ b/tests/regression_tests/collision_track/case_7_all_parameters_used_together/inputs_true.dat @@ -38,9 +38,9 @@ eigenvalue - 100 + 80 5 - 1 + 4 -2.0 -2.0 -2.0 2.0 2.0 2.0 @@ -56,7 +56,7 @@ 1 11 U238 U235 H1 U234 100000.0 - 300 + 100 1 diff --git a/tests/regression_tests/collision_track/case_7_all_parameters_used_together/results_true.dat b/tests/regression_tests/collision_track/case_7_all_parameters_used_together/results_true.dat deleted file mode 100644 index d4d1d1e5ad4..00000000000 --- a/tests/regression_tests/collision_track/case_7_all_parameters_used_together/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_8_2threads/inputs_true.dat b/tests/regression_tests/collision_track/case_8_2threads/inputs_true.dat deleted file mode 100644 index 514932c1a69..00000000000 --- a/tests/regression_tests/collision_track/case_8_2threads/inputs_true.dat +++ /dev/null @@ -1,57 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - eigenvalue - 100 - 5 - 1 - - - -2.0 -2.0 -2.0 2.0 2.0 2.0 - - - true - - - - 200 - - 1 - - diff --git a/tests/regression_tests/collision_track/case_8_2threads/results_true.dat b/tests/regression_tests/collision_track/case_8_2threads/results_true.dat deleted file mode 100644 index d4d1d1e5ad4..00000000000 --- a/tests/regression_tests/collision_track/case_8_2threads/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/test.py b/tests/regression_tests/collision_track/test.py index 00e1e3de410..2f9de9ac00e 100644 --- a/tests/regression_tests/collision_track/test.py +++ b/tests/regression_tests/collision_track/test.py @@ -59,28 +59,13 @@ """ -import os - import openmc -import openmc.lib import pytest from tests.testing_harness import CollisionTrackTestHarness from tests.regression_tests import config -@pytest.fixture(scope="function") -def two_threads(monkeypatch): - """Set the number of OMP threads to 2 for the test.""" - monkeypatch.setenv("OMP_NUM_THREADS", "2") - - -@pytest.fixture(scope="function") -def single_process(monkeypatch): - """Set the number of MPI process to 1 for the test.""" - monkeypatch.setitem(config, "mpi_np", "1") - - @pytest.fixture(scope="module") def model_1(): """Cylindrical core contained in a first box which is contained in a larger box. @@ -181,9 +166,9 @@ def model_1(): # ============================================================================= model.settings = openmc.Settings() - model.settings.particles = 100 + model.settings.particles = 80 model.settings.batches = 5 - model.settings.inactive = 1 + model.settings.inactive = 4 model.settings.seed = 1 bounds = [ @@ -203,19 +188,19 @@ def model_1(): @pytest.mark.parametrize( "folder, model_name, parameter", - [("case_1_Reactions", "model_1", {"max_collisions": 300, "reactions": ["(n,fission)", 101]}), + [("case_1_Reactions", "model_1", {"max_collisions": 100, "reactions": ["(n,fission)", 101]}), ("case_2_Cell_ID", "model_1", { - "max_collisions": 300, "cell_ids": [22]}), + "max_collisions": 100, "cell_ids": [22]}), ("case_3_Material_ID", "model_1", { - "max_collisions": 300, "material_ids": [1]}), + "max_collisions": 100, "material_ids": [1]}), ("case_4_Nuclide_ID", "model_1", { - "max_collisions": 300, "nuclides": ["O16", "U235"]}), + "max_collisions": 100, "nuclides": ["O16", "U235"]}), ("case_5_Universe_ID", "model_1", { - "max_collisions": 300, "cell_ids": [22], "universe_ids": [77]}), + "max_collisions": 100, "cell_ids": [22], "universe_ids": [77]}), ("case_6_deposited_energy_threshold", "model_1", { - "max_collisions": 300, "deposited_E_threshold": 5.5e5}), + "max_collisions": 100, "deposited_E_threshold": 5.5e5}), ("case_7_all_parameters_used_together", "model_1", { - "max_collisions": 300, + "max_collisions": 100, "reactions": ["elastic", 18, "(n,disappear)"], "material_ids": [1, 11], "universe_ids": [77], @@ -235,21 +220,3 @@ def test_collision_track_several_cases( "statepoint.5.h5", model=model, workdir=folder ) harness.main() - - -@pytest.mark.skipif(config["event"], reason="Results from history-based mode.") -def test_collision_track_2threads(model_1, two_threads, single_process): - # This test checks that the `max_collisions` setting is honored: - # no collisions beyond the specified limit should be recorded. - # - # For the result to be reproducible, the number of threads and - # the transport mode (history vs. event) must remain fixed. - assert os.environ["OMP_NUM_THREADS"] == "2" - assert config["mpi_np"] == "1" - model_1.settings.collision_track = { - "max_collisions": 200 - } - harness = CollisionTrackTestHarness( - "statepoint.5.h5", model=model_1, workdir="case_8_2threads" - ) - harness.main() diff --git a/tests/regression_tests/mg_temperature/build_2g.py b/tests/regression_tests/mg_temperature/build_2g.py index 1fb7234499b..bca87364e4b 100644 --- a/tests/regression_tests/mg_temperature/build_2g.py +++ b/tests/regression_tests/mg_temperature/build_2g.py @@ -227,7 +227,7 @@ def analytical_solution_2g_therm(xsmin, xsmax=None, wgt=1.0): L = np.array([sa[0] + ss12, 0.0, -ss12, sa[1]]).reshape(2, 2) Q = np.array([nsf[0], nsf[1], 0.0, 0.0]).reshape(2, 2) arr = np.linalg.inv(L).dot(Q) - return np.amax(np.linalg.eigvals(arr)) + return np.amax(np.linalg.eigvals(arr).real) def build_inf_model(xsnames, xslibname, temperature, tempmethod='nearest'): diff --git a/tests/regression_tests/random_ray_auto_convert/infinite_medium/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert/infinite_medium/inputs_true.dat index 86d5ec4abd5..81f8c98cac7 100644 --- a/tests/regression_tests/random_ray_auto_convert/infinite_medium/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert/infinite_medium/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert/material_wise/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert/material_wise/inputs_true.dat index 86d5ec4abd5..81f8c98cac7 100644 --- a/tests/regression_tests/random_ray_auto_convert/material_wise/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert/material_wise/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert/stochastic_slab/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert/stochastic_slab/inputs_true.dat index 86d5ec4abd5..81f8c98cac7 100644 --- a/tests/regression_tests/random_ray_auto_convert/stochastic_slab/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert/stochastic_slab/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat index b00935ef38a..48b0e8256f8 100644 --- a/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat index 472406fa882..be738c3e05d 100644 --- a/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat index 472406fa882..be738c3e05d 100644 --- a/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/model/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/model/inputs_true.dat index 15981f7fa5f..d2d8289ff2e 100644 --- a/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/model/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/model/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/user/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/user/inputs_true.dat index 86d5ec4abd5..81f8c98cac7 100644 --- a/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/user/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/user/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/model/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/model/inputs_true.dat index 15981f7fa5f..d2d8289ff2e 100644 --- a/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/model/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/model/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/user/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/user/inputs_true.dat index 86d5ec4abd5..81f8c98cac7 100644 --- a/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/user/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/user/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert_temperature/infinite_medium/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_temperature/infinite_medium/inputs_true.dat index c60e6a04199..c49020d558a 100644 --- a/tests/regression_tests/random_ray_auto_convert_temperature/infinite_medium/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_temperature/infinite_medium/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert_temperature/material_wise/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_temperature/material_wise/inputs_true.dat index c60e6a04199..c49020d558a 100644 --- a/tests/regression_tests/random_ray_auto_convert_temperature/material_wise/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_temperature/material_wise/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_auto_convert_temperature/stochastic_slab/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_temperature/stochastic_slab/inputs_true.dat index c60e6a04199..c49020d558a 100644 --- a/tests/regression_tests/random_ray_auto_convert_temperature/stochastic_slab/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_temperature/stochastic_slab/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/regression_tests/random_ray_diagonal_stabilization/inputs_true.dat b/tests/regression_tests/random_ray_diagonal_stabilization/inputs_true.dat index 11100e88e12..0ea8c017760 100644 --- a/tests/regression_tests/random_ray_diagonal_stabilization/inputs_true.dat +++ b/tests/regression_tests/random_ray_diagonal_stabilization/inputs_true.dat @@ -2,17 +2,17 @@ mgxs.h5 - + - + - + - + - + diff --git a/tests/testing_harness.py b/tests/testing_harness.py index 1ad91b7a89f..8d156bd6472 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -546,19 +546,19 @@ def __init__(self, statepoint_name, model=None, inputs_true=None, workdir=None): def _test_output_created(self): """Make sure collision_track.h5 has also been created.""" - super()._test_output_created() if self._model.settings.collision_track: assert os.path.exists( "collision_track.h5" ), "collision_track file has not been created." - def _compare_output(self): + def _compare_results(self): """Compare collision_track.h5 files.""" if self._model.settings.collision_track: collision_track_true = self._return_collision_track_data( "collision_track_true.h5") collision_track_test = self._return_collision_track_data( "collision_track.h5") + assert collision_track_true.shape == collision_track_test.shape np.testing.assert_allclose( collision_track_true, collision_track_test, rtol=1e-07) @@ -582,15 +582,18 @@ def build_inputs(self): def _overwrite_results(self): """Also add the 'collision_track.h5' file during overwriting.""" - super()._overwrite_results() if os.path.exists("collision_track.h5"): shutil.copyfile("collision_track.h5", "collision_track_true.h5") + def _write_results(self, results_string): + # The result file for this test are written by the OpenMC executable itself + pass + @staticmethod def _return_collision_track_data(filepath): """ - Read a collision_track file and return a sorted array composed - of flatten arrays of collision information. + Read a collision_track file and return a sorted array composed of + flattened collision records. Parameters ---------- @@ -600,42 +603,58 @@ def _return_collision_track_data(filepath): Returns ------- data : np.array - Sorted array composed of flatten arrays of collision_track data for - each collision information + Sorted array composed of flattened collision-track records. """ - data = [] - keys = [] - - # Read source file source = openmc.read_collision_track_file(filepath) - for src in source: - r = src['r'] - u = src['u'] - e = src['E'] - de = src['dE'] - time = src['time'] - wgt = src['wgt'] - delayed_group = src['delayed_group'] - cell_id = src['cell_id'] - nuclide_id = src['nuclide_id'] - material_id = src['material_id'] - universe_id = src['universe_id'] - n_collision = src['n_collision'] - event_mt = src['event_mt'] - key = ( - f"{r[0]:.10e} {r[1]:.10e} {r[2]:.10e} {u[0]:.10e} {u[1]:.10e} {u[2]:.10e}" - f"{e:.10e} {de:.10e} {time:.10e} {wgt:.10e} {event_mt} {delayed_group} {cell_id}" - f"{nuclide_id} {material_id} {universe_id} {n_collision} " - ) - keys.append(key) - values = [*r, *u, e, de, time, wgt, event_mt, - delayed_group, cell_id, nuclide_id, material_id, - universe_id, n_collision] - assert len(values) == 17 - data.append(values) - - data = np.array(data) - keys = np.array(keys) - sorted_idx = np.argsort(keys, kind='stable') + columns = [ + source['r']['x'], + source['r']['y'], + source['r']['z'], + source['u']['x'], + source['u']['y'], + source['u']['z'], + source['E'], + source['dE'], + source['time'], + source['wgt'], + source['event_mt'], + source['delayed_group'], + source['cell_id'], + source['nuclide_id'], + source['material_id'], + source['universe_id'], + source['n_collision'], + source['particle'], + source['parent_id'], + source['progeny_id'], + ] + data = np.column_stack(columns) + + # Sort by the complete record, prioritizing stable integer identifiers + # before floating-point fields. This removes dependence on the order in + # which threads append otherwise reproducible collision records. + sort_columns = [ + source['parent_id'], + source['progeny_id'], + source['n_collision'], + source['particle'], + source['cell_id'], + source['material_id'], + source['universe_id'], + source['nuclide_id'], + source['event_mt'], + source['delayed_group'], + source['r']['x'], + source['r']['y'], + source['r']['z'], + source['u']['x'], + source['u']['y'], + source['u']['z'], + source['E'], + source['dE'], + source['time'], + source['wgt'], + ] + sorted_idx = np.lexsort(tuple(reversed(sort_columns))) return data[sorted_idx] diff --git a/tests/unit_tests/test_collision_track.py b/tests/unit_tests/test_collision_track.py index 25344a3bd0c..96676e2cd5c 100644 --- a/tests/unit_tests/test_collision_track.py +++ b/tests/unit_tests/test_collision_track.py @@ -34,6 +34,7 @@ def geometry(): {"max_collisions": 200, "mcpl": True} ], + ids=str ) def test_xml_serialization(parameter, run_in_tmpdir): """Check that the different use cases can be written and read in XML.""" @@ -45,7 +46,7 @@ def test_xml_serialization(parameter, run_in_tmpdir): assert read_settings.collision_track == parameter -@pytest.fixture(scope="module") +@pytest.fixture def model(): """Simple hydrogen sphere divided in two hemispheres by a z-plane to form 2 cells.""" @@ -127,3 +128,43 @@ def test_format_similarity(run_in_tmpdir, model): np.testing.assert_allclose(data_h5, data_mcpl, rtol=1e-05) # tolerance not that low due to the strings that is saved in MCPL, # not enough precision! + + +def test_photon_particles(run_in_tmpdir, model): + """Test that the collision track can be used to track photon particles.""" + model.settings.collision_track = {"max_collisions": 200, "cell_ids": [1, 2]} + + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box(*model.geometry.bounding_box), + energy=openmc.stats.delta_function(1e5), + particle='photon' + ) + model.run() + + with h5py.File("collision_track.h5", "r") as f: + source = f["collision_track_bank"] + + assert len(source) < 200 + + allowed_particles = (openmc.ParticleType.PHOTON, openmc.ParticleType.ELECTRON) + + for point in source: + particle_type = openmc.ParticleType(point['particle']) + assert particle_type in allowed_particles + + if particle_type == openmc.ParticleType.ELECTRON: + assert point['nuclide_id'] == 0 + + +def test_collision_track_two_threads(model, run_in_tmpdir): + # This test checks that the `max_collisions` setting is honored: + # no collisions beyond the specified limit should be recorded. + # + # The exact set of events in the capped bank is not reproducible with + # multiple threads because the bank stores whichever thread appends first + # until capacity is reached. + model.settings.collision_track = {"max_collisions": 200} + model.run(threads=2, particles=500) + + collision_track = openmc.read_collision_track_hdf5("collision_track.h5") + assert len(collision_track) == 200 diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py index 26529e6ce96..0b1937facd9 100644 --- a/tests/unit_tests/test_deplete_microxs.py +++ b/tests/unit_tests/test_deplete_microxs.py @@ -179,6 +179,66 @@ def capture_run(**kwargs): assert ef.values[0] == pytest.approx(energies[0]) assert ef.values[-1] == pytest.approx(energies[-1]) + +def _simple_model(): + model = openmc.Model() + mat = openmc.Material(components={'H1': 1.0, 'H2': 1.0}, + density=5.0, density_units='g/cm3') + sphere = openmc.Sphere(r=10.0, boundary_type='vacuum') + cell = openmc.Cell(region=-sphere, fill=mat) + model.geometry = openmc.Geometry([cell]) + model.settings.particles = 100 + model.settings.batches = 5 + model.settings.run_mode = 'fixed source' + return model, mat + + +def test_hybrid_tally_defaults_to_all_nuclides(run_in_tmpdir): + energies = [0., 0.625, 2.0e7] + kwargs = { + 'nuclides': ['H1', 'H2'], + 'reactions': ['(n,2n)', '(n,gamma)'], + 'energies': energies, + 'reaction_rate_mode': 'flux', + 'chain_file': CHAIN_FILE, + } + + model, mat = _simple_model() + default_fluxes, default_micros = get_microxs_and_flux( + model, [mat], reaction_rate_opts={'reactions': ['(n,2n)']}, **kwargs + ) + + model, mat = _simple_model() + explicit_fluxes, explicit_micros = get_microxs_and_flux( + model, [mat], + reaction_rate_opts={ + 'nuclides': ['H1', 'H2'], + 'reactions': ['(n,2n)'] + }, + **kwargs + ) + + np.testing.assert_allclose(default_fluxes[0], explicit_fluxes[0]) + np.testing.assert_allclose(default_micros[0].data, explicit_micros[0].data) + assert default_micros[0].nuclides == explicit_micros[0].nuclides + assert default_micros[0].reactions == explicit_micros[0].reactions + + +def test_flux_mode_returns_one_group_flux(run_in_tmpdir): + model, mat = _simple_model() + fluxes, micros = get_microxs_and_flux( + model, [mat], + nuclides=['H1'], + reactions=['(n,2n)'], + energies=[0., 0.625, 2.0e7], + reaction_rate_mode='flux', + chain_file=CHAIN_FILE, + ) + + assert fluxes[0].shape == (1,) + assert micros[0].data.shape == (1, 1, 1) + assert fluxes[0][0] > 0.0 + # --------------------------------------------------------------------------- # Tests for MicroXS.merge() # --------------------------------------------------------------------------- diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 6926bc11b18..f62306e834d 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -895,22 +895,23 @@ def test_load_nuclide(lib_init): openmc.lib.load_nuclide('Pu3') +class LegacySlicePlot: + origin = (0.0, 0.0, 0.0) + width = 1.26 + height = 1.26 + basis = 'xy' + h_res = 3 + v_res = 3 + level = -1 + + def test_id_map(lib_init): expected_ids = np.array([[(3, 0, 3), (2, 0, 2), (3, 0, 3)], [(2, 0, 2), (1, 0, 1), (2, 0, 2)], [(3, 0, 3), (2, 0, 2), (3, 0, 3)]], dtype='int32') - # create a plot object - s = openmc.lib.plot._PlotBase() - s.width = 1.26 - s.height = 1.26 - s.v_res = 3 - s.h_res = 3 - s.origin = (0.0, 0.0, 0.0) - s.basis = 'xy' - s.level = -1 - - ids = openmc.lib.plot.id_map(s) + with pytest.warns(FutureWarning, match="deprecated"): + ids = openmc.lib.id_map(LegacySlicePlot()) assert np.array_equal(expected_ids, ids) @@ -920,17 +921,8 @@ def test_property_map(lib_init): [ (293.6, 6.55), (293.6, 10.29769), (293.6, 6.55)], [(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)]], dtype='float') - # create a plot object - s = openmc.lib.plot._PlotBase() - s.width = 1.26 - s.height = 1.26 - s.v_res = 3 - s.h_res = 3 - s.origin = (0.0, 0.0, 0.0) - s.basis = 'xy' - s.level = -1 - - properties = openmc.lib.plot.property_map(s) + with pytest.warns(FutureWarning, match="deprecated"): + properties = openmc.lib.property_map(LegacySlicePlot()) assert np.allclose(expected_properties, properties, atol=1e-04) diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 9234b2d2721..028a83b0bdc 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -1038,3 +1038,32 @@ def test_id_map_to_rgb(): ) # Check that overlap region is green assert np.allclose(rgb_overlap[5:, 5:], [0.0, 1.0, 0.0]) + + +def test_convert_to_multigroup_preserves_material_names(run_in_tmpdir): + """convert_to_multigroup leaves the user's material names unchanged and keys + the MGXS library by a unique sanitised name + id, so distinct materials that + share a name do not collapse to a single cross section.""" + a = openmc.Material(name="Steel Plate #1") + a.add_element("Fe", 1.0) + a.set_density("g/cm3", 7.9) + b = openmc.Material(name="Steel Plate #1") # same name, distinct material + b.add_element("Fe", 1.0) + b.set_density("g/cm3", 7.9) + + s1 = openmc.Sphere(r=1.0) + s2 = openmc.Sphere(r=2.0, boundary_type="vacuum") + c1 = openmc.Cell(fill=a, region=-s1) + c2 = openmc.Cell(fill=b, region=+s1 & -s2) + model = openmc.Model(openmc.Geometry([c1, c2]), openmc.Materials([a, b])) + + # Pre-create the library so MGXS generation (and transport) is skipped. + Path("mgxs.h5").touch() + model.convert_to_multigroup(method="material_wise", mgxs_path="mgxs.h5") + + # User names are preserved, not sanitised or de-duplicated. + assert [m.name for m in model.materials] == ["Steel Plate #1", "Steel Plate #1"] + # Each material reads a unique, sanitised library entry (name + id). + macro = [m._macroscopic for m in model.materials] + assert macro == [f"Steel_Plate__1_{a.id}", f"Steel_Plate__1_{b.id}"] + assert len(set(macro)) == 2 diff --git a/tests/unit_tests/test_slice_data.py b/tests/unit_tests/test_slice_data.py new file mode 100644 index 00000000000..cc5fb047473 --- /dev/null +++ b/tests/unit_tests/test_slice_data.py @@ -0,0 +1,168 @@ +import numpy as np +import openmc +from openmc.examples import pwr_pin_cell + + +def test_slice_data_basic(run_in_tmpdir): + """Test basic slice_data functionality.""" + model = pwr_pin_cell() + geom_data, prop_data = model.slice_data( + origin=(0, 0, 0), + width=(1.0, 1.0), + pixels=(100, 100), + basis='xy' + ) + + # Without filter, should have 3 fields + assert geom_data.shape == (100, 100, 3) + assert geom_data.dtype == np.int32 + assert prop_data.shape == (100, 100, 2) + assert prop_data.dtype == np.float64 + + # Check we have valid geometry + assert np.any(geom_data[:, :, 0] >= 0) # Valid cell IDs + assert np.any(prop_data[:, :, 0] > 0) # Valid temperatures + + +def test_slice_data_no_properties(run_in_tmpdir): + """Test slice_data without property data.""" + model = pwr_pin_cell() + geom_data, prop_data = model.slice_data( + origin=(0, 0, 0), + width=(1.0, 1.0), + pixels=(50, 50), + include_properties=False + ) + + # Without filter, should have 3 fields + assert geom_data.shape == (50, 50, 3) + assert prop_data is None + + +def test_slice_data_with_filter(run_in_tmpdir): + """Test slice_data with a cell filter.""" + model = pwr_pin_cell() + cell_ids = [c.id for c in model.geometry.get_all_cells().values()] + cell_filter = openmc.CellFilter(cell_ids) + + geom_data, _ = model.slice_data( + origin=(0, 0, 0), + width=(1.0, 1.0), + pixels=(50, 50), + filter=cell_filter, + include_properties=False + ) + + # With filter, should have 4 fields + assert geom_data.shape == (50, 50, 4) + + # Filter bin index should be populated where cells exist + filter_bins = geom_data[:, :, 3] + valid_cells = geom_data[:, :, 0] >= 0 + assert np.any(filter_bins[valid_cells] >= 0) + + +def test_slice_data_overlaps(run_in_tmpdir): + """Test slice_data with overlap detection.""" + model = pwr_pin_cell() + geom_data, _ = model.slice_data( + origin=(0, 0, 0), + width=(1.0, 1.0), + pixels=(50, 50), + show_overlaps=True, + include_properties=False + ) + + # Without filter, should have 3 fields + assert geom_data.shape == (50, 50, 3) + # Check for overlap markers (-3) if any exist + # Note: This test may pass without finding overlaps if geometry is correct + + +def test_slice_data_overlaps_with_filter(run_in_tmpdir): + """Test that overlaps don't overwrite filter bin data.""" + model = pwr_pin_cell() + cell_ids = [c.id for c in model.geometry.get_all_cells().values()] + cell_filter = openmc.CellFilter(cell_ids) + + geom_data, _ = model.slice_data( + origin=(0, 0, 0), + width=(1.0, 1.0), + pixels=(50, 50), + filter=cell_filter, + show_overlaps=True, + include_properties=False + ) + + assert geom_data.shape == (50, 50, 4) + + # If any overlaps exist, verify filter bin is still valid (not -3) + overlap_pixels = geom_data[:, :, 0] == -3 + if np.any(overlap_pixels): + # Filter bins at overlap locations should NOT be -3 + filter_bins_at_overlaps = geom_data[overlap_pixels, 3] + assert not np.all(filter_bins_at_overlaps == -3), \ + "Filter bins should be preserved even where overlaps are detected" + + +def test_slice_data_different_bases(run_in_tmpdir): + """Test slice_data with different basis planes.""" + model = pwr_pin_cell() + + for basis in ['xy', 'xz', 'yz']: + geom_data, prop_data = model.slice_data( + origin=(0, 0, 0), + width=(1.0, 1.0), + pixels=(25, 25), + basis=basis + ) + + assert geom_data.shape == (25, 25, 3) + assert prop_data.shape == (25, 25, 2) + + +def test_slice_data_oriented_spans(run_in_tmpdir): + """Test slice_data with oriented span vectors.""" + model = pwr_pin_cell() + + geom_data, prop_data = model.slice_data( + origin=(0, 0, 0), + u_span=(1.0, 0.0, 0.0), + v_span=(0.0, 0.0, 1.0), + pixels=(25, 25) + ) + + assert geom_data.shape == (25, 25, 3) + assert prop_data.shape == (25, 25, 2) + + +def test_slice_data_level(run_in_tmpdir): + """Test slice_data with specific universe level.""" + model = pwr_pin_cell() + geom_data, _ = model.slice_data( + origin=(0, 0, 0), + width=(1.0, 1.0), + pixels=(50, 50), + level=0, # Root universe only + include_properties=False + ) + + assert geom_data.shape == (50, 50, 3) + + +def test_id_map_reverted(run_in_tmpdir): + """Test that id_map returns 3D array without filter support.""" + model = pwr_pin_cell() + id_data = model.id_map( + origin=(0, 0, 0), + width=(1.0, 1.0), + pixels=(50, 50), + basis='xy' + ) + + # Should have 3 fields (cell_id, cell_instance, material_id) + assert id_data.shape == (50, 50, 3) + assert id_data.dtype == np.int32 + + # Check valid data + assert np.any(id_data[:, :, 0] >= 0) # Valid cell IDs