diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 911654d318f..cbc73041644 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -125,6 +125,9 @@ int openmc_nuclide_name(int index, const char** name); int openmc_plot_geometry(); int openmc_id_map(const void* slice, int32_t* data_out); int openmc_property_map(const void* slice, double* data_out); +int openmc_slice_plot(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); 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/plot.h b/include/openmc/plot.h index 3813101c65c..b011d217525 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -17,6 +17,8 @@ #include "openmc/particle.h" #include "openmc/position.h" #include "openmc/random_lcg.h" +#include "openmc/tallies/filter.h" +#include "openmc/tallies/filter_match.h" #include "openmc/xml_interface.h" namespace openmc { @@ -147,28 +149,48 @@ 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 - xt::xtensor data_; //!< 2D array of cell & material ids + xt::xtensor data_; //!< 3D array of cell ID, cell instance, + //!< and material ID }; 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 xt::xtensor 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 + xt::xtensor + id_data_; //!< [v_res, h_res, 3 or 4]: cell, instance, mat, [filter_bin] + xt::xtensor + property_data_; //!< [v_res, h_res, 2]: temperature, density + bool include_filter_; //!< Whether filter bin index is included +}; + //=============================================================================== // Plot class //=============================================================================== @@ -176,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,8 +210,10 @@ class SlicePlotBase { // Members public: Position origin_; //!< Plot origin in geometry - Position width_; //!< Plot width in geometry - PlotBasis basis_; //!< Plot basis (XY/XZ/YZ) + Position u_span_; //!< Full-width span vector in geometry + Position v_span_; //!< Full-height span vector in geometry + Position width_; //!< Axis-aligned plot width in geometry + PlotBasis basis_; //!< Plot basis (XY/XZ/YZ) for axis-aligned slices array pixels_; //!< Plot size in pixels bool slice_color_overlaps_; //!< Show overlapping cells? int slice_level_ {-1}; //!< Plot universe level @@ -197,60 +221,55 @@ class SlicePlotBase { }; 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 + Position u_step = u_span_ / static_cast(width); + Position 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); @@ -259,7 +278,7 @@ 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)) { data.set_overlap(y, x); diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py index 90af80d5b76..32aef4f0698 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -58,6 +58,10 @@ class _PlotBase(Structure): ----------------- origin_ : openmc.lib.plot._Position A position defining the origin of the plot. + u_span_ : openmc.lib.plot._Position + Full-width span vector defining the plot's horizontal axis. + v_span_ : openmc.lib.plot._Position + Full-height span vector defining the plot's vertical axis. width_ : openmc.lib.plot._Position The width of the plot along the x, y, and z axes, respectively basis_ : c_int @@ -88,6 +92,8 @@ class _PlotBase(Structure): The universe level for the plot (default: -1 -> all universes shown) """ _fields_ = [('origin_', _Position), + ('u_span_', _Position), + ('v_span_', _Position), ('width_', _Position), ('basis_', c_int), ('pixels_', 3*c_size_t), @@ -98,6 +104,30 @@ def __init__(self): self.level_ = -1 self.basis_ = 1 self.color_overlaps_ = False + self._update_spans() + + def _update_spans(self): + if self.basis_ == 1: + self.u_span_.x = self.width_.x + self.u_span_.y = 0.0 + self.u_span_.z = 0.0 + self.v_span_.x = 0.0 + self.v_span_.y = self.width_.y + self.v_span_.z = 0.0 + elif self.basis_ == 2: + self.u_span_.x = self.width_.x + self.u_span_.y = 0.0 + self.u_span_.z = 0.0 + self.v_span_.x = 0.0 + self.v_span_.y = 0.0 + self.v_span_.z = self.width_.y + elif self.basis_ == 3: + self.u_span_.x = 0.0 + self.u_span_.y = self.width_.x + self.u_span_.z = 0.0 + self.v_span_.x = 0.0 + self.v_span_.y = 0.0 + self.v_span_.z = self.width_.y @property def origin(self): @@ -116,6 +146,7 @@ def width(self): @width.setter def width(self, width): self.width_.x = width + self._update_spans() @property def height(self): @@ -124,6 +155,7 @@ def height(self): @height.setter def height(self, height): self.width_.y = height + self._update_spans() @property def basis(self): @@ -150,6 +182,7 @@ def basis(self, basis): self.basis_ = 2 elif basis == 'yz': self.basis_ = 3 + self._update_spans() return if isinstance(basis, int): @@ -157,6 +190,7 @@ def basis(self, basis): if basis not in valid_bases: raise ValueError(f"{basis} is not a valid plot basis.") self.basis_ = basis + self._update_spans() return raise ValueError(f"{basis} of type {type(basis)} is an invalid plot basis") @@ -231,8 +265,7 @@ def id_map(plot): contains, in order, cell IDs, cell instances, and material IDs. """ - img_data = np.zeros((plot.v_res, plot.h_res, 3), - dtype=np.dtype('int32')) + img_data = np.zeros((plot.v_res, plot.h_res, 3), dtype=np.int32) _dll.openmc_id_map(plot, img_data.ctypes.data_as(POINTER(c_int32))) return img_data @@ -263,6 +296,155 @@ def property_map(plot): _dll.openmc_property_map(plot, prop_data.ctypes.data_as(POINTER(c_double))) return prop_data + +_dll.openmc_slice_plot.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, # color_overlaps + c_int, # level + c_int32, # filter_index + POINTER(c_int32), # geom_data + POINTER(c_double), # property_data (can be None) +] +_dll.openmc_slice_plot.restype = c_int +_dll.openmc_slice_plot.errcheck = _error_handler + + +def slice_plot(origin, width=None, basis='xy', u_span=None, v_span=None, + pixels=None, color_overlaps=False, level=-1, filter=None, + include_properties=True): + """Generate a 2D raster of geometry and property data for plotting. + + Parameters + ---------- + 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] + color_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 + 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): + basis = basis.lower() + if basis not in basis_map: + raise ValueError(f"{basis} is not a valid plot basis.") + 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.") + 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_plot( + origin_arr, + u_span_arr, + v_span_arr, + pixels_arr, + color_overlaps, + level, + filter_index, + geom_data.ctypes.data_as(POINTER(c_int32)), + prop_ptr + ) + + return geom_data, property_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 a512c4d380c..4790359810e 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1102,7 +1102,134 @@ def id_map( init_kwargs.setdefault('args', ['-c']) with openmc.lib.TemporarySession(self, **init_kwargs): - return openmc.lib.id_map(plot_obj) + ids = openmc.lib.id_map(plot_obj) + + return ids + + def slice_plot( + 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, + color_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. + color_overlaps : bool, optional + Whether to 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 + + 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.") + 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 + filter_ids = {f.id for t in self.tallies for f in t.filters} + original_length = len(self.tallies) + if filter is not None and filter.id not in filter_ids: + temp_tally = openmc.Tally() + temp_tally.filters = [filter] + temp_tally.scores = ['flux'] + self.tallies.append(temp_tally) + + with openmc.lib.TemporarySession(self, **init_kwargs): + geom_data, property_data = openmc.lib.slice_plot( + origin=origin, + width=width, + basis=basis, + u_span=u_span, + v_span=v_span, + pixels=pixels, + color_overlaps=color_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( @@ -1697,10 +1824,10 @@ 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, + correction: str | None, + directory: PathLike, ) -> openmc.mgxs.Library: """ Automatically generate a multi-group cross section libray from a model diff --git a/src/plot.cpp b/src/plot.cpp index e3b1e84c80e..607beb84415 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -34,6 +34,7 @@ #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/string_utils.h" +#include "openmc/tallies/filter.h" namespace openmc { @@ -45,10 +46,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) { @@ -65,7 +68,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) xt::view(data_, y, x, xt::all()) = 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())); @@ -1011,6 +1097,8 @@ void Plot::create_voxel() const 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_; @@ -1937,7 +2025,7 @@ extern "C" int openmc_property_map(const void* plot, double* data_out) 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; } @@ -1953,6 +2041,70 @@ extern "C" int openmc_property_map(const void* plot, double* data_out) return 0; } +extern "C" int openmc_slice_plot(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 + Position u_span_pos {u_span[0], u_span[1], u_span[2]}; + Position 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.width_ = Position {u_norm, v_norm, 0.0}; + plot_params.basis_ = SlicePlotBase::PlotBasis::xy; + plot_params.pixels_[0] = pixels[0]; + plot_params.pixels_[1] = pixels[1]; + plot_params.slice_color_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/unit_tests/test_slice_plot.py b/tests/unit_tests/test_slice_plot.py new file mode 100644 index 00000000000..8b1919c4158 --- /dev/null +++ b/tests/unit_tests/test_slice_plot.py @@ -0,0 +1,168 @@ +import numpy as np +import openmc +from openmc.examples import pwr_pin_cell + + +def test_slice_plot_basic(run_in_tmpdir): + """Test basic slice_plot functionality.""" + model = pwr_pin_cell() + geom_data, prop_data = model.slice_plot( + 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_plot_no_properties(run_in_tmpdir): + """Test slice_plot without property data.""" + model = pwr_pin_cell() + geom_data, prop_data = model.slice_plot( + 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_plot_with_filter(run_in_tmpdir): + """Test slice_plot 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_plot( + 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_plot_overlaps(run_in_tmpdir): + """Test slice_plot with overlap detection.""" + model = pwr_pin_cell() + geom_data, _ = model.slice_plot( + origin=(0, 0, 0), + width=(1.0, 1.0), + pixels=(50, 50), + color_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_plot_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_plot( + origin=(0, 0, 0), + width=(1.0, 1.0), + pixels=(50, 50), + filter=cell_filter, + color_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_plot_different_bases(run_in_tmpdir): + """Test slice_plot with different basis planes.""" + model = pwr_pin_cell() + + for basis in ['xy', 'xz', 'yz']: + geom_data, prop_data = model.slice_plot( + 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_plot_oriented_spans(run_in_tmpdir): + """Test slice_plot with oriented span vectors.""" + model = pwr_pin_cell() + + geom_data, prop_data = model.slice_plot( + 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_plot_level(run_in_tmpdir): + """Test slice_plot with specific universe level.""" + model = pwr_pin_cell() + geom_data, _ = model.slice_plot( + 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