—

Interactive Reaction Path Diagrams for CH₄ plasma chemistry.

This example uses ipywidgets to create interactive displays of reaction path diagrams from Cantera simulations.

Inspired by https://cantera.org/3.1/examples/python/kinetics/interactive_path_diagram.html.

Requires: cantera >= 3.0.0, matplotlib >= 2.0, ipywidgets, graphviz, scipy

Import the required libraries.

import cantera as ct
import graphviz
import matplotlib.pyplot as plt
import numpy as np

from bloc.paths import get_mechanism_path

print(f"Using Cantera version: {ct.__version__}")

# Determine if we're running in a Jupyter Notebook. If so, we can enable the interactive
# diagrams. Otherwise, just draw output for a single set of inputs.
try:
    from IPython.core.getipython import get_ipython

    config = get_ipython()
    if config is None:
        raise ImportError("console")

except (ImportError, AttributeError):
    is_interactive = False
else:
    is_interactive = True

if is_interactive:
    from IPython.display import display
    from ipywidgets import interact, widgets

print(
    "Running in interactive mode"
    if is_interactive
    else "Running in non-interactive mode"
)

Create the gas object and set the initial conditions.

# Load the mechanism.
mechanism = get_mechanism_path("Caltech.yaml")
gas = ct.Solution(mechanism)

# Set the initial temperature, pressure and composition.
gas.TPX = 300, ct.one_atm, "CH4:1.0"
gas.HPX = gas.enthalpy_mass + 12e6, ct.one_atm, "CH4:1.0"

# Set the residence time for the reactor network.
residence_time = 1e0  # s

# Create an ideal gas constant pressure reactor and set solver tolerances.
reactor = ct.IdealGasConstPressureMoleReactor(gas, energy="on")
reactor_network = ct.ReactorNet([reactor])
precon = ct.AdaptivePreconditioner()  # Create the preconditioner
reactor_network.preconditioner = precon  # Add it to the network

# Plot options.
xlim = (1e-12, residence_time)

# Wanted species.
wanted_species = "C2H2"

# Store the reactions and reaction equations.
gas_reactions = gas.reactions()
gas_reaction_equations = gas.reaction_equations()

Store time, pressure, temperature and mole fractions.

# Store the states in a SolutionArray.
states = ct.SolutionArray(gas, 1, extra={"t": [0.0]})

# Run the simulation
steps = 0
time_steps = np.logspace(-6, np.log10(residence_time), 100)
for time_step in time_steps:
    reactor_network.advance(time_step)
    states.append(gas.state, t=reactor_network.time)
    print(f"Step {steps} at time {reactor_network.time} s")
    steps += 1


time = states.t  # type: ignore

Plot temperature evolution with time.

fig, ax = plt.subplots()
ax.plot(
    time,
    states.T,
)
ax.set_ylabel("Temperature [K]")
ax.set_title("Temperature evolution")
ax.set_xlabel("Time [s]")
ax.set_xscale("log")

Plot species evolution with time.

fig, ax = plt.subplots()

species_to_plot = set(["CH4", "H2", "C2H2", "C2H4", "C2H6", "H", wanted_species])
for species in species_to_plot:
    X = states(species).X
    ax.plot(time, X, label=species)

    x = time[np.argmax(X)]
    y = X[np.argmax(X)]
    if x < 1e-11:
        x = 1e-11
        y = X[np.where(time >= 1e-11)[0][0]]


ax.legend()
ax.set_xlabel("Time [s]")
ax.set_ylabel("Mole fraction [%]")
ax.set_title("Mole fraction evolution with time.")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlim(xlim)
ax.set_ylim(1e-6, 1)

Interactive reaction path diagram.

When executed as a Jupyter Notebook, the plotted time step, threshold and element can be changed using the slider provided by IPyWidgets.

def plot_reaction_path_diagrams(plot_step, threshold, details, element):
    P = states.P[plot_step]
    T = states.T[plot_step]
    X = states.X[plot_step]
    gas.TPX = T, P, X

    diagram = ct.ReactionPathDiagram(gas, element)
    diagram.threshold = threshold
    diagram.title = f"time = {time[plot_step]:.2g} s"

    diagram.show_details = details
    if is_interactive:
        graph = graphviz.Source(diagram.get_dot())
        display(graph)
    else:
        graph = graphviz.Source(diagram.get_dot(), format="svg")
        return graph


if is_interactive:
    interact(
        plot_reaction_path_diagrams,
        plot_step=widgets.IntSlider(value=100, min=0, max=steps - 1, step=10),
        threshold=widgets.FloatSlider(value=0.1, min=0.001, max=0.4, step=0.01),
        details=widgets.ToggleButton(),
        element=widgets.Dropdown(
            options=gas.element_names,
            value="C",
            description="Element",
            disabled=False,
        ),
    )
    diagram = ""
else:
    # For non-interactive use, just draw the diagram for a specified time step.
    diagram = plot_reaction_path_diagrams(
        plot_step=520, threshold=0.1, details=False, element="C"
    )


class PlotGraphviz:
    # See https://stackoverflow.com/questions/65008861/capturing-graphviz-figures-in-sphinx-gallery.
    def __init__(self, dot_string):
        self.dot_string = dot_string

    def _repr_html_(self):
        return graphviz.Source(self.dot_string)._repr_image_svg_xml()


PlotGraphviz(str(diagram))