{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Preconditioned integration of a coupled reactor network\n\nThis example demonstrates the effect of including cross-reactor Jacobian terms from\nwalls and flow devices in the sparse preconditioner matrix.\n\nThe network represents the sectors of a simplified annular combustor with non-uniform\nfuel/air mixing. A single methane supply and a single air supply feed each sector\nthrough a pair of mass flow controllers, with the relative flow split set to give\neach sector a distinct equivalence ratio \u2014 as would happen if the fuel injectors\naround the combustor were running at slightly different operating points. All\nsectors vent to a shared exhaust and adjacent sectors are separated by thin metal\nliner walls with high heat transfer. The strong wall coupling drives all sectors\ntoward a nearly uniform temperature while the equivalence-ratio variation keeps the\nspecies distributions distinct.\n\nFour integration strategies are compared:\n\n- Preconditioned, with wall and flow-device Jacobian terms included (default)\n- Preconditioned, with wall terms excluded (``skip-walls`` option)\n- Preconditioned, with flow-device terms excluded (``skip-flow-devices`` option)\n- Direct integration without preconditioning\n\nBoth the wall coupling and the flow-device coupling introduce significant\ncross-reactor entries in the Jacobian. Omitting either set of terms substantially\nincreases the number of Krylov iterations and integrator step rejections; in this\nnetwork, dropping the wall terms is enough to make the preconditioned solve slower\nthan direct integration with no preconditioning at all.\n\nRequires: cantera >= 4.0.0, matplotlib >= 2.0\n\n.. tags:: Python, combustion, reactor network, preconditioner\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import cantera as ct\nimport matplotlib.pyplot as plt\nfrom time import perf_counter\n\nplt.rcParams[\"figure.constrained_layout.use\"] = True" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulation Setup\n\nEach sector is an :ct:`IdealGasMoleReactor`. A single fuel supply (methane) and\nair supply feed each sector through a pair of :ct:`MassFlowController`\nobjects whose rates are chosen to produce the prescribed equivalence ratio at the\nspecified total flow rate. Each sector vents to a shared exhaust reservoir through\na :ct:`Valve`. Adjacent sectors are coupled by a :ct:`Wall` with a high\nheat-transfer coefficient and a small but non-zero expansion-rate coefficient,\nmodeling thin metal partitions between the sectors.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "END_TIME = 1.0\nN_OUTPUT_STEPS = 240\nWALL_U = 5.0e5 # W/m\u00b2/K; very tight thermal coupling through the liner\nWALL_K = 1.0e-4 # m/(Pa\u00b7s); modest mechanical compliance of the liner\nK_EXHAUST = 1.0e-4 # kg/s/Pa; exhaust relief valve\nP_SUPPLY = 5.0 * ct.one_atm\nP_EXHAUST = 3.5 * ct.one_atm\nSECTOR_VOLUME = 0.05 # m\u00b3\nSECTOR_MASS_FLOW = 0.1 # kg/s\nT_SUPPLY = 450.0\nFUEL = \"CH4:1\"\nOXIDIZER = \"O2:1, N2:3.76\"\nEQUIVALENCE_RATIOS = [0.70, 0.85, 1.00, 1.15, 1.30, 1.45]\n\ngas = ct.Solution(\"gri30.yaml\", transport_model=\"none\")\n\ndef make_network(skip_walls=False, skip_flow_devices=False, preconditioner=True):\n \"\"\"Build the annular-combustor network for a given preconditioner configuration.\"\"\"\n gas.TPX = T_SUPPLY, P_SUPPLY, FUEL\n fuel_supply = ct.Reservoir(gas, name=\"fuel\")\n fuel_species = list(gas.mole_fraction_dict().keys())\n\n gas.TPX = T_SUPPLY, P_SUPPLY, OXIDIZER\n air_supply = ct.Reservoir(gas, name=\"air\")\n\n reactors = []\n for phi in EQUIVALENCE_RATIOS:\n gas.set_equivalence_ratio(phi, FUEL, OXIDIZER)\n # Split the prescribed total mass flow between the fuel and air streams to\n # produce the requested equivalence ratio at the sector inlet.\n mass_fuel = sum(gas[fuel_species].Y)\n mdot_fuel = SECTOR_MASS_FLOW * mass_fuel\n mdot_air = SECTOR_MASS_FLOW - mdot_fuel\n\n # Initialize each sector at the hot adiabatic equilibrium of its own fuel-air\n # mixture so that integration starts from an already-burning state and the\n # interesting dynamics are the cross-coupled approach to steady operation.\n gas.TP = T_SUPPLY, P_SUPPLY\n gas.equilibrate(\"HP\")\n sector = ct.IdealGasMoleReactor(gas, name=f\"Reactor @ \u03c6={phi:.2f}\")\n sector.volume = SECTOR_VOLUME\n\n ct.MassFlowController(fuel_supply, sector, mdot=mdot_fuel, name=\"Fuel MFC\")\n ct.MassFlowController(air_supply, sector, mdot=mdot_air, name=\"Air MFC\")\n reactors.append(sector)\n\n gas.TPX = 1500.0, P_EXHAUST, \"N2:1\"\n exhaust = ct.Reservoir(gas, name=\"exhaust\")\n for sector in reactors:\n ct.Valve(sector, exhaust, K=K_EXHAUST)\n\n for left, right in zip(reactors, reactors[1:]):\n ct.Wall(left, right, U=WALL_U, K=WALL_K)\n ct.Wall(reactors[-1], reactors[0], U=WALL_U, K=WALL_K)\n\n net = ct.ReactorNet(reactors)\n if preconditioner:\n settings = {\"skip-falloff\": True, \"skip-third-bodies\": True}\n if skip_walls:\n settings[\"skip-walls\"] = True\n if skip_flow_devices:\n settings[\"skip-flow-devices\"] = True\n net.derivative_settings = settings\n net.preconditioner = ct.AdaptivePreconditioner()\n return net, reactors\n\n\ndef integrate_network(**kwargs):\n \"\"\"Integrate the network and record per-sector trajectories and solver stats.\"\"\"\n net, reactors = make_network(**kwargs)\n co2 = reactors[0].phase.species_index(\"CO2\")\n co = reactors[0].phase.species_index(\"CO\")\n\n times = []\n states = [ct.SolutionArray(gas, extra=[\"t\"]) for _ in reactors]\n\n tic = perf_counter()\n for n in range(N_OUTPUT_STEPS):\n net.advance((n + 1) * END_TIME / N_OUTPUT_STEPS)\n times.append(net.time)\n for i, r in enumerate(reactors):\n states[i].append(state=r.phase.state, t=net.time)\n elapsed = perf_counter() - tic\n\n stats = net.solver_stats\n return {\n \"net\": net,\n \"elapsed\": elapsed,\n \"steps\": stats[\"steps\"],\n \"rhs_evals\": stats[\"rhs_evals\"],\n \"err_test_fails\": stats.get(\"err_test_fails\", 0),\n \"lin_iters\": stats.get(\"lin_iters\", 0),\n \"lin_conv_fails\": stats.get(\"lin_conv_fails\", 0),\n \"nonlinear_conv_fails\": stats.get(\"nonlinear_conv_fails\", 0),\n \"prec_evals\": stats.get(\"prec_evals\", 0),\n \"times\": times,\n \"states\": states,\n }\n\n\ndef print_stats(label, stats):\n print(f\"\\n{label}\")\n print(\"-\" * len(label))\n for key in (\n \"elapsed\",\n \"steps\",\n \"rhs_evals\",\n \"err_test_fails\",\n \"lin_iters\",\n \"lin_conv_fails\",\n \"nonlinear_conv_fails\",\n \"prec_evals\",\n ):\n value = stats[key]\n if isinstance(value, float):\n print(f\" {key:22s} {value:12.6g}\")\n else:\n print(f\" {key:22s} {value:12d}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Preconditioned integration with all cross-reactor terms included\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "full = integrate_network()\nprint_stats(\"Preconditioned, walls + flow devices included\", full)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Network structure\n\nHere we show the reactor network structure. To make the diagram more readable, we\nmodify the layout to use a left-to-right orientation and arrange the sectors in a\nsingle column. The wall-velocity arrows are omitted to avoid cluttering the diagram.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "diagram = full[\"net\"].draw(\n graph_attr={\"rankdir\": \"LR\", \"nodesep\": \"1.0\", \"ranksep\": \"1.5\"},\n show_wall_velocity=False,\n)\nsectors = full[\"net\"].reactors\nwith diagram.subgraph() as column:\n column.attr(rank=\"same\")\n for r in sectors:\n column.node(r.name)\n # Invisible edges fix the order along the column; without them graphviz can\n # rotate the ring of sectors arbitrarily.\n for upper, lower in zip(sectors, sectors[1:]):\n column.edge(upper.name, lower.name, style=\"invis\")\n# Drop the local subgraph handle so the Sphinx-Gallery scraper does not render\n# it as a second standalone diagram.\ndel column\ndiagram" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Preconditioned integration with wall terms skipped\n\nWithout the wall cross-derivatives, the preconditioner misses the strong thermal\ncoupling between sectors. The Krylov solver requires many more iterations per step\nand generates frequent convergence failures.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "skip_walls = integrate_network(skip_walls=True)\nprint_stats(\"Preconditioned, wall terms skipped\", skip_walls)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Preconditioned integration with flow-device terms skipped\n\nThe mass flow controllers and exhaust valves also produce cross-reactor entries\nthrough the inlet enthalpy and pressure-driven flow derivatives. Dropping those\nentries is less damaging than dropping the wall terms but still markedly degrades\nthe preconditioner.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "skip_flow = integrate_network(skip_flow_devices=True)\nprint_stats(\"Preconditioned, flow-device terms skipped\", skip_flow)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Direct integration without preconditioning\n\nThe reference dense direct solve is faster than the skip-walls configuration but\nslower than the full preconditioner, showing that for this kind of tightly coupled\nnetwork a preconditioner missing key cross-reactor terms can be worse than no\npreconditioner at all.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "direct = integrate_network(preconditioner=False)\nprint_stats(\"Direct integration (no preconditioner)\", direct)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solver cost comparison\n\nThe ratios summarize how each configuration compares to the fully preconditioned\nsolve. The full preconditioner is the fastest and the wall cross-terms are the\nmost important: skipping them produces the slowest configuration of the four,\nwhile skipping the flow-device terms gives up most of the remaining preconditioner\nbenefit.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(\"\\nIntegration cost relative to full preconditioner\")\nprint(\"------------------------------------------------\")\nheader = f\"{'configuration':38s} {'elapsed':>10s} {'steps':>8s} {'lin_iters':>10s}\"\nprint(header)\nprint(\"-\" * len(header))\nfor label, stats in (\n (\"preconditioned, full\", full),\n (\"preconditioned, skip walls\", skip_walls),\n (\"preconditioned, skip flow devices\", skip_flow),\n (\"direct (no preconditioner)\", direct),\n):\n e = stats[\"elapsed\"] / full[\"elapsed\"]\n s = stats[\"steps\"] / full[\"steps\"]\n li_denom = full[\"lin_iters\"] or 1\n li = stats[\"lin_iters\"] / li_denom\n print(f\"{label:38s} {e:10.2f} {s:8.2f} {li:10.2f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Trajectories\n\nAll four configurations produce essentially the same results; the plot shows\nthe results from the fully preconditioned solve. The left panel shows the\nstrong wall coupling pulling the sectors to a nearly common temperature that\nasymptotically rises as the network approaches steady operation. The right\npanel shows that the chemistry stays distinct between sectors: CO\u2082 peaks\nnear the stoichiometric point ($\\phi = 1$) and falls off on either side,\nwhile CO is essentially absent in the lean sectors and accumulates monotonically\nin the rich sectors.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, (ax_T, ax_X) = plt.subplots(1, 2, figsize=(13, 4))\nphi_labels = [rf\"$\\phi = {phi}$\" for phi in EQUIVALENCE_RATIOS]\nfor i in range(len(EQUIVALENCE_RATIOS)):\n ax_T.plot(full[\"states\"][i].t, full[\"states\"][i].T, label=phi_labels[i])\nax_T.set(xlabel=\"time [s]\", ylabel=\"temperature [K]\")\nax_T.legend(fontsize=8, ncol=2)\n\nfor species in (\"CO2\", \"CO\", \"H2O\", \"O2\"):\n X = [states(species).X[-1] for states in full[\"states\"]]\n ax_X.plot(EQUIVALENCE_RATIOS, X, \"o-\", label=f\"{species} mole fraction\")\nax_X.set(xlabel=r\"equivalence ratio, $\\phi$\", ylabel=\"mole fraction, $X$\")\nax_X.legend(fontsize=8)\nplt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.5" } }, "nbformat": 4, "nbformat_minor": 0 }