[Examples] Revise preconditioned integration example

- Modify formatting for Sphinx Gallery
- Use a larger mechanism file from the example_data repo to better
  demonstrate the benefits of preconditioning

samples/python/reactors/preconditioned_integration.py

@@ -3,7 +3,7 @@ Acceleration of reactor integration using a sparse preconditioned solver
 ========================================================================
 
 Ideal gas, constant-pressure, adiabatic kinetics simulation that compares preconditioned
-and non-preconditioned integration of nDodecane.
+and non-preconditioned integration of n-hexane.
 
 Requires: cantera >= 3.0.0, matplotlib >= 2.0
 
@@ -12,26 +12,26 @@ Requires: cantera >= 3.0.0, matplotlib >= 2.0
 import cantera as ct
 import numpy as np
 import matplotlib.pyplot as plt
+plt.rcParams['figure.constrained_layout.use'] = True
 from timeit import default_timer
 
-def integrate_reactor(n_reactors=1, preconditioner=True):
-    """
-        Compare the integrations of a preconditioned reactor network and a
-        non-preconditioned reactor network. Increase the number of reactors in the
-        network with the keyword argument `n_reactors` to see greater differences in
-        performance.
-    """
-    # Use a reduced n-dodecane mechanism with PAH formation pathways
-    gas = ct.Solution('nDodecane_Reitz.yaml', 'nDodecane_IG')
-    # Create multiple reactors and set initial contents
+# %%
+# Simulation setup
+# ----------------
+#
+# Create a reactor network for simulating the constant pressure ignition of a
+# stoichiometric n-hexane/air mixture, with or without the use of the preconditioned
+# solver.
+def integrate_reactor(preconditioner=True):
+    # Use a detailed n-hexane mechanism with 1268 species
+    gas = ct.Solution('example_data/n-hexane-NUIG-2015.yaml')
     gas.TP = 1000, ct.one_atm
-    gas.set_equivalence_ratio(1, 'c12h26', 'n2:3.76, o2:1.0')
-    reactors = [ct.IdealGasConstPressureMoleReactor(gas) for n in range(n_reactors)]
+    gas.set_equivalence_ratio(1, 'NC6H14', 'N2:3.76, O2:1.0')
+    reactor = ct.IdealGasConstPressureMoleReactor(gas)
     # set volume for reactors
-    for r in reactors:
-        r.volume = 0.1
+    reactor.volume = 0.1
     # Create reactor network
-    sim = ct.ReactorNet(reactors)
+    sim = ct.ReactorNet([reactor])
     # Add preconditioner
     if preconditioner:
         sim.derivative_settings = {"skip-third-bodies":True, "skip-falloff":True}
@@ -40,10 +40,10 @@ def integrate_reactor(n_reactors=1, preconditioner=True):
     # Advance to steady state
     integ_time = default_timer()
     # solution array for state data
-    states = ct.SolutionArray(reactors[0].thermo, extra=['time'])
+    states = ct.SolutionArray(reactor.thermo, extra=['time'])
     # advance to steady state manually
     while (sim.time < 0.1):
-        states.append(reactors[0].thermo.state, time=sim.time)
+        states.append(reactor.thermo.state, time=sim.time)
         sim.step()
     integ_time = default_timer() - integ_time
     # Return time to integrate
@@ -53,35 +53,41 @@ def integrate_reactor(n_reactors=1, preconditioner=True):
         print(f"Non-preconditioned Integration Time: {integ_time:f}")
     # Get and output solver stats
     for key, value in sim.solver_stats.items():
-        print(key, value)
+        print(f"{key:>24s}: {value}")
     print("\n")
     # return some variables for plotting
-    return states.time, states.T, states('CO2').Y, states('C12H26').Y
+    return states.time, states.T, states('CO2').Y, states('NC6H14').Y
 
-if __name__ == "__main__":
-    # integrate both to steady state
-    timep, Tp, CO2p, C12H26p = integrate_reactor(preconditioner=True)
-    timenp, Tnp, CO2np, C12H26np  = integrate_reactor(preconditioner=False)
-    # plot selected state variables
-    fig, axs = plt.subplots(1, 3)
-    fig.tight_layout()
-    ax1, ax2, ax3 = axs
-    # temperature plot
-    ax1.set_xlabel("Time")
-    ax1.set_ylabel("Temperature")
-    ax1.plot(timenp, Tnp, linewidth=2)
-    ax1.plot(timep, Tp, linewidth=2, linestyle=":")
-    ax1.legend(["Normal", "Preconditioned"])
-    # CO2 plot
-    ax2.set_xlabel("Time")
-    ax2.set_ylabel("CO2")
-    ax2.plot(timenp, CO2np, linewidth=2)
-    ax2.plot(timep, CO2p, linewidth=2, linestyle=":")
-    ax2.legend(["Normal", "Preconditioned"])
-    # C12H26 plot
-    ax3.set_xlabel("Time")
-    ax3.set_ylabel("C12H26")
-    ax3.plot(timenp, C12H26np, linewidth=2)
-    ax3.plot(timep, C12H26p, linewidth=2, linestyle=":")
-    ax3.legend(["Normal", "Preconditioned"])
-    plt.show()
+# %%
+# Integrate with sparse, preconditioned solver
+# --------------------------------------------
+timep, Tp, CO2p, NC6H14p = integrate_reactor(preconditioner=True)
+
+# %%
+# Integrate with direct linear solver
+# -----------------------------------
+timenp, Tnp, CO2np, NC6H14np  = integrate_reactor(preconditioner=False)
+
+# %%
+# Plot selected state variables
+# -----------------------------
+fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(5, 8))
+# temperature plot
+ax1.set_xlabel("Time")
+ax1.set_ylabel("Temperature")
+ax1.plot(timenp, Tnp, linewidth=2)
+ax1.plot(timep, Tp, linewidth=2, linestyle=":")
+ax1.legend(["Normal", "Preconditioned"])
+# CO2 plot
+ax2.set_xlabel("Time")
+ax2.set_ylabel("CO2")
+ax2.plot(timenp, CO2np, linewidth=2)
+ax2.plot(timep, CO2p, linewidth=2, linestyle=":")
+ax2.legend(["Normal", "Preconditioned"])
+# C12H26 plot
+ax3.set_xlabel("Time")
+ax3.set_ylabel("NC6H14")
+ax3.plot(timenp, NC6H14np, linewidth=2)
+ax3.plot(timep, NC6H14p, linewidth=2, linestyle=":")
+ax3.legend(["Normal", "Preconditioned"])
+plt.show()