Implement surface flux tallies (#3742)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Author: GuySten

docs/source/io_formats/settings.rst

@@ -1234,6 +1234,23 @@ attributes/sub-elements:
           are not eligible to store any particles when using ``cell``, ``cellfrom``
           or ``cellto`` attributes. It is recommended to use surface IDs instead.
 
+------------------------------------
+``<surface_grazing_cutoff>`` Element
+------------------------------------
+
+The ``<surface_grazing_cutoff>`` element specifies the surface flux cosine cutoff.
+
+  *Default*: 0.001
+
+-----------------------------------
+``<surface_grazing_ratio>`` Element
+-----------------------------------
+
+The ``<surface_grazing_ratio>`` element specifies the surface flux cosine 
+substitution ratio.
+
+  *Default*: 0.5
+
 ------------------------------
 ``<survival_biasing>`` Element
 ------------------------------

docs/source/methods/tallies.rst

@@ -205,7 +205,71 @@ had a collision at every event. Thus, for tallies with outgoing-energy filters
 or for tallies of scattering moments (which require the scattering cosine of
 the change-in-angle), we must use an analog estimator.
 
-.. TODO: Add description of surface current tallies
+-----------------------------------
+Surface-Integrated Flux and Current
+-----------------------------------
+
+Surface tallies allow you to measure particle behavior as they cross specific
+boundaries in your geometry. Unlike volume tallies, which integrate over a
+volumetric region, surface tallies capture the current or flux passing through a
+surface. Surface tallies are estimated using an analog estimator.
+
+Current Score
+-------------
+
+When tallying the current across a surface, we simply count the weight of
+particles that cross the surface of interest:
+
+
+.. math::
+    :label: analog-current-estimator
+
+    J = \frac{1}{W} \sum_{i \in S} w_i.
+
+where :math:`J` is the area-integrated current passing through surface
+:math:`S`, :math:`W` is the total starting weight of the particles, and
+:math:`w_i` is the weight of the particle as it crosses the surface :math:`S`.
+
+Flux Score
+----------
+
+When tallying flux over a surface, we use the relationship between current and
+flux:
+
+
+.. math::
+    :label: surface-flux-estimator
+
+    \phi_S = \frac{1}{W} \sum_{i \in S} \frac{w_i}{|\mu|}.
+
+where :math:`\phi_S` is the area-integrated flux over surface :math:`S`,
+:math:`W` is the total starting weight of the particles, :math:`w_i` is the
+weight of the particle as it crosses the surface :math:`S` and :math:`\mu` is
+the cosine of angle between the particle direction and the surface normal.
+
+This equation diverges when the particle crossing the surface is nearly parallel
+to it (that is, as :math:`\mu` approaches zero). To remove this divergence,
+OpenMC scores:
+
+.. math::
+    :label: modified-surface-flux-estimator
+
+    \phi_S = \frac{1}{W} \sum_{i \in S} w_i f(\mu).
+
+and the function :math:`f` is defined by:
+
+.. math::
+     f(\mu) = \begin{cases}
+      \frac{1}{|\mu|} & |\mu| > \mu_\text{cut} \\
+      \frac{1}{c\mu_\text{cut}} & |\mu| \le \mu_\text{cut}
+     \end{cases}
+
+where :math:`\mu_\text{cut}` is the grazing cosine cutoff and :math:`c` is the
+cosine substitution ratio. The parameters :math:`\mu_\text{cut}` and :math:`c`
+can be set by the user via the :attr:`openmc.Settings.surface_grazing_cutoff`
+and :attr:`openmc.Settings.surface_grazing_ratio` attributes, respectively. The
+default values for these parameters are 0.001 and 0.5 as recommended by
+`Favorite, Thomas, and Booth <https://doi.org/10.13182/NSE09-72>`_.
 
 .. _tallies_statistics:
 

docs/source/usersguide/tallies.rst

@@ -261,18 +261,21 @@ The following tables show all valid scores:
     +----------------------+---------------------------------------------------+
     |Score                 | Description                                       |
     +======================+===================================================+
-    |current               |Used in combination with a meshsurface filter:     |
+    |current               |It may not be used in conjunction with any other   |
+    |                      |score except flux.                                 |
+    |                      |                                                   |
+    |                      |When used in combination with a meshsurface filter:|
     |                      |Partial currents on the boundaries of each cell in |
-    |                      |a mesh. It may not be used in conjunction with any |
-    |                      |other score. Only energy and mesh filters may be   |
-    |                      |used.                                              |
-    |                      |Used in combination with a surface filter:         |
+    |                      |a mesh.                                            |
+    |                      |                                                   |
+    |                      |When used in combination with a surface filter:    |
     |                      |Net currents on any surface previously defined in  |
     |                      |the geometry. It may be used along with any other  |
     |                      |filter, except meshsurface filters.                |
     |                      |Surfaces can alternatively be defined with cell    |
     |                      |from and cell filters thereby resulting in tallying|
     |                      |partial currents.                                  |
+    |                      |                                                   |
     |                      |Units are particles per source particle.           |
     +----------------------+---------------------------------------------------+
     |events                |Number of scoring events. Units are events per     |

include/openmc/settings.h

@@ -181,6 +181,8 @@ extern int64_t ssw_cell_id;       //!< Cell id for the surface source
                                   //!< write setting
 extern SSWCellType ssw_cell_type; //!< Type of option for the cell
                                   //!< argument of surface source write
+extern double surface_grazing_cutoff; //!< surface flux cosine cutoff
+extern double surface_grazing_ratio;  //!< surface flux substitution ratio
 extern TemperatureMethod
   temperature_method; //!< method for choosing temperatures
 extern double

include/openmc/tallies/tally_scoring.h

@@ -101,11 +101,19 @@ void score_tracklength_tally(Particle& p, double distance);
 //! \param total_distance The distance in [cm] traveled by the particle
 void score_timed_tracklength_tally(Particle& p, double total_distance);
 
-//! Score surface or mesh-surface tallies for particle currents.
+//! Score mesh-surface tallies for particle currents.
 //
 //! \param p The particle being tracked
 //! \param tallies A vector of the indices of the tallies to score to
-void score_surface_tally(Particle& p, const vector<int>& tallies);
+void score_meshsurface_tally(Particle& p, const vector<int>& tallies);
+
+//! Score surface tallies for particle currents.
+//
+//! \param p The particle being tracked
+//! \param tallies A vector of the indices of the tallies to score to
+//! \param surf The surface being crossed
+void score_surface_tally(
+  Particle& p, const vector<int>& tallies, const Surface& surf);
 
 //! Score the pulse-height tally
 //! This is triggered at the end of every particle history

openmc/settings.py

@@ -289,6 +289,14 @@ class Settings:
         :cellto: Cell ID used to determine if particles crossing identified
                  surfaces are to be banked. Particles going to this declared
                  cell will be banked (int)
+    surface_grazing_cutoff : float
+        Surface flux cosine cutoff. If not specified, the default value is
+        0.001. For more information, see the surface tally section in the theory
+        manual.
+    surface_grazing_ratio : float
+        Surface flux cosine substitution ratio. If not specified, the default
+        value is 0.5. For more information, see the surface tally section in the
+        theory manual.
     survival_biasing : bool
         Indicate whether survival biasing is to be used
     tabular_legendre : dict
@@ -399,6 +407,8 @@ def __init__(self, **kwargs):
         self._uniform_source_sampling = None
         self._seed = None
         self._stride = None
+        self._surface_grazing_cutoff = None
+        self._surface_grazing_ratio = None
         self._survival_biasing = None
         self._free_gas_threshold = None
 
@@ -710,6 +720,27 @@ def stride(self, stride: int):
         cv.check_greater_than('random number generator stride', stride, 0)
         self._stride = stride
 
+    @property
+    def surface_grazing_cutoff(self) -> float:
+        return self._surface_grazing_cutoff
+
+    @surface_grazing_cutoff.setter
+    def surface_grazing_cutoff(self, surface_grazing_cutoff: float):
+        cv.check_type('surface grazing cutoff', surface_grazing_cutoff, float)
+        cv.check_greater_than('surface grazing cutoff', surface_grazing_cutoff, 0.0)
+        cv.check_less_than('surface grazing cutoff', surface_grazing_cutoff, 1.0)
+        self._surface_grazing_cutoff = surface_grazing_cutoff
+
+    @property
+    def surface_grazing_ratio(self) -> float:
+        return self._surface_grazing_ratio
+
+    @surface_grazing_ratio.setter
+    def surface_grazing_ratio(self, surface_grazing_ratio: float):
+        cv.check_type('surface grazing ratio', surface_grazing_ratio, float)
+        cv.check_greater_than('surface grazing ratio', surface_grazing_ratio, 0.0)
+        self._surface_grazing_ratio = surface_grazing_ratio
+
     @property
     def survival_biasing(self) -> bool:
         return self._survival_biasing
@@ -1625,6 +1656,16 @@ def _create_stride_subelement(self, root):
             element = ET.SubElement(root, "stride")
             element.text = str(self._stride)
 
+    def _create_surface_grazing_cutoff_subelement(self, root):
+        if self._surface_grazing_cutoff is not None:
+            element = ET.SubElement(root, "surface_grazing_cutoff")
+            element.text = str(self._surface_grazing_cutoff)
+
+    def _create_surface_grazing_ratio_subelement(self, root):
+        if self._surface_grazing_ratio is not None:
+            element = ET.SubElement(root, "surface_grazing_ratio")
+            element.text = str(self._surface_grazing_ratio)
+
     def _create_survival_biasing_subelement(self, root):
         if self._survival_biasing is not None:
             element = ET.SubElement(root, "survival_biasing")
@@ -2128,6 +2169,16 @@ def _stride_from_xml_element(self, root):
         if text is not None:
             self.stride = int(text)
 
+    def _surface_grazing_cutoff_from_xml_element(self, root):
+        text = get_text(root, 'surface_grazing_cutoff')
+        if text is not None:
+            self.surface_grazing_cutoff = float(text)
+
+    def _surface_grazing_ratio_from_xml_element(self, root):
+        text = get_text(root, 'surface_grazing_ratio')
+        if text is not None:
+            self.surface_grazing_ratio = float(text)
+
     def _survival_biasing_from_xml_element(self, root):
         text = get_text(root, 'survival_biasing')
         if text is not None:
@@ -2435,6 +2486,8 @@ def to_xml_element(self, mesh_memo=None):
         self._create_ptables_subelement(element)
         self._create_seed_subelement(element)
         self._create_stride_subelement(element)
+        self._create_surface_grazing_cutoff_subelement(element)
+        self._create_surface_grazing_ratio_subelement(element)
         self._create_survival_biasing_subelement(element)
         self._create_cutoff_subelement(element)
         self._create_entropy_mesh_subelement(element, mesh_memo)
@@ -2549,6 +2602,8 @@ def from_xml_element(cls, elem, meshes=None):
         settings._ptables_from_xml_element(elem)
         settings._seed_from_xml_element(elem)
         settings._stride_from_xml_element(elem)
+        settings._surface_grazing_cutoff_from_xml_element(elem)
+        settings._surface_grazing_ratio_from_xml_element(elem)
         settings._survival_biasing_from_xml_element(elem)
         settings._cutoff_from_xml_element(elem)
         settings._entropy_mesh_from_xml_element(elem, meshes)

src/finalize.cpp

@@ -123,6 +123,8 @@ int openmc_finalize()
   settings::restart_run = false;
   settings::run_CE = true;
   settings::run_mode = RunMode::UNSET;
+  settings::surface_grazing_cutoff = 0.001;
+  settings::surface_grazing_ratio = 0.5;
   settings::solver_type = SolverType::MONTE_CARLO;
   settings::source_latest = false;
   settings::source_rejection_fraction = 0.05;

src/particle.cpp

@@ -302,6 +302,8 @@ void Particle::event_cross_surface()
   surface() = boundary().surface();
   n_coord() = boundary().coord_level();
 
+  const auto& surf {*model::surfaces[surface_index()].get()};
+
   if (boundary().lattice_translation()[0] != 0 ||
       boundary().lattice_translation()[1] != 0 ||
       boundary().lattice_translation()[2] != 0) {
@@ -312,15 +314,14 @@ void Particle::event_cross_surface()
     event() = TallyEvent::LATTICE;
   } else {
     // Particle crosses surface
-    const auto& surf {model::surfaces[surface_index()].get()};
     // If BC, add particle to surface source before crossing surface
-    if (surf->surf_source_ && surf->bc_) {
-      add_surf_source_to_bank(*this, *surf);
+    if (surf.surf_source_ && surf.bc_) {
+      add_surf_source_to_bank(*this, surf);
     }
-    this->cross_surface(*surf);
+    this->cross_surface(surf);
     // If no BC, add particle to surface source after crossing surface
-    if (surf->surf_source_ && !surf->bc_) {
-      add_surf_source_to_bank(*this, *surf);
+    if (surf.surf_source_ && !surf.bc_) {
+      add_surf_source_to_bank(*this, surf);
     }
     if (settings::weight_window_checkpoint_surface) {
       apply_weight_windows(*this);
@@ -329,7 +330,7 @@ void Particle::event_cross_surface()
   }
   // Score cell to cell partial currents
   if (!model::active_surface_tallies.empty()) {
-    score_surface_tally(*this, model::active_surface_tallies);
+    score_surface_tally(*this, model::active_surface_tallies, surf);
   }
 }
 
@@ -346,7 +347,7 @@ void Particle::event_collide()
   // pre-collision direction to figure out what mesh surfaces were crossed
 
   if (!model::active_meshsurf_tallies.empty())
-    score_surface_tally(*this, model::active_meshsurf_tallies);
+    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
 
   // Clear surface component
   surface() = SURFACE_NONE;
@@ -648,7 +649,7 @@ void Particle::cross_vacuum_bc(const Surface& surf)
     // physically moving the particle forward slightly
 
     r() += TINY_BIT * u();
-    score_surface_tally(*this, model::active_meshsurf_tallies);
+    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
   }
 
   // Score to global leakage tally
@@ -680,13 +681,13 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
   // with a mesh boundary
 
   if (!model::active_surface_tallies.empty()) {
-    score_surface_tally(*this, model::active_surface_tallies);
+    score_surface_tally(*this, model::active_surface_tallies, surf);
   }
 
   if (!model::active_meshsurf_tallies.empty()) {
     Position r {this->r()};
     this->r() -= TINY_BIT * u();
-    score_surface_tally(*this, model::active_meshsurf_tallies);
+    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
     this->r() = r;
   }
 
@@ -736,7 +737,7 @@ void Particle::cross_periodic_bc(
   if (!model::active_meshsurf_tallies.empty()) {
     Position r {this->r()};
     this->r() -= TINY_BIT * u();
-    score_surface_tally(*this, model::active_meshsurf_tallies);
+    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
     this->r() = r;
   }
 

src/settings.cpp

@@ -136,6 +136,8 @@ int64_t ssw_max_particles;
 int64_t ssw_max_files;
 int64_t ssw_cell_id {C_NONE};
 SSWCellType ssw_cell_type {SSWCellType::None};
+double surface_grazing_cutoff {0.001};
+double surface_grazing_ratio {0.5};
 TemperatureMethod temperature_method {TemperatureMethod::NEAREST};
 double temperature_tolerance {10.0};
 double temperature_default {293.6};
@@ -678,6 +680,14 @@ void read_settings_xml(pugi::xml_node root)
     free_gas_threshold = std::stod(get_node_value(root, "free_gas_threshold"));
   }
 
+  // Surface grazing
+  if (check_for_node(root, "surface_grazing_cutoff"))
+    surface_grazing_cutoff =
+      std::stod(get_node_value(root, "surface_grazing_cutoff"));
+  if (check_for_node(root, "surface_grazing_ratio"))
+    surface_grazing_ratio =
+      std::stod(get_node_value(root, "surface_grazing_ratio"));
+
   // Survival biasing
   if (check_for_node(root, "survival_biasing")) {
     survival_biasing = get_node_value_bool(root, "survival_biasing");

src/tallies/filter_surface.cpp

@@ -52,11 +52,7 @@ void SurfaceFilter::get_all_bins(
   auto search = map_.find(p.surface_index());
   if (search != map_.end()) {
     match.bins_.push_back(search->second);
-    if (p.surface() < 0) {
-      match.weights_.push_back(-1.0);
-    } else {
-      match.weights_.push_back(1.0);
-    }
+    match.weights_.push_back(1.0);
   }
 }