Add TimedMeshFilter

CMakeLists.txt

@@ -423,6 +423,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_sptl_legendre.cpp
   src/tallies/filter_surface.cpp
   src/tallies/filter_time.cpp
+  src/tallies/filter_timed_mesh.cpp
   src/tallies/filter_universe.cpp
   src/tallies/filter_zernike.cpp
   src/tallies/tally.cpp

include/openmc/tallies/filter.h

@@ -43,6 +43,7 @@ enum class FilterType {
   SPATIAL_LEGENDRE,
   SURFACE,
   TIME,
+  TIMED_MESH,
   UNIVERSE,
   ZERNIKE,
   ZERNIKE_RADIAL

include/openmc/tallies/filter_timed_mesh.h

@@ -0,0 +1,74 @@
+#ifndef OPENMC_TALLIES_FILTER_TIMED_MESH_H
+#define OPENMC_TALLIES_FILTER_TIMED_MESH_H
+
+#include <cstdint>
+
+#include "openmc/constants.h"
+#include "openmc/position.h"
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Indexes the time-location of particle events to a time gridn and a regular
+//  mesh. For tracklength tallies, it will produce multiple valid bins and the
+//  bin weight will correspond to the fraction of the track length that lies in
+//! that bin.
+//==============================================================================
+
+class TimedMeshFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~TimedMeshFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "timedmesh"; }
+  FilterType type() const override { return FilterType::TIMED_MESH; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  int32_t mesh() const { return mesh_; }
+
+  void set_mesh(int32_t mesh);
+
+  void set_translation(const Position& translation);
+
+  void set_translation(const double translation[3]);
+
+  const Position& translation() const { return translation_; }
+
+  bool translated() const { return translated_; }
+
+  const vector<double>& time_grid() const { return time_grid_; }
+
+  void set_time_grid(gsl::span<const double> time_grid);
+
+  void reset_bins();
+
+protected:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  int32_t mesh_; //!< Index of the mesh
+  int mesh_n_bins_;
+  bool translated_ {false}; //!< Whether or not the filter is translated
+  Position translation_ {0.0, 0.0, 0.0}; //!< Filter translation
+  vector<double> time_grid_ {0.0, INFTY};
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_TIMED_MESH_H

openmc/filter.py

@@ -25,7 +25,7 @@
     'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
     'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
-    'collision', 'time'
+    'collision', 'time', 'timedmesh'
 )
 
 _CURRENT_NAMES = (
@@ -1119,6 +1119,157 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
         return pd.concat([df, pd.DataFrame(filter_dict)])
 
 
+class TimedMeshFilter(Filter):
+    """Bins tally event locations by mesh elements and time grids.
+
+    Parameters
+    ----------
+    mesh : openmc.MeshBase
+        The mesh object that events will be tallied onto
+    time_grid : iterable of float
+        A list of values for which each successive pair constitutes a range of
+        time in [s] for a single time bin
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    mesh : openmc.MeshBase
+        The mesh object that events will be tallied onto
+    time_grid : numpy.ndarray
+        An array of values for which each successive pair constitutes a range of
+        time in [s] for a single time bin
+    id : int
+        Unique identifier for the filter
+    translation : Iterable of float
+        This array specifies a vector that is used to translate (shift) the mesh
+        for this filter
+    bins : list of tuple
+        A list of mesh indices for each filter bin, e.g. [(1, 1, 1), (2, 1, 1),
+        ...]
+    num_bins : Integral
+        The number of filter bins
+
+    """
+
+    def __init__(self, mesh, time_grid, filter_id=None):
+        self._time_grid = np.array([0.0, np.inf])  # Temporary
+        self.mesh = mesh
+        self.time_grid = time_grid
+        self.id = filter_id
+        self._translation = None
+
+    def __hash__(self):
+        string = type(self).__name__ + '\n'
+        string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id)
+        string += '{: <16}=\t{}\n'.format('\tTime grid', self.time_grid)
+        return hash(string)
+
+    def __repr__(self):
+        string = type(self).__name__ + '\n'
+        string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id)
+        string += '{: <16}=\t{}\n'.format('\tTime grid', self.time_grid)
+        string += '{: <16}=\t{}\n'.format('\tID', self.id)
+        string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation)
+        return string
+
+    @classmethod
+    def from_hdf5(cls, group, **kwargs):
+        if group['type'][()].decode() != cls.short_name.lower():
+            raise ValueError("Expected HDF5 data for filter type '"
+                             + cls.short_name.lower() + "' but got '"
+                             + group['type'][()].decode() + " instead")
+
+        if 'meshes' not in kwargs:
+            raise ValueError(cls.__name__ + " requires a 'meshes' keyword "
+                             "argument.")
+
+        mesh_id = group['mesh_bins'][()]
+        mesh_obj = kwargs['meshes'][mesh_id]
+        filter_id = int(group.name.split('/')[-1].lstrip('filter '))
+
+        time_grid = group['time_bins'][()]
+
+        out = cls(mesh_obj, time_grid, filter_id=filter_id)
+
+        translation = group.get('translation')
+        if translation:
+            out.translation = translation[()]
+
+        return out
+
+    @property
+    def mesh(self):
+        return self._mesh
+
+    @property
+    def time_grid(self):
+        return self._time_grid
+
+    @mesh.setter
+    def mesh(self, mesh):
+        cv.check_type('filter mesh', mesh, openmc.MeshBase)
+        self._mesh = mesh
+        self._reset_bins()
+
+    @time_grid.setter
+    def time_grid(self, time_grid):
+        self._time_grid = time_grid
+        self._reset_bins()
+
+    def _reset_bins(self):
+        mesh = self._mesh
+        time_grid = self._time_grid
+
+        # Mesh bins
+        if isinstance(mesh, openmc.UnstructuredMesh):
+            if mesh.has_statepoint_data:
+                mesh_bins = list(range(len(mesh.volumes)))
+            else:
+                mesh_bins = []
+        else:
+            mesh_bins = list(mesh.indices)
+
+        # Reset bins
+        self.bins = []
+        for i in range(len(self._time_grid) - 1):
+            for j in range(len(mesh_bins)):
+                new_bin = ()
+                new_bin += ([float(time_grid[i]), float(time_grid[i+1])],)
+                new_bin += (mesh_bins[j],)
+                self.bins.append(new_bin)
+
+    @property
+    def shape(self):
+        return (self.time_grid - 1,) + self.mesh.dimension
+
+    @property
+    def translation(self):
+        return self._translation
+
+    @translation.setter
+    def translation(self, t):
+        cv.check_type('mesh filter translation', t, Iterable, Real)
+        cv.check_length('mesh filter translation', t, 3)
+        self._translation = np.asarray(t)
+
+    def can_merge(self, other):
+        return False
+
+    def to_xml_element(self):
+        element = super().to_xml_element()
+        element[0].text = None
+
+        time_bins = ET.SubElement(element[0], 'time_bins')
+        mesh_bins = ET.SubElement(element[0], 'mesh_bins')
+
+        time_bins.text = ' '.join(str(x) for x in self.time_grid)
+        mesh_bins.text = str(self.mesh.id)
+        if self.translation is not None:
+            element.set('translation', ' '.join(map(str, self.translation)))
+        return element
+
+
 class CollisionFilter(Filter):
     """Bins tally events based on the number of collisions.
 

openmc/tallies.py

@@ -3174,7 +3174,8 @@ def _create_mesh_subelements(self, root_element, memo=None):
         already_written = memo if memo else set()
         for tally in self:
             for f in tally.filters:
-                if isinstance(f, openmc.MeshFilter):
+                if isinstance(f, openmc.MeshFilter) or\
+                isinstance(f, openmc.TimedMeshFilter):
                     if f.mesh.id in already_written:
                         continue
                     if len(f.mesh.name) > 0:

src/particle.cpp

@@ -44,6 +44,15 @@ namespace openmc {
 
 double Particle::speed() const
 {
+  // Multigroup speed?
+  if (!settings::run_CE) {
+    auto& macro_xs = data::mg.macro_xs_[this->material()];
+    int macro_t = this->mg_xs_cache().t;
+    int macro_a = macro_xs.get_angle_index(this->u());
+    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr,
+                   nullptr, nullptr, macro_t, macro_a);
+  }
+
   // Determine mass in eV/c^2
   double mass;
   switch (this->type()) {
@@ -244,6 +253,9 @@ void Particle::event_advance()
     double push_back_distance = speed() * dt;
     this->move_distance(-push_back_distance);
     hit_time_boundary = true;
+
+    // Reduce the distance traveled for tallying
+    distance -= push_back_distance;
   }
 
   // Score track-length tallies

src/physics_mg.cpp

@@ -127,6 +127,7 @@ void create_fission_sites(Particle& p)
     SourceSite site;
     site.r = p.r();
     site.particle = ParticleType::neutron;
+    site.time = p.time();
     site.wgt = 1. / weight;
     site.parent_id = p.id();
     site.progeny_id = p.n_progeny()++;

src/tallies/filter.cpp

@@ -33,6 +33,7 @@
 #include "openmc/tallies/filter_sptl_legendre.h"
 #include "openmc/tallies/filter_surface.h"
 #include "openmc/tallies/filter_time.h"
+#include "openmc/tallies/filter_timed_mesh.h"
 #include "openmc/tallies/filter_universe.h"
 #include "openmc/tallies/filter_zernike.h"
 #include "openmc/xml_interface.h"
@@ -148,6 +149,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<SphericalHarmonicsFilter>(id);
   } else if (type == "time") {
     return Filter::create<TimeFilter>(id);
+  } else if (type == "timedmesh") {
+    return Filter::create<TimedMeshFilter>(id);
   } else if (type == "universe") {
     return Filter::create<UniverseFilter>(id);
   } else if (type == "zernike") {