[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

@@ -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

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