Added an RMS script (using iJulia) for flux diagram and ROP

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 `<name>.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
+}

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/<you>/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 `<name>.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`.