Figure of Merit implementation (#3363)

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

docs/source/methods/tallies.rst

@@ -387,6 +387,33 @@ of this is that the longer you run a simulation, the better you know your
 results. Therefore, by running a simulation long enough, it is possible to
 reduce the stochastic uncertainty to arbitrarily low levels.
 
+Figure of Merit
++++++++++++++++
+
+The figure of merit (FOM) is an indicator that accounts for both the statistical
+uncertainty and the execution time and represents how much information is
+obtained per unit time in the simulation. The FOM is defined as
+
+.. math::
+    :label: figure_of_merit
+
+    FOM = \frac{1}{r^2 t},
+
+where :math:`t` is the total execution time and :math:`r` is the relative error
+defined as
+
+.. math::
+    :label: relative_error
+
+    r = \frac{s_\bar{X}}{\bar{x}}.
+
+Based on this definition, one can see that a higher FOM is desirable. The FOM is
+useful as a comparative tool. For example, if a variance reduction technique is
+being applied to a simulation, the FOM with variance reduction can be compared
+to the FOM without variance reduction to ascertain whether the reduction in
+variance outweighs the potential increase in execution time (e.g., due to
+particle splitting).
+
 Confidence Intervals
 ++++++++++++++++++++
 

openmc/tallies.py

@@ -95,6 +95,10 @@ class Tally(IDManagerMixin):
         An array containing the sample mean for each bin
     std_dev : numpy.ndarray
         An array containing the sample standard deviation for each bin
+    figure_of_merit : numpy.ndarray
+        An array containing the figure of merit for each bin
+
+        .. versionadded:: 0.15.3
     derived : bool
         Whether or not the tally is derived from one or more other tallies
     sparse : bool
@@ -127,6 +131,7 @@ class Tally(IDManagerMixin):
         self._sum_sq = None
         self._mean = None
         self._std_dev = None
+        self._simulation_time = None
         self._with_batch_statistics = False
         self._derived = False
         self._sparse = False
@@ -385,6 +390,9 @@ class Tally(IDManagerMixin):
                 self._sum = sps.lil_matrix(self._sum.flatten(), self._sum.shape)
                 self._sum_sq = sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape)
 
+            # Read simulation time (needed for figure of merit)
+            self._simulation_time = f["runtime"]["simulation"][()]
+
         # Indicate that Tally results have been read
         self._results_read = True
 
@@ -462,6 +470,16 @@ class Tally(IDManagerMixin):
         else:
             return self._std_dev
 
+    @property
+    def figure_of_merit(self):
+        mean = self.mean
+        std_dev = self.std_dev
+        fom = np.zeros_like(mean)
+        nonzero = np.abs(mean) > 0
+        fom[nonzero] = 1.0 / (
+            (std_dev[nonzero] / mean[nonzero])**2 * self._simulation_time)
+        return fom
+
     @property
     def with_batch_statistics(self):
         return self._with_batch_statistics

tests/unit_tests/test_tallies.py

@@ -1,5 +1,5 @@
 import numpy as np
-
+import pytest
 import openmc
 
 
@@ -96,6 +96,22 @@ def test_tally_equivalence():
     assert tally_a == tally_b
 
 
+def test_figure_of_merit(sphere_model, run_in_tmpdir):
+    # Run model with a few simple tally scores
+    tally = openmc.Tally()
+    tally.scores = ['total', 'absorption', 'scatter']
+    sphere_model.tallies = [tally]
+    sp_path = sphere_model.run(apply_tally_results=True)
+
+    # Get execution time and relative error
+    with openmc.StatePoint(sp_path) as sp:
+        time = sp.runtime['simulation']
+    rel_err = tally.std_dev / tally.mean
+
+    # Check that figure of merit is calculated correctly
+    assert tally.figure_of_merit == pytest.approx(1 / (rel_err**2 * time))
+
+
 def test_tally_application(sphere_model, run_in_tmpdir):
     # Create a tally with most possible gizmos
     tally = openmc.Tally(name='test tally')