Update example formatting for better sphinx-gallery rendering

doc/sphinx/ctutils.py

@@ -6,7 +6,6 @@ class ResetArgv:
     wants_plot = {
         "adiabatic.py",
         "premixed_counterflow_twin_flame.py",
-        "piston.py",
         "reactor1.py",
         "reactor2.py",
         "sensitivity1.py",

samples/python/kinetics/blowers_masel.py

@@ -5,10 +5,6 @@
 A simple example to demonstrate the difference between Blowers-Masel
 reaction and elementary reaction.
 
-The first two reactions have the same reaction equations with Arrhenius and
-Blowers-Masel rate parameters, respectively. The Blowers-Masel parameters are the same
-as the Arrhenius parameters with an additional value, bond energy.
-
 First we show that the forward rate constants of the first 2 different reactions are
 different because of the different rate expression, then we print the forward rate
 constants for reaction 2 and reaction 3 to show that even two reactions that have the
@@ -27,19 +23,24 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
-#Create an elementary reaction O+H2<=>H+OH
-r1 = ct.Reaction({"O":1, "H2":1}, {"H":1, "OH":1},
-                 ct.ArrheniusRate(3.87e1, 2.7, 6260*1000*4.184))
+# %%
+# Create an elementary reaction:
+r1 = ct.Reaction(equation="O + H2 <=> H + OH",
+                 rate=ct.ArrheniusRate(3.87e1, 2.7, 6260*1000*4.184))
 
-#Create a Blowers-Masel reaction O+H2<=>H+OH
-r2 = ct.Reaction({"O":1, "H2":1}, {"H":1, "OH":1},
-                 ct.BlowersMaselRate(3.87e1, 2.7, 6260*1000*4.184, 1e9))
+# %%
+# Create a Blowers-Masel reaction:
+r2 = ct.Reaction(equation="O + H2 <=> H + OH",
+                 rate=ct.BlowersMaselRate(3.87e1, 2.7, 6260*1000*4.184, 1e9))
 
-#Create a Blowers-Masel reaction with same parameters with r2
-#reaction equation is H+CH4<=>CH3+H2
-r3 = ct.Reaction({"H":1, "CH4":1}, {"CH3":1, "H2":1},
-                 ct.BlowersMaselRate(3.87e1, 2.7, 6260*1000*4.184, 1e9))
+# %%
+# Create a Blowers-Masel reaction with same parameters as ``r2``, but for different
+# reactants and products
+r3 = ct.Reaction(equation="H + CH4 <=> CH3 + H2",
+                 rate=ct.BlowersMaselRate(3.87e1, 2.7, 6260*1000*4.184, 1e9))
 
+# %%
+# Evaluate these reaction rates as part of a gas mixture
 gas = ct.Solution(thermo='ideal-gas', kinetics='gas',
                   species=ct.Solution('gri30.yaml').species(), reactions=[r1, r2, r3])
 
@@ -49,33 +50,47 @@
 r2_rc = gas.forward_rate_constants[1]
 r3_rc = gas.forward_rate_constants[2]
 
-print("The first and second reactions have same reaction equation,"
-      " but they have different reaction types, so the forward rate"
-      " constant of the first reaction is {0:.3f} kmol/(m^3.s),"
-      " the forward rate constant of the second reaction is {1:.3f} kmol/(m^3.s).".format(r1_rc, r2_rc))
-
-print("The rate parameters of second and the third reactions are same,"
-      " but the forward rate constant of second reaction is {0:.3f} kmol/(m^3.s),"
-      " the forward rate constant of the third reaction is"
-      " {1:.3f} kmol/(m^3.s).".format(r2_rc, r3_rc))
+# %%
+# The first two reactions have the same reaction equations with Arrhenius and
+# Blowers-Masel rate parameters, respectively. The Blowers-Masel parameters are the same
+# as the Arrhenius parameters with an additional value, bond energy. The forward rate
+# constants for these two reactions are different because of the difference in the rate
+# expression:
+print(f"{r1_rc=:.3f} kmol/m³/s")
+print(f"{r2_rc=:.3f} kmol/m³/s")
+
+# %%
+# The rate parameters of second and the third reactions are same, but their reactants
+# and products are different. Since the enthalpy of reaction is used in the
+# Blowers-Masel rate expression, the rate constants are different:
+print(f"{r2_rc=:.3f} kmol/m³/s")
+print(f"{r3_rc=:.3f} kmol/m³/s")
+
+# %%
+# Effect of temperature on reaction rates
+# ---------------------------------------
+# Compare the reaction forward rate constant change of Blowers-Masel reaction and
+# elementary reaction with respect to the temperature.
 
-# Comparing the reaction forward rate constant change of
-# Blowers-Masel reaction and elementary reaction with
-# respect to the temperature.
 r1_kf = []
 r2_kf = []
 T_range = np.arange(300, 3500, 100)
 for temp in T_range:
     gas.TP = temp, ct.one_atm
     r1_kf.append(gas.forward_rate_constants[0])
     r2_kf.append(gas.forward_rate_constants[1])
-plt.plot(T_range, r1_kf, label='Reaction 1 (Elementary)')
-plt.plot(T_range, r2_kf, label='Reaction 2 (Blowers-Masel)')
-plt.xlabel("Temperature(K)")
-plt.ylabel(r"Forward Rate Constant (kmol/(m$^3\cdot$ s))")
-plt.title("Comparison of $k_f$ vs. Temperature For Reaction 1 and 2",y=1.1)
-plt.legend()
+
+fig, ax = plt.subplots()
+ax.plot(T_range, r1_kf, label='Reaction 1 (Elementary)')
+ax.plot(T_range, r2_kf, label='Reaction 2 (Blowers-Masel)')
+ax.set(xlabel="Temperature(K)", ylabel=r"Forward Rate Constant (kmol/(m$^3\cdot$ s))")
+ax.set_title("Comparison of $k_f$ vs. Temperature For Reaction 1 and 2")
+ax.legend()
 plt.savefig("kf_to_T.png")
+plt.show()
+
+# %%
+# Plot the activation energy change of reaction 2 with respect to the enthalpy change
 
 # This is the function to change the enthalpy of a species
 # so that the enthalpy change of reactions involving this species can be changed
@@ -98,8 +113,6 @@ def change_species_enthalpy(gas, species_name, dH):
     gas.modify_species(index, species)
     return gas.delta_enthalpy[1]
 
-# Plot the activation energy change of reaction 2 with respect to the
-# enthalpy change
 E0 = gas.reaction(1).rate.activation_energy
 upper_limit_enthalpy = 5 * E0
 lower_limit_enthalpy = -5 * E0
@@ -111,10 +124,9 @@ def change_species_enthalpy(gas, species_name, dH):
     gas.reaction(1).rate.delta_enthalpy = delta_enthalpy
     Ea_list.append(gas.reaction(1).rate.activation_energy)
 
-plt.figure()
-plt.plot(deltaH_list, Ea_list)
-plt.xlabel("Enthalpy Change (J/kmol)")
-plt.ylabel("Activation Energy Change (J/kmol)")
-plt.title(r"$E_a$ vs. $\Delta H$ For O+H2<=>H+OH", y=1.1)
+fig, ax = plt.subplots()
+ax.plot(deltaH_list, Ea_list)
+ax.set(xlabel="Enthalpy Change (J/kmol)", ylabel="Activation Energy Change (J/kmol)")
+ax.set_title(r"$E_a$ vs. $\Delta H$ For $\mathrm{ O+H_2 \rightleftharpoons H+OH }$")
 plt.savefig("Ea_to_H.png")
 plt.show()

samples/python/kinetics/custom_reactions.py

@@ -25,8 +25,13 @@
 species = gas0.species()
 reactions = gas0.reactions()
 
-# construct reactions based on CustomRate: replace two reactions with equivalent
-# custom versions (uses Func1 functors defined by Python lambda functions)
+# %%
+# Define custom reaction rates
+# ----------------------------
+
+# %%
+# Construct reactions based on `CustomRate`: replace two reactions with equivalent
+# custom versions (uses `Func1` functors defined by Python lambda functions)
 reactions[2] = ct.Reaction(
     equation='H2 + O <=> H + OH',
     rate=lambda T: 38.7 * T**2.7 * exp(-3150.1542797022735/T))
@@ -37,8 +42,9 @@
 gas1a = ct.Solution(thermo='ideal-gas', kinetics='gas',
                     species=species, reactions=reactions)
 
-# construct reactions based on CustomRate: replace two reactions with equivalent
-# custom versions (uses Func1 functors defined by C++ class Arrhenius1)
+# %%
+# construct reactions based on `CustomRate`: replace two reactions with equivalent
+# custom versions (uses `Func1` functors defined by C++ class :ct:`Arrhenius1`)
 reactions[2] = ct.Reaction(
     equation='H2 + O <=> H + OH',
     rate=ct.Func1("Arrhenius", [38.7, 2.7, 3150.1542797022735]))
@@ -49,7 +55,8 @@
 gas1b = ct.Solution(thermo='ideal-gas', kinetics='gas',
                     species=species, reactions=reactions)
 
-# construct reactions based on ExtensibleRate: replace two reactions with equivalent
+# %%
+# Construct reactions based on `ExtensibleRate`: replace two reactions with equivalent
 # extensible versions
 class ExtensibleArrheniusData(ct.ExtensibleRateData):
     __slots__ = ("T",)
@@ -107,6 +114,7 @@ def eval(self, data):
 gas2 = ct.Solution(thermo="ideal-gas", kinetics="gas",
                    species=species, reactions=reactions)
 
+# %%
 # construct test case - simulate ignition
 
 def ignition(gas, dT=0):
@@ -124,7 +132,9 @@ def ignition(gas, dT=0):
 
     return t2 - t1, net.solver_stats['steps']
 
-# output results
+# %%
+# Run simulations and output results
+# ----------------------------------
 
 repeat = 100
 print(f"Average time of {repeat} simulation runs for '{mech}' ({fuel})")

samples/python/kinetics/diamond_cvd.py

@@ -17,11 +17,12 @@
 """
 
 import csv
+import matplotlib.pyplot as plt
+import pandas as pd
 import cantera as ct
 
-print('\n******  CVD Diamond Example  ******\n')
-
-# import the model for the diamond (100) surface and the adjacent bulk phases
+# %%
+# Import the model for the diamond (100) surface and the adjacent bulk phases
 d = ct.Interface("diamond.yaml", "diamond_100")
 g = d.adjacent["gas"]
 dbulk = d.adjacent["diamond"]
@@ -37,6 +38,8 @@
 
 xh0 = x[ih]
 
+# %%
+# Calculate growth rate as a function of H mole fraction in the gas
 with open('diamond.csv', 'w', newline='') as f:
     writer = csv.writer(f)
     writer.writerow(['H mole Fraction', 'Growth Rate (microns/hour)'] +
@@ -56,20 +59,18 @@
     print('H concentration, growth rate, and surface coverages '
           'written to file diamond.csv')
 
-try:
-    import matplotlib.pyplot as plt
-    import pandas as pd
-    data = pd.read_csv('diamond.csv')
-
-    data.plot(x="H mole Fraction", y="Growth Rate (microns/hour)", legend=False)
-    plt.xlabel('H Mole Fraction')
-    plt.ylabel('Growth Rate (microns/hr)')
-    plt.show()
-
-    names = [name for name in data.columns if not name.startswith(('H mole', 'Growth'))]
-    data.plot(x='H mole Fraction', y=names, legend=True)
-    plt.xlabel('H Mole Fraction')
-    plt.ylabel('Coverage')
-    plt.show()
-except ImportError:
-    print("Install matplotlib and pandas to plot the outputs")
+# %%
+# Plot the results
+data = pd.read_csv('diamond.csv')
+
+data.plot(x="H mole Fraction", y="Growth Rate (microns/hour)", legend=False)
+plt.xlabel('H Mole Fraction')
+plt.ylabel('Growth Rate (microns/hr)')
+plt.show()
+
+# %%
+names = [name for name in data.columns if not name.startswith(('H mole', 'Growth'))]
+data.plot(x='H mole Fraction', y=names, legend=True)
+plt.xlabel('H Mole Fraction')
+plt.ylabel('Coverage')
+plt.show()

samples/python/kinetics/extract_submechanism.py

@@ -23,7 +23,8 @@
 all_species = ct.Species.list_from_file(input_file)
 species = []
 
-# Filter species
+# %%
+# Filter species.
 for S in all_species:
     comp = S.composition
     if 'C' in comp and 'H' in comp:
@@ -41,7 +42,8 @@
 species_names = {S.name for S in species}
 print('Species: {0}'.format(', '.join(S.name for S in species)))
 
-# Filter reactions, keeping only those that only involve the selected species
+# %%
+# Filter reactions, keeping only those that only involve the selected species.
 ref_phase = ct.Solution(thermo='ideal-gas', kinetics='gas', species=all_species)
 all_reactions = ct.Reaction.list_from_file(input_file, ref_phase)
 reactions = []
@@ -64,10 +66,13 @@
                    transport_model="mixture-averaged",
                    species=species, reactions=reactions)
 
-# Save the resulting mechanism for later use
+# %%
+# Save the resulting mechanism for later use.
 gas2.update_user_header({"description": "CO-H2 submechanism extracted from GRI-Mech 3.0"})
 gas2.write_yaml("gri30-CO-H2-submech.yaml", header=True)
 
+# %%
+# Simulate a flame with each mechanism.
 def solve_flame(gas):
     gas.TPX = 373, 0.05*ct.one_atm, 'H2:0.4, CO:0.6, O2:1, N2:3.76'
 
@@ -81,7 +86,6 @@ def solve_flame(gas):
     sim.solve(0, auto=True)
     return sim
 
-
 t1 = default_timer()
 sim1 = solve_flame(gas1)
 t2 = default_timer()
@@ -90,6 +94,9 @@ def solve_flame(gas):
 t3 = default_timer()
 print('Solved with CO/H2 submechanism in {0:.2f} seconds'.format(t3-t2))
 
+# %%
+# Plot the results to show that the submechanism gives the same results as the original.
+
 plt.plot(sim1.grid, sim1.heat_release_rate,
          lw=2, label='GRI 3.0')
 

samples/python/kinetics/jet_stirred_reactor.py

@@ -41,7 +41,7 @@
 
 f, ax = plt.subplots(1, 3)
 plt.subplots_adjust(wspace=0.6)
-colours = ["xkcd:grey",'xkcd:purple']
+colors = ["xkcd:grey",'xkcd:purple']
 file = 'example_data/ammonia-CO-H2-Alzueta-2023.yaml'
 models = {'Original': 'baseline', 'LMR-R': 'linear-Burke'}
 
@@ -108,9 +108,9 @@ def getTemperatureDependence(gas, inputs):
     tempDependence = getTemperatureDependence(gas,inputs)
     ax[0].plot(tempDependence.index,
                np.subtract(tempDependence['temperature'],tempDependence.index),
-               color=colours[k],label=m)
-    ax[1].plot(tempDependence.index, tempDependence['O2']*100, color=colours[k])
-    ax[2].plot(tempDependence.index, tempDependence['H2']*100, color=colours[k])
+               color=colors[k],label=m)
+    ax[1].plot(tempDependence.index, tempDependence['O2']*100, color=colors[k])
+    ax[2].plot(tempDependence.index, tempDependence['H2']*100, color=colors[k])
 ax[0].plot(inputs['data']['T_range'], inputs['data']['deltaT'], 'o', fillstyle='none',
            color='k', label="Sabia et al.")
 ax[1].plot(inputs['data']['T_range'], inputs['data']['X_O2'], 'o', fillstyle='none',