{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Timestep Controls for 1D Flame Solvers\n\nThis example demonstrates four solver controls that affect pseudo-transient\nconvergence in one-dimensional flame calculations:\n\n1. ``time_step_growth_strategy``, which controls when the timestep grows after a\n successful transient step.\n2. ``time_step_growth_factor``, which controls how much the timestep grows when\n growth is triggered.\n3. ``time_step_regrid``, which controls how many times the solver may refine the\n grid and retry after timestepping is exhausted.\n4. ``max_time_step_count``, which controls the number of transient steps allowed\n before a timestep attempt gives up.\n\nThe first comparison uses a burner-stabilized premixed flame. The following two\ncomparisons use a high-pressure counterflow diffusion flame, where solver convergence\ncan depend on both regridding and the transient-step budget.\n\nRequires: cantera >= 4.0; matplotlib >= 2.0\n\n.. tags:: Python, combustion, 1D flow, premixed flame, diffusion flame,\n strained flame, plotting\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import annotations\n\nimport time\n\nimport matplotlib.pyplot as plt\nplt.rcParams['figure.constrained_layout.use'] = True\nimport numpy as np\n\nimport cantera as ct" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Timestep Growth Strategies\n\nThe one-dimensional solver falls back to pseudo-transient timestepping when the Newton\nsolver cannot converge directly to steady state. After a successful transient step,\nCantera can either grow the timestep by a fixed factor or apply one of several\nheuristics that wait for evidence that the larger step is helpful.\n\nThis first problem uses a premixed methane / air burner-stabilized flame. The physical\nsolution is the same for each case, so the differences in the table are from the\ntimestep-growth rule. We test several combinations of growth strategies and growth\nfactors.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "growth_strategies = (\n (\"fixed-growth\", 1.5), # default growth strategy and factor\n (\"fixed-growth\", 2.0),\n (\"steady-norm\", 2.0),\n (\"transient-residual\", 2.0),\n (\"residual-ratio\", 2.0),\n (\"newton-iterations\", 2.0),\n (\"steady-norm\", 6.0),\n (\"newton-iterations\", 6.0),\n)\n\ndef make_growth_flame() -> ct.BurnerFlame:\n gas = ct.Solution(\"gri30.yaml\")\n gas.set_equivalence_ratio(0.8, \"CH4\", {\"O2\": 1.0, \"N2\": 3.76})\n gas.TP = 300.0, ct.one_atm\n\n flame = ct.BurnerFlame(gas, width=0.02)\n flame.burner.mdot = 0.09\n flame.set_refine_criteria(ratio=3.0, slope=0.25, curve=0.35, prune=0.05)\n return flame\n\n\ndef run_growth_strategy(strategy: str, growth_factor: float) -> dict[str, object]:\n flame = make_growth_flame()\n timesteps = []\n flame.set_time_step_callback(\n lambda dt: timesteps.append(float(dt)) or 0\n )\n flame.time_step_growth_factor = growth_factor\n flame.time_step_growth_strategy = strategy\n\n tic = time.perf_counter()\n print(\"** Running growth strategy:\", strategy)\n flame.solve(loglevel=1, refine_grid=True)\n elapsed = time.perf_counter() - tic\n\n return {\n \"strategy\": f\"{strategy} (growth factor = {growth_factor:.1f})\",\n \"n_steps\": len(timesteps),\n \"dt_sum\": float(np.sum(timesteps)),\n \"dt_min\": float(np.min(timesteps)),\n \"dt_max\": float(np.max(timesteps)),\n \"jacobians\": int(sum(flame.jacobian_count_stats)),\n \"T_max\": float(np.max(flame.T)),\n \"wall_time\": elapsed,\n \"dts\": np.asarray(timesteps, dtype=float),\n }\n\n\ngrowth_results = [run_growth_strategy(strategy, growth_factor)\n for strategy, growth_factor in growth_strategies]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The timestep history is collected with ``set_time_step_callback``.\nThe columns below report the number of timesteps, the total transient\ntime covered by those steps, the largest timestep, the number of\nJacobian evaluations, and the maximum temperature of the converged flame.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "growth_header = (\n f\"{'Strategy':<42} {'Steps':>8} {'sum(dt) [s]':>12} \"\n f\"{'max(dt) [s]':>12} {'Jac':>5} {'T_max [K]':>12} {'Runtime [s]':>9}\"\n)\nprint(growth_header)\nprint(\"-\" * len(growth_header))\nfor row in growth_results:\n print(\n f\"{row['strategy']:<42} {row['n_steps']:>8d} \"\n f\"{row['dt_sum']:>12.4e} {row['dt_max']:>12.4e} \"\n f\"{row['jacobians']:>5d} {row['T_max']:>12.4f} {row['wall_time']:>9.2f}\"\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In short, ``fixed-growth`` always applies ``time_step_growth_factor`` after a\nsuccessful step. ``steady-norm`` and ``transient-residual`` require the\ncorresponding residual norm to decrease. ``residual-ratio`` scales the growth\nusing the transient residual improvement, capped by\n``time_step_growth_factor``. ``newton-iterations`` grows the step only after a\nlow-iteration Newton solve. Setting ``time_step_growth_factor = 1.0`` disables\ngrowth.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Timestep Histories\n\nThe timestep histories show how the strategy changes the sequence of transient\nsteps even when all cases reach the same steady fixed-temperature solution.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cmap = plt.get_cmap('Dark2')\n\nfig, ax = plt.subplots()\nfor i,row in enumerate(growth_results):\n steps = np.arange(1, row[\"n_steps\"] + 1)\n ax.semilogy(\n steps + 0.1*i, row[\"dts\"], \"o-\", linewidth=1.6, markersize=4,\n color=cmap(i), label=row[\"strategy\"]\n )\nax.set_xlabel(\"timestep number\")\nax.set_ylabel(\"dt [s]\")\nax.grid(True, which=\"both\", alpha=0.3)\nax.legend(fontsize=8)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Regridding After Timestep Failure\n\nThe next comparison uses a high-pressure hydrogen / oxygen counterflow\ndiffusion flame. The initial grid has only 20 points across a 30 mm domain, which is\ntoo coarse for the steady-state solver to converge. The ``time_step_regrid`` option\nlets the solver refine the grid and retry when a timestepping sequence has reached\n``max_time_step_count``.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "diffusion_width = 30e-3\ndiffusion_initial_points = 20\ndiffusion_fuel_mdot = 0.3\ndiffusion_oxidizer_mdot_factor = 3.0\n\n\ndef make_regrid_flame(\n pressure: float,\n regrid_max: int,\n max_time_step_count: int = 200,\n growth_strategy: str = \"fixed-growth\",\n growth_factor: float = 1.5,\n) -> ct.CounterflowDiffusionFlame:\n gas = ct.Solution(\"h2o2.yaml\")\n flame = ct.CounterflowDiffusionFlame(\n gas, grid=np.linspace(0.0, diffusion_width, diffusion_initial_points)\n )\n flame.max_time_step_count = max_time_step_count\n flame.set_refine_criteria(ratio=2.0, slope=0.06, curve=0.08, prune=0.02)\n flame.time_step_regrid = regrid_max\n flame.time_step_growth_strategy = growth_strategy\n flame.time_step_growth_factor = growth_factor\n\n flame.P = pressure\n flame.fuel_inlet.X = \"H2:1.0\"\n flame.fuel_inlet.T = 800.0\n flame.oxidizer_inlet.X = \"O2:1.0\"\n flame.oxidizer_inlet.T = 711.0\n\n rho_fuel = flame.fuel_inlet.phase.density\n rho_oxidizer = flame.oxidizer_inlet.phase.density\n flame.fuel_inlet.mdot = diffusion_fuel_mdot\n flame.oxidizer_inlet.mdot = (\n diffusion_fuel_mdot / rho_fuel * rho_oxidizer * diffusion_oxidizer_mdot_factor\n )\n return flame\n\n\ndef run_regrid_case(case: dict[str, object]) -> dict[str, object]:\n flame = make_regrid_flame(\n pressure=case[\"pressure\"],\n regrid_max=case[\"regrid_max\"],\n max_time_step_count=case[\"max_time_step_count\"],\n growth_strategy=case[\"growth_strategy\"],\n growth_factor=case[\"growth_factor\"],\n )\n flame.clear_stats()\n\n tic = time.perf_counter()\n try:\n flame.solve(loglevel=0, auto=False)\n success = True\n except ct.CanteraError as err:\n success = False\n elapsed = time.perf_counter() - tic\n\n return {\n \"label\": case[\"label\"],\n \"plot_label\": case[\"plot_label\"],\n \"pressure\": case[\"pressure\"],\n \"success\": success,\n \"grid_points\": len(flame.grid),\n \"timesteps\": int(sum(flame.time_step_stats)),\n \"jacobians\": int(sum(flame.jacobian_count_stats)),\n \"peak_temperature\": float(np.max(flame.T)) if success else None,\n \"grid\": np.array(flame.grid, copy=True),\n \"temperature\": np.array(flame.T, copy=True),\n \"wall_time\": elapsed,\n \"max_time_step_count\": case[\"max_time_step_count\"],\n \"regrid_max\": case[\"regrid_max\"],\n \"growth_strategy\": case[\"growth_strategy\"],\n \"growth_factor\": case[\"growth_factor\"],\n }\n\n\ndef print_regrid_rows(results: list[dict[str, object]], label_width: int = 28) -> None:\n header = (\n f\"{'Mode':<{label_width}} {'Success':>7} {'Grid':>6} \"\n f\"{'Steps':>7} {'Jac':>5} {'Tmax [K]':>10} {'Wall [s]':>9}\"\n )\n print(header)\n print(\"-\" * len(header))\n for row in results:\n peak_temperature = (\n f\"{row['peak_temperature']:>10.1f}\"\n if row[\"peak_temperature\"] is not None else f\"{'--':>10}\"\n )\n print(\n f\"{row['label']:<{label_width}} {str(row['success']):>7} \"\n f\"{row['grid_points']:>6d} {row['timesteps']:>7d} \"\n f\"{row['jacobians']:>5d} {peak_temperature} {row['wall_time']:>9.2f}\"\n )\n\n\nregrid_pressure = 7e6\nregrid_cases = (\n {\n \"label\": \"time_step_regrid = 0\",\n \"plot_label\": \"0\",\n \"pressure\": regrid_pressure,\n \"regrid_max\": 0,\n \"max_time_step_count\": 200,\n \"growth_strategy\": \"fixed-growth\",\n \"growth_factor\": 1.5,\n },\n {\n \"label\": \"time_step_regrid = 1\",\n \"plot_label\": \"1\",\n \"pressure\": regrid_pressure,\n \"regrid_max\": 1,\n \"max_time_step_count\": 200,\n \"growth_strategy\": \"fixed-growth\",\n \"growth_factor\": 1.5,\n },\n {\n \"label\": \"time_step_regrid = 3\",\n \"plot_label\": \"3\",\n \"pressure\": regrid_pressure,\n \"regrid_max\": 3,\n \"max_time_step_count\": 200,\n \"growth_strategy\": \"fixed-growth\",\n \"growth_factor\": 1.5,\n },\n)\n\nregrid_results = [run_regrid_case(case) for case in regrid_cases]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here the growth-factor and growth-strategy settings are held at their defaults. Only\nthe number of regrid retries changes.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(\"Regrid retries after timestep failure\")\nprint(\n f\"Case: H2/O2 counterflow diffusion flame at {regrid_pressure / 1e6:.1f} MPa \"\n f\"on an initial {diffusion_initial_points}-point grid\"\n)\nprint(\"Varied option: time_step_regrid\")\nprint()\nprint_regrid_rows(regrid_results, label_width=43)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "``time_step_regrid = 0`` disables retry-after-regrid. Positive values allow up to that\nmany grid-refinement retries after timestep failure. If the refinement criteria do not\nchange the grid, the retry path aborts because there is no new discretization to try.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Regrid Retry Comparison\n\nThe first two cases exit before reaching a steady solution. With three regrid\nretries, the solver has enough opportunities to add points and continue.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "colors = [\"#ae2012\" if not row[\"success\"] else \"#0a9396\" for row in regrid_results]\nfig, ax = plt.subplots()\nbars = ax.bar(\n [row[\"plot_label\"] for row in regrid_results],\n [row[\"timesteps\"] for row in regrid_results],\n color=colors,\n)\nax.set_xlabel(\"time_step_regrid\")\nax.set_ylabel(\"timesteps before exit\")\nax.grid(True, axis=\"y\", alpha=0.3)\nfor bar, row in zip(bars, regrid_results):\n status = \"ok\" if row[\"success\"] else \"fail\"\n ax.text(\n bar.get_x() + bar.get_width() / 2,\n bar.get_height() + 10,\n f\"grid={row['grid_points']}\\n{status}\",\n ha=\"center\",\n va=\"bottom\",\n fontsize=8,\n )\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Flame Profile\n\nThe successful 7 MPa case is not just a solver-status change: after regridding,\nthe solver finds a steady-state solution for a resolved diffusion flame with a hot\nreaction zone in the counterflow domain.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "regrid_result = next(\n row for row in regrid_results if row[\"success\"] and row[\"regrid_max\"] == 3\n)\n\nfig, ax = plt.subplots()\nax.plot(\n 1e3 * regrid_result[\"grid\"],\n regrid_result[\"temperature\"],\n \".-\",\n)\nax.set_xlabel(\"distance from fuel inlet [mm]\")\nax.set_ylabel(\"temperature [K]\")\nax.grid(True, alpha=0.3)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A More Difficult Pressure\n\nAt 15 MPa, the same coarse counterflow problem is more demanding. Regridding is\nenabled for all cases in this section, but a 200-step timestep budget is still\ninsufficient. The first comparison shows the effect of increasing\n``max_time_step_count`` while leaving the growth strategy fixed.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "difficult_pressure = 15e6\ndifficult_budget_cases = (\n {\n \"label\": \"max_time_step_count = 200\",\n \"plot_label\": \"200\",\n \"pressure\": difficult_pressure,\n \"max_time_step_count\": 200,\n \"regrid_max\": 3,\n \"growth_strategy\": \"fixed-growth\",\n \"growth_factor\": 1.5,\n },\n {\n \"label\": \"max_time_step_count = 1250\",\n \"plot_label\": \"1250\",\n \"pressure\": difficult_pressure,\n \"max_time_step_count\": 1250,\n \"regrid_max\": 3,\n \"growth_strategy\": \"fixed-growth\",\n \"growth_factor\": 1.5,\n },\n)\n\ndifficult_budget_results = [run_regrid_case(case) for case in difficult_budget_cases]\n\nprint(\"Solving with a larger timestep budget\")\nprint(\n f\"Case: H2/O2 counterflow diffusion flame at \"\n f\"{difficult_pressure / 1e6:.1f} MPa\"\n)\nprint(\"Fixed settings: time_step_regrid = 3, time_step_growth_strategy = fixed-growth\")\nprint()\nprint_regrid_rows(difficult_budget_results, label_width=28)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once the timestep budget is large enough enable solution of the steady-state problem,\nthe growth strategy can still affect the work required to reach the converged steady\nsolution. The fixed-growth row below reuses the successful result from the previous\ntable, and the remaining rows rerun the same problem with adaptive growth strategies.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "timestep_budget_result = next(\n row for row in difficult_budget_results if row[\"max_time_step_count\"] == 1250\n)\ndifficult_strategy_cases = (\n {\n \"label\": \"residual-ratio\",\n \"plot_label\": \"ratio\",\n \"pressure\": difficult_pressure,\n \"max_time_step_count\": 1250,\n \"regrid_max\": 3,\n \"growth_strategy\": \"residual-ratio\",\n \"growth_factor\": 1.5,\n },\n {\n \"label\": \"newton-iterations\",\n \"plot_label\": \"newton\",\n \"pressure\": difficult_pressure,\n \"max_time_step_count\": 1250,\n \"regrid_max\": 3,\n \"growth_strategy\": \"newton-iterations\",\n \"growth_factor\": 1.5,\n },\n)\n\ndifficult_strategy_results = [\n {\n **timestep_budget_result,\n \"label\": \"fixed-growth\",\n \"plot_label\": \"fixed\",\n },\n *[run_regrid_case(case) for case in difficult_strategy_cases],\n]\n\nprint(\"Growth-strategy sensitivity after recovery\")\nprint(\n \"Fixed settings: max_time_step_count = 1250, time_step_regrid = 3, \"\n \"time_step_growth_factor = 1.5\"\n)\nprint()\nprint_regrid_rows(difficult_strategy_results, label_width=18)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "``max_time_step_count`` limits the number of transient steps before a timestep attempt\ngives up. Larger values can allow harder cases to converge, but solving these cases\nwill be slow. On this 15 MPa case, adaptive growth strategies can reduce wall time\nonce the timestep budget is large enough.\n\n" ] } ], "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.6" } }, "nbformat": 4, "nbformat_minor": 0 }