Use energy-dependent fission yields when constructing depletion matrix

docs/source/pythonapi/deplete.rst

@@ -6,7 +6,7 @@
 
 .. module:: openmc.deplete
 
-Several functions are provided that implement different time-integration
+Several classes are provided that implement different time-integration
 algorithms for depletion calculations, which are described in detail in Colin
 Josey's thesis, `Development and analysis of high order neutron
 transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
@@ -25,7 +25,7 @@ transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
     SICELIIntegrator
     SILEQIIntegrator
 
-Each of these functions expects a "transport operator" to be passed. An operator
+Each of these classes expects a "transport operator" to be passed. An operator
 specific to OpenMC is available using the following class:
 
 .. autosummary::
@@ -65,6 +65,8 @@ for a depletion chain:
    DecayTuple
    Nuclide
    ReactionTuple
+   FissionYieldDistribution
+   FissionYield
 
 The following classes are used during a depletion simulation and store auxiliary
 data, such as number densities and reaction rates for each material.
@@ -75,13 +77,25 @@ data, such as number densities and reaction rates for each material.
    :template: myclass.rst
 
    AtomNumber
-   ChainFissionHelper
-   DirectReactionRateHelper
    OperatorResult
    ReactionRates
    Results
    ResultsList
 
+The following classes are used to help the :class:`openmc.deplete.Operator`
+compute quantities like effective fission yields, reaction rates, and
+total system energy.
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myclass.rst
+
+   helpers.AveragedFissionYieldHelper
+   helpers.ChainFissionHelper
+   helpers.ConstantFissionYieldHelper
+   helpers.DirectReactionRateHelper
+   helpers.FissionYieldCutoffHelper
 
 The following classes are abstract classes that can be used to extend the
 :mod:`openmc.deplete` capabilities:
@@ -91,8 +105,10 @@ The following classes are abstract classes that can be used to extend the
    :nosignatures:
    :template: myclass.rst
 
-   ReactionRateHelper
    EnergyHelper
+   FissionYieldHelper
+   ReactionRateHelper
+   TalliedFissionYieldHelper
    TransportOperator
 
 Custom integrators can be developed by subclassing from the following abstract

openmc/deplete/_matrix_funcs.py

@@ -1,77 +1,78 @@
 """Functions to form the special matrix for depletion"""
 
 
-def celi_f1(chain, rates):
-    return (5 / 12 * chain.form_matrix(rates[0])
-            + 1 / 12 * chain.form_matrix(rates[1]))
+def celi_f1(chain, rates, fission_yields=None):
+    return (5 / 12 * chain.form_matrix(rates[0], fission_yields)
+            + 1 / 12 * chain.form_matrix(rates[1], fission_yields))
 
 
-def celi_f2(chain, rates):
-    return (1 / 12 * chain.form_matrix(rates[0])
-            + 5 / 12 * chain.form_matrix(rates[1]))
+def celi_f2(chain, rates, fission_yields=None):
+    return (1 / 12 * chain.form_matrix(rates[0], fission_yields)
+            + 5 / 12 * chain.form_matrix(rates[1], fission_yields))
 
 
-def cf4_f1(chain, rates):
-    return 1 / 2 * chain.form_matrix(rates)
+def cf4_f1(chain, rates, fission_yields=None):
+    return 1 / 2 * chain.form_matrix(rates, fission_yields)
 
 
-def cf4_f2(chain, rates):
-    return -1 / 2 * chain.form_matrix(rates[0]) + chain.form_matrix(rates[1])
+def cf4_f2(chain, rates, fission_yields=None):
+    return (-1 / 2 * chain.form_matrix(rates[0], fission_yields)
+            + chain.form_matrix(rates[1], fission_yields))
 
 
-def cf4_f3(chain, rates):
-    return (1 / 4 * chain.form_matrix(rates[0])
-            + 1 / 6 * chain.form_matrix(rates[1])
-            + 1 / 6 * chain.form_matrix(rates[2])
-            - 1 / 12 * chain.form_matrix(rates[3]))
+def cf4_f3(chain, rates, fission_yields=None):
+    return (1 / 4 * chain.form_matrix(rates[0], fission_yields)
+            + 1 / 6 * chain.form_matrix(rates[1], fission_yields)
+            + 1 / 6 * chain.form_matrix(rates[2], fission_yields)
+            - 1 / 12 * chain.form_matrix(rates[3], fission_yields))
 
 
-def cf4_f4(chain, rates):
-    return (-1 / 12 * chain.form_matrix(rates[0])
-            + 1 / 6 * chain.form_matrix(rates[1])
-            + 1 / 6 * chain.form_matrix(rates[2])
-            + 1 / 4 * chain.form_matrix(rates[3]))
+def cf4_f4(chain, rates, fission_yields=None):
+    return (-1 / 12 * chain.form_matrix(rates[0], fission_yields)
+            + 1 / 6 * chain.form_matrix(rates[1], fission_yields)
+            + 1 / 6 * chain.form_matrix(rates[2], fission_yields)
+            + 1 / 4 * chain.form_matrix(rates[3], fission_yields))
 
 
-def rk4_f1(chain, rates):
-    return 1 / 2 * chain.form_matrix(rates)
+def rk4_f1(chain, rates, fission_yields=None):
+    return 1 / 2 * chain.form_matrix(rates, fission_yields)
 
 
-def rk4_f4(chain, rates):
-    return (1 / 6 * chain.form_matrix(rates[0])
-            + 1 / 3 * chain.form_matrix(rates[1])
-            + 1 / 3 * chain.form_matrix(rates[2])
-            + 1 / 6 * chain.form_matrix(rates[3]))
+def rk4_f4(chain, rates, fission_yields=None):
+    return (1 / 6 * chain.form_matrix(rates[0], fission_yields)
+            + 1 / 3 * chain.form_matrix(rates[1], fission_yields)
+            + 1 / 3 * chain.form_matrix(rates[2], fission_yields)
+            + 1 / 6 * chain.form_matrix(rates[3], fission_yields))
 
 
-def leqi_f1(chain, inputs):
-    f1 = chain.form_matrix(inputs[0])
-    f2 = chain.form_matrix(inputs[1])
+def leqi_f1(chain, inputs, fission_yields):
+    f1 = chain.form_matrix(inputs[0], fission_yields)
+    f2 = chain.form_matrix(inputs[1], fission_yields)
     dt_l, dt = inputs[2], inputs[3]
     return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2
 
 
-def leqi_f2(chain, inputs):
-    f1 = chain.form_matrix(inputs[0])
-    f2 = chain.form_matrix(inputs[1])
+def leqi_f2(chain, inputs, fission_yields=None):
+    f1 = chain.form_matrix(inputs[0], fission_yields)
+    f2 = chain.form_matrix(inputs[1], fission_yields)
     dt_l, dt = inputs[2], inputs[3]
     return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2
 
 
-def leqi_f3(chain, inputs):
-    f1 = chain.form_matrix(inputs[0])
-    f2 = chain.form_matrix(inputs[1])
-    f3 = chain.form_matrix(inputs[2])
+def leqi_f3(chain, inputs, fission_yields=None):
+    f1 = chain.form_matrix(inputs[0], fission_yields)
+    f2 = chain.form_matrix(inputs[1], fission_yields)
+    f3 = chain.form_matrix(inputs[2], fission_yields)
     dt_l, dt = inputs[3], inputs[4]
     return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1
             + (dt ** 2 + 6 * dt * dt_l + 5 * dt_l ** 2)
             / (12 * dt_l * (dt + dt_l)) * f2 + dt_l / (12 * (dt + dt_l)) * f3)
 
 
-def leqi_f4(chain, inputs):
-    f1 = chain.form_matrix(inputs[0])
-    f2 = chain.form_matrix(inputs[1])
-    f3 = chain.form_matrix(inputs[2])
+def leqi_f4(chain, inputs, fission_yields=None):
+    f1 = chain.form_matrix(inputs[0], fission_yields)
+    f2 = chain.form_matrix(inputs[1], fission_yields)
+    f3 = chain.form_matrix(inputs[2], fission_yields)
     dt_l, dt = inputs[3], inputs[4]
     return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1
             + (dt ** 2 + 2 * dt * dt_l + dt_l ** 2)