Merge develop into ww_edge_cases to resolve conflict

.github/workflows/ci.yml

@@ -177,7 +177,7 @@ jobs:
 
       - name: Setup tmate debug session
         continue-on-error: true
-        if: ${{ contains(env.COMMIT_MESSAGE, '[gha-debug]') }}
+        if: ${{ failure() && contains(env.COMMIT_MESSAGE, '[gha-debug]') }}
         uses: mxschmitt/action-tmate@v3
         timeout-minutes: 10
 

docs/source/io_formats/collision_track.rst

@@ -10,7 +10,7 @@ may also be written after each batch when multiple files are requested
 (``collision_track.N.h5``) or when the run is performed in parallel. The file
 contains the information needed to reconstruct each recorded collision.
 
-The current revision of the collision track file format is 1.1.
+The current revision of the collision track file format is 1.2.
 
 **/**
 
@@ -33,7 +33,7 @@ The current revision of the collision track file format is 1.1.
               - ``event_mt`` (*int*) -- ENDF MT number identifying the reaction.
               - ``delayed_group`` (*int*) -- Delayed neutron group index (non-zero for delayed events).
               - ``cell_id`` (*int*) -- ID of the cell in which the collision occurred.
-              - ``nuclide_id`` (*int*) -- ZA identifier of the nuclide (ZZZAAAM format).
+              - ``nuclide_id`` (*int*) -- PDG number of the nuclide (100ZZZAAAM).
               - ``material_id`` (*int*) -- ID of the material containing the collision site.
               - ``universe_id`` (*int*) -- ID of the universe containing the collision site.
               - ``n_collision`` (*int*) -- Collision counter for the particle history.

docs/source/io_formats/settings.rst

@@ -98,6 +98,11 @@ sub-elements:
     A list of strings representing the nuclide, to define specific
     define specific target nuclide collisions to be banked.
 
+    .. note::
+      Electron and positron collision-track events are not associated with
+      a specific nuclide. If a ``nuclides`` entry is specified, these events
+      are omitted.
+
     *Default*: None
 
   :reactions:
@@ -606,30 +611,30 @@ found in the :ref:`random ray user guide <random_ray>`.
     *Default*: None
 
   :adjoint_source:
-    Specifies an adjoint fixed source for adjoint transport simulations, and 
-    follows the format for :ref:`source_element`. The distributions which make 
-    up the adjoint source are subject to the same restrictions as forward 
+    Specifies an adjoint fixed source for adjoint transport simulations, and
+    follows the format for :ref:`source_element`. The distributions which make
+    up the adjoint source are subject to the same restrictions as forward
     fixed sources in Random Ray mode.
 
     *Default*: None
-  
+
   :adjoint:
-    Specifies whether to perform adjoint transport. The default is 'False', 
+    Specifies whether to perform adjoint transport. The default is 'False',
     corresponding to forward transport.
 
     *Default*: None
-  
+
   :volume_estimator:
-    Specifies choice of volume estimator for the random ray solver. Options 
+    Specifies choice of volume estimator for the random ray solver. Options
     are 'naive', 'simulation_averaged', or 'hybrid'. The default is 'hybrid'.
 
     *Default*: None
 
   :volume_normalized_flux_tallies:
-    Specifies whether to normalize flux tallies by volume (bool). The 
-    default is 'False'. When enabled, flux tallies will be reported in units 
-    of cm/cm^3. When disabled, flux tallies will be reported in units of cm 
-    (i.e., total distance traveled by neutrons in the spatial tally 
+    Specifies whether to normalize flux tallies by volume (bool). The
+    default is 'False'. When enabled, flux tallies will be reported in units
+    of cm/cm^3. When disabled, flux tallies will be reported in units of cm
+    (i.e., total distance traveled by neutrons in the spatial tally
     region).
 
     *Default*: None
@@ -1757,11 +1762,11 @@ mesh-based weight windows.
         The ratio of the lower to upper weight window bounds.
 
         *Default*: 5.0
-    
+
     For FW-CADIS:
 
       :targets:
-        A sequence of IDs corresponding to the tallies which cover phase 
+        A sequence of IDs corresponding to the tallies which cover phase
         space regions of interest for local variance reduction.
 
         *Default*: None

docs/source/usersguide/settings.rst

@@ -792,6 +792,11 @@ collision_track.h5 file at the end of the simulation. The file contains
 300 recorded collisions that occurred in materials with IDs 1 or 2, involving
 fission or (n,2n) reactions on the nuclides U-238 or O-16, within cells
 with IDs 5 and 12.
+
+.. note::
+   Electron and positron collision-track events are not associated with a
+   specific nuclide. If a ``nuclides`` entry is specified, these events are omitted.
+
 The file can be read using :func:`openmc.read_collision_track_file`.
 The example below shows how to extract the data from the collision_track
 feature and displays the fields stored in the file:

include/openmc/capi.h

@@ -123,8 +123,13 @@ int openmc_new_filter(const char* type, int32_t* index);
 int openmc_next_batch(int* status);
 int openmc_nuclide_name(int index, const char** name);
 int openmc_plot_geometry();
+// Deprecated; use openmc_slice_data.
 int openmc_id_map(const void* slice, int32_t* data_out);
+// Deprecated; use openmc_slice_data.
 int openmc_property_map(const void* slice, double* data_out);
+int openmc_slice_data(const double origin[3], const double u_span[3],
+  const double v_span[3], const size_t pixels[2], bool show_overlaps, int level,
+  int32_t filter_index, int32_t* geom_data, double* property_data);
 int openmc_get_plot_index(int32_t id, int32_t* index);
 int openmc_plot_get_id(int32_t index, int32_t* id);
 int openmc_plot_set_id(int32_t index, int32_t id);

include/openmc/constants.h

@@ -35,7 +35,7 @@ constexpr array<int, 2> VERSION_VOXEL {2, 0};
 constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
 constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
 constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
-constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 1};
+constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 2};
 
 // ============================================================================
 // ADJUSTABLE PARAMETERS

include/openmc/nuclide.h

@@ -84,6 +84,9 @@ class Nuclide {
   double collapse_rate(int MT, double temperature, span<const double> energy,
     span<const double> flux) const;
 
+  //! Return a ParticleType object representing this nuclide
+  ParticleType particle_type() const { return {Z_, A_, metastable_}; }
+
   //============================================================================
   // Data members
   std::string name_; //!< Name of nuclide, e.g. "U235"

include/openmc/plot.h

@@ -18,6 +18,8 @@
 #include "openmc/position.h"
 #include "openmc/random_lcg.h"
 #include "openmc/ray.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/tallies/filter_match.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -148,10 +150,11 @@ class PlottableInterface {
 
 struct IdData {
   // Constructor
-  IdData(size_t h_res, size_t v_res);
+  IdData(size_t h_res, size_t v_res, bool include_filter = false);
 
   // Methods
-  void set_value(size_t y, size_t x, const GeometryState& p, int level);
+  void set_value(size_t y, size_t x, const Particle& p, int level,
+    Filter* filter = nullptr, FilterMatch* match = nullptr);
   void set_overlap(size_t y, size_t x);
 
   // Members
@@ -160,24 +163,42 @@ struct IdData {
 
 struct PropertyData {
   // Constructor
-  PropertyData(size_t h_res, size_t v_res);
+  PropertyData(size_t h_res, size_t v_res, bool include_filter = false);
 
   // Methods
-  void set_value(size_t y, size_t x, const GeometryState& p, int level);
+  void set_value(size_t y, size_t x, const Particle& p, int level,
+    Filter* filter = nullptr, FilterMatch* match = nullptr);
   void set_overlap(size_t y, size_t x);
 
   // Members
   tensor::Tensor<double> data_; //!< 2D array of temperature & density data
 };
 
+struct RasterData {
+  // Constructor
+  RasterData(size_t h_res, size_t v_res, bool include_filter = false);
+
+  // Methods
+  void set_value(size_t y, size_t x, const Particle& p, int level,
+    Filter* filter = nullptr, FilterMatch* match = nullptr);
+  void set_overlap(size_t y, size_t x);
+
+  // Members
+  tensor::Tensor<int32_t>
+    id_data_; //!< [v_res, h_res, 3 or 4]: cell, instance, mat, [filter_bin]
+  tensor::Tensor<double>
+    property_data_;     //!< [v_res, h_res, 2]: temperature, density
+  bool include_filter_; //!< Whether filter bin index is included
+};
+
 //===============================================================================
 // Plot class
 //===============================================================================
 
 class SlicePlotBase {
 public:
   template<class T>
-  T get_map() const;
+  T get_map(int32_t filter_index = -1) const;
 
   enum class PlotBasis { xy = 1, xz = 2, yz = 3 };
 
@@ -188,70 +209,65 @@ class SlicePlotBase {
 
   // Members
 public:
-  Position origin_;           //!< Plot origin in geometry
-  Position width_;            //!< Plot width in geometry
-  PlotBasis basis_;           //!< Plot basis (XY/XZ/YZ)
-  array<size_t, 3> pixels_;   //!< Plot size in pixels
-  bool slice_color_overlaps_; //!< Show overlapping cells?
-  int slice_level_ {-1};      //!< Plot universe level
+  Position origin_;         //!< Plot origin in geometry
+  Direction u_span_;        //!< Full-width span vector in geometry
+  Direction v_span_;        //!< Full-height span vector in geometry
+  array<size_t, 3> pixels_; //!< Plot size in pixels
+  bool show_overlaps_;      //!< Show overlapping cells?
+  int slice_level_ {-1};    //!< Plot universe level
 private:
 };
 
 template<class T>
-T SlicePlotBase::get_map() const
+T SlicePlotBase::get_map(int32_t filter_index) const
 {
 
   size_t width = pixels_[0];
   size_t height = pixels_[1];
 
-  // get pixel size
-  double in_pixel = (width_[0]) / static_cast<double>(width);
-  double out_pixel = (width_[1]) / static_cast<double>(height);
+  // Determine if filter is being used
+  bool include_filter = (filter_index >= 0);
+  Filter* filter = nullptr;
+  if (include_filter) {
+    filter = model::tally_filters[filter_index].get();
+  }
 
   // size data array
-  T data(width, height);
-
-  // setup basis indices and initial position centered on pixel
-  int in_i, out_i;
-  Position xyz = origin_;
-  switch (basis_) {
-  case PlotBasis::xy:
-    in_i = 0;
-    out_i = 1;
-    break;
-  case PlotBasis::xz:
-    in_i = 0;
-    out_i = 2;
-    break;
-  case PlotBasis::yz:
-    in_i = 1;
-    out_i = 2;
-    break;
-  default:
-    UNREACHABLE();
-  }
+  T data(width, height, include_filter);
 
-  // set initial position
-  xyz[in_i] = origin_[in_i] - width_[0] / 2. + in_pixel / 2.;
-  xyz[out_i] = origin_[out_i] + width_[1] / 2. - out_pixel / 2.;
+  // compute pixel steps and top-left pixel center
+  Direction u_step = u_span_ / static_cast<double>(width);
+  Direction v_step = v_span_ / static_cast<double>(height);
+
+  Position start =
+    origin_ - 0.5 * u_span_ + 0.5 * v_span_ + 0.5 * u_step - 0.5 * v_step;
+
+  // Validate that span vectors define a valid plane
+  Position cross = u_span_.cross(v_span_);
+  if (cross.norm() == 0.0) {
+    fatal_error("Slice span vectors are invalid (zero area).");
+  }
 
-  // arbitrary direction
-  Direction dir = {1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.0};
+  // Use an arbitrary direction that is not aligned with any coordinate axis.
+  // The direction has no physical meaning for plotting but is used by
+  // Surface::sense() to break ties when a pixel is coincident with a surface.
+  Direction dir = {1.0 / std::sqrt(2.0), 1.0 / std::sqrt(2.0), 0.0};
 
 #pragma omp parallel
   {
-    GeometryState p;
-    p.r() = xyz;
+    Particle p;
+    p.r() = start;
     p.u() = dir;
     p.coord(0).universe() = model::root_universe;
     int level = slice_level_;
     int j {};
+    FilterMatch match;
 
 #pragma omp for
     for (int y = 0; y < height; y++) {
-      p.r()[out_i] = xyz[out_i] - out_pixel * y;
+      Position row = start - v_step * static_cast<double>(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<double>(x);
         p.n_coord() = 1;
         // local variables
         bool found_cell = exhaustive_find_cell(p);
@@ -260,9 +276,9 @@ T SlicePlotBase::get_map() const
           j = level;
         }
         if (found_cell) {
-          data.set_value(y, x, p, j);
+          data.set_value(y, x, p, j, filter, &match);
         }
-        if (slice_color_overlaps_ && check_cell_overlap(p, false)) {
+        if (show_overlaps_ && check_cell_overlap(p, false)) {
           data.set_overlap(y, x);
         }
       } // inner for
@@ -297,6 +313,8 @@ class Plot : public PlottableInterface, public SlicePlotBase {
   void print_info() const override;
 
   PlotType type_;                 //!< Plot type (Slice/Voxel)
+  Position width_;                //!< Axis-aligned width from plot.xml
+  PlotBasis basis_;               //!< Basis from plot.xml for slice plots
   int meshlines_width_;           //!< Width of lines added to the plot
   int index_meshlines_mesh_ {-1}; //!< Index of the mesh to draw on the plot
   RGBColor meshlines_color_;      //!< Color of meshlines on the plot

openmc/deplete/independent_operator.py

@@ -339,8 +339,12 @@ def get_material_rates(self, mat_index, nuc_index, react_index):
 
             for i_nuc in nuc_index:
                 nuc = self.nuc_ind_map[i_nuc]
+                if nuc not in xs._index_nuc:
+                    continue
                 for i_rx in react_index:
                     rx = self.rx_ind_map[i_rx]
+                    if rx not in xs._index_rx:
+                        continue
 
                     # Determine reaction rate by multiplying xs in [b] by flux
                     # in [n-cm/src] to give [(reactions/src)*b-cm/atom]