diff --git a/AUTHORS b/AUTHORS index 1f3cbf71b19..792a71a898e 100644 --- a/AUTHORS +++ b/AUTHORS @@ -60,6 +60,7 @@ Ingmar Schoegl (@ischoegl), Louisiana State University Santosh Shanbhogue (@santoshshanbhogue), Massachusetts Institute of Technology Travis Sikes (@tsikes) Harsh Sinha (@sin-ha) +Sai Krishna Sirumalla (@skrsna), Northeastern University David Sondak Raymond Speth (@speth), Massachusetts Institute of Technology Su Sun (@ssun30), Northeastern University diff --git a/doc/sphinx/conf.py b/doc/sphinx/conf.py index 32d48185704..2080889af88 100644 --- a/doc/sphinx/conf.py +++ b/doc/sphinx/conf.py @@ -14,6 +14,7 @@ import sys, os, re from pathlib import Path from sphinx_gallery.sorting import ExplicitOrder +from sphinx_gallery.scrapers import figure_rst # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the @@ -51,6 +52,41 @@ 'sphinx_copybutton', ] + +class GraphvizScraper(): + """ + Capture Graphviz objects that are assigned to variables in the global namespace. + """ + def __init__(self): + # IDs of graphviz objects that have already been seen and processed + self.processed = set() + + def __repr__(self): + return 'GraphvizScraper' + + def __call__(self, block, block_vars, gallery_conf): + import graphviz + # We use a list to collect references to image names + image_names = list() + + # The `image_path_iterator` is created by Sphinx-Gallery, it will yield + # a path to a file name that adheres to Sphinx-Gallery naming convention. + image_path_iterator = block_vars['image_path_iterator'] + + # Define a list of our already-created figure objects. + for obj in block_vars["example_globals"].values(): + if isinstance(obj, graphviz.Source) and id(obj) not in self.processed: + self.processed.add(id(obj)) + image_path = Path(next(image_path_iterator)).with_suffix(".svg") + obj.format = "svg" + obj.render(image_path.with_suffix("")) + image_names.append(image_path) + + # Use the `figure_rst` helper function to generate the reST for this + # code block's figures. + return figure_rst(image_names, gallery_conf['src_dir']) + + sphinx_gallery_conf = { 'filename_pattern': '\.py', 'example_extensions': {'.py', '.cpp', '.h', '.c', '.f', '.f90', '.m'}, @@ -59,6 +95,7 @@ 'image_srcset': ["2x"], 'remove_config_comments': True, 'ignore_repr_types': r'matplotlib\.(text|axes|legend)', + 'image_scrapers': ('matplotlib', GraphvizScraper()), 'examples_dirs': [ '../samples/python/', '../samples/cxx/', diff --git a/samples/python/reactors/interactive_path_diagram.py b/samples/python/reactors/interactive_path_diagram.py new file mode 100644 index 00000000000..915b25b7ba3 --- /dev/null +++ b/samples/python/reactors/interactive_path_diagram.py @@ -0,0 +1,284 @@ +""" +Interactive Reaction Path Diagrams +================================== + +This example uses ``ipywidgets`` to create interactive displays of reaction path +diagrams from Cantera simulations. + +Requires: cantera >= 3.0.0, matplotlib >= 2.0, ipywidgets, graphviz, scipy + +.. tags:: Python, combustion, reactor network, plotting, reaction path analysis + +.. tip:: + To try the interactive features, download the Jupyter notebook version of this + example: :download:`interactive_path_diagram.ipynb`. +""" + +# %% +import numpy as np +from scipy import integrate +import graphviz +import os +from matplotlib import pyplot as plt +from collections import defaultdict +import cantera as ct + +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 import get_ipython + if "IPKernelApp" not in get_ipython().config: + raise ImportError("console") + if "VSCODE_PID" in os.environ: + raise ImportError("vscode") +except (ImportError, AttributeError): + is_interactive = False +else: + is_interactive = True + +if is_interactive: + from IPython.display import display + from matplotlib_inline.backend_inline import set_matplotlib_formats + set_matplotlib_formats('pdf', 'svg') + from ipywidgets import widgets, interact + + +# %% +# When using Cantera, the first thing you usually need is an object representing some +# phase of matter. Here, we'll create a gas mixture using GRI-Mech: + +gas = ct.Solution("gri30.yaml") + +# %% +# Use Shock tube ignition delay measurement conditions corresponding to the experiments +# by Spadaccini and Colket [1]_. +# +# * CH₄-C₂H₆-O₂-Ar (3.29%-0.21%-7%-89.5%) +# * :math:`\phi` = 1.045 +# * P = 6.1 - 7.6 atm +# * T = 1356 - 1688 K + +# Set temperature, pressure, and composition +gas.TPX = 1550.0, 6.5 * ct.one_atm, "CH4:3.29, C2H6:0.21, O2:7 , Ar:89.5" + +# %% +# Residence time is close to ignition delay reported by Spadaccini and Colket (1994). + +residence_time = 1e-3 + +# %% +# Create a batch reactor object and set solver tolerances + +reactor = ct.IdealGasConstPressureReactor(gas, energy="on") +reactor_network = ct.ReactorNet([reactor]) +reactor_network.atol = 1e-12 +reactor_network.rtol = 1e-12 + +# %% +# Store time, pressure, temperature and mole fractions + +profiles = defaultdict(list) +time = 0 +steps = 0 +while time < residence_time: + profiles["time"].append(time) + profiles["pressure"].append(gas.P) + profiles["temperature"].append(gas.T) + profiles["mole_fractions"].append(gas.X) + time = reactor_network.step() + steps += 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 = profiles["pressure"][plot_step] + T = profiles["temperature"][plot_step] + X = profiles["mole_fractions"][plot_step] + time = profiles["time"][plot_step] + gas.TPX = T, P, X + + diagram = ct.ReactionPathDiagram(gas, element) + diagram.threshold = threshold + diagram.title = f"time = {time:.2g} s" + + diagram.show_details = details + graph = graphviz.Source(diagram.get_dot()) + if is_interactive: + display(graph) + else: + 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, + ), + ) +else: + # For non-interactive use, just draw the diagram for a specified time step + diagram = plot_reaction_path_diagrams( + plot_step=100, + threshold=0.1, + details=False, + element="C" + ) + +# %% +# Interactive plot of instantaneous fluxes +# ---------------------------------------- +# +# Find reactions containing the species of interest, C₂H₆ in this case. + +C2H6_stoichiometry = np.zeros_like(gas.reactions()) +for i, r in enumerate(gas.reactions()): + C2H6_moles = r.products.get("C2H6", 0) - r.reactants.get("C2H6", 0) + C2H6_stoichiometry[i] = C2H6_moles +C2H6_reaction_indices = C2H6_stoichiometry.nonzero()[0] + +# %% +# The following cell calculates net rates of progress of reactions containing the +# species of interest and stores them. + +profiles["C2H6_production_rates"] = [] +for i in range(len(profiles["time"])): + X = profiles["mole_fractions"][i] + t = profiles["time"][i] + T = profiles["temperature"][i] + P = profiles["pressure"][i] + gas.TPX = (T, P, X) + C2H6_production_rates = ( + gas.net_rates_of_progress + * C2H6_stoichiometry # [kmol/m^3/s] + * gas.volume_mass # Specific volume [m^3/kg]. + ) # overall, mol/s/g (g total in reactor, same basis as N_atoms_in_fuel) + + profiles["C2H6_production_rates"].append( + C2H6_production_rates[C2H6_reaction_indices] + ) + +# Create the instantaneous flux plot. When executed as a Jupyter Notebook, the threshold +# for annotating of reaction strings can be changed using the slider provided by +# IPyWidgets. + +plt.rcParams["figure.constrained_layout.use"] = True + +def plot_instantaneous_fluxes(profiles, annotation_cutoff): + profiles = profiles + fig = plt.figure(figsize=(6, 6)) + plt.plot(profiles["time"], np.array(profiles["C2H6_production_rates"])) + + for i, C2H6_production_rate in enumerate( + np.array(profiles["C2H6_production_rates"]).T + ): + peak_index = abs(C2H6_production_rate).argmax() + peak_time = profiles["time"][peak_index] + peak_C2H6_production = C2H6_production_rate[peak_index] + reaction_string = gas.reaction_equations(C2H6_reaction_indices)[i] + + if abs(peak_C2H6_production) > annotation_cutoff: + plt.annotate( + reaction_string.replace("<=>", "="), + xy=(peak_time, peak_C2H6_production), + xytext=( + peak_time * 2, + ( + peak_C2H6_production + + 0.003 + * (peak_C2H6_production / abs(peak_C2H6_production)) + * (abs(peak_C2H6_production) > 0.005) + * (peak_C2H6_production < 0.06) + ), + ), + arrowprops=dict( + arrowstyle="->", + color="black", + relpos=(0, 0.6), + linewidth=2, + ), + horizontalalignment="left", + ) + + plt.xlabel("Time (s)") + plt.ylabel("Net rates of C2H6 production") + plt.show() + +if is_interactive: + interact( + plot_instantaneous_fluxes, + annotation_cutoff=widgets.FloatSlider(value=0.1, min=1e-2, max=4, steps=10), + profiles=widgets.fixed(profiles) + ) +else: + plot_instantaneous_fluxes(annotation_cutoff=0.1, profiles=profiles) + +# %% +# Interactive plot of integrated fluxes +# ------------------------------------- +# +# When executed as a Jupyter Notebook, the threshold for annotating of reaction strings +# can be changed using the slider provided by iPyWidgets + +# Integrate fluxes over time +integrated_fluxes = integrate.cumulative_trapezoid( + np.array(profiles["C2H6_production_rates"]), + np.array(profiles["time"]), + axis=0, + initial=0, +) + +def plot_integrated_fluxes(profiles, integrated_fluxes, annotation_cutoff): + profiles = profiles + integrated_fluxes = integrated_fluxes + fig = plt.figure(figsize=(6, 6)) + plt.plot(profiles["time"], integrated_fluxes) + final_time = profiles["time"][-1] + for i, C2H6_production in enumerate(integrated_fluxes.T): + total_C2H6_production = C2H6_production[-1] + reaction_string = gas.reaction_equations(C2H6_reaction_indices)[i] + + if abs(total_C2H6_production) > annotation_cutoff: + plt.text(final_time * 1.06, total_C2H6_production, reaction_string, + fontsize=8) + + plt.xlabel("Time (s)") + plt.ylabel("Integrated net rate of progress") + plt.title("Cumulative C₂H₆ formation") + plt.show() + +if is_interactive: + interact( + plot_integrated_fluxes, + annotation_cutoff=widgets.FloatLogSlider( + value=1e-5, min=-5, max=-4, base=10, step=0.1 + ), + profiles=widgets.fixed(profiles), + integrated_fluxes=widgets.fixed(integrated_fluxes) + ) +else: + plot_integrated_fluxes( + profiles=profiles, + integrated_fluxes=integrated_fluxes, + annotation_cutoff=1e-5 + ) + +# %% +# References +# ---------- +# .. [1] L. J. Spadaccini and M. B. Colket (1994). "Ignition delay characteristics of +# methane fuels", *Progress in Energy and Combustion Science,* 20:5, 431-460. +# Prog. Energy Combust. Sci. 20, 431. +# https://doi.org/10.1016/0360-1285(94)90011-6.