[Reactor] Add documentation for advance_to_steady_state

Also, update tutorials to use this feature.

doc/sphinx/reactors.rst

@@ -366,7 +366,7 @@ Time Integration
 Cantera provides an ODE solver for solving the stiff equations of reacting
 systems. If installed in combination with SUNDIALS, their optimized solver is
 used. Starting off the current state of the system, it can be advanced in time
-by two methods:
+by one of the following methods:
 
 - ``step()``: The step method computes the state of the system at the a priori
   unspecified time `t_{\rm new}`. The time `t_{\rm new}` is internally computed
@@ -382,6 +382,13 @@ by two methods:
   Internally, several ``step()`` calls are typically performed to reach the
   accurate state at time `t_{\rm new}`.
 
+- ``advance_to_steady_state(max_steps, residual_threshold, atol,
+  write_residuals)`` [Python interface only]: If the steady state solution of a
+  reactor network is of interest, this method can be used. Internally, the
+  steady state is approached by time stepping. The network is considered to be
+  at steady state if the feature-scaled residual of the state vector is below a
+  given threshold value (which by default is 10 times the time step rtol).
+
 The use of the ``advance`` method in a loop has the advantage that it produces
 results corresponding to a predefined time series. These are associated with a
 predefined memory consumption and well comparable between simulation runs with

interfaces/cython/cantera/examples/reactors/mix1.py

@@ -59,17 +59,8 @@ outlet = ct.Valve(mixer, downstream, K=10.0)
 sim = ct.ReactorNet([mixer])
 
 # Since the mixer is a reactor, we need to integrate in time to reach steady
-# state. A few residence times should be enough.
-print('{0:>14s} {1:>14s} {2:>14s}  {3:>14s}  {4:>14s}'.format(
-      't [s]', 'T [K]', 'h [J/kg]', 'P [Pa]', 'X_CH4'))
-
-t = 0.0
-for n in range(30):
-    tres = mixer.mass/(mfc1.mdot(t) + mfc2.mdot(t))
-    t += 0.5*tres
-    sim.advance(t)
-    print('{0:14.5g} {1:14.5g} {2:14.5g}  {3:14.5g}  {4:14.5g}'.format(
-          t, mixer.T, mixer.thermo.h, mixer.thermo.P, mixer.thermo['CH4'].X[0]))
+# state
+sim.advance_to_steady_state()
 
 # view the state of the gas in the mixer
 print(mixer.thermo.report())

interfaces/cython/cantera/examples/reactors/pfr.py

@@ -125,19 +125,8 @@ for n in range(n_steps):
     gas2.TDY = r2.thermo.TDY
     upstream.syncState()
     # integrate the reactor forward in time until steady state is reached
-    sim2.set_initial_time(0)  # forces reinitialization
-    time = 0
-    all_done = False
-    # determine steady state from H2 mole fraction
-    X_H2_previous = r2.thermo['H2'].X
-    while not all_done:
-        time += dt
-        sim2.advance(time)
-        if np.abs(r2.thermo['H2'].X - X_H2_previous) < 1.e-10:
-            # check whether surface coverages are in steady state.
-            all_done = True
-        else:
-            X_H2_previous = r2.thermo['H2'].X
+    sim2.reinitialize()
+    sim2.advance_to_steady_state()
     # compute velocity and transform into time
     u2[n] = mass_flow_rate2 / area / r2.thermo.density
     t_r2[n] = r2.mass / mass_flow_rate2  # residence time in this reactor

interfaces/cython/cantera/examples/reactors/surf_pfr.py

@@ -110,33 +110,8 @@ for n in range(NReactors):
     # Set the state of the reservoir to match that of the previous reactor
     gas.TDY = r.thermo.TDY
     upstream.syncState()
-
-    time = 0
-    all_done = False
-    sim.set_initial_time(0) # forces reinitialization
-    while not all_done:
-        time += dt
-        sim.advance(time)
-
-        if time > 10 * dt:
-            # check whether surface coverages are in steady state. This will be
-            # the case if the creation and destruction rates for a surface (but
-            # not gas) species are equal.
-            all_done = True
-
-            # Note: netProduction = creation - destruction. By supplying the
-            # surface object as an argument, only the values for the surface
-            # species are returned by these methods
-            sdot = surf.get_net_production_rates(surf)
-            cdot = surf.get_creation_rates(surf)
-            ddot = surf.get_destruction_rates(surf)
-
-            for ks in range(surf.n_species):
-                ratio = abs(sdot[ks]/(cdot[ks] + ddot[ks]))
-                if ratio > 1.0e-9:
-                    all_done = False
-                    break
-
+    sim.reinitialize()
+    sim.advance_to_steady_state()
     dist = n * rlen * 1.0e3   # distance in mm
 
     if not n % 10: