From c547b57e94c4a16b69c20d982f8f3388d91c94a1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 24 Apr 2025 15:03:48 -0500 Subject: [PATCH 01/13] Add filter_get_bins interface (not efficient for plotter) --- openmc/lib/filter.py | 49 +++++++++++++++++++++++++++- src/tallies/filter.cpp | 72 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 120 insertions(+), 1 deletion(-) diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index dc81c3487b9..9e64e0731a3 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -1,4 +1,4 @@ -from collections.abc import Mapping +from collections.abc import Mapping, Sequence from ctypes import c_int, c_int32, c_double, c_char_p, POINTER, \ create_string_buffer, c_size_t from weakref import WeakValueDictionary @@ -14,6 +14,7 @@ from .error import _error_handler from .material import Material from .mesh import _get_mesh +from .plot import _Position __all__ = [ @@ -69,6 +70,11 @@ _dll.openmc_filter_set_id.argtypes = [c_int32, c_int32] _dll.openmc_filter_set_id.restype = c_int _dll.openmc_filter_set_id.errcheck = _error_handler +_dll.openmc_filter_get_plot_bins.argtypes = [ + c_int32, _Position, _Position, c_int, POINTER(c_int), POINTER(c_int32) +] +_dll.openmc_filter_get_plot_bins.restype = c_int +_dll.openmc_filter_get_plot_bins.errcheck = _error_handler _dll.openmc_get_filter_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_filter_index.restype = c_int _dll.openmc_get_filter_index.errcheck = _error_handler @@ -155,6 +161,7 @@ _dll.openmc_zernike_filter_set_order.errcheck = _error_handler _dll.tally_filters_size.restype = c_size_t + class Filter(_FortranObjectWithID): __instances = WeakValueDictionary() @@ -203,6 +210,46 @@ def n_bins(self): _dll.openmc_filter_get_num_bins(self._index, n) return n.value + def get_plot_bins( + self, + origin: Sequence[float], + width: Sequence[float], + basis: str, + pixels: Sequence[int] + ) -> np.ndarray: + """Get filter bin indices for a rasterized plot. + + .. versionadded:: 0.15.3 + + Parameters + ---------- + origin : iterable of float + Origin of the plotting view. Should have length 3. + width : iterable of float + Width of the plotting view. Should have length 2. + basis : {'xy', 'xz', 'yz'} + Plotting basis. + pixels : iterable of int + Number of pixels in each direction. Should have length 2. + + Returns + ------- + 2D numpy array with mesh bin indices corresponding to each pixel within + the plotting view. + + """ + origin = _Position(*origin) + width = _Position(*width) + basis = {'xy': 1, 'xz': 2, 'yz': 3}[basis] + pixel_array = (c_int*2)(*pixels) + img_data = np.zeros((pixels[1], pixels[0]), dtype=np.dtype('int32')) + + _dll.openmc_filter_get_plot_bins( + self._index, origin, width, basis, pixel_array, + img_data.ctypes.data_as(POINTER(c_int32)) + ) + return img_data + class EnergyFilter(Filter): filter_type = 'energy' diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 16e4c055021..acfd0fb556e 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -10,6 +10,8 @@ #include "openmc/capi.h" #include "openmc/constants.h" // for MAX_LINE_LEN; #include "openmc/error.h" +#include "openmc/geometry.h" +#include "openmc/position.h" #include "openmc/tallies/filter_azimuthal.h" #include "openmc/tallies/filter_cell.h" #include "openmc/tallies/filter_cell_instance.h" @@ -251,6 +253,76 @@ extern "C" int openmc_filter_get_num_bins(int32_t index, int* n_bins) return 0; } +extern "C" int openmc_filter_get_plot_bins(int32_t index, Position origin, + Position width, int basis, int* pixels, int32_t* data) +{ + if (int err = verify_filter(index)) + return err; + const auto& filter = model::tally_filters[index].get(); + + int pixel_width = pixels[0]; + int pixel_height = pixels[1]; + + // get pixel size + double in_pixel = (width[0]) / static_cast(pixel_width); + double out_pixel = (width[1]) / static_cast(pixel_height); + + // setup basis indices and initial position centered on pixel + int in_i, out_i; + Position xyz = origin; + enum class PlotBasis { xy = 1, xz = 2, yz = 3 }; + PlotBasis basis_enum = static_cast(basis); + switch (basis_enum) { + 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(); + } + + // 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.; + +#pragma omp parallel + { + Particle p; + p.r() = xyz; + p.u() = {1.0, 0.0, 0.0}; + p.coord(0).universe() = model::root_universe; + FilterMatch match; + +#pragma omp for + for (int y = 0; y < pixel_height; y++) { + p.r()[out_i] = xyz[out_i] - out_pixel * y; + for (int x = 0; x < pixel_width; x++) { + p.r()[in_i] = xyz[in_i] + in_pixel * x; + p.n_coord() = 1; + bool found_cell = exhaustive_find_cell(p); + filter->get_all_bins(p, TallyEstimator::COLLISION, match); + if (match.bins_.empty()) { + data[pixel_width * y + x] = -1; + } else { + data[pixel_width * y + x] = match.bins_[0]; + } + match.bins_.clear(); + match.weights_.clear(); + } + } + } + + return 0; +} + extern "C" int openmc_get_filter_index(int32_t id, int32_t* index) { auto it = model::filter_map.find(id); From 7f668bc858933dad2535512ffa2b08d090574c0e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 25 Jul 2025 17:06:11 +0700 Subject: [PATCH 02/13] Combine filter index search with id_map --- include/openmc/plot.h | 25 +++++++--- openmc/lib/plot.py | 18 ++++--- openmc/model/model.py | 21 +++++++- src/plot.cpp | 108 ++++++++++++++++++++++++++---------------- 4 files changed, 115 insertions(+), 57 deletions(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 02ff2848ded..8ee7ed781d5 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 { @@ -148,11 +150,13 @@ struct IdData { IdData(size_t h_res, size_t v_res); // 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, FilterMatch* match); void set_overlap(size_t y, size_t x); // Members - xt::xtensor data_; //!< 2D array of cell & material ids + xt::xtensor data_; //!< 4D array of cell ID, cell instance, + //!< material ID, and filter index }; struct PropertyData { @@ -160,7 +164,8 @@ struct PropertyData { PropertyData(size_t h_res, size_t v_res); // 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, FilterMatch* match); void set_overlap(size_t y, size_t x); // Members @@ -174,7 +179,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 }; @@ -195,12 +200,17 @@ 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]; + Filter* filter = nullptr; + if (filter_index > 0) { + filter = &model::tally_filters[filter_index].get(); + } + // get pixel size double in_pixel = (width_[0]) / static_cast(width); double out_pixel = (width_[1]) / static_cast(height); @@ -237,12 +247,13 @@ T SlicePlotBase::get_map() const #pragma omp parallel { - GeometryState p; + Particle p; p.r() = xyz; 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++) { @@ -257,7 +268,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 68f61821c57..db9874ec0dc 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -204,12 +204,12 @@ def __repr__(self): return '\n'.join(out_str) -_dll.openmc_id_map.argtypes = [POINTER(_PlotBase), POINTER(c_int32)] +_dll.openmc_id_map.argtypes = [POINTER(_PlotBase), c_int32, POINTER(c_int32)] _dll.openmc_id_map.restype = c_int _dll.openmc_id_map.errcheck = _error_handler -def id_map(plot): +def id_map(plot, filter=None): """ Generate a 2-D map of cell and material IDs. Used for in-memory image generation. @@ -218,18 +218,22 @@ def id_map(plot): ---------- plot : openmc.lib.plot._PlotBase Object describing the slice of the model to be generated + filter : openmc.Filter, optional + If provided, the information for each pixel also includes an index in + the filter corresponding to the pixel position. Returns ------- id_map : numpy.ndarray - A NumPy array with shape (vertical pixels, horizontal pixels, 3) of + A NumPy array with shape (vertical pixels, horizontal pixels, 4) of OpenMC property ids with dtype int32. The last dimension of the array - contains, in order, cell IDs, cell instances, and material IDs. + contains, in order, cell IDs, cell instances, material IDs, and filter + indices. """ - 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))) + filter_id = filter.id if filter is not None else -1 + img_data = np.zeros((plot.v_res, plot.h_res, 4), dtype=np.int32) + _dll.openmc_id_map(plot, filter_id, img_data.ctypes.data_as(POINTER(c_int32))) return img_data diff --git a/openmc/model/model.py b/openmc/model/model.py index 82fe87ca5f4..d97d1e49a80 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1038,6 +1038,7 @@ def id_map( pixels: int | Sequence[int] = 40000, basis: str = 'xy', color_overlaps: bool = False, + filter: openmc.Filter | None = None, **init_kwargs ) -> np.ndarray: """Generate an ID map for domains based on the plot parameters @@ -1070,6 +1071,9 @@ def id_map( 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. + filter : openmc.Filter, optional + If provided, the information for each pixel also includes an index + in the filter corresponding to the pixel position. **init_kwargs Keyword arguments passed to :meth:`Model.init_lib`. @@ -1101,8 +1105,23 @@ def id_map( 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): - return openmc.lib.id_map(plot_obj) + ids = openmc.lib.id_map(plot_obj) + + # If filter was temporarily added, remove it + if len(self.tallies) > original_length: + self.tallies.pop() + return ids @add_plot_params def plot( diff --git a/src/plot.cpp b/src/plot.cpp index daeb0bde7c6..fee9847fadf 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 { @@ -48,7 +49,8 @@ constexpr int32_t OVERLAP {-3}; IdData::IdData(size_t h_res, size_t v_res) : 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) { @@ -70,6 +72,18 @@ void IdData::set_value(size_t y, size_t x, const GeometryState& p, int level) Material* m = model::materials.at(p.material()).get(); data_(y, x, 2) = m->id_; } + + // set filter index + if (filter) { + filter->get_all_bins(p, TallyEstimator::COLLISION, *match); + if (match->bins_.empty()) { + data_(y, x, 3) = -1; + } else { + data_(y, x, 3) = match->bins_[0]; + } + match->bins_.clear(); + match->weights_.clear(); + } } void IdData::set_overlap(size_t y, size_t x) @@ -81,8 +95,8 @@ PropertyData::PropertyData(size_t h_res, size_t v_res) : 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; @@ -189,9 +203,10 @@ void Plot::print_info() const void read_plots_xml() { - // Check if plots.xml exists; this is only necessary when the plot runmode is - // initiated. Otherwise, we want to read plots.xml because it may be called - // later via the API. In that case, its ok for a plots.xml to not exist + // Check if plots.xml exists; this is only necessary when the plot runmode + // is initiated. Otherwise, we want to read plots.xml because it may be + // called later via the API. In that case, its ok for a plots.xml to not + // exist std::string filename = settings::path_input + "plots.xml"; if (!file_exists(filename) && settings::run_mode == RunMode::PLOTTING) { fatal_error(fmt::format("Plots XML file '{}' does not exist!", filename)); @@ -929,14 +944,15 @@ void Plot::draw_mesh_lines(ImageData& data) const } /* outputs a binary file that can be input into silomesh for 3D geometry - * visualization. It works the same way as create_image by dragging a particle - * across the geometry for the specified number of voxels. The first 3 int's in - * the binary are the number of x, y, and z voxels. The next 3 double's are - * the widths of the voxels in the x, y, and z directions. The next 3 double's - * are the x, y, and z coordinates of the lower left point. Finally the binary - * is filled with entries of four int's each. Each 'row' in the binary contains - * four int's: 3 for x,y,z position and 1 for cell or material id. For 1 - * million voxels this produces a file of approximately 15MB. + * visualization. It works the same way as create_image by dragging a + * particle across the geometry for the specified number of voxels. The first + * 3 int's in the binary are the number of x, y, and z voxels. The next 3 + * double's are the widths of the voxels in the x, y, and z directions. The + * next 3 double's are the x, y, and z coordinates of the lower left point. + * Finally the binary is filled with entries of four int's each. Each 'row' in + * the binary contains four int's: 3 for x,y,z position and 1 for cell or + * material id. For 1 million voxels this produces a file of approximately + * 15MB. */ void Plot::create_voxel() const { @@ -1182,8 +1198,9 @@ bool WireframeRayTracePlot::trackstack_equivalent( return false; // Check if neighboring cells are different - // if (track1[t1_i ? t1_i - 1 : 0].id != track2[t2_i ? t2_i - 1 : 0].id) - // return false; if (track1[t1_i < track1.size() - 1 ? t1_i + 1 : t1_i + // if (track1[t1_i ? t1_i - 1 : 0].id != track2[t2_i ? t2_i - 1 : + // 0].id) return false; if (track1[t1_i < track1.size() - 1 ? t1_i + 1 + // : t1_i // ].id != // track2[t2_i < track2.size() - 1 ? t2_i + 1 : t2_i].id) return // false; @@ -1243,8 +1260,8 @@ ImageData WireframeRayTracePlot::create_image() const size_t height = pixels()[1]; ImageData data({width, height}, not_found_); - // This array marks where the initial wireframe was drawn. We convolve it with - // a filter that gets adjusted with the wireframe thickness in order to + // This array marks where the initial wireframe was drawn. We convolve it + // with a filter that gets adjusted with the wireframe thickness in order to // thicken the lines. xt::xtensor wireframe_initial({width, height}, 0); @@ -1278,9 +1295,9 @@ ImageData WireframeRayTracePlot::create_image() const for (int iter = 0; iter <= pixels()[1] / n_threads; iter++) { // Save bottom line of current work chunk to compare against later. This - // used to be inside the below if block, but it causes a spurious line to - // be drawn at the bottom of the image. Not sure why, but moving it here - // fixes things. + // used to be inside the below if block, but it causes a spurious line + // to be drawn at the bottom of the image. Not sure why, but moving it + // here fixes things. if (tid == n_threads - 1) old_segments = this_line_segments[n_threads - 1]; @@ -1305,8 +1322,8 @@ ImageData WireframeRayTracePlot::create_image() const // There must be at least two cell intersections to color, front and // back of the cell. Maybe an infinitely thick cell could be present - // with no back, but why would you want to color that? It's easier to - // just skip that edge case and not even color it. + // with no back, but why would you want to color that? It's easier + // to just skip that edge case and not even color it. if (segments.size() <= 1) continue; @@ -1341,8 +1358,8 @@ ImageData WireframeRayTracePlot::create_image() const } } // end "if" vert in correct range - // We require a barrier before comparing vertical neighbors' intersection - // stacks. i.e. all threads must be done with their line. + // We require a barrier before comparing vertical neighbors' + // intersection stacks. i.e. all threads must be done with their line. #pragma omp barrier // Now that the horizontal line has finished rendering, we can fill in @@ -1366,8 +1383,8 @@ ImageData WireframeRayTracePlot::create_image() const } } - // We need another barrier to ensure threads don't proceed to modify their - // intersection stacks on that horizontal line while others are + // We need another barrier to ensure threads don't proceed to modify + // their intersection stacks on that horizontal line while others are // potentially still working on the above. #pragma omp barrier vert += n_threads; @@ -1633,7 +1650,8 @@ void Ray::trace() // radiation transport model. // // After phase one is done, we can starting tracing from cell to cell within - // the model. This step can use neighbor lists to accelerate the ray tracing. + // the model. This step can use neighbor lists to accelerate the ray + // tracing. // Attempt to initialize the particle. We may have to enter a loop to move // it up to the edge of the model. @@ -1731,13 +1749,13 @@ void Ray::trace() inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10); // Call the specialized logic for this type of ray. Note that we do not - // call this if the advance distance is very small. Unfortunately, it seems - // darn near impossible to get the particle advanced to the model boundary - // and through it without sometimes accidentally calling on_intersection - // twice. This incorrectly shades the region as occluded when it might not - // actually be. By screening out intersection distances smaller than a - // threshold 10x larger than the scoot distance used to advance up to the - // model boundary, we can avoid that situation. + // call this if the advance distance is very small. Unfortunately, it + // seems darn near impossible to get the particle advanced to the model + // boundary and through it without sometimes accidentally calling + // on_intersection twice. This incorrectly shades the region as occluded + // when it might not actually be. By screening out intersection distances + // smaller than a threshold 10x larger than the scoot distance used to + // advance up to the model boundary, we can avoid that situation. if (call_on_intersection) { on_intersection(); if (stop_) @@ -1800,8 +1818,8 @@ void PhongRay::on_intersection() // TODO // Not sure what can cause a surface token to be invalid here, although it // sometimes happens for a few pixels. It's very very rare, so proceed by - // coloring the pixel with the overlap color. It seems to happen only for a - // few pixels on the outer boundary of a hex lattice. + // coloring the pixel with the overlap color. It seems to happen only for + // a few pixels on the outer boundary of a hex lattice. // // We cannot detect it in the outer loop, and it only matters here, so // that's why the error handling is a little different than for a lost @@ -1872,9 +1890,9 @@ void PhongRay::on_intersection() compute_distance(); } else { - // If it's not facing the light, we color with the diffuse contribution, so - // next we check if we're going to occlude the last reflected surface. if - // so, color by the diffuse contribution instead + // If it's not facing the light, we color with the diffuse contribution, + // so next we check if we're going to occlude the last reflected surface. + // if so, color by the diffuse contribution instead if (orig_hit_id_ == -1) fatal_error("somehow a ray got reflected but not original ID set?"); @@ -1885,7 +1903,8 @@ void PhongRay::on_intersection() } } -extern "C" int openmc_id_map(const void* plot, int32_t* data_out) +extern "C" int openmc_id_map( + const void* plot, int32_t filter_index, int32_t* data_out) { auto plt = reinterpret_cast(plot); @@ -1894,11 +1913,16 @@ extern "C" int openmc_id_map(const void* plot, int32_t* data_out) return OPENMC_E_INVALID_ARGUMENT; } + if (filter_index >= 0) { + if (int err = verify_filter(filter_index)) + return err; + } + if (plt->slice_color_overlaps_ && model::overlap_check_count.size() == 0) { model::overlap_check_count.resize(model::cells.size()); } - auto ids = plt->get_map(); + auto ids = plt->get_map(filter_index); // write id data to array std::copy(ids.data_.begin(), ids.data_.end(), data_out); From 2287926ba062f3d2d626e266d87f1329e1161fd7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 31 Jan 2026 16:51:17 -0600 Subject: [PATCH 03/13] Implement new openmc_raster_plot function --- include/openmc/capi.h | 3 + include/openmc/plot.h | 40 +++++--- openmc/lib/plot.py | 109 ++++++++++++++++++-- openmc/model/model.py | 95 +++++++++++++++-- src/plot.cpp | 234 +++++++++++++++++++++++++++++++++++++----- 5 files changed, 421 insertions(+), 60 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 533ebb65a49..820612186a9 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_raster_plot(const double origin[3], const double width[2], int basis, + const size_t pixels[2], bool color_overlaps, int level, int32_t filter_index, + int32_t* geom_data, double* property_data); int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x, int* nx, double** grid_y, int* ny, double** grid_z, int* nz); int openmc_rectilinear_mesh_set_grid(int32_t index, const double* grid_x, diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 8ee7ed781d5..512979c0e18 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -150,26 +150,41 @@ struct IdData { IdData(size_t h_res, size_t v_res); // Methods - void set_value(size_t y, size_t x, const Particle& p, int level, - Filter* filter, FilterMatch* match); + void set_value(size_t y, size_t x, const Particle& p, int level); void set_overlap(size_t y, size_t x); // Members - xt::xtensor data_; //!< 4D array of cell ID, cell instance, - //!< material ID, and filter index + 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); + // Methods + void set_value(size_t y, size_t x, const Particle& p, int level); + 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); + // Methods void set_value(size_t y, size_t x, const Particle& p, int level, Filter* filter, FilterMatch* match); void set_overlap(size_t y, size_t x); // Members - xt::xtensor data_; //!< 2D array of temperature & density data + 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 }; //=============================================================================== @@ -179,7 +194,10 @@ struct PropertyData { class SlicePlotBase { public: template - T get_map(int32_t filter_index = -1) const; + T get_map() const; + + // Specialized method for RasterData that takes filter_index + RasterData get_raster_map(int32_t filter_index) const; enum class PlotBasis { xy = 1, xz = 2, yz = 3 }; @@ -200,17 +218,12 @@ class SlicePlotBase { }; template -T SlicePlotBase::get_map(int32_t filter_index) const +T SlicePlotBase::get_map() const { size_t width = pixels_[0]; size_t height = pixels_[1]; - Filter* filter = nullptr; - if (filter_index > 0) { - filter = &model::tally_filters[filter_index].get(); - } - // get pixel size double in_pixel = (width_[0]) / static_cast(width); double out_pixel = (width_[1]) / static_cast(height); @@ -253,7 +266,6 @@ T SlicePlotBase::get_map(int32_t filter_index) const 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++) { @@ -268,7 +280,7 @@ T SlicePlotBase::get_map(int32_t filter_index) const j = level; } if (found_cell) { - data.set_value(y, x, p, j, filter, &match); + data.set_value(y, x, p, j); } 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 db9874ec0dc..c142ed164fc 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -204,12 +204,12 @@ def __repr__(self): return '\n'.join(out_str) -_dll.openmc_id_map.argtypes = [POINTER(_PlotBase), c_int32, POINTER(c_int32)] +_dll.openmc_id_map.argtypes = [POINTER(_PlotBase), POINTER(c_int32)] _dll.openmc_id_map.restype = c_int _dll.openmc_id_map.errcheck = _error_handler -def id_map(plot, filter=None): +def id_map(plot): """ Generate a 2-D map of cell and material IDs. Used for in-memory image generation. @@ -218,22 +218,17 @@ def id_map(plot, filter=None): ---------- plot : openmc.lib.plot._PlotBase Object describing the slice of the model to be generated - filter : openmc.Filter, optional - If provided, the information for each pixel also includes an index in - the filter corresponding to the pixel position. Returns ------- id_map : numpy.ndarray - A NumPy array with shape (vertical pixels, horizontal pixels, 4) of + 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, material IDs, and filter - indices. + contains, in order, cell IDs, cell instances, and material IDs. """ - filter_id = filter.id if filter is not None else -1 - img_data = np.zeros((plot.v_res, plot.h_res, 4), dtype=np.int32) - _dll.openmc_id_map(plot, filter_id, img_data.ctypes.data_as(POINTER(c_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 @@ -262,3 +257,95 @@ def property_map(plot): prop_data = np.zeros((plot.v_res, plot.h_res, 2)) _dll.openmc_property_map(plot, prop_data.ctypes.data_as(POINTER(c_double))) return prop_data + + +_dll.openmc_raster_plot.argtypes = [ + POINTER(c_double * 3), # origin + POINTER(c_double * 2), # width + c_int, # basis + 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_raster_plot.restype = c_int +_dll.openmc_raster_plot.errcheck = _error_handler + + +def raster_plot(origin, width, basis, pixels, 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] + basis : {'xy', 'xz', 'yz'} or int + Plot basis + 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.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 + """ + # Convert basis string to int + basis_map = {'xy': 1, 'xz': 2, 'yz': 3} + if isinstance(basis, str): + basis = basis_map[basis.lower()] + + # Prepare ctypes arrays + origin_arr = (c_double * 3)(*origin) + width_arr = (c_double * 2)(*width) + 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_raster_plot( + origin_arr, + width_arr, + basis, + pixels_arr, + color_overlaps, + level, + filter_index, + geom_data.ctypes.data_as(POINTER(c_int32)), + prop_ptr + ) + + return geom_data, property_data diff --git a/openmc/model/model.py b/openmc/model/model.py index d97d1e49a80..4fdfca11f6d 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1038,7 +1038,6 @@ def id_map( pixels: int | Sequence[int] = 40000, basis: str = 'xy', color_overlaps: bool = False, - filter: openmc.Filter | None = None, **init_kwargs ) -> np.ndarray: """Generate an ID map for domains based on the plot parameters @@ -1071,9 +1070,6 @@ def id_map( 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. - filter : openmc.Filter, optional - If provided, the information for each pixel also includes an index - in the filter corresponding to the pixel position. **init_kwargs Keyword arguments passed to :meth:`Model.init_lib`. @@ -1105,6 +1101,83 @@ def id_map( init_kwargs.setdefault('output', False) init_kwargs.setdefault('args', ['-c']) + with openmc.lib.TemporarySession(self, **init_kwargs): + ids = openmc.lib.id_map(plot_obj) + + return ids + + def raster_plot( + self, + origin: Sequence[float] | None = None, + width: Sequence[float] | None = None, + pixels: int | Sequence[int] = 40000, + basis: str = 'xy', + 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. + 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 + + 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} @@ -1116,12 +1189,16 @@ def id_map( self.tallies.append(temp_tally) with openmc.lib.TemporarySession(self, **init_kwargs): - ids = openmc.lib.id_map(plot_obj) + geom_data, property_data = openmc.lib.raster_plot( + origin, width, basis, pixels, color_overlaps, level, filter, + include_properties + ) # If filter was temporarily added, remove it if len(self.tallies) > original_length: self.tallies.pop() - return ids + + return geom_data, property_data @add_plot_params def plot( @@ -1716,10 +1793,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 fee9847fadf..655c873a185 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -49,8 +49,7 @@ constexpr int32_t OVERLAP {-3}; IdData::IdData(size_t h_res, size_t v_res) : data_({v_res, h_res, 3}, NOT_FOUND) {} -void IdData::set_value(size_t y, size_t x, const Particle& p, int level, - Filter* filter, FilterMatch* match) +void IdData::set_value(size_t y, size_t x, const Particle& p, int level) { // set cell data if (p.n_coord() <= level) { @@ -67,23 +66,10 @@ void IdData::set_value(size_t y, size_t x, const Particle& 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_; } - - // set filter index - if (filter) { - filter->get_all_bins(p, TallyEstimator::COLLISION, *match); - if (match->bins_.empty()) { - data_(y, x, 3) = -1; - } else { - data_(y, x, 3) = match->bins_[0]; - } - match->bins_.clear(); - match->weights_.clear(); - } } void IdData::set_overlap(size_t y, size_t x) @@ -95,8 +81,7 @@ PropertyData::PropertyData(size_t h_res, size_t v_res) : data_({v_res, h_res, 2}, NOT_FOUND) {} -void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level, - Filter* filter, FilterMatch* match) +void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level) { Cell* c = model::cells.at(p.lowest_coord().cell()).get(); data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN; @@ -111,6 +96,158 @@ 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; +} + +//============================================================================== +// SlicePlotBase::get_raster_map implementation +//============================================================================== + +RasterData SlicePlotBase::get_raster_map(int32_t filter_index) const +{ + size_t width = pixels_[0]; + size_t height = pixels_[1]; + + bool include_filter = (filter_index >= 0); + Filter* filter = nullptr; + if (include_filter) { + filter = model::tally_filters[filter_index].get(); + } + + // get pixel size + double in_pixel = (width_[0]) / static_cast(width); + double out_pixel = (width_[1]) / static_cast(height); + + // size data array + RasterData data(width, height, include_filter); + + // 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(); + } + + // 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.; + + // arbitrary direction + Direction dir = {1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.0}; + +#pragma omp parallel + { + Particle p; + p.r() = xyz; + 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; + for (int x = 0; x < width; x++) { + p.r()[in_i] = xyz[in_i] + in_pixel * x; + p.n_coord() = 1; + // local variables + bool found_cell = exhaustive_find_cell(p); + j = p.n_coord() - 1; + if (level >= 0) { + j = level; + } + if (found_cell) { + data.set_value(y, x, p, j, filter, &match); + } + if (slice_color_overlaps_ && check_cell_overlap(p, false)) { + data.set_overlap(y, x); + } + } // inner for + } + } + + return data; +} + //============================================================================== // Global variables //============================================================================== @@ -1903,8 +2040,7 @@ void PhongRay::on_intersection() } } -extern "C" int openmc_id_map( - const void* plot, int32_t filter_index, int32_t* data_out) +extern "C" int openmc_id_map(const void* plot, int32_t* data_out) { auto plt = reinterpret_cast(plot); @@ -1913,16 +2049,11 @@ extern "C" int openmc_id_map( return OPENMC_E_INVALID_ARGUMENT; } - if (filter_index >= 0) { - if (int err = verify_filter(filter_index)) - return err; - } - if (plt->slice_color_overlaps_ && model::overlap_check_count.size() == 0) { model::overlap_check_count.resize(model::cells.size()); } - auto ids = plt->get_map(filter_index); + auto ids = plt->get_map(); // write id data to array std::copy(ids.data_.begin(), ids.data_.end(), data_out); @@ -1935,7 +2066,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; } @@ -1951,4 +2082,55 @@ extern "C" int openmc_property_map(const void* plot, double* data_out) return 0; } +extern "C" int openmc_raster_plot(const double origin[3], const double width[2], + int basis, const size_t pixels[2], bool color_overlaps, int level, + int32_t filter_index, int32_t* geom_data, double* property_data) +{ + // Validate basis + if (basis < 1 || basis > 3) { + set_errmsg("Invalid basis value. Must be 1 (xy), 2 (xz), or 3 (yz)."); + 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_raster_map logic + SlicePlotBase plot_params; + plot_params.origin_ = Position {origin[0], origin[1], origin[2]}; + plot_params.width_ = Position {width[0], width[1], 0.0}; + plot_params.basis_ = static_cast(basis); + 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_raster_map to generate data + auto data = plot_params.get_raster_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; +} + } // namespace openmc From b12b92e4700500a43939dc3e3ea7f4dff85fc5ca Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 31 Jan 2026 17:06:58 -0600 Subject: [PATCH 04/13] Avoid duplicated code in get_raster_map --- include/openmc/plot.h | 33 +++++++++------ src/plot.cpp | 99 ++++--------------------------------------- 2 files changed, 29 insertions(+), 103 deletions(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 512979c0e18..8524d2d5f42 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -147,10 +147,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 Particle& 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,10 +161,11 @@ 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 Particle& 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 @@ -172,11 +174,11 @@ struct PropertyData { struct RasterData { // Constructor - RasterData(size_t h_res, size_t v_res, bool include_filter); + 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, FilterMatch* match); + Filter* filter = nullptr, FilterMatch* match = nullptr); void set_overlap(size_t y, size_t x); // Members @@ -194,10 +196,7 @@ struct RasterData { class SlicePlotBase { public: template - T get_map() const; - - // Specialized method for RasterData that takes filter_index - RasterData get_raster_map(int32_t filter_index) const; + T get_map(int32_t filter_index = -1) const; enum class PlotBasis { xy = 1, xz = 2, yz = 3 }; @@ -218,18 +217,25 @@ 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]; + // 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(); + } + // get pixel size double in_pixel = (width_[0]) / static_cast(width); double out_pixel = (width_[1]) / static_cast(height); // size data array - T data(width, height); + T data(width, height, include_filter); // setup basis indices and initial position centered on pixel int in_i, out_i; @@ -266,6 +272,7 @@ T SlicePlotBase::get_map() const 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++) { @@ -280,7 +287,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/src/plot.cpp b/src/plot.cpp index 655c873a185..5778b7918dd 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -46,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 Particle& 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) { @@ -77,11 +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 Particle& 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; @@ -164,90 +167,6 @@ void RasterData::set_overlap(size_t y, size_t x) property_data_(y, x, 1) = OVERLAP; } -//============================================================================== -// SlicePlotBase::get_raster_map implementation -//============================================================================== - -RasterData SlicePlotBase::get_raster_map(int32_t filter_index) const -{ - size_t width = pixels_[0]; - size_t height = pixels_[1]; - - bool include_filter = (filter_index >= 0); - Filter* filter = nullptr; - if (include_filter) { - filter = model::tally_filters[filter_index].get(); - } - - // get pixel size - double in_pixel = (width_[0]) / static_cast(width); - double out_pixel = (width_[1]) / static_cast(height); - - // size data array - RasterData data(width, height, include_filter); - - // 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(); - } - - // 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.; - - // arbitrary direction - Direction dir = {1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.0}; - -#pragma omp parallel - { - Particle p; - p.r() = xyz; - 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; - for (int x = 0; x < width; x++) { - p.r()[in_i] = xyz[in_i] + in_pixel * x; - p.n_coord() = 1; - // local variables - bool found_cell = exhaustive_find_cell(p); - j = p.n_coord() - 1; - if (level >= 0) { - j = level; - } - if (found_cell) { - data.set_value(y, x, p, j, filter, &match); - } - if (slice_color_overlaps_ && check_cell_overlap(p, false)) { - data.set_overlap(y, x); - } - } // inner for - } - } - - return data; -} - //============================================================================== // Global variables //============================================================================== @@ -2114,8 +2033,8 @@ extern "C" int openmc_raster_plot(const double origin[3], const double width[2], plot_params.slice_color_overlaps_ = color_overlaps; plot_params.slice_level_ = level; - // Use get_raster_map to generate data - auto data = plot_params.get_raster_map(filter_index); + // 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); From e86c060d0b416f48318ae71236c13843a5440936 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 31 Jan 2026 17:24:55 -0600 Subject: [PATCH 05/13] Remove get_plot_bins for Filter --- openmc/lib/filter.py | 49 +--------------------------- src/tallies/filter.cpp | 72 ------------------------------------------ 2 files changed, 1 insertion(+), 120 deletions(-) diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 9e64e0731a3..dc81c3487b9 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -1,4 +1,4 @@ -from collections.abc import Mapping, Sequence +from collections.abc import Mapping from ctypes import c_int, c_int32, c_double, c_char_p, POINTER, \ create_string_buffer, c_size_t from weakref import WeakValueDictionary @@ -14,7 +14,6 @@ from .error import _error_handler from .material import Material from .mesh import _get_mesh -from .plot import _Position __all__ = [ @@ -70,11 +69,6 @@ _dll.openmc_filter_set_id.argtypes = [c_int32, c_int32] _dll.openmc_filter_set_id.restype = c_int _dll.openmc_filter_set_id.errcheck = _error_handler -_dll.openmc_filter_get_plot_bins.argtypes = [ - c_int32, _Position, _Position, c_int, POINTER(c_int), POINTER(c_int32) -] -_dll.openmc_filter_get_plot_bins.restype = c_int -_dll.openmc_filter_get_plot_bins.errcheck = _error_handler _dll.openmc_get_filter_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_filter_index.restype = c_int _dll.openmc_get_filter_index.errcheck = _error_handler @@ -161,7 +155,6 @@ _dll.openmc_zernike_filter_set_order.errcheck = _error_handler _dll.tally_filters_size.restype = c_size_t - class Filter(_FortranObjectWithID): __instances = WeakValueDictionary() @@ -210,46 +203,6 @@ def n_bins(self): _dll.openmc_filter_get_num_bins(self._index, n) return n.value - def get_plot_bins( - self, - origin: Sequence[float], - width: Sequence[float], - basis: str, - pixels: Sequence[int] - ) -> np.ndarray: - """Get filter bin indices for a rasterized plot. - - .. versionadded:: 0.15.3 - - Parameters - ---------- - origin : iterable of float - Origin of the plotting view. Should have length 3. - width : iterable of float - Width of the plotting view. Should have length 2. - basis : {'xy', 'xz', 'yz'} - Plotting basis. - pixels : iterable of int - Number of pixels in each direction. Should have length 2. - - Returns - ------- - 2D numpy array with mesh bin indices corresponding to each pixel within - the plotting view. - - """ - origin = _Position(*origin) - width = _Position(*width) - basis = {'xy': 1, 'xz': 2, 'yz': 3}[basis] - pixel_array = (c_int*2)(*pixels) - img_data = np.zeros((pixels[1], pixels[0]), dtype=np.dtype('int32')) - - _dll.openmc_filter_get_plot_bins( - self._index, origin, width, basis, pixel_array, - img_data.ctypes.data_as(POINTER(c_int32)) - ) - return img_data - class EnergyFilter(Filter): filter_type = 'energy' diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index acfd0fb556e..16e4c055021 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -10,8 +10,6 @@ #include "openmc/capi.h" #include "openmc/constants.h" // for MAX_LINE_LEN; #include "openmc/error.h" -#include "openmc/geometry.h" -#include "openmc/position.h" #include "openmc/tallies/filter_azimuthal.h" #include "openmc/tallies/filter_cell.h" #include "openmc/tallies/filter_cell_instance.h" @@ -253,76 +251,6 @@ extern "C" int openmc_filter_get_num_bins(int32_t index, int* n_bins) return 0; } -extern "C" int openmc_filter_get_plot_bins(int32_t index, Position origin, - Position width, int basis, int* pixels, int32_t* data) -{ - if (int err = verify_filter(index)) - return err; - const auto& filter = model::tally_filters[index].get(); - - int pixel_width = pixels[0]; - int pixel_height = pixels[1]; - - // get pixel size - double in_pixel = (width[0]) / static_cast(pixel_width); - double out_pixel = (width[1]) / static_cast(pixel_height); - - // setup basis indices and initial position centered on pixel - int in_i, out_i; - Position xyz = origin; - enum class PlotBasis { xy = 1, xz = 2, yz = 3 }; - PlotBasis basis_enum = static_cast(basis); - switch (basis_enum) { - 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(); - } - - // 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.; - -#pragma omp parallel - { - Particle p; - p.r() = xyz; - p.u() = {1.0, 0.0, 0.0}; - p.coord(0).universe() = model::root_universe; - FilterMatch match; - -#pragma omp for - for (int y = 0; y < pixel_height; y++) { - p.r()[out_i] = xyz[out_i] - out_pixel * y; - for (int x = 0; x < pixel_width; x++) { - p.r()[in_i] = xyz[in_i] + in_pixel * x; - p.n_coord() = 1; - bool found_cell = exhaustive_find_cell(p); - filter->get_all_bins(p, TallyEstimator::COLLISION, match); - if (match.bins_.empty()) { - data[pixel_width * y + x] = -1; - } else { - data[pixel_width * y + x] = match.bins_[0]; - } - match.bins_.clear(); - match.weights_.clear(); - } - } - } - - return 0; -} - extern "C" int openmc_get_filter_index(int32_t id, int32_t* index) { auto it = model::filter_map.find(id); From 2e81cb07950a3c964be0ec4cedc5211a60a12540 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 31 Jan 2026 17:30:43 -0600 Subject: [PATCH 06/13] Add tests for raster_plot --- tests/unit_tests/test_raster_plot.py | 153 +++++++++++++++++++++++++++ 1 file changed, 153 insertions(+) create mode 100644 tests/unit_tests/test_raster_plot.py diff --git a/tests/unit_tests/test_raster_plot.py b/tests/unit_tests/test_raster_plot.py new file mode 100644 index 00000000000..46e162fc584 --- /dev/null +++ b/tests/unit_tests/test_raster_plot.py @@ -0,0 +1,153 @@ +import numpy as np +import openmc +from openmc.examples import pwr_pin_cell + + +def test_raster_plot_basic(run_in_tmpdir): + """Test basic raster_plot functionality.""" + model = pwr_pin_cell() + geom_data, prop_data = model.raster_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_raster_plot_no_properties(run_in_tmpdir): + """Test raster_plot without property data.""" + model = pwr_pin_cell() + geom_data, prop_data = model.raster_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_raster_plot_with_filter(run_in_tmpdir): + """Test raster_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.raster_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_raster_plot_overlaps(run_in_tmpdir): + """Test raster_plot with overlap detection.""" + model = pwr_pin_cell() + geom_data, _ = model.raster_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_raster_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.raster_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_raster_plot_different_bases(run_in_tmpdir): + """Test raster_plot with different basis planes.""" + model = pwr_pin_cell() + + for basis in ['xy', 'xz', 'yz']: + geom_data, prop_data = model.raster_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_raster_plot_level(run_in_tmpdir): + """Test raster_plot with specific universe level.""" + model = pwr_pin_cell() + geom_data, _ = model.raster_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 From d80198bb6d3abb9772bc537ba6907ec89c19732d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 2 Feb 2026 11:38:31 -0600 Subject: [PATCH 07/13] Small doc fix --- openmc/lib/plot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py index c142ed164fc..f324746410a 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -293,7 +293,7 @@ def raster_plot(origin, width, basis, pixels, color_overlaps=False, level=-1, Whether to detect overlapping cells level : int, optional Universe level (-1 for deepest) - filter : openmc.Filter, optional + filter : openmc.lib.Filter, optional Filter for bin index lookup include_properties : bool, optional Whether to compute temperature/density From 3db1e542e8d4f46b0d92094efcb63b97385981ca Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 3 Feb 2026 09:52:58 -0600 Subject: [PATCH 08/13] Initial shot at arbitrary orientation slice plots --- include/openmc/capi.h | 6 +- include/openmc/plot.h | 51 ++++++----------- openmc/lib/plot.py | 128 +++++++++++++++++++++++++++++++++++++----- src/plot.cpp | 48 +++++++++++++--- 4 files changed, 176 insertions(+), 57 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 820612186a9..e17acb642ce 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -125,9 +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_raster_plot(const double origin[3], const double width[2], int basis, - const size_t pixels[2], bool color_overlaps, int level, int32_t filter_index, - int32_t* geom_data, double* property_data); +int openmc_raster_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_rectilinear_mesh_get_grid(int32_t index, double** grid_x, int* nx, double** grid_y, int* ny, double** grid_z, int* nz); int openmc_rectilinear_mesh_set_grid(int32_t index, const double* grid_x, diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 8524d2d5f42..e042392b5ac 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -208,8 +208,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 @@ -230,44 +232,27 @@ T SlicePlotBase::get_map(int32_t filter_index) const filter = model::tally_filters[filter_index].get(); } - // get pixel size - double in_pixel = (width_[0]) / static_cast(width); - double out_pixel = (width_[1]) / static_cast(height); - // size data array T data(width, height, include_filter); - // 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(); - } + // compute pixel steps and top-left pixel center + Position u_step = u_span_ / static_cast(width); + Position v_step = v_span_ / static_cast(height); - // 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.; + Position start = + origin_ - 0.5 * u_span_ + 0.5 * v_span_ + 0.5 * u_step - 0.5 * v_step; - // arbitrary direction - Direction dir = {1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.0}; + Direction dir = u_span_.cross(v_span_); + double dir_norm = dir.norm(); + if (dir_norm == 0.0) { + fatal_error("Slice span vectors are invalid (zero area)."); + } + dir /= dir_norm; #pragma omp parallel { Particle p; - p.r() = xyz; + p.r() = start; p.u() = dir; p.coord(0).universe() = model::root_universe; int level = slice_level_; @@ -276,9 +261,9 @@ T SlicePlotBase::get_map(int32_t filter_index) const #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); diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py index f324746410a..c37280931ae 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -54,6 +54,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 @@ -84,6 +88,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), @@ -94,6 +100,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): @@ -112,6 +142,7 @@ def width(self): @width.setter def width(self, width): self.width_.x = width + self._update_spans() @property def height(self): @@ -120,6 +151,7 @@ def height(self): @height.setter def height(self, height): self.width_.y = height + self._update_spans() @property def basis(self): @@ -146,6 +178,7 @@ def basis(self, basis): self.basis_ = 2 elif basis == 'yz': self.basis_ = 3 + self._update_spans() return if isinstance(basis, int): @@ -153,6 +186,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") @@ -261,8 +295,8 @@ def property_map(plot): _dll.openmc_raster_plot.argtypes = [ POINTER(c_double * 3), # origin - POINTER(c_double * 2), # width - c_int, # basis + 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 @@ -274,8 +308,10 @@ def property_map(plot): _dll.openmc_raster_plot.errcheck = _error_handler -def raster_plot(origin, width, basis, pixels, color_overlaps=False, level=-1, - filter=None, include_properties=True): +def raster_plot(origin, width=None, basis='xy', pixels=None, + color_overlaps=False, level=-1, filter=None, + include_properties=True, u_span=None, v_span=None, + axes=None): """ Generate a 2D raster of geometry and property data for plotting. @@ -284,9 +320,10 @@ def raster_plot(origin, width, basis, pixels, color_overlaps=False, level=-1, origin : sequence of float Center position of the plot [x, y, z] width : sequence of float - Width of the plot [horizontal, vertical] + Width of the plot [horizontal, vertical]. Mutually exclusive with + u_span/v_span/axes. basis : {'xy', 'xz', 'yz'} or int - Plot basis + Plot basis. Ignored if u_span/v_span/axes are provided. pixels : sequence of int Number of pixels [horizontal, vertical] color_overlaps : bool, optional @@ -297,6 +334,14 @@ def raster_plot(origin, width, basis, pixels, color_overlaps=False, level=-1, Filter for bin index lookup include_properties : bool, optional Whether to compute temperature/density + 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. + axes : tuple of (u_span, v_span), optional + Convenience alternative to pass both span vectors. Returns ------- @@ -308,14 +353,71 @@ def raster_plot(origin, width, basis, pixels, color_overlaps=False, level=-1, Array of shape (v_res, h_res, 2) with float64 dtype containing [temperature, density], or None if include_properties=False """ - # Convert basis string to int - basis_map = {'xy': 1, 'xz': 2, 'yz': 3} - if isinstance(basis, str): - basis = basis_map[basis.lower()] + if pixels is None: + raise ValueError("pixels must be specified.") + if len(pixels) != 2: + raise ValueError("pixels must be a length-2 sequence.") + + if axes is not None: + if u_span is not None or v_span is not None: + raise ValueError("axes is mutually exclusive with u_span/v_span.") + if len(axes) != 2: + raise ValueError("axes must be a length-2 iterable of span vectors.") + u_span, v_span = axes + + 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/axes.") + + 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) - width_arr = (c_double * 2)(*width) + 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 @@ -338,8 +440,8 @@ def raster_plot(origin, width, basis, pixels, color_overlaps=False, level=-1, _dll.openmc_raster_plot( origin_arr, - width_arr, - basis, + u_span_arr, + v_span_arr, pixels_arr, color_overlaps, level, diff --git a/src/plot.cpp b/src/plot.cpp index 5778b7918dd..5be6de33816 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -495,6 +495,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())); @@ -1057,6 +1073,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_; @@ -2001,13 +2019,25 @@ extern "C" int openmc_property_map(const void* plot, double* data_out) return 0; } -extern "C" int openmc_raster_plot(const double origin[3], const double width[2], - int basis, const size_t pixels[2], bool color_overlaps, int level, - int32_t filter_index, int32_t* geom_data, double* property_data) +extern "C" int openmc_raster_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 basis - if (basis < 1 || basis > 3) { - set_errmsg("Invalid basis value. Must be 1 (xy), 2 (xz), or 3 (yz)."); + // 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; } @@ -2026,8 +2056,10 @@ extern "C" int openmc_raster_plot(const double origin[3], const double width[2], // Create a temporary SlicePlotBase object to reuse get_raster_map logic SlicePlotBase plot_params; plot_params.origin_ = Position {origin[0], origin[1], origin[2]}; - plot_params.width_ = Position {width[0], width[1], 0.0}; - plot_params.basis_ = static_cast(basis); + 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; From adf7ad3ef086b06a85e86107e59f67bf42f2f80a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 3 Feb 2026 23:59:16 -0600 Subject: [PATCH 09/13] Remove 'axes' and reorder arguments in lib.raster_plot --- openmc/lib/plot.py | 35 ++++++++++++----------------------- 1 file changed, 12 insertions(+), 23 deletions(-) diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py index c37280931ae..094a6b10025 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -308,12 +308,10 @@ def property_map(plot): _dll.openmc_raster_plot.errcheck = _error_handler -def raster_plot(origin, width=None, basis='xy', pixels=None, - color_overlaps=False, level=-1, filter=None, - include_properties=True, u_span=None, v_span=None, - axes=None): - """ - Generate a 2D raster of geometry and property data for plotting. +def raster_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 ---------- @@ -321,9 +319,15 @@ def raster_plot(origin, width=None, basis='xy', pixels=None, 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/axes. + u_span/v_span. basis : {'xy', 'xz', 'yz'} or int - Plot basis. Ignored if u_span/v_span/axes are provided. + 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 @@ -334,14 +338,6 @@ def raster_plot(origin, width=None, basis='xy', pixels=None, Filter for bin index lookup include_properties : bool, optional Whether to compute temperature/density - 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. - axes : tuple of (u_span, v_span), optional - Convenience alternative to pass both span vectors. Returns ------- @@ -358,13 +354,6 @@ def raster_plot(origin, width=None, basis='xy', pixels=None, if len(pixels) != 2: raise ValueError("pixels must be a length-2 sequence.") - if axes is not None: - if u_span is not None or v_span is not None: - raise ValueError("axes is mutually exclusive with u_span/v_span.") - if len(axes) != 2: - raise ValueError("axes must be a length-2 iterable of span vectors.") - u_span, v_span = axes - 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/axes.") From fceb98ccdefb79a1e573f68b63c2e321d9628f39 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 4 Feb 2026 00:13:28 -0600 Subject: [PATCH 10/13] Support u_span and v_span in Model.raster_plot --- openmc/lib/plot.py | 3 +-- openmc/model/model.py | 39 +++++++++++++++++++++++++--- tests/unit_tests/test_raster_plot.py | 15 +++++++++++ 3 files changed, 51 insertions(+), 6 deletions(-) diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py index 094a6b10025..c5b049a033a 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -338,7 +338,6 @@ def raster_plot(origin, width=None, basis='xy', u_span=None, v_span=None, Filter for bin index lookup include_properties : bool, optional Whether to compute temperature/density - Returns ------- geom_data : numpy.ndarray @@ -355,7 +354,7 @@ def raster_plot(origin, width=None, basis='xy', u_span=None, v_span=None, 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/axes.") + 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: diff --git a/openmc/model/model.py b/openmc/model/model.py index 4fdfca11f6d..a6e011d05c3 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1112,6 +1112,8 @@ def raster_plot( 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, @@ -1144,6 +1146,12 @@ def raster_plot( 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 @@ -1170,8 +1178,23 @@ def raster_plot( """ 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.") + + 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 @@ -1190,8 +1213,16 @@ def raster_plot( with openmc.lib.TemporarySession(self, **init_kwargs): geom_data, property_data = openmc.lib.raster_plot( - origin, width, basis, pixels, color_overlaps, level, filter, - include_properties + 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 diff --git a/tests/unit_tests/test_raster_plot.py b/tests/unit_tests/test_raster_plot.py index 46e162fc584..f09eba3c5db 100644 --- a/tests/unit_tests/test_raster_plot.py +++ b/tests/unit_tests/test_raster_plot.py @@ -121,6 +121,21 @@ def test_raster_plot_different_bases(run_in_tmpdir): assert prop_data.shape == (25, 25, 2) +def test_raster_plot_oriented_spans(run_in_tmpdir): + """Test raster_plot with oriented span vectors.""" + model = pwr_pin_cell() + + geom_data, prop_data = model.raster_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_raster_plot_level(run_in_tmpdir): """Test raster_plot with specific universe level.""" model = pwr_pin_cell() From ac0e900b7a60114d12d1bbff244f7ca8b6cb82ca Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 8 Feb 2026 14:28:22 -0600 Subject: [PATCH 11/13] Rename raster -> slice for consistency --- include/openmc/capi.h | 2 +- openmc/lib/plot.py | 10 ++-- openmc/model/model.py | 4 +- src/plot.cpp | 9 ++-- ...test_raster_plot.py => test_slice_plot.py} | 46 +++++++++---------- 5 files changed, 35 insertions(+), 36 deletions(-) rename tests/unit_tests/{test_raster_plot.py => test_slice_plot.py} (77%) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index e17acb642ce..6a361914da9 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -125,7 +125,7 @@ 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_raster_plot(const double origin[3], const double u_span[3], +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_rectilinear_mesh_get_grid(int32_t index, double** grid_x, int* nx, diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py index c5b049a033a..2784974efdf 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -293,7 +293,7 @@ def property_map(plot): return prop_data -_dll.openmc_raster_plot.argtypes = [ +_dll.openmc_slice_plot.argtypes = [ POINTER(c_double * 3), # origin POINTER(c_double * 3), # u_span POINTER(c_double * 3), # v_span @@ -304,11 +304,11 @@ def property_map(plot): POINTER(c_int32), # geom_data POINTER(c_double), # property_data (can be None) ] -_dll.openmc_raster_plot.restype = c_int -_dll.openmc_raster_plot.errcheck = _error_handler +_dll.openmc_slice_plot.restype = c_int +_dll.openmc_slice_plot.errcheck = _error_handler -def raster_plot(origin, width=None, basis='xy', u_span=None, v_span=None, +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. @@ -426,7 +426,7 @@ def raster_plot(origin, width=None, basis='xy', u_span=None, v_span=None, property_data = None prop_ptr = None - _dll.openmc_raster_plot( + _dll.openmc_slice_plot( origin_arr, u_span_arr, v_span_arr, diff --git a/openmc/model/model.py b/openmc/model/model.py index a6e011d05c3..749772e2be3 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1106,7 +1106,7 @@ def id_map( return ids - def raster_plot( + def slice_plot( self, origin: Sequence[float] | None = None, width: Sequence[float] | None = None, @@ -1212,7 +1212,7 @@ def raster_plot( self.tallies.append(temp_tally) with openmc.lib.TemporarySession(self, **init_kwargs): - geom_data, property_data = openmc.lib.raster_plot( + geom_data, property_data = openmc.lib.slice_plot( origin=origin, width=width, basis=basis, diff --git a/src/plot.cpp b/src/plot.cpp index 5be6de33816..befd6035c18 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -2019,10 +2019,9 @@ extern "C" int openmc_property_map(const void* plot, double* data_out) return 0; } -extern "C" int openmc_raster_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) +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]}; @@ -2053,7 +2052,7 @@ extern "C" int openmc_raster_plot(const double origin[3], } try { - // Create a temporary SlicePlotBase object to reuse get_raster_map logic + // 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; diff --git a/tests/unit_tests/test_raster_plot.py b/tests/unit_tests/test_slice_plot.py similarity index 77% rename from tests/unit_tests/test_raster_plot.py rename to tests/unit_tests/test_slice_plot.py index f09eba3c5db..8b1919c4158 100644 --- a/tests/unit_tests/test_raster_plot.py +++ b/tests/unit_tests/test_slice_plot.py @@ -3,10 +3,10 @@ from openmc.examples import pwr_pin_cell -def test_raster_plot_basic(run_in_tmpdir): - """Test basic raster_plot functionality.""" +def test_slice_plot_basic(run_in_tmpdir): + """Test basic slice_plot functionality.""" model = pwr_pin_cell() - geom_data, prop_data = model.raster_plot( + geom_data, prop_data = model.slice_plot( origin=(0, 0, 0), width=(1.0, 1.0), pixels=(100, 100), @@ -24,10 +24,10 @@ def test_raster_plot_basic(run_in_tmpdir): assert np.any(prop_data[:, :, 0] > 0) # Valid temperatures -def test_raster_plot_no_properties(run_in_tmpdir): - """Test raster_plot without property data.""" +def test_slice_plot_no_properties(run_in_tmpdir): + """Test slice_plot without property data.""" model = pwr_pin_cell() - geom_data, prop_data = model.raster_plot( + geom_data, prop_data = model.slice_plot( origin=(0, 0, 0), width=(1.0, 1.0), pixels=(50, 50), @@ -39,13 +39,13 @@ def test_raster_plot_no_properties(run_in_tmpdir): assert prop_data is None -def test_raster_plot_with_filter(run_in_tmpdir): - """Test raster_plot with a cell filter.""" +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.raster_plot( + geom_data, _ = model.slice_plot( origin=(0, 0, 0), width=(1.0, 1.0), pixels=(50, 50), @@ -62,10 +62,10 @@ def test_raster_plot_with_filter(run_in_tmpdir): assert np.any(filter_bins[valid_cells] >= 0) -def test_raster_plot_overlaps(run_in_tmpdir): - """Test raster_plot with overlap detection.""" +def test_slice_plot_overlaps(run_in_tmpdir): + """Test slice_plot with overlap detection.""" model = pwr_pin_cell() - geom_data, _ = model.raster_plot( + geom_data, _ = model.slice_plot( origin=(0, 0, 0), width=(1.0, 1.0), pixels=(50, 50), @@ -79,13 +79,13 @@ def test_raster_plot_overlaps(run_in_tmpdir): # Note: This test may pass without finding overlaps if geometry is correct -def test_raster_plot_overlaps_with_filter(run_in_tmpdir): +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.raster_plot( + geom_data, _ = model.slice_plot( origin=(0, 0, 0), width=(1.0, 1.0), pixels=(50, 50), @@ -105,12 +105,12 @@ def test_raster_plot_overlaps_with_filter(run_in_tmpdir): "Filter bins should be preserved even where overlaps are detected" -def test_raster_plot_different_bases(run_in_tmpdir): - """Test raster_plot with different basis planes.""" +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.raster_plot( + geom_data, prop_data = model.slice_plot( origin=(0, 0, 0), width=(1.0, 1.0), pixels=(25, 25), @@ -121,11 +121,11 @@ def test_raster_plot_different_bases(run_in_tmpdir): assert prop_data.shape == (25, 25, 2) -def test_raster_plot_oriented_spans(run_in_tmpdir): - """Test raster_plot with oriented span vectors.""" +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.raster_plot( + 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), @@ -136,10 +136,10 @@ def test_raster_plot_oriented_spans(run_in_tmpdir): assert prop_data.shape == (25, 25, 2) -def test_raster_plot_level(run_in_tmpdir): - """Test raster_plot with specific universe level.""" +def test_slice_plot_level(run_in_tmpdir): + """Test slice_plot with specific universe level.""" model = pwr_pin_cell() - geom_data, _ = model.raster_plot( + geom_data, _ = model.slice_plot( origin=(0, 0, 0), width=(1.0, 1.0), pixels=(50, 50), From 6533034cdf11ec5ffde363ddca695425cf050de0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 16 Feb 2026 10:11:07 -0600 Subject: [PATCH 12/13] Restore original (arbitrary) direction for slice plotting --- include/openmc/plot.h | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 05921080696..b011d217525 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -244,12 +244,16 @@ T SlicePlotBase::get_map(int32_t filter_index) const Position start = origin_ - 0.5 * u_span_ + 0.5 * v_span_ + 0.5 * u_step - 0.5 * v_step; - Direction dir = u_span_.cross(v_span_); - double dir_norm = dir.norm(); - if (dir_norm == 0.0) { + // 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)."); } - dir /= dir_norm; + + // 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 { From 4824efac824d03579603ef2bd1cad06611b39815 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 16 Feb 2026 10:14:58 -0600 Subject: [PATCH 13/13] Revert changes to comments --- src/plot.cpp | 78 +++++++++++++++++++++++++--------------------------- 1 file changed, 37 insertions(+), 41 deletions(-) diff --git a/src/plot.cpp b/src/plot.cpp index 6cbbe6a2dcb..607beb84415 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -259,10 +259,9 @@ void Plot::print_info() const void read_plots_xml() { - // Check if plots.xml exists; this is only necessary when the plot runmode - // is initiated. Otherwise, we want to read plots.xml because it may be - // called later via the API. In that case, its ok for a plots.xml to not - // exist + // Check if plots.xml exists; this is only necessary when the plot runmode is + // initiated. Otherwise, we want to read plots.xml because it may be called + // later via the API. In that case, its ok for a plots.xml to not exist std::string filename = settings::path_input + "plots.xml"; if (!file_exists(filename) && settings::run_mode == RunMode::PLOTTING) { fatal_error(fmt::format("Plots XML file '{}' does not exist!", filename)); @@ -1042,15 +1041,14 @@ void Plot::draw_mesh_lines(ImageData& data) const } /* outputs a binary file that can be input into silomesh for 3D geometry - * visualization. It works the same way as create_image by dragging a - * particle across the geometry for the specified number of voxels. The first - * 3 int's in the binary are the number of x, y, and z voxels. The next 3 - * double's are the widths of the voxels in the x, y, and z directions. The - * next 3 double's are the x, y, and z coordinates of the lower left point. - * Finally the binary is filled with entries of four int's each. Each 'row' in - * the binary contains four int's: 3 for x,y,z position and 1 for cell or - * material id. For 1 million voxels this produces a file of approximately - * 15MB. + * visualization. It works the same way as create_image by dragging a particle + * across the geometry for the specified number of voxels. The first 3 int's in + * the binary are the number of x, y, and z voxels. The next 3 double's are + * the widths of the voxels in the x, y, and z directions. The next 3 double's + * are the x, y, and z coordinates of the lower left point. Finally the binary + * is filled with entries of four int's each. Each 'row' in the binary contains + * four int's: 3 for x,y,z position and 1 for cell or material id. For 1 + * million voxels this produces a file of approximately 15MB. */ void Plot::create_voxel() const { @@ -1298,9 +1296,8 @@ bool WireframeRayTracePlot::trackstack_equivalent( return false; // Check if neighboring cells are different - // if (track1[t1_i ? t1_i - 1 : 0].id != track2[t2_i ? t2_i - 1 : - // 0].id) return false; if (track1[t1_i < track1.size() - 1 ? t1_i + 1 - // : t1_i + // if (track1[t1_i ? t1_i - 1 : 0].id != track2[t2_i ? t2_i - 1 : 0].id) + // return false; if (track1[t1_i < track1.size() - 1 ? t1_i + 1 : t1_i // ].id != // track2[t2_i < track2.size() - 1 ? t2_i + 1 : t2_i].id) return // false; @@ -1360,8 +1357,8 @@ ImageData WireframeRayTracePlot::create_image() const size_t height = pixels()[1]; ImageData data({width, height}, not_found_); - // This array marks where the initial wireframe was drawn. We convolve it - // with a filter that gets adjusted with the wireframe thickness in order to + // This array marks where the initial wireframe was drawn. We convolve it with + // a filter that gets adjusted with the wireframe thickness in order to // thicken the lines. xt::xtensor wireframe_initial({width, height}, 0); @@ -1395,9 +1392,9 @@ ImageData WireframeRayTracePlot::create_image() const for (int iter = 0; iter <= pixels()[1] / n_threads; iter++) { // Save bottom line of current work chunk to compare against later. This - // used to be inside the below if block, but it causes a spurious line - // to be drawn at the bottom of the image. Not sure why, but moving it - // here fixes things. + // used to be inside the below if block, but it causes a spurious line to + // be drawn at the bottom of the image. Not sure why, but moving it here + // fixes things. if (tid == n_threads - 1) old_segments = this_line_segments[n_threads - 1]; @@ -1422,8 +1419,8 @@ ImageData WireframeRayTracePlot::create_image() const // There must be at least two cell intersections to color, front and // back of the cell. Maybe an infinitely thick cell could be present - // with no back, but why would you want to color that? It's easier - // to just skip that edge case and not even color it. + // with no back, but why would you want to color that? It's easier to + // just skip that edge case and not even color it. if (segments.size() <= 1) continue; @@ -1458,8 +1455,8 @@ ImageData WireframeRayTracePlot::create_image() const } } // end "if" vert in correct range - // We require a barrier before comparing vertical neighbors' - // intersection stacks. i.e. all threads must be done with their line. + // We require a barrier before comparing vertical neighbors' intersection + // stacks. i.e. all threads must be done with their line. #pragma omp barrier // Now that the horizontal line has finished rendering, we can fill in @@ -1483,8 +1480,8 @@ ImageData WireframeRayTracePlot::create_image() const } } - // We need another barrier to ensure threads don't proceed to modify - // their intersection stacks on that horizontal line while others are + // We need another barrier to ensure threads don't proceed to modify their + // intersection stacks on that horizontal line while others are // potentially still working on the above. #pragma omp barrier vert += n_threads; @@ -1750,8 +1747,7 @@ void Ray::trace() // radiation transport model. // // After phase one is done, we can starting tracing from cell to cell within - // the model. This step can use neighbor lists to accelerate the ray - // tracing. + // the model. This step can use neighbor lists to accelerate the ray tracing. // Attempt to initialize the particle. We may have to enter a loop to move // it up to the edge of the model. @@ -1849,13 +1845,13 @@ void Ray::trace() inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10); // Call the specialized logic for this type of ray. Note that we do not - // call this if the advance distance is very small. Unfortunately, it - // seems darn near impossible to get the particle advanced to the model - // boundary and through it without sometimes accidentally calling - // on_intersection twice. This incorrectly shades the region as occluded - // when it might not actually be. By screening out intersection distances - // smaller than a threshold 10x larger than the scoot distance used to - // advance up to the model boundary, we can avoid that situation. + // call this if the advance distance is very small. Unfortunately, it seems + // darn near impossible to get the particle advanced to the model boundary + // and through it without sometimes accidentally calling on_intersection + // twice. This incorrectly shades the region as occluded when it might not + // actually be. By screening out intersection distances smaller than a + // threshold 10x larger than the scoot distance used to advance up to the + // model boundary, we can avoid that situation. if (call_on_intersection) { on_intersection(); if (stop_) @@ -1918,8 +1914,8 @@ void PhongRay::on_intersection() // TODO // Not sure what can cause a surface token to be invalid here, although it // sometimes happens for a few pixels. It's very very rare, so proceed by - // coloring the pixel with the overlap color. It seems to happen only for - // a few pixels on the outer boundary of a hex lattice. + // coloring the pixel with the overlap color. It seems to happen only for a + // few pixels on the outer boundary of a hex lattice. // // We cannot detect it in the outer loop, and it only matters here, so // that's why the error handling is a little different than for a lost @@ -1990,9 +1986,9 @@ void PhongRay::on_intersection() compute_distance(); } else { - // If it's not facing the light, we color with the diffuse contribution, - // so next we check if we're going to occlude the last reflected surface. - // if so, color by the diffuse contribution instead + // If it's not facing the light, we color with the diffuse contribution, so + // next we check if we're going to occlude the last reflected surface. if + // so, color by the diffuse contribution instead if (orig_hit_id_ == -1) fatal_error("somehow a ray got reflected but not original ID set?");