Merge 7f497ff04f into 61c8a59cff

CMakeLists.txt

@@ -382,6 +382,7 @@ list(APPEND libopenmc_SOURCES
   src/ifp.cpp
   src/initialize.cpp
   src/lattice.cpp
+  src/majorant.cpp
   src/material.cpp
   src/math_functions.cpp
   src/mcpl_interface.cpp

docs/source/io_formats/statepoint.rst

@@ -28,6 +28,10 @@ The current version of the statepoint file format is 18.2.
              'continuous-energy' or 'multi-group'.
            - **run_mode** (*char[]*) -- Run mode used, either 'eigenvalue' or
              'fixed source'.
+           - **photon_transport** (*bool*) -- Whether photon transport was enabled
+             or not.
+           - **delta_tracking** (*bool*) -- Whether delta tracking was enabled
+             or not.
            - **n_particles** (*int8_t*) -- Number of particles used per generation.
            - **n_batches** (*int*) -- Number of batches to simulate.
            - **current_batch** (*int*) -- The number of batches already simulated.
@@ -178,6 +182,8 @@ All values are given in seconds and are measured on the master process.
              allocating arrays, etc.
            - **reading cross sections** (*double*) -- Time spent loading cross
              section libraries (this is a subset of initialization).
+           - **build majorant** (*double*) -- Time spent constructing the majorant
+             cross sections. This is only saved if running with delta tracking.
            - **simulation** (*double*) -- Time spent between initialization and
              finalization.
            - **transport** (*double*) -- Time spent transporting particles.

docs/source/methods/neutron_physics.rst

@@ -52,6 +52,111 @@ the formula usually used to calculate the distance to next collision is
 
     \ell = -\frac{\ln \xi}{\Sigma_t}
 
+.. _surface_tracking:
+
+~~~~~~~~~~~~~~~~~~~~~~~~~
+Surface Tracking
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The development of Equation :eq:`sample-distance-2` requires the assumption
+that the medium under consideration is homogeneous. To accomodate heterogeneous
+geometries, a resampling scheme is used. First, the distance to the next
+collision is sampled with Equation :eq:`sample-distance-2`. Then, the distance
+to the nearest surface from the particle position along its current trajectory
+is computed (discussed in the :ref:`methods_geometry` section). If the distance
+to the nearest surface is smaller than the distance to the next collision, the
+sampled distance is not statistically valid. The particle is moved to the
+surface and is considered to be contained by the next geometric region. Cross
+sections are recomputed, and a new distance to the next collision is sampled
+with Equation :eq:`sample-distance-2`. This process repeats until the distance
+to the next collision is smaller than the distance to the nearest surface,
+which is when a collision is accepted. This procedure is known as surface
+tracking.
+
+Surface tracking is quite efficient when used in problems with short mean free
+paths relative to the size of individual regions in the problem geometry.
+Surface tracking also admits the use of the track length estimator (discussed
+in the :ref:`methods_tallies` section). In problems with long mean free paths
+relative to the size of geometry regions, surface tracking will require a large
+number of surface distance calculations per collision. The cost of finding the
+nearest surface is also non-trivial for problems that contain many geometric
+regions at the same cell level (e.g. TRISO-fueled fission reactors).
+
+.. _delta_tracking:
+
+~~~~~~~~~~~~~~~~~~~~~~~~~
+Delta Tracking
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The disadvantages of surface tracking for certain classes of problems motivates
+the development of alternative approaches which do not require distance
+to surface checks. Delta tracking (also known by Woodcock tracking,
+delta scattering, and null scattering) is one approach to
+avoid surface geometry queries [Woodcock]_. In delta tracking, the domain is
+homogenized to a majorant cross section (:math:`\Sigma_{maj}(E)`) which is
+computed as the maximum total cross section over the entire problem
+
+.. math::
+    :label: majorant-xs-1
+
+    \Sigma_{maj}(E) = \max_{\mathbf{r}}\left(\Sigma_{t}(\mathbf{r}, E)\right).
+
+Particles move through this homogenized problem by sampling a distance to the
+next collision with the majorant cross section instead of the total cross section
+
+.. math::
+    :label: sample-distance-maj
+
+    \ell = -\frac{\ln \xi}{\Sigma_{maj}(E)}.
+
+To recover the spatial heterogeneity present in the original simulation, the
+delta tracking method formally defines the majorant cross section to be the
+sum of the total cross section and a spatially-varying fictitious delta
+scattering cross sections :math:`\Sigma_{\delta}(\mathbf{r}, E)`
+
+.. math::
+    :label: majorant-xs-2
+
+    \Sigma_{maj}(E) = \Sigma_{t}(\mathbf{r}, E) + \Sigma_{\delta}(\mathbf{r},
+    E).
+
+At the collision point computed from :eq:`sample-distance-maj` one of
+two reactions could occur. The first is a real collision, which is processed
+as usual. The second is known as a delta scatter event (also referred to
+as a virtual or null collision), which is the reaction type associated with
+:math:`\Sigma_{\delta}(\mathbf{r}, E)`. A rejection sampling test is used to
+determine which collision type occurs; a random number :math:`\xi` on the
+interval :math:`[0,1)` is drawn and used to check
+
+.. math::
+    :label: delta-real-collision
+
+    \xi < \frac{\Sigma_t (\mathbf{r}, E)}{\Sigma_{maj} (E)}.
+
+If the condition above is true, the collision is accepted as real. If the condition
+is false, a delta scatter event has occurred and the particle continues along
+its trajectory with the same energy and direction. Boundary conditions
+are applied by testing the distance to the nearest external boundary and
+comparing this to the distance sampled with Equation :eq:`sample-distance-maj`.
+If the distance to the nearest boundary is less than the sampled distance to
+the next collision, the particle crosses the external boundary.
+
+Delta tracking is advantageous as it only requires point location checks to
+determine the total cross section at each collision point to test
+:eq:`delta-real-collision`. This allows for the use of continuously-varying
+material properties and avoids computationally expensive distance-to-nearest-surface
+calculations. Problems which contain small regions with large total
+cross sections (such as burnable absorbers) will have majorant cross sections
+several orders of magnitude larger than the total cross section over the majority
+of the domain [Leppänen]_. This decreases the number of real collisions, and
+therefore the effectiveness of delta tracking. Material discontinuities are not
+considered in delta tracking, which prohibits the use of track length
+estimators for quantities restricted to material/geometric subdomains (such as
+reaction rates) and forces the use of the higher-variance collision
+estimator (discussed in detail in the :ref:`methods_tallies` section). When compared
+with surface tracking, delta tracking often performs better in problems where
+the particle mean free path is larger than the distance between surfaces.
+
 ----------------------------------------------------
 :math:`(n,\gamma)` and Other Disappearance Reactions
 ----------------------------------------------------
@@ -1631,6 +1736,13 @@ the unresolved range to get the actual cross sections. Lastly, the total cross
 section is calculated as the sum of the elastic, fission, capture, and inelastic
 cross sections.
 
+Unresolved resonance probability tables pose a challenge when computing a majorant
+cross section for :ref:`delta_tracking`. OpenMC implements a conservative approach:
+the maximum total cross section is computed over all bands, which is then
+interpolated to the corresponding energy using either linear or logarithmic
+interpolation. This ensures the majorant bounds the total cross section at the
+cost of increasing the number of delta scatters in the unresolved range.
+
 -----------------------------
 Variance Reduction Techniques
 -----------------------------
@@ -1735,6 +1847,10 @@ types.
 .. [Gelbard] Ely M. Gelbard, "Epithermal Scattering in VIM," FRA-TM-123, Argonne
    National Laboratory (1979).
 
+.. [Leppänen] J. Leppänen. "Performance of Woodcock Delta-Tracking in Lattice
+   Physics Applications using the Serpent Monte Carlo Reactor Physics Burnup
+   Calculation Code", *Annals of Nuclear Energy*, 37:715-722, 2010.
+
 .. [Squires] G. L. Squires, *Introduction to the Theory of Thermal Neutron
    Scattering*, Cambridge University Press (1978).
 
@@ -1742,6 +1858,11 @@ types.
    Neutrons*, North-Holland Publishing Co., Amsterdam (1966). **Note:** This
    book can be obtained for free from the OECD_.
 
+.. [Woodcock] E.R. Woodcock, T. Murphy, P.J. Hemmings, and T.C. Longworth.
+   "Techniques used in the GEM Code for Monte Carlo Neutronics Calculations
+   in Reactors and other Systems of Complex Geometry", ANL-7050,
+   Argonne National Laboratory (1965).
+
 .. |sab| replace:: S(:math:`\alpha,\beta,T`)
 
 .. _SIGMA1 method: https://doi.org/10.13182/NSE76-1

include/openmc/event.h

@@ -87,6 +87,15 @@ void free_event_queues(void);
 //! \param buffer_idx The particle's actual index in the particle buffer
 void dispatch_xs_event(int64_t buffer_idx);
 
+//! Execute the death event for all particles
+//
+//! \param n_particles The number of particles in the particle buffer
+void process_death_events(int64_t n_particles);
+
+//==============================================================================
+// Surface tracking
+//==============================================================================
+
 //! Execute the initialization event for all particles
 //
 //! \param n_particles The number of particles in the particle buffer
@@ -107,11 +116,6 @@ void process_surface_crossing_events();
 //! Execute the collision event for all particles in this event's buffer
 void process_collision_events();
 
-//! Execute the death event for all particles
-//
-//! \param n_particles The number of particles in the particle buffer
-void process_death_events(int64_t n_particles);
-
 //! Process event queues until all are empty. Each iteration processes the
 //! longest queue first to maximize vectorization efficiency.
 void process_transport_events();
@@ -125,6 +129,43 @@ void process_transport_events();
 void process_init_secondary_events(int64_t n_particles, int64_t offset,
   const SharedArray<SourceSite>& shared_secondary_bank);
 
+//==============================================================================
+// Delta tracking
+//==============================================================================
+
+//! Specialization of process_init_events() for delta tracking.
+//
+//! \param n_particles The number of particles in the particle buffer
+//! \param source_offset The offset index in the source bank to use
+void process_delta_init_events(int64_t n_particles, int64_t source_offset);
+
+//! Specialization of process_calculate_xs_events() for delta tracking.
+//
+//! \param queue A reference to the desired XS lookup queue
+void process_delta_calculate_xs_events(SharedArray<EventQueueItem>& queue);
+
+//! Execute the delta advance particle event for all particles in this advance
+//! buffer
+void process_delta_advance_particle_events();
+
+//! Specialization of process_surface_crossing_events() for delta tracking.
+void process_delta_surface_crossing_events();
+
+//! Execute the delta tracking collision event for all particles in this event's
+//! buffer
+void process_delta_collision_events();
+
+//! Specialization of process_transport_events() for delta tracking.
+void process_delta_transport_events();
+
+//! Specialization of process_init_secondary_events() for delta tracking.
+//
+//! \param n_particles The number of particles to initialize
+//! \param offset The offset index in the shared secondary bank
+//! \param shared_secondary_bank The shared secondary bank to read from
+void process_delta_init_secondary_events(int64_t n_particles, int64_t offset,
+  const SharedArray<SourceSite>& shared_secondary_bank);
+
 } // namespace openmc
 
 #endif // OPENMC_EVENT_H

include/openmc/geometry.h

@@ -53,6 +53,8 @@ namespace model {
 extern int root_universe;      //!< Index of root universe
 extern "C" int n_coord_levels; //!< Number of CSG coordinate levels
 
+extern vector<int> boundary_surfaces; //!< The surfaces with boundary conditions
+
 extern vector<int64_t> overlap_check_count;
 
 // Overlap data structures get cleared every slice_data run
@@ -107,10 +109,21 @@ void cross_lattice(
 
 //==============================================================================
 //! Find the next boundary a particle will intersect.
+//! \param p A geometry state to compute distances with.
+//! \return Boundary information corresponding to the nearest surface.
 //==============================================================================
 
 BoundaryInfo distance_to_boundary(GeometryState& p);
 
+//==============================================================================
+//! Find the next external boundary a particle will intersect.
+//!
+//! \param p A geometry state to compute distances with.
+//! \return Boundary information corresponding to the nearest external boundary.
+//==============================================================================
+
+BoundaryInfo distance_to_external_boundary(GeometryState& p);
+
 } // namespace openmc
 
 #endif // OPENMC_GEOMETRY_H

include/openmc/majorant.h

@@ -0,0 +1,312 @@
+//! \file majorant.h
+//! Majorant cross section type
+
+#ifndef OPENMC_MAJORANT_H
+#define OPENMC_MAJORANT_H
+
+#include <vector>
+
+#include "openmc/nuclide.h"
+
+namespace openmc {
+
+class NeutronMajorant;
+class PhotonMajorant;
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace data {
+extern std::unique_ptr<NeutronMajorant> n_majorant;
+extern std::unique_ptr<PhotonMajorant> p_majorant;
+} // namespace data
+
+//==============================================================================
+
+class Majorant {
+public:
+  //----------------------------------------------------------------------------
+  // Constructor
+
+  Majorant(int i_universe);
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  //! A function to unionize particle energy grids.
+  virtual void compute_unionized_grid() = 0;
+
+  //! A function to populate the majorant cross section.
+  void compute_majorant();
+
+protected:
+  //----------------------------------------------------------------------------
+  // Protected Methods
+
+  //! Virtual function for the minimum transport energy in the majorant.
+  //!
+  //! \return The minimum transport energy associated with the majorant
+  virtual double min_transport_energy() const = 0;
+
+  //! Virtual function for the maximum transport energy in the majorant.
+  //!
+  //! \return The maximum transport energy associated with the majorant
+  virtual double max_transport_energy() const = 0;
+
+  //! Compute a per-material macroscopic majorant cross section in units
+  //! of [cm^-1]
+  //!
+  //! \param[in] mat The material to compute the majorant cross section of
+  //! \param[in] to_grid The grid points to evaluate the majorant at in [eV]
+  //! \param[out] mat_maj The array to write the macroscopic majorant to.
+  //!   The resulting cross section has units of [cm^-1]
+  virtual void fill_material_maj_xs(int i_material, double max_density_mult,
+    const std::vector<double>& to_grid, std::vector<double>& mat_maj) const = 0;
+
+  //! Post-processes the energy grid by calling std::sort(), std::unique().
+  //! This also removes all energy values below the transport minimum and
+  //! above the transport maximum for a given particle type.
+  //!
+  //! \param[out] grid The energy grid to post-process. This is performed
+  //!   in-place.
+  void post_process_grid(Nuclide::EnergyGrid& grid) const;
+
+  //----------------------------------------------------------------------------
+  // Protected data members
+
+  //! Index into the universe array for the universe which this majorant uses
+  //! to fetch material properties.
+  int maj_universe_ = C_NONE;
+  //!< A vector of materials contained in maj_universe_.
+  std::vector<int> contained_materials_;
+  //! A map of each material index and the corresponding maximum density
+  //! multiplier applied to that material by a cell.
+  std::unordered_map<int, double> max_density_mult_;
+  //! The unionized energy grid.
+  Nuclide::EnergyGrid grid_;
+  //! Macroscopic majorant cross sections at each energy point in grid_.
+  std::vector<double> xs_;
+}; // class Majorant
+
+//==============================================================================
+
+class NeutronMajorant : public Majorant {
+
+public:
+  //----------------------------------------------------------------------------
+  // Constructors
+
+  NeutronMajorant(int i_universe);
+
+  //----------------------------------------------------------------------------
+  // Public Methods
+
+  virtual void compute_unionized_grid() override final;
+
+  //! Calculate the macroscopic majorant cross section.
+  //!
+  //! \param[in] energy The energy to compute the cross section at in [eV]
+  //! \return The neutron majorant cross section at the given energy in [cm^-1]
+  double calculate_neutron_xs(double energy) const;
+
+  //----------------------------------------------------------------------------
+  // Data members
+
+  //! A dilation factor to ensure floating-point round off and inexact majorant
+  //! URR and S(a,b) cross sections don't bias results.
+  constexpr static double safety_factor_ {1.01};
+
+protected:
+  //----------------------------------------------------------------------------
+  // Protected Methods
+
+  //! Minimum neutron transport energy.
+  //
+  //! \return The minimum transport energy associated with the majorant [eV]
+  virtual double min_transport_energy() const override
+  {
+    return data::energy_min[i_neutron_];
+  }
+
+  //! Maximum neutron transport energy.
+  //
+  //! \return The maximum transport energy associated with the majorant [eV]
+  virtual double max_transport_energy() const override
+  {
+    return data::energy_max[i_neutron_];
+  }
+
+  //! Compute a per-material macroscopic majorant cross section.
+  //!
+  //! \param[in] i_material Index into the materials array.
+  //! \param[in] to_grid The grid points to evaluate the majorant at in [eV].
+  //! \param[out] mat_maj The array to write the macroscopic majorant to.
+  //!   The resulting cross section has units of [cm^-1].
+  virtual void fill_material_maj_xs(int i_material, double max_density_mult,
+    const std::vector<double>& to_grid,
+    std::vector<double>& mat_maj) const override;
+
+private:
+  //----------------------------------------------------------------------------
+  // Private Methods
+
+  //! Compute the maximum smooth microscopic total cross section.
+  //!
+  //! \param[in] i_nuclide Index into the nuclides array.
+  //! \param[in] energy The energy to evaluate the cross section at in [eV].
+  //! \return The maximum smooth cross section in [barn].
+  double calculate_max_smooth_xs(int i_nuclide, double energy) const;
+
+  //! Compute the maximum microscopic total URR cross section.
+  //!
+  //! \param[in] i_nuclide Index into the nuclides array.
+  //! \param[in] energy The energy to evaluate the cross section at in [eV].
+  //! \param[in] smooth_xs The smooth total cross section in units of [eV] to
+  //!   use (if needed by the ptable).
+  //! \return The maximum URR total cross section in [barn].
+  double calculate_max_urr_xs(
+    int i_nuclide, double energy, double smooth_xs) const;
+
+  //! Compute the maximum microscopic S(a,b) total cross section.
+  //!
+  //! \param[in] energy The energy to evaluate the cross section at in [eV].
+  //! \param[in] i_sab The index into the thermal scattering table array for
+  //!   this nuclide.
+  //! \param[in] sab_frac The fraction of the bound cross section to use vs the
+  //!   free gas cross section.
+  //! \param[in] nuc The nuclide to compute the microscopic total cross section
+  //!   of.
+  //! \return The maximum S(a,b) total cross section in [barn].
+  double calculate_max_sab_tot_xs(
+    int i_nuclide, int i_sab, double sab_frac, double energy) const;
+
+  //! Get the grid index for energy interpolation.
+  //!
+  //! \param[in] energy The energy to evaluate the cross section at in [eV].
+  //! \param[in] grid The energy grid to search for an energy grid index.
+  //! \return The grid index.
+  int get_i_grid(double energy, const Nuclide::EnergyGrid& grid) const;
+
+  //----------------------------------------------------------------------------
+  // Private data members
+
+  static constexpr int i_neutron_ = ParticleType::neutron().transport_index();
+}; // class NeutronMajorant
+
+//==============================================================================
+
+class PhotonMajorant : public Majorant {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors
+
+  PhotonMajorant(int i_universe);
+
+  //----------------------------------------------------------------------------
+  // Public Methods
+
+  virtual void compute_unionized_grid() override final;
+
+  //! Calculate the macroscopic majorant cross section.
+  //!
+  //! \param[in] energy The energy to compute the cross section at in [eV]
+  //! \return The photon majorant cross section at the given energy in [cm^-1]
+  double calculate_photon_xs(double energy) const;
+
+  //----------------------------------------------------------------------------
+  // Data members
+
+  //! A dilation factor to catch interpolation error.
+  constexpr static double safety_factor_ {1.01};
+
+protected:
+  //----------------------------------------------------------------------------
+  // Protected Methods
+
+  //! Minimum photon transport energy.
+  //
+  //! \return The log of the minimum transport energy associated with the
+  //! majorant
+  virtual double min_transport_energy() const override
+  {
+    return std::log(data::energy_min[i_photon_]);
+  }
+
+  //! Maximum photon transport energy.
+  //
+  //! \return The log of the maximum transport energy associated with the
+  //! majorant
+  virtual double max_transport_energy() const override
+  {
+    return std::log(data::energy_max[i_photon_]);
+  }
+
+  //! Compute a per-material macroscopic majorant cross section.
+  //!
+  //! \param[in] i_material Index into the materials array
+  //! \param[in] to_grid The grid points to evaluate the majorant at in [eV]
+  //! \param[out] mat_maj The array to write the macroscopic majorant to. The
+  //!   resulting cross section has units of [cm^-1]
+  virtual void fill_material_maj_xs(int i_material, double max_density_mult,
+    const std::vector<double>& to_grid,
+    std::vector<double>& mat_maj) const override;
+
+private:
+  //----------------------------------------------------------------------------
+  // Private Methods
+
+  //! Compute the maximum smooth microscopic total cross section.
+  //!
+  //! \param[in] i_element Index in the elements array.
+  //! \param[in] energy The energy to evaluate the cross section at in [eV].
+  //! \return The maximum microscopic total cross section in [barn].
+  double calculate_elem_tot_xs(int i_element, double log_energy) const;
+
+  //! Get the grid index for energy interpolation. This is templated due to the
+  //! use of tensor::Tensor<double> in PhotonInteraction.
+  //!
+  //! \param[in] energy The energy to evaluate the cross section at in [eV].
+  //! \param[in] grid The energy grid to search for an energy grid index.
+  //! \return The grid index.
+  template<typename T>
+  int get_i_grid(double log_energy, const T& energy_grid) const
+  {
+    int n_grid = energy_grid.size();
+    int i_grid;
+    if (log_energy <= energy_grid[0]) {
+      i_grid = 0;
+    } else if (log_energy > energy_grid[n_grid - 1]) {
+      i_grid = n_grid - 2;
+    } else {
+      // We use upper_bound_index here because sometimes photons are created
+      // with energies that exactly match a grid point
+      i_grid =
+        upper_bound_index(energy_grid.cbegin(), energy_grid.cend(), log_energy);
+    }
+
+    // check for case where two energy points are the same
+    if (energy_grid[i_grid] == energy_grid[i_grid + 1])
+      ++i_grid;
+
+    return i_grid;
+  }
+
+  //----------------------------------------------------------------------------
+  // Private data members
+
+  static constexpr int i_photon_ = ParticleType::photon().transport_index();
+}; // class PhotonMajorant
+
+//==============================================================================
+// Static functions
+//==============================================================================
+
+//! A function to create majorant cross sections.
+void create_majorants();
+
+//! A function to reset majorant cross sections.
+void reset_majorants();
+} // namespace openmc
+
+#endif // OPENMC_MAJORANT_H

include/openmc/nuclide.h

@@ -33,6 +33,10 @@ public:
   // Types, aliases
   using EmissionMode = ReactionProduct::EmissionMode;
   struct EnergyGrid {
+    // Init method.
+    void init();
+
+    // Data members.
     vector<int> grid_index;
     vector<double> energy;
   };

include/openmc/particle.h

@@ -69,6 +69,7 @@ public:
   // Coarse-grained particle events
   void event_calculate_xs();
   void event_advance();
+  void event_delta_advance();
   void event_cross_surface();
   void event_collide();
   void event_revive_from_secondary(const SourceSite& site);
@@ -111,6 +112,16 @@ public:
   virtual void mark_as_lost(const char* message) override;
   using GeometryState::mark_as_lost;
 
+  //! Update the cached majorant cross section for this particle.
+  void update_majorant();
+
+  //! Check if the majorant is valid. If its not, the particle is killed and
+  //! a restart file is written.
+  //
+  //! \return true if the majorant is invalid, false if not. If true,
+  //!   the particle is also killed.
+  bool kill_invalid_maj();
+
   //! create a particle restart HDF5 file
   void write_restart() const;
 

include/openmc/particle_data.h

@@ -534,7 +534,7 @@ private:
   bool write_track_ {false};
 
   uint64_t seeds_[N_STREAMS];
-  int stream_;
+  int stream_ {STREAM_TRACKING};
 
   vector<SourceSite> local_secondary_bank_;
 
@@ -573,6 +573,11 @@ private:
 
   int64_t n_progeny_ {0};
 
+  //! Flag to indicate whether or not delta tracking is active.
+  bool delta_tracking_ {false};
+  //! Most recent value of the majorant cross section.
+  double majorant_ {0.0};
+
 public:
   //----------------------------------------------------------------------------
   // Constructors
@@ -765,6 +770,14 @@ public:
   // Number of progeny produced by this particle
   int64_t& n_progeny() { return n_progeny_; }
 
+  //! Gets whether the particle is being tracked with delta tracking or not.
+  bool& delta_tracking() { return delta_tracking_; }
+  const bool& delta_tracking() const { return delta_tracking_; }
+
+  //! Gets the majorant cross section.
+  double& majorant() { return majorant_; }
+  const double& majorant() const { return majorant_; }
+
   //! Gets the pointer to the particle's current PRN seed
   uint64_t* current_seed() { return seeds_ + stream_; }
   const uint64_t* current_seed() const { return seeds_ + stream_; }

include/openmc/settings.h

@@ -67,6 +67,7 @@ extern bool create_delayed_neutrons; //!< create delayed fission neutrons?
 extern "C" bool cmfd_run;            //!< is a CMFD run?
 extern bool
   delayed_photon_scaling;   //!< Scale fission photon yield to include delayed
+extern bool delta_tracking; //!< use delta tracking
 extern "C" bool entropy_on; //!< calculate Shannon entropy?
 extern "C" bool
   event_based;      //!< use event-based mode (instead of history-based)

include/openmc/simulation.h

@@ -118,6 +118,9 @@ void transport_history_based();
 //! secondary bank
 void transport_history_based_shared_secondary();
 
+//! Simulate a single particle history from birth to death using delta tracking
+void transport_delta_history_based_single_particle(Particle& p);
+
 //! Simulate all particle histories using event-based parallelism
 void transport_event_based();
 

include/openmc/timer.h

@@ -21,6 +21,7 @@ extern Timer time_finalize;
 extern Timer time_inactive;
 extern Timer time_initialize;
 extern Timer time_read_xs;
+extern Timer time_build_majorant;
 extern Timer time_statepoint;
 extern Timer time_tallies;
 extern Timer time_total;

openmc/examples.py

@@ -6,6 +6,98 @@ import openmc
 
 PINCELL_PITCH = 1.26 # cm
 
+def delta_tracking_lattice(
+    run_photon: bool = False,
+    boundary_type: str = 'reflective',
+    densities: list[float] | None = None) -> openmc.Model:
+    """Create a simple PWR-style lattice for testing delta tracking.
+
+    Parameters
+    ----------
+    run_photon : bool, optional
+        If coupled neutron-photon transport should be run or not.
+    boundary_type : str, optional
+        The boundary to apply to the outer surface of the infinite lattice.
+    densities : list of float, optional
+        The distributed cell densities to apply to cell 1 (the fuel) in the
+        infinite lattice.
+
+    Returns
+    -------
+    model : openmc.Model
+        A PWR-style 2x2 infinite assembly model
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    DELTA_PIN_RADIUS = 0.4
+
+    # Create some simple materials. UO2 fuel for the inner cylinder in the pin,
+    # and water for the remainder of the domain.
+    uo2 = openmc.Material(name='UO2')
+    uo2.set_density('g/cm3', 10.0)
+    uo2.add_nuclide('U235', 1.0)
+    uo2.add_nuclide('O16', 2.0)
+    water = openmc.Material(name='light water')
+    water.add_nuclide('H1', 2.0)
+    water.add_nuclide('O16', 1.0)
+    water.set_density('g/cm3', 1.0)
+    water.add_s_alpha_beta('c_H_in_H2O')
+    model.materials.extend([uo2, water])
+
+    # Create the geometry, starting with the fuel pincell.
+    cyl = openmc.ZCylinder(r=DELTA_PIN_RADIUS)
+    pin = openmc.model.pin([cyl], [uo2, water])
+
+    # Create a 2x2 lattice to allow for distributed properties.
+    lattice = openmc.RectLattice()
+    lattice.lower_left = (-PINCELL_PITCH, -PINCELL_PITCH)
+    lattice.pitch = (PINCELL_PITCH, PINCELL_PITCH)
+    lattice.universes = [[pin, pin],
+                         [pin, pin]]
+    box = openmc.model.RectangularPrism(
+        2.0 * PINCELL_PITCH, 2.0 * PINCELL_PITCH,
+        origin=(0.0, 0.0),
+        boundary_type=boundary_type
+    )
+
+    # Set distributed densities if required.
+    if densities != None:
+        pin.cells[1].density = densities
+
+    # Finally, save the geometry.
+    model.geometry = openmc.Geometry([openmc.Cell(fill=lattice, region=-box)])
+    model.geometry.merge_surfaces = True
+
+    # Add a mesh tally for the neutron total reaction rate.
+    msh = openmc.RegularMesh()
+    msh.lower_left = (-PINCELL_PITCH, -PINCELL_PITCH)
+    msh.upper_right = (PINCELL_PITCH, PINCELL_PITCH)
+    msh.dimension = (2, 2)
+    t = openmc.Tally()
+    t.filters = [openmc.ParticleFilter(bins='neutron'), openmc.MeshFilter(mesh=msh)]
+    t.scores = ['total']
+    t.estimator = 'collision'
+    model.tallies.append(t)
+
+    # If photon transport is required, add a mesh tally for the photon total reaction rate.
+    if run_photon:
+        t = openmc.Tally()
+        t.filters = [openmc.ParticleFilter(bins='photon'), openmc.MeshFilter(mesh=msh)]
+        t.scores = ['total']
+        t.estimator = 'collision'
+        model.tallies.append(t)
+
+    # Set some simulation settings.
+    model.settings.batches = 10
+    model.settings.inactive = 5
+    model.settings.particles = 1000
+    model.settings.delta_tracking = True
+    model.settings.photon_transport = run_photon
+
+    return model
+
 def pwr_pin_cell() -> openmc.Model:
     """Create a PWR pin-cell model.
 
@@ -968,7 +1060,7 @@ def random_ray_lattice(second_temp = False) -> openmc.Model:
         Whether or not the cross sections should contain two temperature datapoints.
         The first data point is the C5G7 cross sections, which corresponds to a temperature
         of 294 K. The second data point is the C5G7 cross sections multiplied by 1/2,
-        which corresponds to a temperature of 3934 K. This temperature dependence is
+        which corresponds to a temperature of 394 K. This temperature dependence is
         fictitious; it is used for testing temperature feedback in the random ray solver.
 
     Returns
@@ -1314,17 +1406,17 @@ def random_ray_three_region_cube() -> openmc.Model:
 def random_ray_three_region_cube_with_detectors() -> openmc.Model:
     """Create a three region cube model with two external tally regions.
 
-    This is an adaptation of the simple monoenergetic problem of a cube with 
-    three concentric cubic regions. The innermost region is near void (with 
-    Sigma_t around 10^-5) and contains an external isotropic source term, the 
-    middle region is a mild scatterer (with Sigma_t around 10^-3), and the 
-    outer region of the cube is a scatterer and absorber (with Sigma_t around 
+    This is an adaptation of the simple monoenergetic problem of a cube with
+    three concentric cubic regions. The innermost region is near void (with
+    Sigma_t around 10^-5) and contains an external isotropic source term, the
+    middle region is a mild scatterer (with Sigma_t around 10^-3), and the
+    outer region of the cube is a scatterer and absorber (with Sigma_t around
     1).
 
-    Two cubic "detector" regions are found outside this geometry, one along the 
-    y-axis near z=0, and the other in the upper right corner of the system. 
-    The size of each detector is scaled to be equal to that of the source 
-    region. The model returned by this function contains cell tallies on each 
+    Two cubic "detector" regions are found outside this geometry, one along the
+    y-axis near z=0, and the other in the upper right corner of the system.
+    The size of each detector is scaled to be equal to that of the source
+    region. The model returned by this function contains cell tallies on each
     detector.
 
     Returns
@@ -1498,29 +1590,29 @@ def random_ray_three_region_cube_with_detectors() -> openmc.Model:
         fill=absorber_mat,
         region=detector2_region
     )
-    
+
     external_x = (
-        +x_high & +y_low & +z_low & -x_outer & 
+        +x_high & +y_low & +z_low & -x_outer &
         ((-y_outer & -z_high) | (-y_high & +z_high & -z_outer))
     )
     external_y = (
-        +y_high & -y_outer & 
+        +y_high & -y_outer &
         (
-            (+detector1_right & -x_high & +z_low & -z_outer) | 
-            (-detector1_right & +x_low & +detector1_top & -z_outer) | 
+            (+detector1_right & -x_high & +z_low & -z_outer) |
+            (-detector1_right & +x_low & +detector1_top & -z_outer) |
             (+x_high & -x_outer & +z_low & -z_high)
         )
     )
     external_z = (
-        +x_low & +y_low & +z_high & -z_outer & 
+        +x_low & +y_low & +z_high & -z_outer &
         ((-y_outer & -x_high) | (-y_high & +x_high & -x_outer))
     )
-    external_cell = openmc.Cell(fill=cavity_mat, 
-                                region=(external_x | external_y | external_z), 
+    external_cell = openmc.Cell(fill=cavity_mat,
+                                region=(external_x | external_y | external_z),
                                 name='outside cube')
 
     root = openmc.Universe(
-        name='root universe', 
+        name='root universe',
         cells=[cube_domain, detector1, detector2, external_cell]
     )
 
@@ -1604,8 +1696,8 @@ def random_ray_three_region_cube_with_detectors() -> openmc.Model:
     source_tally.estimator = estimator
 
     # Instantiate a Tallies collection and export to XML
-    tallies = openmc.Tallies([detector1_tally, 
-                              detector2_tally, 
+    tallies = openmc.Tallies([detector1_tally,
+                              detector2_tally,
                               absorber_tally,
                               cavity_tally,
                               source_tally])

openmc/settings.py

@@ -94,6 +94,10 @@ class Settings:
         release of delayed photons.
 
         .. versionadded:: 0.12
+    delta_tracking : bool
+        Whether transport should be performed with delta tracking or not.
+
+        .. versionadded:: 0.15.4
     electron_treatment : {'led', 'ttb'}
         Whether to deposit all energy from electrons locally ('led') or create
         secondary bremsstrahlung photons ('ttb').
@@ -241,7 +245,7 @@ class Settings:
             also tends to dampen the convergence rate of the solver, thus requiring
             more iterations to converge.
         :adjoint_source:
-            Source object used to define localized adjoint source/detector response 
+            Source object used to define localized adjoint source/detector response
             function.
 
         .. versionadded:: 0.15.0
@@ -482,7 +486,7 @@ class Settings:
         self._delayed_photon_scaling = None
         self._material_cell_offsets = None
         self._log_grid_bins = None
-
+        self._delta_tracking = None
         self._event_based = None
         self._max_particles_in_flight = None
         self._max_particle_events = None
@@ -1240,6 +1244,15 @@ class Settings:
         cv.check_type('delayed photon scaling', value, bool)
         self._delayed_photon_scaling = value
 
+    @property
+    def delta_tracking(self):
+        return self._delta_tracking
+
+    @delta_tracking.setter
+    def delta_tracking(self, value):
+        cv.check_type('delta_tracking', value, bool)
+        self._delta_tracking = value
+
     @property
     def material_cell_offsets(self) -> bool:
         return self._material_cell_offsets
@@ -2031,11 +2044,11 @@ class Settings:
                         path = f"./mesh[@id='{mesh.id}']"
                         if root.find(path) is None:
                             root.append(mesh.to_xml_element())
-                            if mesh_memo is not None:    
+                            if mesh_memo is not None:
                                 mesh_memo.add(mesh.id)
                 elif key == 'adjoint_source':
                     subelement = ET.SubElement(element, 'adjoint_source')
-                    # Check that all entries are valid SourceBase instances, in case 
+                    # Check that all entries are valid SourceBase instances, in case
                     # the random_ray setter was not used to populate dict entries.
                     if not isinstance(value, MutableSequence):
                         value = [value]
@@ -2062,6 +2075,11 @@ class Settings:
             element = ET.SubElement(root, "free_gas_threshold")
             element.text = str(self._free_gas_threshold)
 
+    def _create_delta_tracking_subelement(self, root):
+        if self._delta_tracking:
+            elem = ET.SubElement(root, "delta_tracking")
+            elem.text = str(self._delta_tracking).lower()
+
     def _eigenvalue_from_xml_element(self, root):
         elem = root.find('eigenvalue')
         if elem is not None:
@@ -2560,6 +2578,11 @@ class Settings:
         if text is not None:
             self.free_gas_threshold = float(text)
 
+    def _delta_tracking_from_xml_element(self, root):
+        text = get_text(root, 'delta_tracking')
+        if text is not None:
+            self.delta_tracking = text in ('true', '1')
+
     def to_xml_element(self, mesh_memo=None):
         """Create a 'settings' element to be written to an XML file.
 
@@ -2637,6 +2660,7 @@ class Settings:
         self._create_use_decay_photons_subelement(element)
         self._create_source_rejection_fraction_subelement(element)
         self._create_free_gas_threshold_subelement(element)
+        self._create_delta_tracking_subelement(element)
 
         # Clean the indentation in the file to be user-readable
         clean_indentation(element)
@@ -2755,6 +2779,7 @@ class Settings:
         settings._use_decay_photons_from_xml_element(elem)
         settings._source_rejection_fraction_from_xml_element(elem)
         settings._free_gas_threshold_from_xml_element(elem)
+        settings._delta_tracking_from_xml_element(elem)
 
         return settings
 

openmc/statepoint.py

@@ -97,6 +97,8 @@ class StatePoint:
         Working directory for simulation
     photon_transport : bool
         Indicate whether photon transport is active
+    delta_tracking : bool
+        Indicate whether delta tracking is active
     run_mode : str
         Simulation run mode, e.g. 'eigenvalue'
     runtime : dict
@@ -350,6 +352,10 @@ class StatePoint:
     def photon_transport(self):
         return self._f.attrs['photon_transport'] > 0
 
+    @property
+    def delta_tracking(self):
+        return self._f.attrs['delta_tracking'] > 0
+
     @property
     def run_mode(self):
         return self._f['run_mode'][()].decode()

src/boundary_condition.cpp

@@ -138,6 +138,13 @@ void TranslationalPeriodicBC::handle_particle(
   auto new_r = p.r() + translation_;
   int new_surface = p.surface() > 0 ? j_surf_ + 1 : -(j_surf_ + 1);
 
+  // Need to nudge the particle back into the domain when using delta tracking
+  // as event_delta_advance() does not go right up to the surface to fix
+  // tunneling.
+  if (p.delta_tracking()) {
+    new_r += FP_REL_PRECISION * p.u();
+  }
+
   // Handle the effects of the surface albedo on the particle's weight.
   BoundaryCondition::handle_albedo(p, surf);
 

src/eigenvalue.cpp

@@ -49,9 +49,15 @@ void calculate_generation_keff()
   const auto& gt = simulation::global_tallies;
 
   // Get keff for this generation by subtracting off the starting value
-  simulation::keff_generation =
-    gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) -
-    simulation::keff_generation;
+  if (settings::delta_tracking) {
+    simulation::keff_generation =
+      gt(GlobalTally::K_COLLISION, TallyResult::VALUE) -
+      simulation::keff_generation;
+  } else {
+    simulation::keff_generation =
+      gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) -
+      simulation::keff_generation;
+  }
 
   double keff_reduced;
 #ifdef OPENMC_MPI
@@ -478,30 +484,39 @@ int openmc_get_keff(double* k_combined)
   // exceptions and an expression specifically derived for the combination of
   // two estimators (vice three) should be used instead.
 
-  // First we will identify if there are any matching estimators
+  // If delta tracking is enabled, use the collision and absorption estimators
+  // only. Otherwise, we will identify if there are any matching estimators. If
+  // none match, all three estimates are used.
   int i, j;
   bool use_three = false;
-  if ((std::abs(kv[0] - kv[1]) / kv[0] < FP_REL_PRECISION) &&
-      (std::abs(cov(0, 0) - cov(1, 1)) / cov(0, 0) < FP_REL_PRECISION)) {
-    // 0 and 1 match, so only use 0 and 2 in our comparisons
-    i = 0;
-    j = 2;
-
-  } else if ((std::abs(kv[0] - kv[2]) / kv[0] < FP_REL_PRECISION) &&
-             (std::abs(cov(0, 0) - cov(2, 2)) / cov(0, 0) < FP_REL_PRECISION)) {
-    // 0 and 2 match, so only use 0 and 1 in our comparisons
+  if (settings::delta_tracking) {
     i = 0;
     j = 1;
-
-  } else if ((std::abs(kv[1] - kv[2]) / kv[1] < FP_REL_PRECISION) &&
-             (std::abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) < FP_REL_PRECISION)) {
-    // 1 and 2 match, so only use 0 and 1 in our comparisons
-    i = 0;
-    j = 1;
-
   } else {
-    // No two estimators match, so set boolean to use all three estimators.
-    use_three = true;
+    if ((std::abs(kv[0] - kv[1]) / kv[0] < FP_REL_PRECISION) &&
+        (std::abs(cov(0, 0) - cov(1, 1)) / cov(0, 0) < FP_REL_PRECISION)) {
+      // 0 and 1 match, so only use 0 and 2 in our comparisons
+      i = 0;
+      j = 2;
+
+    } else if ((std::abs(kv[0] - kv[2]) / kv[0] < FP_REL_PRECISION) &&
+               (std::abs(cov(0, 0) - cov(2, 2)) / cov(0, 0) <
+                 FP_REL_PRECISION)) {
+      // 0 and 2 match, so only use 0 and 1 in our comparisons
+      i = 0;
+      j = 1;
+
+    } else if ((std::abs(kv[1] - kv[2]) / kv[1] < FP_REL_PRECISION) &&
+               (std::abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) <
+                 FP_REL_PRECISION)) {
+      // 1 and 2 match, so only use 0 and 1 in our comparisons
+      i = 0;
+      j = 1;
+
+    } else {
+      // No two estimators match, so set boolean to use all three estimators.
+      use_three = true;
+    }
   }
 
   if (use_three) {
@@ -569,16 +584,17 @@ int openmc_get_keff(double* k_combined)
     // Urbatsch, but are simpler than for the three estimators case since the
     // block matrices of the three estimator equations reduces to scalars here
 
-    // Store the commonly used term
-    double f = kv[i] - kv[j];
-    double g = cov(i, i) + cov(j, j) - 2.0 * cov(i, j);
+    // Store common terms.
+    const double f = cov(i, i) + cov(j, j) - 2.0 * cov(i, j);
+    const double g = cov(i, i) - cov(i, j);
+    const double h = kv[i] - kv[j];
 
-    // Calculate combined estimate of k-effective
-    k_combined[0] = kv[i] - (cov(i, i) - cov(i, j)) / g * f;
+    // Eq. 36 in LA-12658-MS
+    k_combined[0] = kv[i] - g * h / f;
 
-    // Calculate standard deviation of combined estimate
-    k_combined[1] = (cov(i, i) * cov(j, j) - cov(i, j) * cov(i, j)) *
-                    (g + n * f * f) / (n * (n - 2) * g * g);
+    // This is re-derived from the \Sigma matrix based on Eq. 40 in LA-12658-MS
+    k_combined[1] =
+      (cov(i, i) - g * g / f) * ((n - 1.0) / n + h * h / f) / (n - 2.0);
     k_combined[1] = std::sqrt(k_combined[1]);
   }
   return 0;

src/event.cpp

@@ -62,6 +62,21 @@ void dispatch_xs_event(int64_t buffer_idx)
   }
 }
 
+void process_death_events(int64_t n_particles)
+{
+  simulation::time_event_death.start();
+#pragma omp parallel for schedule(runtime)
+  for (int64_t i = 0; i < n_particles; i++) {
+    Particle& p = simulation::particles[i];
+    p.event_death();
+  }
+  simulation::time_event_death.stop();
+}
+
+//==============================================================================
+// Functions for surface tracking
+//==============================================================================
+
 void process_init_events(int64_t n_particles, int64_t source_offset)
 {
   simulation::time_event_init.start();
@@ -169,17 +184,6 @@ void process_collision_events()
   simulation::time_event_collision.stop();
 }
 
-void process_death_events(int64_t n_particles)
-{
-  simulation::time_event_death.start();
-#pragma omp parallel for schedule(runtime)
-  for (int64_t i = 0; i < n_particles; i++) {
-    Particle& p = simulation::particles[i];
-    p.event_death();
-  }
-  simulation::time_event_death.stop();
-}
-
 void process_transport_events()
 {
   while (true) {
@@ -221,4 +225,189 @@ void process_init_secondary_events(int64_t n_particles, int64_t offset,
   simulation::time_event_init.stop();
 }
 
+//==============================================================================
+// Functions for delta tracking
+//==============================================================================
+
+void process_delta_init_events(int64_t n_particles, int64_t source_offset)
+{
+  simulation::time_event_init.start();
+#pragma omp parallel for schedule(runtime)
+  for (int64_t i = 0; i < n_particles; i++) {
+    simulation::particles[i].delta_tracking() = true;
+    initialize_particle_track(
+      simulation::particles[i], source_offset + i + 1, false);
+
+    simulation::advance_particle_queue[i] = {simulation::particles[i], i};
+  }
+  simulation::advance_particle_queue.resize(n_particles);
+  simulation::time_event_init.stop();
+}
+
+void process_delta_calculate_xs_events(SharedArray<EventQueueItem>& queue)
+{
+  simulation::time_event_calculate_xs.start();
+
+  // TODO: If using C++17, we could perform a parallel sort of the queue by
+  // particle type, material type, and then energy, in order to improve cache
+  // locality and reduce thread divergence on GPU. However, the parallel
+  // algorithms typically require linking against an additional library (Intel
+  // TBB). Prior to C++17, std::sort is a serial only operation, which in this
+  // case makes it too slow to be practical for most test problems.
+  //
+  // std::sort(std::execution::par_unseq, queue.data(), queue.data() +
+  // queue.size());
+
+  int64_t offset = simulation::collision_queue.size();
+
+  simulation::collision_queue.resize(offset + queue.size());
+
+#pragma omp parallel for schedule(runtime)
+  for (int64_t i = 0; i < queue.size(); i++) {
+    Particle* p = &simulation::particles[queue[i].idx];
+    p->event_calculate_xs();
+
+    // After executing a calculate_xs event in delta tracking, particles will
+    // always require a collision event. Therefore, we don't need to use
+    // the protected enqueuing function.
+    simulation::collision_queue[offset + i] = queue[i];
+  }
+
+  queue.resize(0);
+
+  simulation::time_event_calculate_xs.stop();
+}
+
+void process_delta_advance_particle_events()
+{
+  simulation::time_event_advance_particle.start();
+
+#pragma omp parallel for schedule(runtime)
+  for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {
+    int64_t buffer_idx = simulation::advance_particle_queue[i].idx;
+    Particle& p = simulation::particles[buffer_idx];
+    p.event_delta_advance();
+    if (!p.alive())
+      continue;
+
+    if (p.type() == ParticleType::electron() ||
+        p.type() == ParticleType::positron()) {
+      // Electrons / positrons collide in place and don't require cross section
+      // calculations. Can append to the collision queue directly.
+      simulation::collision_queue.thread_safe_append({p, buffer_idx});
+      continue;
+    }
+
+    if (p.collision_distance() < p.boundary().distance()) {
+      // We need to compute cross sections prior to processing a collision.
+      dispatch_xs_event(buffer_idx);
+    } else {
+      simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx});
+    }
+  }
+
+  simulation::advance_particle_queue.resize(0);
+
+  simulation::time_event_advance_particle.stop();
+}
+
+void process_delta_surface_crossing_events()
+{
+  simulation::time_event_surface_crossing.start();
+
+#pragma omp parallel for schedule(runtime)
+  for (int64_t i = 0; i < simulation::surface_crossing_queue.size(); i++) {
+    int64_t buffer_idx = simulation::surface_crossing_queue[i].idx;
+    Particle& p = simulation::particles[buffer_idx];
+    p.event_cross_surface();
+    p.event_check_limit_and_revive();
+    if (p.alive())
+      simulation::advance_particle_queue.thread_safe_append({p, buffer_idx});
+  }
+
+  simulation::surface_crossing_queue.resize(0);
+
+  simulation::time_event_surface_crossing.stop();
+}
+
+void process_delta_collision_events()
+{
+  simulation::time_event_collision.start();
+
+#pragma omp parallel for schedule(runtime)
+  for (int64_t i = 0; i < simulation::collision_queue.size(); i++) {
+    int64_t buffer_idx = simulation::collision_queue[i].idx;
+    Particle& p = simulation::particles[buffer_idx];
+
+    if (p.type() == ParticleType::electron() ||
+        p.type() == ParticleType::positron()) {
+      p.event_collide();
+    } else {
+      if (p.kill_invalid_maj()) {
+        continue;
+      }
+
+      if (prn(p.current_seed()) < (p.macro_xs().total / p.majorant())) {
+        // Real collision, need to process the collision prior to enqueuing an
+        // advance event.
+        p.event_collide();
+      }
+    }
+
+    p.event_check_limit_and_revive();
+    if (p.alive()) {
+      simulation::advance_particle_queue.thread_safe_append({p, buffer_idx});
+    }
+  }
+
+  simulation::collision_queue.resize(0);
+
+  simulation::time_event_collision.stop();
+}
+
+void process_delta_transport_events()
+{
+  while (true) {
+    int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
+      simulation::calculate_nonfuel_xs_queue.size(),
+      simulation::advance_particle_queue.size(),
+      simulation::surface_crossing_queue.size(),
+      simulation::collision_queue.size()});
+
+    if (max == 0) {
+      break;
+    } else if (max == simulation::calculate_fuel_xs_queue.size()) {
+      process_delta_calculate_xs_events(simulation::calculate_fuel_xs_queue);
+    } else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
+      process_delta_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
+    } else if (max == simulation::advance_particle_queue.size()) {
+      process_delta_advance_particle_events();
+    } else if (max == simulation::surface_crossing_queue.size()) {
+      process_delta_surface_crossing_events();
+    } else if (max == simulation::collision_queue.size()) {
+      process_delta_collision_events();
+    }
+  }
+}
+
+void process_delta_init_secondary_events(int64_t n_particles, int64_t offset,
+  const SharedArray<SourceSite>& shared_secondary_bank)
+{
+  simulation::time_event_init.start();
+#pragma omp parallel for schedule(runtime)
+  for (int64_t i = 0; i < n_particles; i++) {
+    simulation::particles[i].delta_tracking() = true;
+    initialize_particle_track(simulation::particles[i], offset + i + 1, true);
+    const SourceSite& site = shared_secondary_bank[offset + i];
+    simulation::particles[i].event_revive_from_secondary(site);
+
+    if (simulation::particles[i].alive()) {
+      simulation::particles[i].event_calculate_xs();
+      simulation::advance_particle_queue.thread_safe_append(
+        {simulation::particles[i], i});
+    }
+  }
+  simulation::time_event_init.stop();
+}
+
 } // namespace openmc

src/geometry.cpp

@@ -24,6 +24,8 @@ namespace model {
 int root_universe {-1};
 int n_coord_levels;
 
+vector<int> boundary_surfaces;
+
 vector<int64_t> overlap_check_count;
 
 vector<OverlapKey> overlap_keys;
@@ -477,6 +479,30 @@ BoundaryInfo distance_to_boundary(GeometryState& p)
   return info;
 }
 
+//==============================================================================
+
+BoundaryInfo distance_to_external_boundary(GeometryState& p)
+{
+  BoundaryInfo info;
+
+  info.distance() = INFTY;
+  info.surface() = 0;
+  info.coord_level() = 1;
+  for (auto s_idx : model::boundary_surfaces) {
+    const auto& s = model::surfaces[s_idx];
+    double surf_dist = s->distance(p.r(), p.u(), false);
+    if (surf_dist < info.distance()) {
+      info.distance() = surf_dist;
+      info.surface() = s_idx + 1;
+      if (s->sense(p.r(), p.u())) {
+        info.surface() *= -1;
+      }
+    }
+  }
+
+  return info;
+}
+
 //==============================================================================
 // C API
 //==============================================================================

src/geometry_aux.cpp

@@ -270,6 +270,20 @@ void get_temperatures(
 
 //==============================================================================
 
+void detect_boundary_surfaces()
+{
+  for (int i = 0; i < model::surfaces.size(); i++) {
+    // if the surface has a non-transmission boundary condition,
+    // add it to the list of surfaces to track during delta tracking
+    const auto& s = model::surfaces[i];
+    if (s->bc_) {
+      model::boundary_surfaces.push_back(i);
+    }
+  }
+}
+
+//==============================================================================
+
 void finalize_geometry()
 {
   // Perform some final operations to set up the geometry
@@ -280,6 +294,10 @@ void finalize_geometry()
   // Assign temperatures to cells that don't have temperatures already assigned
   assign_temperatures();
 
+  // Find all boundary surfaces. Used in delta tracking to trace through the
+  // geometry.
+  detect_boundary_surfaces();
+
   // Determine number of nested coordinate levels in the geometry
   model::n_coord_levels = maximum_levels(model::root_universe);
 }

src/majorant.cpp

@@ -0,0 +1,601 @@
+#include <fmt/core.h>
+
+#include "openmc/capi.h"
+#include "openmc/constants.h"
+#include "openmc/geometry.h"
+#include "openmc/interpolate.h"
+#include "openmc/majorant.h"
+#include "openmc/material.h"
+#include "openmc/nuclide.h"
+#include "openmc/photon.h"
+#include "openmc/search.h"
+#include "openmc/settings.h"
+#include "openmc/simulation.h"
+#include "openmc/thermal.h"
+#include "openmc/timer.h"
+#include "openmc/universe.h"
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace data {
+std::unique_ptr<NeutronMajorant> n_majorant;
+std::unique_ptr<PhotonMajorant> p_majorant;
+
+} // namespace data
+
+//==============================================================================
+// Majorant implementation
+//==============================================================================
+
+Majorant::Majorant(int i_universe) : maj_universe_(i_universe)
+{
+  if (maj_universe_ == C_NONE || maj_universe_ >= model::universes.size()) {
+    fatal_error(fmt::format("Invalid majorant universe: {}", maj_universe_));
+  }
+
+  // First, find unique cells contained in this universe.
+  std::unordered_set<int> unique_mat_cells;
+
+  const auto& maj_uni = model::universes[maj_universe_];
+  for (int i_cell : maj_uni->cells_) {
+    const auto& uni_cell = model::cells[i_cell];
+
+    // If the cell is filled with a material, it won't have any sub-cells.
+    if (uni_cell->type_ == Fill::MATERIAL) {
+      if (unique_mat_cells.count(i_cell) == 0) {
+        unique_mat_cells.emplace(i_cell);
+      }
+    } else {
+      // This cell is filled with a universe or lattice. Need to get the list of
+      // cells and cell instances.
+      const auto contained_cells = uni_cell->get_contained_cells();
+      for (const auto& [i_con_cell, contained_instances] : contained_cells) {
+        if (unique_mat_cells.count(i_con_cell) == 0) {
+          unique_mat_cells.emplace(i_con_cell);
+        }
+      }
+    }
+  }
+
+  // Next, find all materials contained in the majorant's universe. This also
+  // obtains the maximum density multiplier applied to that material.
+  std::unordered_set<int> unique_materials;
+  for (int i_cell : unique_mat_cells) {
+    auto& cell = model::cells[i_cell];
+
+    for (int instance = 0; instance < cell->n_instances(); ++instance) {
+      int i_material = cell->material(instance);
+      // Skip over void materials.
+      if (i_material == MATERIAL_VOID) {
+        continue;
+      }
+
+      // Check to see if we've found the contained material yet. If not, add
+      // to the set of materials discovered and add to the map of density
+      // multipliers.
+      if (unique_materials.count(i_material) == 0) {
+        unique_materials.emplace(i_material);
+        max_density_mult_[i_material] = cell->density_mult(instance);
+      } else {
+        // We've found this material already. Need to take the maximum density
+        // multiplier.
+        max_density_mult_.at(i_material) = std::max(
+          max_density_mult_.at(i_material), cell->density_mult(instance));
+      }
+    }
+  }
+
+  // Insert the elements from the set.
+  contained_materials_.assign(unique_materials.begin(), unique_materials.end());
+}
+
+void Majorant::compute_majorant()
+{
+  // Fill with zeros.
+  xs_.resize(grid_.energy.size(), 0.0);
+
+  std::vector<double> material_maj_xs;
+  for (int i_material : contained_materials_) {
+    // Populate the per-material majorant cross section. We pass in
+    // 'material_maj_xs' instead of returning a vector with the per-material
+    // majorant every time to avoid costly reallocations and copy operations,
+    // which have a fairly large impact on the time it takes to build the
+    // majorant. The function 'fill_material_maj_xs(...)' is responsible for
+    // resizing 'material_maj_xs' and populating each value at a given energy
+    // grid point.
+    fill_material_maj_xs(i_material, max_density_mult_.at(i_material),
+      grid_.energy, material_maj_xs);
+
+    // Compute the full majorant by taking the max over each material cross
+    // section.
+    for (int i_energy = 0; i_energy < xs_.size(); ++i_energy) {
+      xs_[i_energy] = std::max(xs_[i_energy], material_maj_xs[i_energy]);
+    }
+  }
+}
+
+void Majorant::post_process_grid(Nuclide::EnergyGrid& grid) const
+{
+  // Fetch the minimum and maximum transport energies from the superclass.
+  const double E_min = min_transport_energy();
+  const double E_max = max_transport_energy();
+
+  std::sort(grid.energy.begin(), grid.energy.end());
+  auto unique_end = std::unique(grid.energy.begin(), grid.energy.end());
+  grid.energy.resize(std::distance(grid.energy.begin(), unique_end));
+
+  // Remove all values below the minimum neutron energy.
+  auto min_it = grid.energy.begin();
+  while (*min_it < E_min) {
+    min_it++;
+  }
+  grid.energy.erase(grid.energy.begin(), min_it + 1);
+
+  // Insert the minimum neutron energy at the beginning.
+  grid.energy.insert(grid.energy.begin(), E_min);
+
+  // Remove all values above the maximum neutron energy.
+  auto max_it = --grid.energy.end();
+  while (*max_it > E_max) {
+    max_it--;
+  }
+  grid.energy.erase(max_it - 1, grid.energy.end());
+
+  // Insert the maximum neutron energy at the end.
+  grid.energy.insert(grid.energy.end(), E_max);
+}
+
+//==============================================================================
+// NeutronMajorant implementation
+//==============================================================================
+
+NeutronMajorant::NeutronMajorant(int i_universe) : Majorant(i_universe) {}
+
+double NeutronMajorant::calculate_neutron_xs(double energy) const
+{
+  const int i_grid = get_i_grid(energy, grid_);
+  return interpolate_lin_lin(grid_.energy[i_grid], grid_.energy[i_grid + 1],
+    xs_[i_grid], xs_[i_grid + 1], energy);
+}
+
+void NeutronMajorant::compute_unionized_grid()
+{
+  // This function generates a unionized cross section grid between smooth cross
+  // sections and URR probability table grids.
+  std::unordered_set<int> processed_nuclides;
+  for (int i_mat : contained_materials_) {
+    const auto& mat = model::materials[i_mat];
+    for (auto i_nuclide : mat->nuclide_) {
+      // Only unionize nuclides we haven't checked yet.
+      if (processed_nuclides.count(i_nuclide) > 0) {
+        continue;
+      }
+
+      const auto& nuclide = data::nuclides[i_nuclide];
+      // ======================================================================
+      // Unionizing the URR energy grids.
+      if (nuclide->urr_present_ && settings::urr_ptables_on) {
+        for (const auto& nuc_urr : nuclide->urr_data_) {
+          grid_.energy.insert(
+            grid_.energy.end(), nuc_urr.energy_.begin(), nuc_urr.energy_.end());
+        }
+      }
+
+      // ======================================================================
+      // Unionize the smooth cross section energy grids.
+      for (const auto& nuc_grid : nuclide->grid_) {
+        grid_.energy.insert(
+          grid_.energy.end(), nuc_grid.energy.begin(), nuc_grid.energy.end());
+      }
+
+      processed_nuclides.insert(i_nuclide);
+    }
+  }
+
+  // Post-process the energy grid now that all points from nuclides are
+  // included. This sorts the energy points, removes duplicates, and removes all
+  // energies exceeding neutron transport bounds.
+  post_process_grid(grid_);
+
+  // Initialize the grid for fast lookups. This only applies to neutrons.
+  grid_.init();
+}
+
+void NeutronMajorant::fill_material_maj_xs(int i_material,
+  double max_density_mult, const std::vector<double>& to_grid,
+  std::vector<double>& mat_maj) const
+{
+  const auto& mat = *model::materials[i_material];
+
+  mat_maj.resize(to_grid.size());
+
+#pragma omp parallel for
+  for (int i_energy = 0; i_energy < to_grid.size(); ++i_energy) {
+    mat_maj[i_energy] = 0.0;
+    const double union_energy = to_grid[i_energy];
+
+    int mat_sab_table_idx = 0;
+    bool check_sab = (mat.thermal_tables_.size() > 0);
+
+    for (int i = 0; i < mat.nuclide_.size(); ++i) {
+      // ======================================================================
+      // CHECK FOR S(A,B) TABLE
+      int i_sab = C_NONE;
+      double sab_frac = 0.0;
+
+      // Check if this nuclide matches one of the S(a,b) tables specified.
+      // This relies on thermal_tables_ being sorted by .index_nuclide
+      if (check_sab) {
+        const auto& sab {mat.thermal_tables_[mat_sab_table_idx]};
+        if (i == sab.index_nuclide) {
+          // Get index in sab_tables
+          i_sab = sab.index_table;
+          sab_frac = sab.fraction;
+
+          // If particle energy is greater than the highest energy for the
+          // S(a,b) table, then don't use the S(a,b) table
+          if (union_energy > data::thermal_scatt[i_sab]->energy_max_) {
+            i_sab = C_NONE;
+          }
+
+          // Increment position in thermal_tables_
+          ++mat_sab_table_idx;
+
+          // Don't check for S(a,b) tables if there are no more left
+          if (mat_sab_table_idx == mat.thermal_tables_.size()) {
+            check_sab = false;
+          }
+        }
+      }
+
+      // ======================================================================
+      // Compute the maximum smooth total cross section. This is either the
+      // free gas cross section at energies larger than the Bragg edge, or
+      // the bound cross section in the thermal scattering region.
+      double micro_smooth_tot_xs = 0.0;
+      if (i_sab >= 0) {
+        // Thermal scattering cross sections using S(a,b) tables.
+        micro_smooth_tot_xs = calculate_max_sab_tot_xs(
+          mat.nuclide_[i], i_sab, sab_frac, union_energy);
+      } else {
+        // Free gas smooth cross section
+        micro_smooth_tot_xs =
+          calculate_max_smooth_xs(mat.nuclide_[i], union_energy);
+      }
+
+      // ======================================================================
+      // Compute the URR cross section. This shouldn't intersect with the
+      // S(a,b) cross section.
+      double micro_urr_xs = calculate_max_urr_xs(
+        mat.nuclide_[i], union_energy, micro_smooth_tot_xs);
+
+      // ======================================================================
+      // Accumulate the macroscopic cross section.
+      mat_maj[i_energy] += std::max(micro_smooth_tot_xs, micro_urr_xs) *
+                           mat.atom_density(i, max_density_mult);
+    }
+  }
+}
+
+double NeutronMajorant::calculate_max_smooth_xs(
+  int i_nuclide, double energy) const
+{
+  const auto& nuc = *data::nuclides[i_nuclide];
+
+  double max_smooth_tot_xs = 0.0;
+  for (int i_temp = 0; i_temp < nuc.kTs_.size(); ++i_temp) {
+    const auto& nuc_grid = nuc.grid_[i_temp];
+    int i_grid = get_i_grid(energy, nuc_grid);
+    auto total = nuc.xs_[i_temp].slice(openmc::tensor::all, 0);
+    double xs = interpolate_lin_lin(nuc_grid.energy[i_grid],
+      nuc_grid.energy[i_grid + 1], total[i_grid], total[i_grid + 1], energy);
+    max_smooth_tot_xs = std::max(max_smooth_tot_xs, xs);
+  }
+
+  return max_smooth_tot_xs;
+}
+
+double NeutronMajorant::calculate_max_urr_xs(
+  int i_nuclide, double energy, double smooth_xs) const
+{
+  // A tolerance on the URR check to make sure we include the URR energy grid
+  // bounds.
+  constexpr double URR_FUZZY_CHECK = 1e-6;
+
+  const auto& nuc = *data::nuclides[i_nuclide];
+  if (!nuc.urr_present_) {
+    return 0.0;
+  }
+
+  double max_urr_xs = 0.0;
+  for (const auto& urr : nuc.urr_data_) {
+    if (!(urr.energy_in_bounds(energy - URR_FUZZY_CHECK) ||
+          urr.energy_in_bounds(energy + URR_FUZZY_CHECK))) {
+      continue;
+    }
+
+    int i_energy;
+    if (energy <= urr.energy_.front()) {
+      i_energy = 0;
+    } else if (energy >= urr.energy_.back()) {
+      i_energy = urr.energy_.size() - 2;
+    } else {
+      i_energy =
+        lower_bound_index(&urr.energy_.front(), &urr.energy_.back(), energy);
+    }
+
+    // Find the maximum URR cross sections for the two bounding energy points.
+    double max_urr_xs_E0 = 0.0;
+    double max_urr_xs_E1 = 0.0;
+    for (int i_cdf = 0; i_cdf < urr.n_cdf(); ++i_cdf) {
+      max_urr_xs_E0 =
+        std::max(max_urr_xs_E0, urr.xs_values_(i_energy, i_cdf).total);
+      max_urr_xs_E1 =
+        std::max(max_urr_xs_E1, urr.xs_values_(i_energy + 1, i_cdf).total);
+    }
+    // Handle the rare case where the points could be negative.
+    max_urr_xs_E0 = std::max(max_urr_xs_E0, 0.0);
+    max_urr_xs_E1 = std::max(max_urr_xs_E1, 0.0);
+
+    // Interpolate the bounding energy points.
+    double interp_urr_xs = 0.0;
+    if (urr.interp_ == Interpolation::lin_lin) {
+      interp_urr_xs = interpolate_lin_lin(urr.energy_[i_energy],
+        urr.energy_[i_energy + 1], max_urr_xs_E0, max_urr_xs_E1, energy);
+    } else if (urr.interp_ == Interpolation::log_log) {
+      interp_urr_xs = interpolate_log_log(urr.energy_[i_energy],
+        urr.energy_[i_energy + 1], max_urr_xs_E0, max_urr_xs_E1, energy);
+    }
+
+    // Multiply by the smooth cross section (after interpolation) if required.
+    if (urr.multiply_smooth_) {
+      interp_urr_xs *= smooth_xs;
+    }
+
+    max_urr_xs = std::max({max_urr_xs, interp_urr_xs, smooth_xs});
+  }
+
+  return max_urr_xs;
+}
+
+double NeutronMajorant::calculate_max_sab_tot_xs(
+  int i_nuclide, int i_sab, double sab_frac, double energy) const
+{
+  const auto& nuc = *data::nuclides[i_nuclide];
+  const auto& thermal = *data::thermal_scatt[i_sab];
+
+  // Loop over the nuclide's temperature grid to ensure we're consistent.
+  double max_sab_total = 0.0;
+  for (int i_nuc_temp = 0; i_nuc_temp < nuc.kTs_.size(); ++i_nuc_temp) {
+    double nuc_kT = nuc.kTs_[i_nuc_temp] * nuc.kTs_[i_nuc_temp];
+
+    // Compute the elastic and inelastic scattering cross sections. The S(a,b)
+    // cross sections are interpolated to match the nuclide temperature point.
+    double thermal_elastic;
+    double thermal_inelastic;
+    const auto& tkTs = thermal.kTs_;
+    if (tkTs.size() > 1) {
+      if (nuc_kT < tkTs.front()) {
+        thermal.data_.front().calculate_xs(
+          energy, &thermal_elastic, &thermal_inelastic);
+      } else if (nuc_kT > tkTs.back()) {
+        thermal.data_.back().calculate_xs(
+          energy, &thermal_elastic, &thermal_inelastic);
+      } else {
+        // Find temperatures that bound the actual temperature
+        int i_sab_temp = 0;
+        while (
+          tkTs[i_sab_temp + 1] < nuc_kT && i_sab_temp + 1 < tkTs.size() - 1) {
+          ++i_sab_temp;
+        }
+        // Interpolate the scattering cross sections to the nuclide temperature
+        // grid point.
+        double T0_elastic, T1_elastic, T0_inelastic, T1_inelastic;
+        thermal.data_[i_sab_temp].calculate_xs(
+          energy, &T0_elastic, &T0_inelastic);
+        thermal.data_[i_sab_temp + 1].calculate_xs(
+          energy, &T1_elastic, &T1_inelastic);
+        thermal_elastic = interpolate_lin_lin(tkTs[i_sab_temp],
+          tkTs[i_sab_temp + 1], T0_elastic, T1_elastic, nuc_kT);
+        thermal_inelastic = interpolate_lin_lin(tkTs[i_sab_temp],
+          tkTs[i_sab_temp + 1], T0_inelastic, T1_inelastic, nuc_kT);
+      }
+    } else {
+      thermal.data_[0].calculate_xs(
+        energy, &thermal_elastic, &thermal_inelastic);
+    }
+
+    // Compute the free gas total and elastic cross sections interpolated on the
+    // majorant grid.
+    const auto& nuc_grid = nuc.grid_[i_nuc_temp];
+    int i_grid = get_i_grid(energy, nuc_grid);
+    const auto& free_tot = nuc.xs_[i_nuc_temp].slice(openmc::tensor::all, 0);
+    const auto& free_ela = nuc.reactions_[0]->xs_[i_nuc_temp].value;
+    double tot_xs =
+      interpolate_lin_lin(nuc_grid.energy[i_grid], nuc_grid.energy[i_grid + 1],
+        free_tot[i_grid], free_tot[i_grid + 1], energy);
+    double ela_xs =
+      interpolate_lin_lin(nuc_grid.energy[i_grid], nuc_grid.energy[i_grid + 1],
+        free_ela[i_grid], free_ela[i_grid + 1], energy);
+
+    double thermal_xs = sab_frac * (thermal_elastic + thermal_inelastic);
+    double sab_corrected_total = tot_xs + thermal_xs - sab_frac * ela_xs;
+    max_sab_total = std::max(sab_corrected_total, max_sab_total);
+  }
+
+  return max_sab_total;
+}
+
+int NeutronMajorant::get_i_grid(
+  double energy, const Nuclide::EnergyGrid& grid) const
+{
+  // Find energy index on energy grid
+  int i_log_union =
+    std::log(energy / data::energy_min[i_neutron_]) / simulation::log_spacing;
+
+  int i_grid;
+  if (energy <= grid.energy.front()) {
+    i_grid = 0;
+  } else if (energy >= grid.energy.back()) {
+    i_grid = grid.energy.size() - 2;
+  } else {
+    // Determine bounding indices based on which equal log-spaced
+    // interval the energy is in
+    int i_low = grid.grid_index[i_log_union];
+    int i_high = grid.grid_index[i_log_union + 1] + 1;
+
+    // This catches the very rare case where floating point comparisons fail.
+    if (i_low >= grid.energy.size() || i_high >= grid.energy.size()) {
+      i_grid = grid.energy.size() - 2;
+    } else {
+      // Perform binary search over reduced range
+      i_grid = i_low + lower_bound_index(
+                         &grid.energy[i_low], &grid.energy[i_high], energy);
+    }
+  }
+
+  // check for rare case where two energy points are the same
+  if (grid.energy[i_grid] == grid.energy[i_grid + 1])
+    ++i_grid;
+
+  return i_grid;
+}
+
+//==============================================================================
+// PhotonMajorant implementation
+//==============================================================================
+
+PhotonMajorant::PhotonMajorant(int i_universe) : Majorant(i_universe) {}
+
+void PhotonMajorant::compute_unionized_grid()
+{
+  // This function generates a unionized cross section grid for all elements.
+  std::unordered_set<int> processed_elements;
+  for (int i_mat : contained_materials_) {
+    const auto& mat = model::materials[i_mat];
+    for (int i = 0; i < mat->nuclide_.size(); ++i) {
+      // Only unionize elements we haven't checked yet.
+      if (processed_elements.count(mat->element_[i]) > 0) {
+        continue;
+      }
+
+      const auto& element = data::elements[mat->element_[i]];
+      grid_.energy.insert(
+        grid_.energy.end(), element->energy_.begin(), element->energy_.end());
+
+      processed_elements.insert(mat->element_[i]);
+    }
+  }
+
+  // Post-process the energy grid now that all points from photon interactions
+  // are included. This sorts the energy points, removes duplicates, and removes
+  // all energies exceeding photon transport bounds.
+  post_process_grid(grid_);
+}
+
+double PhotonMajorant::calculate_photon_xs(double energy) const
+{
+  double log_energy = std::log(energy);
+  int i_grid = get_i_grid<std::vector<double>>(log_energy, grid_.energy);
+
+  // calculate interpolation factor
+  double f = (log_energy - grid_.energy[i_grid]) /
+             (grid_.energy[i_grid + 1] - grid_.energy[i_grid]);
+
+  // interpolate the total cross section
+  return std::exp(xs_[i_grid] + f * (xs_[i_grid + 1] - xs_[i_grid]));
+}
+
+void PhotonMajorant::fill_material_maj_xs(int i_material,
+  double max_density_mult, const std::vector<double>& to_grid,
+  std::vector<double>& mat_maj) const
+{
+  const auto& mat = *model::materials[i_material];
+
+  mat_maj.resize(to_grid.size());
+
+#pragma omp parallel for
+  for (int i_energy = 0; i_energy < to_grid.size(); ++i_energy) {
+    mat_maj[i_energy] = 0.0;
+    const double union_log_energy = to_grid[i_energy];
+
+    for (int i = 0; i < mat.nuclide_.size(); ++i) {
+      const int i_element = mat.element_[i];
+
+      mat_maj[i_energy] += calculate_elem_tot_xs(i_element, union_log_energy) *
+                           mat.atom_density(i, max_density_mult);
+    }
+    mat_maj[i_energy] = std::log(mat_maj[i_energy]);
+  }
+}
+
+double PhotonMajorant::calculate_elem_tot_xs(
+  int i_element, double log_energy) const
+{
+  const auto& elem = *data::elements[i_element];
+  int i_grid = get_i_grid<tensor::Tensor<double>>(log_energy, elem.energy_);
+
+  // calculate interpolation factor
+  double f = (log_energy - elem.energy_(i_grid)) /
+             (elem.energy_(i_grid + 1) - elem.energy_(i_grid));
+
+  // Calculate microscopic coherent cross section
+  double coherent =
+    std::exp(elem.coherent_(i_grid) +
+             f * (elem.coherent_(i_grid + 1) - elem.coherent_(i_grid)));
+
+  // Calculate microscopic incoherent cross section
+  double incoherent =
+    std::exp(elem.incoherent_(i_grid) +
+             f * (elem.incoherent_(i_grid + 1) - elem.incoherent_(i_grid)));
+
+  // Calculate microscopic photoelectric cross section
+  double photoelectric = 0.0;
+  tensor::View<const double> xs_lower = elem.cross_sections_.slice(i_grid);
+  tensor::View<const double> xs_upper = elem.cross_sections_.slice(i_grid + 1);
+
+  for (int i = 0; i < xs_upper.size(); ++i)
+    if (xs_lower(i) != 0)
+      photoelectric += std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i)));
+
+  // Calculate microscopic pair production cross section
+  double pair_production =
+    std::exp(elem.pair_production_total_(i_grid) +
+             f * (elem.pair_production_total_(i_grid + 1) -
+                   elem.pair_production_total_(i_grid)));
+
+  // Calculate microscopic total cross section
+  return coherent + incoherent + photoelectric + pair_production;
+}
+
+//! Create a majorant cross section for photons or neutrons.
+void create_majorants()
+{
+  simulation::time_build_majorant.start();
+
+  write_message("Constructing a neutron majorant cross section");
+  data::n_majorant = std::make_unique<NeutronMajorant>(model::root_universe);
+  data::n_majorant->compute_unionized_grid();
+  data::n_majorant->compute_majorant();
+
+  if (settings::photon_transport) {
+    write_message("Constructing a photon majorant cross section");
+    data::p_majorant = std::make_unique<PhotonMajorant>(model::root_universe);
+    data::p_majorant->compute_unionized_grid();
+    data::p_majorant->compute_majorant();
+  }
+
+  simulation::time_build_majorant.stop();
+}
+
+//! Reset the photon and neutron majorant cross sections.
+void reset_majorants()
+{
+  openmc::data::n_majorant.reset(nullptr);
+  openmc::data::p_majorant.reset(nullptr);
+}
+} // namespace openmc

src/material.cpp

@@ -60,6 +60,10 @@ Material::Material(pugi::xml_node node)
   }
 
   if (check_for_node(node, "cfg")) {
+    if (settings::delta_tracking) {
+      fatal_error("NCrystal materials are presently not supported when running "
+                  "with delta tracking!");
+    }
     auto cfg = get_node_value(node, "cfg");
     write_message(
       5, "NCrystal config string for material #{}: '{}'", this->id(), cfg);

src/nuclide.cpp

@@ -493,7 +493,7 @@ void Nuclide::create_derived(
   }
 }
 
-void Nuclide::init_grid()
+void Nuclide::EnergyGrid::init()
 {
   int neutron = ParticleType::neutron().transport_index();
   double E_min = data::energy_min[neutron];
@@ -506,23 +506,27 @@ void Nuclide::init_grid()
   // Create equally log-spaced energy grid
   auto umesh = tensor::linspace(0.0, M * spacing, M + 1);
 
-  for (auto& grid : grid_) {
-    // Resize array for storing grid indices
-    grid.grid_index.resize(M + 1);
+  grid_index.resize(M + 1);
 
-    // Determine corresponding indices in nuclide grid to energies on
-    // equal-logarithmic grid
-    int j = 0;
-    for (int k = 0; k <= M; ++k) {
-      while (std::log(grid.energy[j + 1] / E_min) <= umesh(k)) {
-        // Ensure that for isotopes where maxval(grid.energy) << E_max that
-        // there are no out-of-bounds issues.
-        if (j + 2 == grid.energy.size())
-          break;
-        ++j;
-      }
-      grid.grid_index[k] = j;
+  // Determine corresponding indices in nuclide grid to energies on
+  // equal-logarithmic grid
+  int j = 0;
+  for (int k = 0; k <= M; ++k) {
+    while (std::log(energy[j + 1] / E_min) <= umesh(k)) {
+      // Ensure that for isotopes where maxval(grid.energy) << E_max that
+      // there are no out-of-bounds issues.
+      if (j + 2 == energy.size())
+        break;
+      ++j;
     }
+    grid_index[k] = j;
+  }
+}
+
+void Nuclide::init_grid()
+{
+  for (auto& grid : grid_) {
+    grid.init();
   }
 }
 

src/output.cpp

@@ -439,6 +439,9 @@ void print_runtime()
   // display time elapsed for various sections
   show_time("Total time for initialization", time_initialize.elapsed());
   show_time("Reading cross sections", time_read_xs.elapsed(), 1);
+  if (settings::delta_tracking) {
+    show_time("Total time building majorants", time_build_majorant.elapsed());
+  }
   show_time("Total time in simulation",
     time_inactive.elapsed() + time_active.elapsed());
   show_time("Time in transport only", time_transport.elapsed(), 1);
@@ -553,9 +556,14 @@ void print_results()
       std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_COLLISION, 0), n);
       fmt::print(" k-effective (Collision)     = {:.5f} +/- {:.5f}\n", mean,
         t_n1 * stdev);
-      std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_TRACKLENGTH, 0), n);
-      fmt::print(" k-effective (Track-length)  = {:.5f} +/- {:.5f}\n", mean,
-        t_n1 * stdev);
+      if (settings::delta_tracking) {
+        fmt::print(" k-effective (Track-length)  = (Delta-tracking enabled)\n");
+      } else {
+        std::tie(mean, stdev) =
+          mean_stdev(&gt(GlobalTally::K_TRACKLENGTH, 0), n);
+        fmt::print(" k-effective (Track-length)  = {:.5f} +/- {:.5f}\n", mean,
+          t_n1 * stdev);
+      }
       std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_ABSORPTION, 0), n);
       fmt::print(" k-effective (Absorption)    = {:.5f} +/- {:.5f}\n", mean,
         t_n1 * stdev);

src/particle.cpp

@@ -15,6 +15,7 @@
 #include "openmc/geometry.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/lattice.h"
+#include "openmc/majorant.h"
 #include "openmc/material.h"
 #include "openmc/message_passing.h"
 #include "openmc/mgxs_interface.h"
@@ -152,6 +153,7 @@ void Particle::from_source(const SourceSite* src)
   material() = C_NONE;
   n_collision() = 0;
   fission() = false;
+  majorant() = 0.0;
   zero_flux_derivs();
   lifetime() = 0.0;
 #ifdef OPENMC_DAGMC_ENABLED
@@ -191,6 +193,10 @@ void Particle::from_source(const SourceSite* src)
   wgt_born() = src->wgt_born;
   wgt_ww_born() = src->wgt_ww_born;
   n_split() = src->n_split;
+
+  if (delta_tracking()) {
+    update_majorant();
+  }
 }
 
 void Particle::event_calculate_xs()
@@ -306,7 +312,8 @@ void Particle::event_advance()
   }
 
   // Score track-length estimate of k-eff
-  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
+  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron() &&
+      !delta_tracking()) {
     keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
   }
 
@@ -321,6 +328,68 @@ void Particle::event_advance()
   }
 }
 
+void Particle::event_delta_advance()
+{
+  if (E() != E_last()) {
+    update_majorant();
+  }
+
+  // Sample distance to next position
+  if (type() == ParticleType::electron() ||
+      type() == ParticleType::positron()) {
+    // Electrons/positrons don't move
+    collision_distance() = 0.0;
+  } else if (majorant() == 0.0) {
+    // For a void majorant (rare but possible for a source in a void),
+    // the collision distance is infinity.
+    collision_distance() = INFINITY;
+  } else {
+    // Sample collision distance based on the majorant for this energy.
+    collision_distance() = -std::log(prn(current_seed())) / majorant();
+  }
+
+  // Update distance to problem boundary. Particles with large majorant
+  // cross sections will tunnel out of the domain if a floating point
+  // tolerance is not specified on the boundary distance calculation.
+  boundary() = distance_to_external_boundary(*this);
+  boundary().distance() -= FP_REL_PRECISION;
+
+  double speed = this->speed();
+  double time_cutoff = settings::time_cutoff[type().transport_index()];
+  double distance_cutoff =
+    (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
+
+  // Move to the external boundary, delta tracking collision site, or time
+  // cutoff distance.
+  double distance =
+    std::min({collision_distance(), boundary().distance(), distance_cutoff});
+  move_distance(distance);
+
+  // Advance particle in time.
+  double dt = distance / speed;
+  time() += dt;
+  lifetime() += dt;
+
+  // Need to locate the particle at the collision site or boundary.
+  for (int j = 0; j < n_coord(); ++j) {
+    coord(j).reset();
+  }
+  if (!exhaustive_find_cell(*this)) {
+    // We've lost this particle.
+    mark_as_lost(fmt::format(
+      "Particle {} could not be located while running delta tracking!", id()));
+    return;
+  }
+
+  // Force re-calculation of material properties at the collision site.
+  material_last() = C_NONE;
+
+  // Set particle weight to zero if it hit the time boundary
+  if (distance == distance_cutoff) {
+    wgt() = 0.0;
+  }
+}
+
 void Particle::event_cross_surface()
 {
   // Saving previous cell data
@@ -385,7 +454,6 @@ void Particle::event_cross_surface()
 
 void Particle::event_collide()
 {
-
   // Score collision estimate of keff
   if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
     keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
@@ -569,7 +637,7 @@ void Particle::event_death()
 #pragma omp atomic
       global_tally_collision += k_collision;
     }
-    if (k_tracklength != 0.0) {
+    if (k_tracklength != 0.0 && !settings::delta_tracking) {
 #pragma omp atomic
       global_tally_tracklength += k_tracklength;
     }
@@ -851,6 +919,30 @@ void Particle::cross_periodic_bc(
   }
 }
 
+void Particle::update_majorant()
+{
+  if (type().is_neutron()) {
+    majorant() = NeutronMajorant::safety_factor_ *
+                 data::n_majorant->calculate_neutron_xs(E());
+  } else if (type().is_photon()) {
+    majorant() = PhotonMajorant::safety_factor_ *
+                 data::p_majorant->calculate_photon_xs(E());
+  }
+}
+
+bool Particle::kill_invalid_maj()
+{
+  if (alive() && (macro_xs().total > majorant())) {
+    mark_as_lost(
+      fmt::format("Ratio of the total cross section ({}) to the majorant "
+                  "cross section ({}) for particle {} ({}) with energy {} is "
+                  "greater than unity!",
+        macro_xs().total, majorant(), id(), type().str(), E()));
+    return true;
+  }
+  return false;
+}
+
 void Particle::mark_as_lost(const char* message)
 {
   // Print warning and write lost particle file

src/settings.cpp

@@ -53,6 +53,7 @@ bool confidence_intervals {false};
 bool create_delayed_neutrons {true};
 bool create_fission_neutrons {true};
 bool delayed_photon_scaling {true};
+bool delta_tracking {false};
 bool entropy_on {false};
 bool event_based {false};
 bool ifp_on {false};
@@ -1229,6 +1230,22 @@ void read_settings_xml(pugi::xml_node root)
     event_based = get_node_value_bool(root, "event_based");
   }
 
+  // Check whether or not to use delta tracking
+  if (check_for_node(root, "delta_tracking")) {
+    delta_tracking = get_node_value_bool(root, "delta_tracking");
+
+    if (temperature_multipole && delta_tracking) {
+      fatal_error(
+        "At present, delta tracking cannot be used with a windowed multipole "
+        "temperature treatment.");
+    }
+
+    if (!run_CE && delta_tracking) {
+      fatal_error("At present, delta tracking can only be used in continuous "
+                  "energy simulations.");
+    }
+  }
+
   // Check whether material cell offsets should be generated
   if (check_for_node(root, "material_cell_offsets")) {
     material_cell_offsets = get_node_value_bool(root, "material_cell_offsets");

src/simulation.cpp

@@ -9,6 +9,7 @@
 #include "openmc/event.h"
 #include "openmc/geometry_aux.h"
 #include "openmc/ifp.h"
+#include "openmc/majorant.h"
 #include "openmc/material.h"
 #include "openmc/message_passing.h"
 #include "openmc/nuclide.h"
@@ -88,6 +89,11 @@ int openmc_simulation_init()
     initialize_data();
   }
 
+  // Create the majorant cross sections for delta tracking.
+  if (settings::delta_tracking) {
+    create_majorants();
+  }
+
   // Determine how much work each process should do
   calculate_work(settings::n_particles);
 
@@ -245,6 +251,9 @@ int openmc_simulation_finalize()
   if (settings::check_overlaps)
     print_overlap_check();
 
+  // Clear majorants as they could change if OpenMC is run again.
+  reset_majorants();
+
   // Reset flags
   simulation::initialized = false;
   return 0;
@@ -633,8 +642,13 @@ void initialize_generation()
       ufs_count_sites();
 
     // Store current value of tracklength k
-    simulation::keff_generation = simulation::global_tallies(
-      GlobalTally::K_TRACKLENGTH, TallyResult::VALUE);
+    if (settings::delta_tracking) {
+      simulation::keff_generation = simulation::global_tallies(
+        GlobalTally::K_COLLISION, TallyResult::VALUE);
+    } else {
+      simulation::keff_generation = simulation::global_tallies(
+        GlobalTally::K_TRACKLENGTH, TallyResult::VALUE);
+    }
   }
 }
 
@@ -765,6 +779,12 @@ void initialize_particle_track(
     write_message("Simulating Particle {}", p.id());
   }
 
+  // Compute the majorant and set the delta tracking flag.
+  if (settings::delta_tracking) {
+    p.delta_tracking() = true;
+    p.update_majorant();
+  }
+
   // Add particle's starting weight to count for normalizing tallies later
   if (!is_secondary) {
 #pragma omp atomic
@@ -838,10 +858,13 @@ void initialize_data()
   // Determine minimum/maximum energy for incident neutron/photon data
   data::energy_max = {INFTY, INFTY, INFTY, INFTY};
   data::energy_min = {0.0, 0.0, 0.0, 0.0};
+  int neutron = ParticleType::neutron().transport_index();
+  int photon = ParticleType::photon().transport_index();
+  int electron = ParticleType::electron().transport_index();
+  int positron = ParticleType::positron().transport_index();
 
   for (const auto& nuc : data::nuclides) {
     if (nuc->grid_.size() >= 1) {
-      int neutron = ParticleType::neutron().transport_index();
       data::energy_min[neutron] =
         std::max(data::energy_min[neutron], nuc->grid_[0].energy.front());
       data::energy_max[neutron] =
@@ -852,7 +875,6 @@ void initialize_data()
   if (settings::photon_transport) {
     for (const auto& elem : data::elements) {
       if (elem->energy_.size() >= 1) {
-        int photon = ParticleType::photon().transport_index();
         int n = elem->energy_.size();
         data::energy_min[photon] =
           std::max(data::energy_min[photon], std::exp(elem->energy_(1)));
@@ -865,9 +887,6 @@ void initialize_data()
       // Determine if minimum/maximum energy for bremsstrahlung is greater/less
       // than the current minimum/maximum
       if (data::ttb_e_grid.size() >= 1) {
-        int photon = ParticleType::photon().transport_index();
-        int electron = ParticleType::electron().transport_index();
-        int positron = ParticleType::positron().transport_index();
         int n_e = data::ttb_e_grid.size();
 
         const std::vector<int> charged = {electron, positron};
@@ -891,7 +910,6 @@ void initialize_data()
     // grid has not been allocated
     if (nuc->grid_.size() > 0) {
       double max_E = nuc->grid_[0].energy.back();
-      int neutron = ParticleType::neutron().transport_index();
       if (max_E == data::energy_max[neutron]) {
         write_message(7, "Maximum neutron transport energy: {} eV for {}",
           data::energy_max[neutron], nuc->name_);
@@ -908,7 +926,6 @@ void initialize_data()
   for (auto& nuc : data::nuclides) {
     nuc->init_grid();
   }
-  int neutron = ParticleType::neutron().transport_index();
   simulation::log_spacing =
     std::log(data::energy_max[neutron] / data::energy_min[neutron]) /
     settings::n_log_bins;
@@ -974,6 +991,40 @@ void transport_history_based_single_particle(Particle& p)
   p.event_death();
 }
 
+void transport_delta_history_based_single_particle(Particle& p)
+{
+  while (p.alive()) {
+    p.event_delta_advance();
+
+    if (p.alive()) {
+      // Electrons and positrons collide in-place, no need to rejection sample.
+      if (p.type() == ParticleType::electron() ||
+          p.type() == ParticleType::positron()) {
+        p.event_collide();
+      }
+
+      if (p.alive() && p.collision_distance() < p.boundary().distance()) {
+        // Collided before hitting an external boundary. Rejection sample the
+        // majorant.
+        p.event_calculate_xs();
+        if (p.kill_invalid_maj()) {
+          break;
+        }
+        if (p.alive() &&
+            (prn(p.current_seed()) < (p.macro_xs().total / p.majorant()))) {
+          p.event_collide();
+        }
+      } else if (p.alive()) {
+        // Crossed an external boundary before colliding.
+        p.event_cross_surface();
+      }
+    }
+
+    p.event_check_limit_and_revive();
+  }
+  p.event_death();
+}
+
 void transport_history_based()
 {
 #pragma omp parallel
@@ -982,7 +1033,11 @@ void transport_history_based()
 #pragma omp for schedule(runtime)
     for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
       initialize_particle_track(p, i_work, false);
-      transport_history_based_single_particle(p);
+      if (settings::delta_tracking) {
+        transport_delta_history_based_single_particle(p);
+      } else {
+        transport_history_based_single_particle(p);
+      }
     }
   }
 }
@@ -1023,7 +1078,11 @@ void transport_history_based_shared_secondary()
 #pragma omp for schedule(runtime)
     for (int64_t i = 1; i <= simulation::work_per_rank; i++) {
       initialize_particle_track(p, i, false);
-      transport_history_based_single_particle(p);
+      if (settings::delta_tracking) {
+        transport_delta_history_based_single_particle(p);
+      } else {
+        transport_history_based_single_particle(p);
+      }
       for (auto& site : p.local_secondary_bank()) {
         thread_bank.push_back(site);
       }
@@ -1083,7 +1142,11 @@ void transport_history_based_shared_secondary()
         initialize_particle_track(p, i, true);
         SourceSite& site = simulation::shared_secondary_bank_read[i - 1];
         p.event_revive_from_secondary(site);
-        transport_history_based_single_particle(p);
+        if (settings::delta_tracking) {
+          transport_delta_history_based_single_particle(p);
+        } else {
+          transport_history_based_single_particle(p);
+        }
         for (auto& secondary_site : p.local_secondary_bank()) {
           thread_bank.push_back(secondary_site);
         }
@@ -1119,8 +1182,13 @@ void transport_event_based()
       std::min(remaining_work, settings::max_particles_in_flight);
 
     // Initialize all particle histories for this subiteration
-    process_init_events(n_particles, source_offset);
-    process_transport_events();
+    if (settings::delta_tracking) {
+      process_delta_init_events(n_particles, source_offset);
+      process_delta_transport_events();
+    } else {
+      process_init_events(n_particles, source_offset);
+      process_transport_events();
+    }
     process_death_events(n_particles);
 
     // Adjust remaining work and source offset variables
@@ -1154,8 +1222,13 @@ void transport_event_based_shared_secondary()
     int64_t n_particles =
       std::min(remaining_work, settings::max_particles_in_flight);
 
-    process_init_events(n_particles, source_offset);
-    process_transport_events();
+    if (settings::delta_tracking) {
+      process_delta_init_events(n_particles, source_offset);
+      process_delta_transport_events();
+    } else {
+      process_init_events(n_particles, source_offset);
+      process_transport_events();
+    }
     process_death_events(n_particles);
 
     collect_event_secondary_banks(n_particles);
@@ -1217,9 +1290,15 @@ void transport_event_based_shared_secondary()
       int64_t n_particles =
         std::min(sec_remaining, settings::max_particles_in_flight);
 
-      process_init_secondary_events(
-        n_particles, sec_offset, simulation::shared_secondary_bank_read);
-      process_transport_events();
+      if (settings::delta_tracking) {
+        process_delta_init_secondary_events(
+          n_particles, sec_offset, simulation::shared_secondary_bank_read);
+        process_delta_transport_events();
+      } else {
+        process_init_secondary_events(
+          n_particles, sec_offset, simulation::shared_secondary_bank_read);
+        process_transport_events();
+      }
       process_death_events(n_particles);
 
       collect_event_secondary_banks(n_particles);

src/state_point.cpp

@@ -114,6 +114,7 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
       break;
     }
     write_attribute(file_id, "photon_transport", settings::photon_transport);
+    write_attribute(file_id, "delta_tracking", settings::delta_tracking);
     write_dataset(file_id, "n_particles", settings::n_particles);
     write_dataset(file_id, "n_batches", settings::n_batches);
 
@@ -313,6 +314,10 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
     hid_t runtime_group = create_group(file_id, "runtime");
     write_dataset(
       runtime_group, "total initialization", time_initialize.elapsed());
+    if (settings::delta_tracking) {
+      write_dataset(
+        runtime_group, "build majorant", time_build_majorant.elapsed());
+    }
     write_dataset(
       runtime_group, "reading cross sections", time_read_xs.elapsed());
     write_dataset(runtime_group, "simulation",

src/tallies/tally.cpp

@@ -350,6 +350,12 @@ Tally::Tally(pugi::xml_node node)
     this->init_triggers(node);
   }
 
+  // If we're running with delta tracking, the default estimator must be set to
+  // collision or analog.
+  if (settings::delta_tracking && estimator_ == TallyEstimator::TRACKLENGTH) {
+    estimator_ = TallyEstimator::COLLISION;
+  }
+
   // =======================================================================
   // SET TALLY ESTIMATOR
 

src/timer.cpp

@@ -16,6 +16,7 @@ Timer time_finalize;
 Timer time_inactive;
 Timer time_initialize;
 Timer time_read_xs;
+Timer time_build_majorant;
 Timer time_statepoint;
 Timer time_tallies;
 Timer time_total;
@@ -75,6 +76,7 @@ void reset_timers()
   simulation::time_finalize.reset();
   simulation::time_inactive.reset();
   simulation::time_initialize.reset();
+  simulation::time_build_majorant.reset();
   simulation::time_read_xs.reset();
   simulation::time_statepoint.reset();
   simulation::time_tallies.reset();

tests/regression_tests/delta_tracking_bcs/periodic/inputs_true.dat

@@ -0,0 +1,60 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="UO2" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="2" name="light water">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <cell id="2" material="2" region="1" universe="1"/>
+    <cell id="3" fill="2" region="2 -3 4 -5" universe="3"/>
+    <lattice id="2">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+1 1 
+1 1 </universes>
+    </lattice>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 0.4"/>
+    <surface id="2" name="minimum x" type="x-plane" boundary="periodic" coeffs="-1.26" periodic_surface_id="3"/>
+    <surface id="3" name="maximum x" type="x-plane" boundary="periodic" coeffs="1.26" periodic_surface_id="2"/>
+    <surface id="4" name="minimum y" type="y-plane" boundary="periodic" coeffs="-1.26" periodic_surface_id="5"/>
+    <surface id="5" name="maximum y" type="y-plane" boundary="periodic" coeffs="1.26" periodic_surface_id="4"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <photon_transport>false</photon_transport>
+    <delta_tracking>true</delta_tracking>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/delta_tracking_bcs/periodic/results_true.dat

@@ -0,0 +1,11 @@
+k-combined:
+1.837808E+00 6.243816E-03
+tally 1:
+1.874300E+01
+7.031590E+01
+1.879200E+01
+7.068119E+01
+1.870100E+01
+7.000963E+01
+1.852700E+01
+6.868985E+01

tests/regression_tests/delta_tracking_bcs/reflective/inputs_true.dat

@@ -0,0 +1,60 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="UO2" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="2" name="light water">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <cell id="2" material="2" region="1" universe="1"/>
+    <cell id="3" fill="2" region="2 -3 4 -5" universe="3"/>
+    <lattice id="2">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+1 1 
+1 1 </universes>
+    </lattice>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 0.4"/>
+    <surface id="2" name="minimum x" type="x-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="3" name="maximum x" type="x-plane" boundary="reflective" coeffs="1.26"/>
+    <surface id="4" name="minimum y" type="y-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="5" name="maximum y" type="y-plane" boundary="reflective" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <photon_transport>false</photon_transport>
+    <delta_tracking>true</delta_tracking>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/delta_tracking_bcs/reflective/results_true.dat

@@ -0,0 +1,11 @@
+k-combined:
+1.850729E+00 4.144286E-03
+tally 1:
+1.828800E+01
+6.694305E+01
+1.823900E+01
+6.658124E+01
+1.894000E+01
+7.184607E+01
+1.852200E+01
+6.869212E+01

tests/regression_tests/delta_tracking_bcs/test.py

@@ -0,0 +1,14 @@
+import openmc
+import pytest
+from openmc.utility_funcs import change_directory
+from openmc.examples import delta_tracking_lattice
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.mark.parametrize("boundary", ['vacuum', 'reflective', 'periodic', 'white'])
+def test_lattice(boundary):
+    with change_directory(boundary):
+        model = delta_tracking_lattice(False, boundary)
+        harness = PyAPITestHarness('statepoint.10.h5', model)
+        harness.main()

tests/regression_tests/delta_tracking_bcs/vacuum/inputs_true.dat

@@ -0,0 +1,60 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="UO2" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="2" name="light water">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <cell id="2" material="2" region="1" universe="1"/>
+    <cell id="3" fill="2" region="2 -3 4 -5" universe="3"/>
+    <lattice id="2">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+1 1 
+1 1 </universes>
+    </lattice>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 0.4"/>
+    <surface id="2" name="minimum x" type="x-plane" boundary="vacuum" coeffs="-1.26"/>
+    <surface id="3" name="maximum x" type="x-plane" boundary="vacuum" coeffs="1.26"/>
+    <surface id="4" name="minimum y" type="y-plane" boundary="vacuum" coeffs="-1.26"/>
+    <surface id="5" name="maximum y" type="y-plane" boundary="vacuum" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <photon_transport>false</photon_transport>
+    <delta_tracking>true</delta_tracking>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/delta_tracking_bcs/vacuum/results_true.dat

@@ -0,0 +1,11 @@
+k-combined:
+6.926354E-02 3.739970E-03
+tally 1:
+9.270000E-01
+1.739570E-01
+1.017000E+00
+2.077250E-01
+8.920000E-01
+1.607780E-01
+9.290000E-01
+1.727250E-01

tests/regression_tests/delta_tracking_bcs/white/inputs_true.dat

@@ -0,0 +1,60 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="UO2" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="2" name="light water">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <cell id="2" material="2" region="1" universe="1"/>
+    <cell id="3" fill="2" region="2 -3 4 -5" universe="3"/>
+    <lattice id="2">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+1 1 
+1 1 </universes>
+    </lattice>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 0.4"/>
+    <surface id="2" name="minimum x" type="x-plane" boundary="white" coeffs="-1.26"/>
+    <surface id="3" name="maximum x" type="x-plane" boundary="white" coeffs="1.26"/>
+    <surface id="4" name="minimum y" type="y-plane" boundary="white" coeffs="-1.26"/>
+    <surface id="5" name="maximum y" type="y-plane" boundary="white" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <photon_transport>false</photon_transport>
+    <delta_tracking>true</delta_tracking>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/delta_tracking_bcs/white/results_true.dat

@@ -0,0 +1,11 @@
+k-combined:
+1.842649E+00 6.248566E-03
+tally 1:
+1.844200E+01
+6.808763E+01
+1.903200E+01
+7.249322E+01
+1.858800E+01
+6.913862E+01
+1.866100E+01
+6.966707E+01

tests/regression_tests/delta_tracking_distribrho/False/inputs_true.dat

@@ -0,0 +1,62 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="UO2" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="2" name="light water">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1">
+      <density>10.0 20.0 10.0 20.0</density>
+    </cell>
+    <cell id="2" material="2" region="1" universe="1"/>
+    <cell id="3" fill="2" region="2 -3 4 -5" universe="3"/>
+    <lattice id="2">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+1 1 
+1 1 </universes>
+    </lattice>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 0.4"/>
+    <surface id="2" name="minimum x" type="x-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="3" name="maximum x" type="x-plane" boundary="reflective" coeffs="1.26"/>
+    <surface id="4" name="minimum y" type="y-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="5" name="maximum y" type="y-plane" boundary="reflective" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <photon_transport>false</photon_transport>
+    <delta_tracking>true</delta_tracking>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/delta_tracking_distribrho/False/results_true.dat

@@ -0,0 +1,11 @@
+k-combined:
+1.860197E+00 2.386186E-03
+tally 1:
+1.633700E+01
+5.350220E+01
+1.886400E+01
+7.123130E+01
+1.617900E+01
+5.238011E+01
+1.833900E+01
+6.734011E+01

tests/regression_tests/delta_tracking_distribrho/True/inputs_true.dat

@@ -0,0 +1,70 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="UO2" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="2" name="light water">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1">
+      <density>10.0 20.0 10.0 20.0</density>
+    </cell>
+    <cell id="2" material="2" region="1" universe="1"/>
+    <cell id="3" fill="2" region="2 -3 4 -5" universe="3"/>
+    <lattice id="2">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+1 1 
+1 1 </universes>
+    </lattice>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 0.4"/>
+    <surface id="2" name="minimum x" type="x-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="3" name="maximum x" type="x-plane" boundary="reflective" coeffs="1.26"/>
+    <surface id="4" name="minimum y" type="y-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="5" name="maximum y" type="y-plane" boundary="reflective" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <photon_transport>true</photon_transport>
+    <delta_tracking>true</delta_tracking>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="3" type="particle">
+      <bins>photon</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+    <tally id="2">
+      <filters>3 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/delta_tracking_distribrho/True/results_true.dat

@@ -0,0 +1,20 @@
+k-combined:
+1.863429E+00 6.240045E-03
+tally 1:
+1.638600E+01
+5.380010E+01
+1.836200E+01
+6.749015E+01
+1.596600E+01
+5.106840E+01
+1.903000E+01
+7.262313E+01
+tally 2:
+6.547989E+01
+8.594162E+02
+1.145708E+02
+2.625502E+03
+6.537182E+01
+8.551704E+02
+1.154863E+02
+2.671966E+03

tests/regression_tests/delta_tracking_distribrho/test.py

@@ -0,0 +1,14 @@
+import openmc
+import pytest
+from openmc.utility_funcs import change_directory
+from openmc.examples import delta_tracking_lattice
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.mark.parametrize("photon", [False, True])
+def test_lattice_checkerboard(photon):
+    with change_directory(str(photon)):
+        model = delta_tracking_lattice(photon, densities=[10.0, 20.0, 10.0, 20.0])
+        harness = PyAPITestHarness('statepoint.10.h5', model)
+        harness.main()

tests/regression_tests/delta_tracking_event/False/inputs_true.dat

@@ -0,0 +1,61 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="UO2" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="2" name="light water">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <cell id="2" material="2" region="1" universe="1"/>
+    <cell id="3" fill="2" region="2 -3 4 -5" universe="3"/>
+    <lattice id="2">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+1 1 
+1 1 </universes>
+    </lattice>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 0.4"/>
+    <surface id="2" name="minimum x" type="x-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="3" name="maximum x" type="x-plane" boundary="reflective" coeffs="1.26"/>
+    <surface id="4" name="minimum y" type="y-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="5" name="maximum y" type="y-plane" boundary="reflective" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <photon_transport>false</photon_transport>
+    <event_based>true</event_based>
+    <delta_tracking>true</delta_tracking>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/delta_tracking_event/False/results_true.dat

@@ -0,0 +1,11 @@
+k-combined:
+1.850729E+00 4.144286E-03
+tally 1:
+1.828800E+01
+6.694305E+01
+1.823900E+01
+6.658124E+01
+1.894000E+01
+7.184607E+01
+1.852200E+01
+6.869212E+01

tests/regression_tests/delta_tracking_event/True/inputs_true.dat

@@ -0,0 +1,69 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="UO2" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="2" name="light water">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <cell id="2" material="2" region="1" universe="1"/>
+    <cell id="3" fill="2" region="2 -3 4 -5" universe="3"/>
+    <lattice id="2">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+1 1 
+1 1 </universes>
+    </lattice>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 0.4"/>
+    <surface id="2" name="minimum x" type="x-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="3" name="maximum x" type="x-plane" boundary="reflective" coeffs="1.26"/>
+    <surface id="4" name="minimum y" type="y-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="5" name="maximum y" type="y-plane" boundary="reflective" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <photon_transport>true</photon_transport>
+    <event_based>true</event_based>
+    <delta_tracking>true</delta_tracking>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="3" type="particle">
+      <bins>photon</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+    <tally id="2">
+      <filters>3 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/delta_tracking_event/True/results_true.dat

@@ -0,0 +1,20 @@
+k-combined:
+1.837786E+00 3.391810E-03
+tally 1:
+1.864400E+01
+6.959490E+01
+1.872400E+01
+7.018429E+01
+1.855700E+01
+6.894280E+01
+1.811100E+01
+6.571011E+01
+tally 2:
+8.578019E+01
+1.473018E+03
+8.729418E+01
+1.526265E+03
+8.810551E+01
+1.553484E+03
+8.716879E+01
+1.520479E+03

tests/regression_tests/delta_tracking_event/test.py

@@ -0,0 +1,15 @@
+import openmc
+import pytest
+from openmc.utility_funcs import change_directory
+from openmc.examples import delta_tracking_lattice
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.mark.parametrize("photon", [False, True])
+def test_lattice(photon):
+    with change_directory(str(photon)):
+        model = delta_tracking_lattice(photon)
+        model.settings.event_based = True
+        harness = PyAPITestHarness('statepoint.10.h5', model)
+        harness.main()

tests/regression_tests/delta_tracking_history/False/inputs_true.dat

@@ -0,0 +1,60 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="UO2" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="2" name="light water">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <cell id="2" material="2" region="1" universe="1"/>
+    <cell id="3" fill="2" region="2 -3 4 -5" universe="3"/>
+    <lattice id="2">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+1 1 
+1 1 </universes>
+    </lattice>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 0.4"/>
+    <surface id="2" name="minimum x" type="x-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="3" name="maximum x" type="x-plane" boundary="reflective" coeffs="1.26"/>
+    <surface id="4" name="minimum y" type="y-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="5" name="maximum y" type="y-plane" boundary="reflective" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <photon_transport>false</photon_transport>
+    <delta_tracking>true</delta_tracking>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/delta_tracking_history/False/results_true.dat

@@ -0,0 +1,11 @@
+k-combined:
+1.850729E+00 4.144286E-03
+tally 1:
+1.828800E+01
+6.694305E+01
+1.823900E+01
+6.658124E+01
+1.894000E+01
+7.184607E+01
+1.852200E+01
+6.869212E+01

tests/regression_tests/delta_tracking_history/True/inputs_true.dat

@@ -0,0 +1,68 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="UO2" depletable="true">
+      <density value="10.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="2" name="light water">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <cell id="2" material="2" region="1" universe="1"/>
+    <cell id="3" fill="2" region="2 -3 4 -5" universe="3"/>
+    <lattice id="2">
+      <pitch>1.26 1.26</pitch>
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <universes>
+1 1 
+1 1 </universes>
+    </lattice>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 0.4"/>
+    <surface id="2" name="minimum x" type="x-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="3" name="maximum x" type="x-plane" boundary="reflective" coeffs="1.26"/>
+    <surface id="4" name="minimum y" type="y-plane" boundary="reflective" coeffs="-1.26"/>
+    <surface id="5" name="maximum y" type="y-plane" boundary="reflective" coeffs="1.26"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <photon_transport>true</photon_transport>
+    <delta_tracking>true</delta_tracking>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-1.26 -1.26</lower_left>
+      <upper_right>1.26 1.26</upper_right>
+    </mesh>
+    <filter id="1" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="2" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="3" type="particle">
+      <bins>photon</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+    <tally id="2">
+      <filters>3 2</filters>
+      <scores>total</scores>
+      <estimator>collision</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/delta_tracking_history/True/results_true.dat

@@ -0,0 +1,20 @@
+k-combined:
+1.837786E+00 3.391810E-03
+tally 1:
+1.864400E+01
+6.959490E+01
+1.872400E+01
+7.018429E+01
+1.855700E+01
+6.894280E+01
+1.811100E+01
+6.571011E+01
+tally 2:
+8.578019E+01
+1.473018E+03
+8.729418E+01
+1.526265E+03
+8.810551E+01
+1.553484E+03
+8.716879E+01
+1.520479E+03

tests/regression_tests/delta_tracking_history/test.py

@@ -0,0 +1,14 @@
+import openmc
+import pytest
+from openmc.utility_funcs import change_directory
+from openmc.examples import delta_tracking_lattice
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.mark.parametrize("photon", [False, True])
+def test_lattice(photon):
+    with change_directory(str(photon)):
+        model = delta_tracking_lattice(photon)
+        harness = PyAPITestHarness('statepoint.10.h5', model)
+        harness.main()

tests/regression_tests/mg_survival_biasing/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-9.889968E-01 9.144186E-03
+9.889968E-01 1.823442E-02

tests/unit_tests/test_settings.py

@@ -94,6 +94,7 @@ def test_export_to_xml(run_in_tmpdir):
     s.max_secondaries = 1_000_000
     s.source_rejection_fraction = 0.01
     s.free_gas_threshold = 800.0
+    s.delta_tracking = True
 
     # Make sure exporting XML works
     s.export_to_xml()
@@ -190,6 +191,7 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.max_secondaries == 1_000_000
     assert s.source_rejection_fraction == 0.01
     assert s.free_gas_threshold == 800.0
+    assert s.delta_tracking == True
 
 
 def test_properties_file_load(tmp_path, mpi_intracomm):