From dfd4b7e671c847cbe2aa6b1e1ed49b0215727990 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Sun, 3 Dec 2017 16:42:00 -0500
Subject: [PATCH] [Python/1D] Automatically increase domain width for free
 flame

If the domain is narrow with respect to the flame width, there can be
significant temperature gradients at the boundary, which lead to either
incorrect flame speeds or solver failures.

When the 'auto' option to FreeFlame.solve is specified, the solver will now
check the gradients at the ends of the domain after each steady-state solve and
increase the width if necessary.

Fixes #385
---
 interfaces/cython/cantera/onedim.py           | 61 +++++++++++++++++++
 interfaces/cython/cantera/test/test_onedim.py | 13 ++++
 2 files changed, 74 insertions(+)

diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py
index e40c786a9..355a4b163 100644
--- a/interfaces/cython/cantera/onedim.py
+++ b/interfaces/cython/cantera/onedim.py
@@ -469,6 +469,67 @@ class FreeFlame(FlameBase):
             self.set_profile(self.gas.species_name(n),
                              locs, [Y0[n], Y0[n], Yeq[n], Yeq[n]])
 
+    def solve(self, loglevel=1, refine_grid=True, auto=False):
+        """
+        Solve the problem.
+
+        :param loglevel:
+            integer flag controlling the amount of diagnostic output. Zero
+            suppresses all output, and 5 produces very verbose output.
+        :param refine_grid:
+            if True, enable grid refinement.
+        :param auto: if True, sequentially execute the different solution stages
+            and attempt to automatically recover from errors. Attempts to first
+            solve on the initial grid with energy enabled. If that does not
+            succeed, a fixed-temperature solution will be tried followed by
+            enabling the energy equation, and then with grid refinement enabled.
+            If non-default tolerances have been specified or multicomponent
+            transport is enabled, an additional solution using these options
+            will be calculated.
+        """
+        if not auto:
+            return super(FreeFlame, self).solve(loglevel, refine_grid, auto)
+
+        # Use a callback function to check that the domain is actually wide
+        # enough to contain the flame after each steady-state solve. If the user
+        # provided a callback, store this so it can called in addition to our
+        # callback, and restored at the end.
+        original_callback = self._steady_callback
+
+        class DomainTooNarrow(Exception): pass
+
+        def check_width(t):
+            T = self.T
+            x = self.grid
+            mRef = (T[-1] - T[0]) / (x[-1] - x[0])
+            mLeft = (T[1] - T[0]) / (x[1] - x[0]) / mRef
+            mRight = (T[-3] - T[-1]) / (x[-3] - x[-1]) / mRef
+
+            # The domain is considered too narrow if gradient at the left or
+            # right edge is significant, compared to the average gradient across
+            # the domain.
+            if mLeft > 0.05 or mRight > 0.05:
+                raise DomainTooNarrow()
+
+            if original_callback:
+                return original_callback(t)
+            else:
+                return 0.0
+
+        self.set_steady_callback(check_width)
+
+        for _ in range(12):
+            try:
+                return super(FreeFlame, self).solve(loglevel, refine_grid, auto)
+            except DomainTooNarrow:
+                self.flame.grid *= 2
+                if loglevel > 0:
+                    print('Expanding domain to accomodate flame thickness. '
+                          'New width: {} m'.format(
+                          self.flame.grid[-1] - self.flame.grid[0]))
+
+        self.set_steady_callback(original_callback)
+
     def get_flame_speed_reaction_sensitivities(self):
         """
         Compute the normalized sensitivities of the laminar flame speed
diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py
index c79540952..99e3df981 100644
--- a/interfaces/cython/cantera/test/test_onedim.py
+++ b/interfaces/cython/cantera/test/test_onedim.py
@@ -157,6 +157,19 @@ class TestFreeFlame(utilities.CanteraTest):
 
         self.assertEqual(self.sim.transport_model, 'Multi')
 
+    def test_auto_width(self):
+        Tin = 300
+        p = ct.one_atm
+        reactants = 'H2:0.65, O2:0.5, AR:2'
+        self.create_sim(p, Tin, reactants, width=0.0001)
+        self.sim.set_refine_criteria(ratio=3, slope=0.3, curve=0.2)
+        self.sim.solve(loglevel=0, refine_grid=True, auto=True)
+
+        self.gas.TPX = Tin, p, reactants
+        self.gas.equilibrate('HP')
+        Tad = self.gas.T
+        self.assertNear(Tad, self.sim.T[-1], 2e-2)
+
     def test_converge_adiabatic(self):
         # Test that the adiabatic flame temperature and species profiles
         # converge to the correct equilibrium values as the grid is refined