Delta tracking

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/methods/neutron_physics.rst

@@ -52,6 +52,106 @@ 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 many geometric regions
+across particle flights, surface tracking will require a large number of surface
+distance calculations and collision distance samples. 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 by adding a fictitious delta scattering cross sections
+:math:`\Sigma_{\delta}(\mathbf{r}, E)` to the medium. The sum of the delta
+scattering cross section and the total cross section is known as the majorant
+cross section (:math:`\Sigma_{maj}(E)`)
+
+.. math::
+    :label: majorant-xs-1
+
+    \Sigma_{maj}(E) = \Sigma_{t}(\mathbf{r}, E) + \Sigma_{\delta}(\mathbf{r},
+    E),
+
+which is computed as the maximum macroscopic total cross section over the
+spatial domain
+
+.. math::
+    :label: majorant-xs-2
+
+    \Sigma_{maj}(E) = \max_{\mathbf{r}}\left(\Sigma_{t}(\mathbf{r}, E)\right).
+
+The delta tracking procedure computes the 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)}.
+
+This is followed by a rejection sampling test to determine if a real
+collision or a delta scattering collision occured at this delta collision
+point. 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 is true, the collision is accepted as real. If the condition
+is false, a delta scatter event has occured 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 compared to surface tracking as it allows
+for continuously-varying material properties due to the use of pointwise
+cross section samples in Equation :eq:`delta-real-collision`. Delta tracking
+often performs better than surface tracking in problems where the particle
+mean free path is larger than the distance between surfaces. Problems
+containing small regions with large total cross sections (such as burnable
+absorbers) will have majorant cross sections several orders of magnitude in
+larger than the total cross section over the majority of the domain
+[Leppänen]_. This decreases the number of accepted collisions, and therefore
+the effectiveness of delta tracking. Material discontinuities are not
+considered in delta tracking, which prohibits the use of track length
+estimators and forces the use of the higher- variance collision estimator
+(discussed in detail in the :ref:`methods_tallies` section).
+
 ----------------------------------------------------
 :math:`(n,\gamma)` and Other Disappearance Reactions
 ----------------------------------------------------
@@ -1631,6 +1731,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,13 +1842,22 @@ 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).
 
 .. [Williams] M. M. R. Williams, *The Slowing Down and Thermalization of
    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

@@ -73,10 +73,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