From 57e3ee8c1aef6b8dff73fd68c870e496c6e8748f Mon Sep 17 00:00:00 2001
From: Thomas Fiala <fiala@td.mw.tum.de>
Date: Mon, 2 Nov 2015 11:33:22 +0100
Subject: [PATCH] [Reactor] Implement the Cython function
 advance_to_steady_state

Resolves #95
Resolves #303
---
 interfaces/cython/cantera/reactor.pyx         | 61 +++++++++++++++++++
 .../cython/cantera/test/test_reactor.py       | 10 +++
 2 files changed, 71 insertions(+)

diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx
index 507f1ed3c..1c97a7279 100644
--- a/interfaces/cython/cantera/reactor.pyx
+++ b/interfaces/cython/cantera/reactor.pyx
@@ -996,6 +996,67 @@ cdef class ReactorNet:
         self.net.getState(&y[0])
         return y
 
+    def advance_to_steady_state(self, int max_steps=10000,
+                                double residual_threshold=0., double atol=0.,
+                                pybool return_residuals=False):
+        r"""
+        Advance the reactor network in time until steady state is reached.
+
+        The steady state is defined by requiring that the state of the system
+        only changes below a certain threshold. The residual is computed using
+        feature scaling:
+
+        .. math:: r = \frac{x(t + \Delta t) - x(t)}{\text{max}(x) + \text{atol}}
+
+        :param max_steps:
+            Maximum number of steps to be taken
+        :param residual_threshold:
+            Threshold below which the feature-scaled residual r should drop such
+            that the network is defines as steady state. By default,
+            residual_threshold is 10 times the solver rtol.
+        :param atol:
+            The smallest expected value of interest. Used for feature scaling.
+            By default, this atol is identical to the solver atol.
+        :param return_residuals:
+            If set to `True`, this function returns the residual time series
+            as a vector with length `max_steps`.
+
+        """
+        # get default tolerances:
+        if not atol:
+            atol = self.rtol
+        if not residual_threshold:
+            residual_threshold = 10. * self.rtol
+        if residual_threshold <= self.rtol:
+            raise Exception('Residual threshold (' + str(residual_threshold) +
+                            ') should be below solver rtol (' +
+                            str(self.rtol) + ')')
+        if return_residuals:
+            residuals = np.empty(max_steps)
+        # check if system is initialized
+        if not self.n_vars:
+            self.reinitialize()
+        max_state_values = self.get_state()  # denominator for feature scaling
+        for step in range(max_steps):
+            previous_state = self.get_state()
+            # take 10 steps (just to increase speed)
+            for n1 in range(10):
+                self.step()
+            state = self.get_state()
+            max_state_values = np.maximum(max_state_values, state)
+            # determine feature_scaled residual
+            residual = np.linalg.norm((state - previous_state)
+                / (max_state_values + atol))
+            if return_residuals:
+                residuals[step] = residual
+            if residual < residual_threshold:
+                break
+        if step == max_steps - 1:
+            raise Exception('Maximum number of steps reached before convergence'
+                            ' below maximum residual')
+        if return_residuals:
+            return residuals[:step + 1]
+
     def __reduce__(self):
         raise NotImplementedError('ReactorNet object is not picklable')
 
diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py
index 4f3ebdbe2..61db7919b 100644
--- a/interfaces/cython/cantera/test/test_reactor.py
+++ b/interfaces/cython/cantera/test/test_reactor.py
@@ -668,6 +668,16 @@ class TestWellStirredReactorIgnition(utilities.CanteraTest):
         t,T = self.integrate(100.0)
         self.assertTrue(T[-1] < 910) # mixture did not ignite
 
+    def test_steady_state(self):
+        self.setup(900.0, 10*ct.one_atm, 1.0, 20.0)
+        residuals = self.net.advance_to_steady_state(return_residuals=True)
+        # test if steady state is reached
+        self.assertTrue(residuals[-1] < 10. * self.net.rtol)
+        # regression test; no external basis for these results
+        self.assertNear(self.combustor.T, 2486.14, 1e-5)
+        self.assertNear(self.combustor.thermo['H2O'].Y[0], 0.103804, 1e-5)
+        self.assertNear(self.combustor.thermo['HO2'].Y[0], 7.71296e-06, 1e-5)
+
 
 class TestConstPressureReactor(utilities.CanteraTest):
     """