From 08960d704752e5a9ee7e25b1b3502b79a407bc63 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Thu, 29 Jun 2023 14:59:07 +0300 Subject: [PATCH] Add RMS notebook for constant-T,P reactor ROP & flux diagrams Julia/ReactionMechanismSimulator notebook: simulate a constant-T,P ideal-gas batch reactor, then plot rate-of-production and a reaction flux diagram. Reworked from the original scratch notebook to be reusable: all inputs (mechanism path, conditions, initial mole fractions, ROP/flux targets and tolerances) live in one User-inputs cell instead of hardcoded personal paths; added section markdown; dropped unused imports and the unused run_sa flag; cleared embedded outputs. Adds README.md (environment setup + usage) and an example/ scaffold for a user-supplied mechanism. Co-Authored-By: Claude Fable 5 --- RMS/Constant T P ideal gas reactor.ipynb | 213 +++++++++++++++++++++++ RMS/README.md | 89 ++++++++++ 2 files changed, 302 insertions(+) create mode 100644 RMS/Constant T P ideal gas reactor.ipynb create mode 100644 RMS/README.md diff --git a/RMS/Constant T P ideal gas reactor.ipynb b/RMS/Constant T P ideal gas reactor.ipynb new file mode 100644 index 0000000..1aaf653 --- /dev/null +++ b/RMS/Constant T P ideal gas reactor.ipynb @@ -0,0 +1,213 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e683ee51", + "metadata": {}, + "source": [ + "# Constant-*T*,*P* ideal-gas reactor — ROP & flux diagram (RMS)\n", + "\n", + "Simulates a homogeneous, isothermal, isobaric (constant *T*, *P*) ideal-gas batch\n", + "reactor with [ReactionMechanismSimulator.jl](https://github.com/ReactionMechanismSimulator/ReactionMechanismSimulator.jl)\n", + "(RMS), then produces:\n", + "\n", + "- a **rate-of-production (ROP)** plot for a chosen species, and\n", + "- a **reaction-flux diagram** at a chosen time.\n", + "\n", + "> **This is a Julia notebook, not Python** — it needs a Julia kernel with RMS\n", + "> installed. See [`README.md`](README.md) for one-time environment setup and for how\n", + "> to obtain a mechanism file. In normal use you only edit the **User inputs** cell;\n", + "> every other cell runs unchanged." + ] + }, + { + "cell_type": "markdown", + "id": "4b6d0280", + "metadata": {}, + "source": [ + "## 1. Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e3758b23", + "metadata": {}, + "outputs": [], + "source": [ + "using ReactionMechanismSimulator\n", + "using DifferentialEquations\n", + "using Sundials # provides the CVODE_BDF stiff ODE solver\n", + "using PyPlot # plotting backend used by plotrops / getfluxdiagram" + ] + }, + { + "cell_type": "markdown", + "id": "f77bfc74", + "metadata": {}, + "source": [ + "## 2. User inputs\n", + "\n", + "Everything you need to change lives in this one cell." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6623b9cb", + "metadata": {}, + "outputs": [], + "source": [ + "# --- Mechanism -------------------------------------------------------------\n", + "# Point `mech` at a mechanism RMS can read (see README → \"Getting a mechanism\"):\n", + "# 1. A native RMS file: mech = \"path/to/chem.rms\" (leave spc_dict = nothing)\n", + "# 2. A Chemkin file + RMG species dictionary:\n", + "# mech = \"path/to/chem_annotated.inp\"\n", + "# spc_dict = \"path/to/species_dictionary.txt\"\n", + "# On first read of a Chemkin file, RMS caches a `.rms` next to it.\n", + "mech = \"path/to/mechanism.rms\" # <- set to your RMS or Chemkin file (see README)\n", + "spc_dict = nothing # set to a species_dictionary.txt path for Chemkin input\n", + "\n", + "# --- Reactor conditions ----------------------------------------------------\n", + "T = 1300.0 # temperature [K]\n", + "P = 126656.0 # pressure [Pa]\n", + "max_time = 0.2 # simulated time horizon [s]\n", + "\n", + "# Initial composition as mole fractions. Keys MUST match species names in the\n", + "# mechanism (RMG-style labels, e.g. \"NH3\", \"O2\", \"Ar\"). They need not sum to 1 —\n", + "# RMS normalizes.\n", + "initial_mole_fractions = Dict(\"NH3\" => 0.08, \"O2\" => 0.06, \"Ar\" => 0.086)\n", + "\n", + "# --- Analysis --------------------------------------------------------------\n", + "rop_species = \"NO\" # species shown in the ROP plot; must exist in the mechanism\n", + "rop_tol = 0.25 # ROP plot: keep reactions contributing > this fraction of the max\n", + "\n", + "flux_time = 0.0015 # time [s] at which to draw the flux diagram\n", + "flux_conc_tol = 1e-12 # min concentration for a species to appear\n", + "flux_rate_tol = 1e-12 # min species rate for a reaction to appear\n", + "\n", + "# --- ODE solver tolerances -------------------------------------------------\n", + "abstol = 1e-20\n", + "reltol = 1e-12" + ] + }, + { + "cell_type": "markdown", + "id": "958605ea", + "metadata": {}, + "source": [ + "## 3. Load the mechanism" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c2aff15", + "metadata": {}, + "outputs": [], + "source": [ + "phase_dict = spc_dict === nothing ? readinput(mech) : readinput(mech, spcdict=spc_dict)\n", + "spcs = phase_dict[\"phase\"][\"Species\"]\n", + "rxns = phase_dict[\"phase\"][\"Reactions\"]\n", + "println(\"Loaded $(length(spcs)) species and $(length(rxns)) reactions.\")" + ] + }, + { + "cell_type": "markdown", + "id": "132dbbd5", + "metadata": {}, + "source": [ + "## 4. Build the ideal-gas phase and constant-*T*,*P* domain" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "370982bc", + "metadata": {}, + "outputs": [], + "source": [ + "ig = IdealGas(spcs, rxns, name=\"gas\")\n", + "\n", + "initialconds = merge(Dict(\"T\" => T, \"P\" => P), initial_mole_fractions)\n", + "domain, y0, p = ConstantTPDomain(phase=ig, initialconds=initialconds)" + ] + }, + { + "cell_type": "markdown", + "id": "8b510f70", + "metadata": {}, + "source": [ + "## 5. Solve the reactor" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c1d68b46", + "metadata": {}, + "outputs": [], + "source": [ + "react = Reactor(domain, y0, (0.0, max_time); p=p)\n", + "sol = solve(react.ode, CVODE_BDF(), abstol=abstol, reltol=reltol)\n", + "bsol = Simulation(sol, domain)" + ] + }, + { + "cell_type": "markdown", + "id": "4c69fc8c", + "metadata": {}, + "source": [ + "## 6. Rate-of-production (ROP) analysis\n", + "\n", + "Plots the net rate of production of `rop_species` and the individual reactions\n", + "contributing more than `rop_tol` of the peak rate." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c16b7c07", + "metadata": {}, + "outputs": [], + "source": [ + "plotrops(bsol, rop_species, tol=rop_tol)" + ] + }, + { + "cell_type": "markdown", + "id": "97d0b049", + "metadata": {}, + "source": [ + "## 7. Reaction-flux diagram\n", + "\n", + "Draws the species/reaction flux network at `flux_time`. Requires Graphviz `dot`\n", + "on the `PATH` (see README)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e532da9a", + "metadata": {}, + "outputs": [], + "source": [ + "getfluxdiagram(bsol, flux_time, concentrationtol=flux_conc_tol, speciesratetolerance=flux_rate_tol)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.10", + "language": "julia", + "name": "julia-1.10" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/RMS/README.md b/RMS/README.md new file mode 100644 index 0000000..f43ba23 --- /dev/null +++ b/RMS/README.md @@ -0,0 +1,89 @@ +# RMS — reactor simulation, ROP & flux diagrams + +[`Constant T P ideal gas reactor.ipynb`](Constant%20T%20P%20ideal%20gas%20reactor.ipynb) +simulates a homogeneous, isothermal, isobaric (constant *T*, *P*) ideal-gas batch +reactor with [ReactionMechanismSimulator.jl](https://github.com/ReactionMechanismSimulator/ReactionMechanismSimulator.jl) +(RMS) and produces a **rate-of-production (ROP)** plot and a **reaction-flux diagram**. + +RMS is a **Julia** package, so this is a Julia notebook — it does *not* run in the +Cantera/Python environment used by the rest of this repo (e.g. `Cantera/ROP`). The two +are complementary: reach for RMS when you want its native `getfluxdiagram` / +`plotrops` and its tight coupling to RMG mechanisms. + +## Environment (one-time setup) + +RMS installation details drift between releases; the +[RMS README](https://github.com/ReactionMechanismSimulator/ReactionMechanismSimulator.jl#installation) +is the authoritative source. The setup below is the common path. + +1. **Install Julia** (1.9+ recommended) — e.g. via [`juliaup`](https://github.com/JuliaLang/juliaup). + Download the installer, review it, then run it (rather than piping straight into a shell): + ```bash + curl -fsSL https://install.julialang.org -o install-julia.sh + sh install-julia.sh + ``` + See the [juliaup README](https://github.com/JuliaLang/juliaup#installation) for Windows + and other install options. + +2. **Install RMS and the packages this notebook uses** from a Julia REPL: + ```julia + using Pkg + Pkg.add(["ReactionMechanismSimulator", "DifferentialEquations", "Sundials", "PyPlot"]) + ``` + +3. **Enable Chemkin input (optional).** To read RMG `chem_annotated.inp` + + `species_dictionary.txt`, RMS calls RMG-Py through `PyCall`, so point PyCall at a + Python that has `rmgpy` installed (e.g. this repo's `rmg_env` conda env): + ```julia + ENV["PYTHON"] = "/home//anaconda3/envs/rmg_env/bin/python" + using Pkg; Pkg.build("PyCall") + ``` + You can skip this if you only ever read native `.rms` files. + +4. **Register the Jupyter Julia kernel** so this notebook opens with a Julia kernel: + ```julia + using Pkg; Pkg.add("IJulia") + using IJulia; installkernel("Julia") + ``` + +5. **Install Graphviz** (`dot`) — RMS renders the flux diagram with it: + ```bash + conda install -c conda-forge graphviz # or: sudo apt install graphviz + ``` + +## Getting a mechanism + +Set `mech` in the notebook to a mechanism file (the default is a `path/to/mechanism.rms` +placeholder). RMS reads either: + +- **A native RMS file** (`.rms` / `.yml`) — no Python needed. Set `mech` to it and leave + `spc_dict = nothing`. +- **An RMG Chemkin pair** — `chem_annotated.inp` + `species_dictionary.txt` from any RMG + job. Set `mech` to the `.inp` and `spc_dict` to the dictionary. On first read RMS caches + a `.rms` next to the `.inp`, which loads faster afterwards. + +Species names in `initial_mole_fractions` and `rop_species` must match the labels used in +your mechanism (RMG-style, e.g. `NH3`, `O2`, `Ar`, `NO`). + +Keep mechanism files outside the repo (or add them to `.gitignore`) — they are large, +run-specific inputs and should not be committed. + +## Running + +1. Open the notebook with a **Julia kernel** (JupyterLab/Notebook, or VS Code with the + Julia extension). +2. Edit the **User inputs** cell — mechanism path, `T`, `P`, `max_time`, initial mole + fractions, `rop_species`, and the flux-diagram time/tolerances. +3. Run all cells top to bottom. + +Expected output: cell 6 shows the ROP plot for `rop_species`; cell 7 renders the flux +diagram at `flux_time`. The first run in a fresh Julia session is slow — RMS and its ODE +dependencies precompile. + +## Notes + +- First-call latency ("time to first plot") is a Julia trait, not an error — subsequent + runs in the same session are fast. +- The stiff solver (`CVODE_BDF` from Sundials) with tight tolerances suits combustion + kinetics; loosen `abstol`/`reltol` if a run is slow and you don't need that precision. +- Flux-diagram output empty or too dense? Adjust `flux_conc_tol` / `flux_rate_tol`.