Adding measure of convergence for tight coupling (#1033)

This PR adds the capability to quantifiably measure the numerical convergence of solving nonlinear systems of equations when considering tight coupling. The current approach within ARMI just iterates on the system of equations using the numCoupledIterations setting. While indispensable for framework development, characterizing the numerical convergence of the nonlinear system of equations is essential to ensuring the mathematical validity of obtained solutions.

Co-authored-by: jakehader <jhader@terrapower.com>

armi/bookkeeping/db/databaseInterface.py

@@ -150,7 +150,7 @@ def interactEveryNode(self, cycle, node):
         - if tight coupling is enabled, the DB will be written in operator.py::Operator::_timeNodeLoop
           via writeDBEveryNode
         """
-        if self.o.cs["numCoupledIterations"]:
+        if self.o.cs["tightCoupling"]:
             # h5 cant handle overwriting so we skip here and write once the tight coupling loop has completed
             return
         self.writeDBEveryNode(cycle, node)

armi/bookkeeping/db/tests/test_databaseInterface.py

@@ -99,8 +99,8 @@ def tearDown(self):
         self.td.__exit__(None, None, None)
 
     def test_interactEveryNodeReturn(self):
-        """test that the DB is NOT written to if cs["numCoupledIterations"] != 0"""
-        self.o.cs["numCoupledIterations"] = 1
+        """test that the DB is NOT written to if cs["tightCoupling"] = True"""
+        self.o.cs["tightCoupling"] = True
         self.dbi.interactEveryNode(0, 0)
         self.assertFalse(self.dbi.database.hasTimeStep(0, 0))
 
@@ -118,11 +118,11 @@ def test_distributable(self):
         self.dbi.interactDistributeState()
         self.assertEqual(self.dbi.distributable(), 4)
 
-    def test_timeNodeLoop_numCoupledIterations(self):
+    def test_timeNodeLoop_tightCoupling(self):
         """test that database is written out after the coupling loop has completed"""
         # clear out interfaces (no need to run physics) but leave database
         self.o.interfaces = [self.dbi]
-        self.o.cs["numCoupledIterations"] = 1
+        self.o.cs["tightCoupling"] = True
         self.assertFalse(self.dbi._db.hasTimeStep(0, 0))
         self.o._timeNodeLoop(0, 0)
         self.assertTrue(self.dbi._db.hasTimeStep(0, 0))

armi/bookkeeping/report/newReportUtils.py

@@ -295,7 +295,7 @@ def _setGeneralSimulationData(core, cs, coreDesignTable):
     coreDesignTable.addRow(["  ", ""])
     coreDesignTable.addRow(["Full Core Model", "{}".format(core.isFullCore)])
     coreDesignTable.addRow(
-        ["Loose Physics Coupling Enabled", "{}".format(bool(cs["looseCoupling"]))]
+        ["Tight Physics Coupling Enabled", "{}".format(bool(cs["tightCoupling"]))]
     )
     coreDesignTable.addRow(["Lattice Physics Enabled for", "{}".format(cs["genXS"])])
     coreDesignTable.addRow(

armi/bookkeeping/report/reportingUtils.py

@@ -41,7 +41,6 @@
 from armi.utils import plotting
 from armi.utils import textProcessors
 from armi.utils import units
-from armi.utils.mathematics import findClosest
 
 
 # Set to prevent the image and text from being too small to read.
@@ -238,6 +237,15 @@ def getInterfaceStackSummary(o):
     return text
 
 
+def writeTightCouplingConvergenceSummary(convergenceSummary):
+    runLog.info("Tight Coupling Convergence Summary: Norm Type = Inf")
+    runLog.info(
+        tabulate.tabulate(
+            convergenceSummary, headers="keys", showindex=True, tablefmt="armi"
+        )
+    )
+
+
 def writeAssemblyMassSummary(r):
     r"""Print out things like Assembly weights to the runLog.
 
@@ -766,8 +774,8 @@ def _setGeneralSimulationData(core, cs, coreDesignTable):
         "Full Core Model", "{}".format(core.isFullCore), coreDesignTable, report.DESIGN
     )
     report.setData(
-        "Loose Physics Coupling Enabled",
-        "{}".format(bool(cs["looseCoupling"])),
+        "Tight Physics Coupling Enabled",
+        "{}".format(bool(cs["tightCoupling"])),
         coreDesignTable,
         report.DESIGN,
     )

armi/interfaces.py

@@ -31,9 +31,13 @@
 from typing import List
 from typing import Dict
 
+import numpy
+from numpy.linalg import norm
+
 from armi import getPluginManagerOrFail, settings, utils
 from armi.utils import textProcessors
 from armi.reactor import parameters
+from armi import runLog
 
 
 class STACK_ORDER:  # pylint: disable=invalid-name, too-few-public-methods
@@ -77,6 +81,156 @@ class STACK_ORDER:  # pylint: disable=invalid-name, too-few-public-methods
     POSTPROCESSING = BOOKKEEPING + 1
 
 
+class TightCoupler:
+    """
+    Data structure that defines tight coupling attributes that are implemented
+    within an Interface and called upon when ``interactCoupled`` is called.
+
+    Parameters
+    ----------
+    param : str
+        The name of a parameter defined in the ARMI Reactor model.
+
+    tolerance : float
+        Defines the allowable error, epsilon, between the current previous
+        parameter value(s) to determine if the selected coupling parameter has
+        been converged.
+
+    maxIters : int
+        Maximum number of tight coupling iterations allowed
+    """
+
+    _SUPPORTED_TYPES = [float, int, list, numpy.ndarray]
+
+    def __init__(self, param, tolerance, maxIters):
+        self.parameter = param
+        self.tolerance = tolerance
+        self.maxIters = maxIters
+        self._numIters = 0
+        self._previousIterationValue = None
+        self.eps = numpy.inf
+
+    def __repr__(self):
+        return f"<{self.__class__.__name__}, Parameter: {self.parameter}, Convergence Criteria: {self.tolerance}, Maximum Coupled Iterations: {self.maxIters}>"
+
+    def storePreviousIterationValue(self, val: _SUPPORTED_TYPES):
+        """
+        Stores the previous iteration value of the given parameter.
+
+        Parameters
+        ----------
+        val : _SUPPORTED_TYPES
+            the value to store. Is commonly equal to interface.getTightCouplingValue()
+
+        Raises
+        ------
+        TypeError
+            Checks the type of the val against ``_SUPPORTED_TYPES`` before storing.
+            If invalid, a TypeError is raised.
+        """
+        if type(val) not in self._SUPPORTED_TYPES:
+            raise TypeError(
+                f"{val} supplied has type {type(val)} which is not supported in {self}. "
+                f"Supported types: {self._SUPPORTED_TYPES}"
+            )
+        self._previousIterationValue = val
+
+    def isConverged(self, val: _SUPPORTED_TYPES) -> bool:
+        """
+        Return boolean indicating if the convergence criteria between the current and previous iteration values are met.
+
+        Parameters
+        ----------
+        val : _SUPPORTED_TYPES
+            the most recent value for computing convergence critera. Is commonly equal to interface.getTightCouplingValue()
+
+        Returns
+        -------
+        boolean
+            True (False) interface is (not) converged
+
+        Notes
+        -----
+        - On convergence, this class is automatically reset to its initial condition to avoid retaining
+          or holding a stale state. Calling this method will increment a counter that when exceeded will
+          clear the state. A warning will be reported if the state is cleared prior to the convergence
+          criteria being met.
+        - For computing convergence of arrays, only up to 2D is allowed. 3D arrays would arise from considering
+          component level parameters. However, converging on component level parameters is not supported at this time.
+
+        Raises
+        ------
+        ValueError
+            If the previous iteration value has not been assigned. The ``storePreviousIterationValue`` method
+            must be called first.
+        RuntimeError
+            Only support calculating norms for up to 2D arrays.
+        """
+        if self._previousIterationValue is None:
+            raise ValueError(
+                f"Cannot check convergence of {self} with no previous iteration value set. "
+                f"Set using `storePreviousIterationValue` first."
+            )
+
+        previous = self._previousIterationValue
+
+        # calculate convergence of val and previous
+        if isinstance(val, (int, float)):
+            self.eps = abs(val - previous)
+        else:
+            dim = self.getListDimension(val)
+            if dim == 1:  # 1D array
+                self.eps = norm(numpy.subtract(val, previous), ord=2)
+            elif dim == 2:  # 2D array
+                epsVec = []
+                for old, new in zip(previous, val):
+                    epsVec.append(norm(numpy.subtract(old, new), ord=2))
+                self.eps = norm(epsVec, ord=numpy.inf)
+            else:
+                raise RuntimeError(
+                    "Currently only support up to 2D arrays for calculating convergence of arrays."
+                )
+
+        # Check if convergence is satisfied. If so, or if reached max number of iters, then
+        # reset the number of iterations
+        converged = self.eps < self.tolerance
+        if converged:
+            self._numIters = 0
+        else:
+            self._numIters += 1
+            if self._numIters == self.maxIters:
+                runLog.warning(
+                    f"Maximum number of iterations for {self.parameter} reached without convergence!"
+                    f"Prescribed convergence criteria is {self.tolerance}."
+                )
+                self._numIters = 0
+
+        return converged
+
+    @staticmethod
+    def getListDimension(listToCheck: list, dim: int = 1) -> int:
+        """return the dimension of a python list
+
+        Parameters
+        ----------
+        listToCheck: list
+            the supplied python list to have its dimension returned
+        dim: int, optional
+            the dimension of the list
+
+        Returns
+        -------
+        dim, int
+            the dimension of the list. Typically 1, 2, or 3 but can be arbitrary order, N.
+        """
+        for v in listToCheck:
+            if isinstance(v, list):
+                dim += 1
+                dim = TightCoupler.getListDimension(v, dim)
+            break
+        return dim
+
+
 class Interface:
     """
     The eponymous Interface between the ARMI Reactor model and modules that operate upon it.
@@ -153,6 +307,7 @@ def __init__(self, r, cs):
         self.cs = settings.getMasterCs() if cs is None else cs
         self.r = r
         self.o = r.o if r else None
+        self.coupler = _setTightCouplerByInterfaceFunction(self, cs)
 
     def __repr__(self):
         return "<Interface {0}>".format(self.name)
@@ -307,7 +462,11 @@ def interactEveryNode(self, cycle, node):
         pass
 
     def interactCoupled(self, iteration):
-        """Called repeatedly at each time node/subcycle when tight physics couping is active."""
+        """Called repeatedly at each time node/subcycle when tight physics coupling is active."""
+        pass
+
+    def getTightCouplingValue(self):
+        """Abstract method to retrieve the value in which tight coupling will converge on."""
         pass
 
     def interactError(self):
@@ -528,6 +687,35 @@ def apply(self, reactor):
         raise NotImplementedError()
 
 
+def _setTightCouplerByInterfaceFunction(interfaceClass, cs):
+    """
+    Return an instance of a ``TightCoupler`` class or ``None``.
+
+    Parameters
+    ----------
+    interfaceClass : Interface
+        Interface class that a ``TightCoupler`` object will be added to.
+
+    cs : Settings
+        Case settings that are parsed to determine if tight coupling is enabled
+        globally and if both a target parameter and convergence criteria defined.
+    """
+    # No tight coupling if there is no function for the Interface defined.
+    if interfaceClass.function is None:
+        return None
+
+    if not cs["tightCoupling"] or (
+        interfaceClass.function not in cs["tightCouplingSettings"]
+    ):
+        return None
+
+    parameter = cs["tightCouplingSettings"][interfaceClass.function]["parameter"]
+    tolerance = cs["tightCouplingSettings"][interfaceClass.function]["convergence"]
+    maxIters = cs["tightCouplingMaxNumIters"]
+
+    return TightCoupler(parameter, tolerance, maxIters)
+
+
 def getActiveInterfaceInfo(cs):
     """
     Return a list containing information for all of the Interface classes that are present.

armi/operators/operator.py

@@ -30,6 +30,7 @@
 import re
 import shutil
 import time
+import collections
 
 from armi import context
 from armi import interfaces
@@ -134,6 +135,7 @@ def __init__(self, cs):
         self._maxBurnSteps = None
         self._powerFractions = None
         self._availabilityFactors = None
+        self._convergenceSummary = None
 
         # Create the welcome headers for the case (case, input, machine, and some basic reactor information)
         reportingUtils.writeWelcomeHeaders(self, cs)
@@ -377,10 +379,13 @@ def _timeNodeLoop(self, cycle, timeNode):
         self.r.p.timeNode = timeNode
         self.interactAllEveryNode(cycle, timeNode)
         # perform tight coupling if requested
-        if self.cs["numCoupledIterations"]:
-            for coupledIteration in range(self.cs["numCoupledIterations"]):
+        if self.couplingIsActive():
+            self._convergenceSummary = collections.defaultdict(list)
+            for coupledIteration in range(self.cs["tightCouplingMaxNumIters"]):
                 self.r.core.p.coupledIteration = coupledIteration + 1
-                self.interactAllCoupled(coupledIteration)
+                converged = self.interactAllCoupled(coupledIteration)
+                if converged:
+                    break
             # database has not yet been written, so we need to write it.
             dbi = self.getInterface("database")
             dbi.writeDBEveryNode(cycle, timeNode)
@@ -617,8 +622,42 @@ def interactAllCoupled(self, coupledIteration):
         ARMI supports tight and loose coupling.
         """
         activeInterfaces = [ii for ii in self.interfaces if ii.enabled()]
+
+        # Store the previous iteration values before calling interactAllCoupled
+        # for each interface.
+        for interface in activeInterfaces:
+            if interface.coupler is not None:
+                interface.coupler.storePreviousIterationValue(
+                    interface.getTightCouplingValue()
+                )
+
         self._interactAll("Coupled", activeInterfaces, coupledIteration)
 
+        return self._checkTightCouplingConvergence(activeInterfaces)
+
+    def _checkTightCouplingConvergence(self, activeInterfaces: list):
+        """check if interfaces are converged
+
+        Parameters
+        ----------
+        activeInterfaces : list
+            the list of active interfaces on the operator
+
+        Notes
+        -----
+        This is split off from self.interactAllCoupled to accomodate testing"""
+        # Summarize the coupled results and the convergence status.
+        converged = []
+        for interface in activeInterfaces:
+            coupler = interface.coupler
+            if coupler is not None:
+                key = f"{interface.name}: {coupler.parameter}"
+                converged.append(coupler.isConverged(interface.getTightCouplingValue()))
+                self._convergenceSummary[key].append(coupler.eps)
+
+        reportingUtils.writeTightCouplingConvergenceSummary(self._convergenceSummary)
+        return all(converged)
+
     def interactAllError(self):
         """Interact when an error is raised by any other interface. Provides a wrap-up option on the way to a crash."""
         for i in self.interfaces:
@@ -648,6 +687,7 @@ def createInterfaces(self):
         armi.interfaces.getActiveInterfaceInfo : Collects the interface classes from relevant
             packages.
         """
+        runLog.header("=========== Creating Interfaces ===========")
         interfaceList = interfaces.getActiveInterfaceInfo(self.cs)
 
         for klass, kwargs in interfaceList:
@@ -1064,4 +1104,4 @@ def setStateToDefault(cs):
 
     def couplingIsActive(self):
         """True if any kind of physics coupling is active."""
-        return self.cs["looseCoupling"] or self.cs["numCoupledIterations"] > 0
+        return self.cs["tightCoupling"]