Skip to content

Timestep regridding for 1D solver#2074

Open
wandadars wants to merge 14 commits into
Cantera:mainfrom
wandadars:timestep_regridding
Open

Timestep regridding for 1D solver#2074
wandadars wants to merge 14 commits into
Cantera:mainfrom
wandadars:timestep_regridding

Conversation

@wandadars
Copy link
Copy Markdown
Contributor

@wandadars wandadars commented Jan 11, 2026

During a discussion on the users group, Ray mentioned that Cantera doesn't perform any regridding if timestepping fails during a solve. This pull request is an attempt to add a capability to the 1D solver that would allow for a regrid to happen if there is a failure in the timestepping. This feature allows for a set number of regridding steps to be taken during a failure of the timestepping. The test case that we were looking at was that of a pressure field being stepped up incrementally on a counterflow diffusion flame domain. The highest pressure possible was increased when the timestep regridding is allowed. This happens when cases have a bad initial condition or a coarse grid and fails right at the start during the timestepping phase.

This is Ray's script that I tweaked a bit and used to test the regridding.

Details
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt

mech_path = 'h2o2.yaml'


# Inlet parameters
fuel_X = "H2:1.0"
oxidizer_X = "O2:1.0"
fuel_T = 800.0
oxidizer_T = 711.0
mdot_fuel = 0.3  # kg/m^2/s
width = 30e-3
P_start = 10e6
P_end = 60e6
n_pressures = 100
pressure_levels = np.linspace(P_start, P_end, n_pressures)

oxidizer_mdot_factor = 3.0
n_points = 50

# Adaptive stepping with backtracking on failures
loglevel = 1

gas = ct.Solution(mech_path)
#gas = ct.Solution(mech_path, 'ohmech-RK')

def make_flame(gas):
    flame = ct.CounterflowDiffusionFlame(gas, grid=np.linspace(0, width, n_points))
    flame.max_time_step_count = 1000
    #flame.set_time_step(1.0e-5, [10,20,30,50])

    flame.set_refine_criteria(ratio=2.0, slope=0.06, curve=0.08, prune=0.02)
    flame.fuel_inlet.mdot = mdot_fuel
    flame.fuel_inlet.X = fuel_X
    flame.fuel_inlet.T = fuel_T
    flame.oxidizer_inlet.X = oxidizer_X
    flame.oxidizer_inlet.T = oxidizer_T
    return flame

def update_inlet_mdot(flame):
    rho_f = flame.fuel_inlet.phase.density
    rho_o = flame.oxidizer_inlet.phase.density
    flame.fuel_inlet.mdot = mdot_fuel
    flame.oxidizer_inlet.mdot = (mdot_fuel / rho_f) * rho_o * oxidizer_mdot_factor

flame = make_flame(gas)
flame.set_time_step_regrid(15)

def plot_diagnostics(flame, pressure, step_index, total_pressures):
    z = flame.grid
    temperature = flame.T
    velocity = flame.velocity
    spread_rate = flame.spread_rate
    heat_release = flame.heat_release_rate
    density = flame.density
    mixture_fraction = flame.mixture_fraction('Bilger')

    fig, axes = plt.subplots(3, 2, figsize=(14, 10), sharex=True)
    axes = axes.ravel()

    line_kwargs = dict(linewidth=2.0, marker='o', markersize=3,
                       markerfacecolor='none', markeredgewidth=0.6)
    axes[0].plot(z, temperature, **line_kwargs)
    axes[0].set_title('Temperature')
    axes[0].set_ylabel('T [K]')

    axes[1].plot(z, heat_release, **line_kwargs)
    axes[1].set_title('Heat Release Rate')
    axes[1].set_ylabel('HRR [W/m^3]')

    axes[2].plot(z, velocity, **line_kwargs)
    axes[2].set_title('Velocity')
    axes[2].set_ylabel('u [m/s]')

    axes[3].plot(z, spread_rate, **line_kwargs)
    axes[3].set_title('Spread Rate')
    axes[3].set_ylabel('V [1/s]')

    axes[4].plot(z, density, **line_kwargs)
    axes[4].set_title('Density')
    axes[4].set_ylabel('rho [kg/m^3]')

    axes[5].plot(z, mixture_fraction, **line_kwargs)
    axes[5].set_title('Bilger Mixture Fraction')
    axes[5].set_ylabel('Z [-]')

    for ax in axes:
        ax.grid(True, alpha=0.3)
    for ax in axes[-2:]:
        ax.set_xlabel('z [m]')

    fig.suptitle(f'Diagnostics at P = {pressure:.3e} Pa '
                 f'({step_index}/{total_pressures})')
    fig.tight_layout(rect=[0, 0.03, 1, 0.95])
    plt.show()
    plt.close(fig)

pressure_outfile = 'pressure_solutions.h5'
total_pressures = len(pressure_levels)
for step_index, pressure in enumerate(pressure_levels, start=1):
    print(f'\n=== Solving at P = {pressure:.3e} Pa ({step_index}/{total_pressures}) ===')
    flame.P = pressure
    flame.fuel_inlet.X = fuel_X
    flame.fuel_inlet.T = fuel_T
    flame.oxidizer_inlet.X = oxidizer_X
    flame.oxidizer_inlet.T = oxidizer_T
    update_inlet_mdot(flame)
    flame.solve(loglevel=loglevel, auto=False)
    flame.save(pressure_outfile, overwrite=True)
    print(f'Saved solution to {pressure_outfile}')
    plot_diagnostics(flame, pressure, step_index, total_pressures)

print(f'Number of debug outputs = {i}')
flame.show_stats()

AI Statement (required)

-Extensive use of generative AI.
Significant portions of code or documentation were generated with AI, including
logic and implementation decisions. All generated code and documentation were
reviewed and understood by the contributor. Examples: Output from agentic coding
tools and/or substantial refactoring by LLMs (web-based or local).

Checklist

  • The pull request includes a clear description of this code change
  • Commit messages have short titles and reference relevant issues
  • Build passes (scons build & scons test) and unit tests address code coverage
  • Style & formatting of contributed code follows contributing guidelines
  • AI Statement is included
  • The pull request is ready for review

@wandadars
Copy link
Copy Markdown
Contributor Author

I don't know how to edit the PR header, but this isn't a WIP anymore.

@speth speth changed the title WIP: Timestep regridding for 1D solver Timestep regridding for 1D solver Feb 3, 2026
@speth
Copy link
Copy Markdown
Member

speth commented Feb 11, 2026

Based on the tests that are failing and timing out, it seems like this is affecting the solution of steady-state reactor networks, which also use the SteadyStateSystem solver, even though those systems don't use regridding.

@codecov
Copy link
Copy Markdown

codecov Bot commented Feb 12, 2026

Codecov Report

❌ Patch coverage is 70.00000% with 33 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.71%. Comparing base (44732e7) to head (b797be0).

Files with missing lines Patch % Lines
src/oneD/Sim1D.cpp 23.52% 12 Missing and 1 partial ⚠️
src/numerics/SteadyStateSystem.cpp 83.58% 10 Missing and 1 partial ⚠️
include/cantera/numerics/SteadyStateSystem.h 60.00% 4 Missing ⚠️
include/cantera/oneD/Sim1D.h 33.33% 3 Missing and 1 partial ⚠️
interfaces/cython/cantera/_onedim.pyx 83.33% 1 Missing ⚠️
Additional details and impacted files
@@           Coverage Diff            @@
##             main    #2074    +/-   ##
========================================
  Coverage   77.71%   77.71%            
========================================
  Files         452      452            
  Lines       53316    53421   +105     
  Branches     8894     8909    +15     
========================================
+ Hits        41432    41516    +84     
- Misses       8870     8888    +18     
- Partials     3014     3017     +3     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Copy link
Copy Markdown
Member

@speth speth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for putting this together, @wandadars. I think there are some useful ideas here. What we need now is something to show that. If there are any very simple cases where this can be exercised, it would be great to add them to the test suite. Besides that, we should add a new example demonstrating these options and cases where they have a beneficial impact on solver convergence (and perhaps when they don't).

Comment thread include/cantera/numerics/SteadyStateSystem.h
Comment thread include/cantera/numerics/SteadyStateSystem.h Outdated
Comment thread include/cantera/numerics/SteadyStateSystem.h Outdated
Comment thread include/cantera/numerics/SteadyStateSystem.h Outdated
Comment thread include/cantera/numerics/SteadyStateSystem.h Outdated
Comment thread src/numerics/SteadyStateSystem.cpp Outdated
Comment thread src/numerics/SteadyStateSystem.cpp Outdated
Comment thread src/oneD/Sim1D.cpp Outdated
Comment thread src/oneD/Sim1D.cpp Outdated
Comment thread test/python/test_onedim.py Outdated
@speth
Copy link
Copy Markdown
Member

speth commented Mar 25, 2026

Thanks for the updates, @wandadars. I just wanted to check, is this ready for another review?

@wandadars
Copy link
Copy Markdown
Contributor Author

Yes @speth I think so. Sorry for the delay.

@speth speth self-assigned this Apr 27, 2026
@speth
Copy link
Copy Markdown
Member

speth commented May 5, 2026

I combined the examples and refactored them significantly for style. After that, I spent some time playing with the different timestep growth strategies implemented. Here is what they do for a more complex diffusion flame than the one in the original example (now including regridding). For all the examples I've tried, it seems like most of the new strategies give more or less the same behavior as the current "fixed-growth" option, with the exception of the "newton-iterations" approach which is often significantly worse in terms of run time:

Named timestep-growth strategies
Case: C2H6 diffusion flame at 300 K fuel, 700 K oxidizer, 5 atm, width = 4 cm
Configured time_step_growth_factor = 2.0

Strategy             Accepted  sum(dt) [s]  max(dt) [s]   Jac     u0 [m/s] Runtime [s]
--------------------------------------------------------------------------------------
fixed-growth              120   2.8707e+00   1.3107e+00   121 0.2947188407      9.03
steady-norm               140   1.0877e+00   1.6384e-01   123 0.2947188407      9.53
transient-residual        120   2.8707e+00   1.3107e+00   121 0.2947188407      9.11
residual-ratio            120   1.1411e+00   3.0358e-01   121 0.2947188407      9.11
newton-iterations         260   1.6747e-01   5.1200e-03   141 0.2947188407     10.78
image

Are there some other cases you saw more dramatic benefits from implementing these alternative time step growth strategies?

@speth speth force-pushed the timestep_regridding branch from 655642e to 7879319 Compare May 5, 2026 15:00
@wandadars
Copy link
Copy Markdown
Contributor Author

IIRC, I think the only one that really seemed to have any improvement was the high pressure case. But I have since changed the approach that I was using to get the high pressure solutions from stepping up solutions and using the last one as the initial condition to using a similarity solution for the initial condition, which had much better performance. So that attempt to squeeze out any speed by using different heuristics was a bit of a moot point for the application that I was using to judge them. I'm all for axing them if they're not showing any improvements for the various cases that Cantera considers.

@wandadars
Copy link
Copy Markdown
Contributor Author

I went back and ran some cases with larger mechanisms like the gri30 one, and it does seem like at least some subset of those heuristics seem to give a performance improvement.

I may be missing some important context here that makes these results not as difference as they seem at face value here. I'll provide the script that generated them for context. If anything the fixed, steady norm, and the newton iterations variants seem to show promising results. The others, not so much.

Details

#!/usr/bin/env python3
"""
One-time benchmark for 1D flame time-step growth heuristics.

This script is intentionally more expensive than the example-sized runner. It
tries to answer whether the named time-step growth heuristics behave
meaningfully differently from fixed growth:

  1. A fast H2/O2 counterflow-premixed flame sweeps growth factors.
  2. A methane/air GRI30 burner flame uses a larger chemistry mechanism so that
    excess Jacobian work shows up clearly in wall time.
    """

from future import annotations

import argparse
import csv
import json
import math
import time
from dataclasses import dataclass
from pathlib import Path
from typing import Callable

import cantera as ct
import matplotlib
import numpy as np

matplotlib.use("Agg")
import matplotlib.pyplot as plt

SCRIPT_DIR = Path(file).resolve().parent
DEFAULT_OUTPUT_DIR = SCRIPT_DIR / "flame_timestep_heuristic_benchmark_outputs"

@DataClass(frozen=True)
class Strategy:
label: str
option: str

@DataClass(frozen=True)
class Scenario:
name: str
title: str
description: str
mechanism: str
make_flame: Callable[[], ct.FlameBase]
solve_flame: Callable[[ct.FlameBase, int], None]
growth_factors: tuple[float, ...]
suites: tuple[str, ...]

STRATEGIES = (
Strategy("fixed-growth", "fixed-growth"),
Strategy("steady-norm", "steady-norm"),
Strategy("transient-residual", "transient-residual"),
Strategy("residual-ratio", "residual-ratio"),
Strategy("newton-iterations", "newton-iterations"),
)

def cantera_datafile(name: str) -> str:
for directory in ct.get_data_directories():
if directory == ".":
continue
candidate = Path(directory) / name
if candidate.exists():
return str(candidate)
return name

H2O2 = cantera_datafile("h2o2.yaml")
GRI30 = cantera_datafile("gri30.yaml")

def make_h2o2_counterflow_premixed() -> ct.CounterflowPremixedFlame:
gas = ct.Solution(H2O2)
gas.TPX = 373.0, 0.05 * ct.one_atm, "H2:1.6, O2:1, AR:7"

flame = ct.CounterflowPremixedFlame(gas, width=0.08)
flame.reactants.mdot = 0.12
flame.products.mdot = 0.06
flame.set_refine_criteria(ratio=3.0, slope=0.12, curve=0.2, prune=0.03)
flame.set_initial_guess()
return flame

def make_gri30_methane_burner() -> ct.BurnerFlame:
gas = ct.Solution(GRI30)
gas.set_equivalence_ratio(0.8, "CH4", {"O2": 1.0, "N2": 3.76})
gas.TP = 300.0, ct.one_atm

flame = ct.BurnerFlame(gas, grid=np.linspace(0.0, 0.03, 40))
flame.burner.mdot = 0.04
flame.set_refine_criteria(ratio=3.0, slope=0.25, curve=0.35, prune=0.05)
return flame

def make_gri30_methane_burner_refined() -> ct.BurnerFlame:
gas = ct.Solution(GRI30)
gas.set_equivalence_ratio(0.8, "CH4", {"O2": 1.0, "N2": 3.76})
gas.TP = 300.0, ct.one_atm

flame = ct.BurnerFlame(gas, grid=np.linspace(0.0, 0.03, 60))
flame.burner.mdot = 0.04
flame.set_refine_criteria(ratio=3.0, slope=0.15, curve=0.2, prune=0.03)
return flame

def make_gri30_twin_premixed() -> ct.CounterflowTwinPremixedFlame:
gas = ct.Solution(GRI30)
gas.set_equivalence_ratio(0.75, "CH4", {"O2": 1.0, "N2": 3.76})
gas.TP = 300.0, ct.one_atm

flame = ct.CounterflowTwinPremixedFlame(gas, width=0.025)
flame.reactants.mdot = gas.density * 2.0
flame.set_refine_criteria(ratio=4.0, slope=0.35, curve=0.45, prune=0.08)
flame.set_initial_guess()
return flame

def solve_auto(flame: ct.FlameBase, loglevel: int) -> None:
flame.solve(loglevel=loglevel, auto=True)

SCENARIOS = (
Scenario(
"h2o2_counterflow_premixed",
"H2/O2 Counterflow Premixed Growth Sweep",
"small mechanism, same flame, increasingly aggressive growth factor",
H2O2,
make_h2o2_counterflow_premixed,
solve_auto,
(2.0, 5.0, 10.0, 20.0, 50.0),
("standard",),
),
Scenario(
"gri30_methane_burner",
"GRI30 Methane Burner Stress Case",
"larger mechanism where extra Jacobian work is visible in wall time",
GRI30,
make_gri30_methane_burner,
solve_auto,
(20.0,),
("standard", "long"),
),
Scenario(
"gri30_methane_burner_sweep",
"GRI30 Methane Burner Growth Sweep",
"same methane burner as the stress case, swept over aggressive growth factors",
GRI30,
make_gri30_methane_burner,
solve_auto,
(2.0, 5.0, 10.0, 20.0),
("long",),
),
Scenario(
"gri30_methane_burner_refined",
"GRI30 Methane Burner Refined-Grid Stress Case",
"larger initial grid and tighter refinement; intended as the heavier proof case",
GRI30,
make_gri30_methane_burner_refined,
solve_auto,
(20.0,),
("long",),
),
Scenario(
"gri30_twin_premixed",
"GRI30 Methane Twin Premixed Stress Case",
"premixed methane/air opposed-flow flame; optional long-running cross-check",
GRI30,
make_gri30_twin_premixed,
solve_auto,
(20.0,),
("long",),
),
)

def sum_stats(values: list[int] | list[float]) -> float:
return float(sum(values)) if values else 0.0

def finite(value: float | None) -> bool:
return value is not None and math.isfinite(value)

def growth_tag(growth_factor: float) -> str:
return f"growth{growth_factor:g}".replace(".", "p")

def run_one(
scenario: Scenario,
strategy: Strategy,
growth_factor: float,
loglevel: int,
) -> tuple[dict[str, object], dict[str, np.ndarray]]:
flame = scenario.make_flame()
flame.time_step_growth_strategy = strategy.option
flame.time_step_growth_factor = growth_factor
flame.time_step_regrid = 3
flame.max_time_step_count = 500

accepted_dts: list[float] = []
flame.set_time_step_callback(lambda dt: accepted_dts.append(float(dt)) or 0)
flame.clear_stats()

tic = time.perf_counter()
success = True
message = ""
try:
    scenario.solve_flame(flame, loglevel)
except ct.CanteraError as err:
    success = False
    message = str(err).splitlines()[0] if str(err) else type(err).__name__
wall_time = time.perf_counter() - tic

dts = np.asarray(accepted_dts, dtype=float)
grid = np.asarray(flame.grid, dtype=float)
temperature = np.asarray(flame.T, dtype=float)

row: dict[str, object] = {
    "scenario": scenario.name,
    "scenario_title": scenario.title,
    "strategy": strategy.label,
    "growth_factor": growth_factor,
    "success": success,
    "accepted_steps": int(dts.size),
    "dt_sum_s": float(np.sum(dts)) if dts.size else None,
    "dt_min_s": float(np.min(dts)) if dts.size else None,
    "dt_median_s": float(np.median(dts)) if dts.size else None,
    "dt_max_s": float(np.max(dts)) if dts.size else None,
    "time_step_stats_total": int(sum_stats(flame.time_step_stats)),
    "jacobians": int(sum_stats(flame.jacobian_count_stats)),
    "evals": int(sum_stats(flame.eval_count_stats)),
    "grid_points": int(grid.size),
    "peak_temperature_K": float(np.max(temperature)) if temperature.size else None,
    "wall_time_s": wall_time,
    "message": message,
}
arrays = {"dts": dts, "grid": grid, "temperature": temperature}
return row, arrays

def row_key(row: dict[str, object]) -> tuple[str, str, float]:
return (str(row["scenario"]), str(row["strategy"]), float(row["growth_factor"]))

def write_csv(path: Path, rows: list[dict[str, object]]) -> None:
fields = [
"scenario",
"scenario_title",
"strategy",
"growth_factor",
"success",
"accepted_steps",
"dt_sum_s",
"dt_min_s",
"dt_median_s",
"dt_max_s",
"time_step_stats_total",
"jacobians",
"evals",
"grid_points",
"peak_temperature_K",
"wall_time_s",
"message",
]
with path.open("w", newline="", encoding="utf-8") as handle:
writer = csv.DictWriter(handle, fieldnames=fields)
writer.writeheader()
writer.writerows(rows)

def format_value(value: object, digits: int = 4) -> str:
if value is None:
return ""
if isinstance(value, float):
if not math.isfinite(value):
return ""
return f"{value:.{digits}g}"
return str(value)

def markdown_table(headers: list[str], rows: list[list[object]]) -> str:
rendered = [[format_value(value) for value in row] for row in rows]
widths = [
max(len(headers[i]), *(len(row[i]) for row in rendered))
for i in range(len(headers))
]
lines = [
"| " + " | ".join(headers[i].ljust(widths[i]) for i in range(len(headers))) + " |",
"| " + " | ".join("-" * widths[i] for i in range(len(headers))) + " |",
]
for row in rendered:
lines.append("| " + " | ".join(row[i].ljust(widths[i]) for i in range(len(row))) + " |")
return "\n".join(lines)

def write_summary_md(path: Path, rows: list[dict[str, object]]) -> None:
table_rows = [
[
row["scenario"],
row["growth_factor"],
row["strategy"],
row["success"],
row["accepted_steps"],
row["dt_max_s"],
row["jacobians"],
row["grid_points"],
row["peak_temperature_K"],
row["wall_time_s"],
]
for row in rows
]
path.write_text(
markdown_table(
[
"Scenario",
"Growth",
"Strategy",
"OK",
"Steps",
"max(dt) [s]",
"Jac",
"Grid",
"Tmax [K]",
"Wall [s]",
],
table_rows,
)
+ "\n",
encoding="utf-8",
)

def write_timestep_history(path: Path, dts: np.ndarray) -> None:
with path.open("w", newline="", encoding="utf-8") as handle:
writer = csv.writer(handle)
writer.writerow(["accepted_step", "dt_s"])
for index, dt in enumerate(dts, start=1):
writer.writerow([index, f"{dt:.16e}"])

def write_temperature_profile(path: Path, grid: np.ndarray, temperature: np.ndarray) -> None:
with path.open("w", newline="", encoding="utf-8") as handle:
writer = csv.writer(handle)
writer.writerow(["grid_m", "temperature_K"])
for z, temp in zip(grid, temperature):
writer.writerow([f"{z:.16e}", f"{temp:.16e}"])

def plot_growth_sweep(
rows: list[dict[str, object]],
scenario_name: str,
output_path: Path,
) -> None:
fig, axes = plt.subplots(2, 2, figsize=(12, 8), constrained_layout=True)
metrics = (
("jacobians", "Jacobian evaluations", "linear"),
("wall_time_s", "Wall time [s]", "linear"),
("accepted_steps", "Accepted pseudo-steps", "linear"),
("dt_max_s", "max accepted dt [s]", "log"),
)
colors = plt.get_cmap("tab10").colors
markers = ("o", "s", "^", "D", "v")
scenario_rows = [row for row in rows if row["scenario"] == scenario_name]
scenario_title = (
str(scenario_rows[0]["scenario_title"]) if scenario_rows else scenario_name
)
growth_factors = sorted({float(row["growth_factor"]) for row in scenario_rows})

for axis, (field, ylabel, scale) in zip(axes.flat, metrics):
    for index, strategy in enumerate(STRATEGIES):
        strategy_rows = sorted(
            [row for row in scenario_rows if row["strategy"] == strategy.label],
            key=lambda row: float(row["growth_factor"]),
        )
        x = [float(row["growth_factor"]) for row in strategy_rows]
        y = [float(row[field]) for row in strategy_rows]
        axis.plot(
            x,
            y,
            marker=markers[index % len(markers)],
            color=colors[index % len(colors)],
            linewidth=1.7,
            label=strategy.label,
        )
    axis.set_xlabel("time_step_growth_factor")
    axis.set_ylabel(ylabel)
    axis.set_xscale("log")
    axis.set_xticks(growth_factors)
    axis.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
    axis.set_yscale(scale)
    axis.grid(True, which="both", alpha=0.3)
axes[0, 0].legend(fontsize=8)
fig.suptitle(f"{scenario_title}: Growth-Factor Sweep")
fig.savefig(output_path, dpi=200)
plt.close(fig)

def plot_strategy_bars(
rows: list[dict[str, object]],
scenario_name: str,
growth_factor: float,
output_path: Path,
) -> None:
selected = [
row for row in rows
if row["scenario"] == scenario_name and float(row["growth_factor"]) == growth_factor
]
selected.sort(key=lambda row: [strategy.label for strategy in STRATEGIES].index(str(row["strategy"])))

labels = [str(row["strategy"]) for row in selected]
x = np.arange(len(labels))
colors = plt.get_cmap("tab10").colors[:len(labels)]
fig, axes = plt.subplots(2, 2, figsize=(12, 8), constrained_layout=True)
metrics = (
    ("wall_time_s", "Wall time [s]", "linear"),
    ("jacobians", "Jacobian evaluations", "linear"),
    ("accepted_steps", "Accepted pseudo-steps", "linear"),
    ("dt_max_s", "max accepted dt [s]", "log"),
)
for axis, (field, ylabel, scale) in zip(axes.flat, metrics):
    y = [float(row[field]) for row in selected]
    bars = axis.bar(x, y, color=colors)
    axis.set_xticks(x, labels, rotation=25, ha="right")
    axis.set_ylabel(ylabel)
    axis.set_yscale(scale)
    axis.grid(True, axis="y", which="both", alpha=0.3)
    for bar, value in zip(bars, y):
        axis.text(
            bar.get_x() + bar.get_width() / 2,
            bar.get_height(),
            format_value(value, digits=3),
            ha="center",
            va="bottom",
            fontsize=8,
        )
title = selected[0]["scenario_title"] if selected else scenario_name
fig.suptitle(f"{title}: growth factor {growth_factor:g}")
fig.savefig(output_path, dpi=200)
plt.close(fig)

def plot_timestep_histories(
rows: list[dict[str, object]],
arrays_by_key: dict[tuple[str, str, float], dict[str, np.ndarray]],
scenario_name: str,
growth_factor: float,
output_path: Path,
) -> None:
selected = [
row for row in rows
if row["scenario"] == scenario_name and float(row["growth_factor"]) == growth_factor
]
selected.sort(key=lambda row: [strategy.label for strategy in STRATEGIES].index(str(row["strategy"])))
colors = plt.get_cmap("tab10").colors
markers = ("o", "s", "^", "D", "v")
all_dts = [
dt
for row in selected
for dt in arrays_by_key[row_key(row)]["dts"]
if dt > 0.0
]

fig, axes = plt.subplots(
    len(selected),
    1,
    figsize=(10, 9),
    sharex=True,
    sharey=True,
    constrained_layout=True,
)
if len(selected) == 1:
    axes = [axes]
for index, (axis, row) in enumerate(zip(axes, selected)):
    dts = arrays_by_key[row_key(row)]["dts"]
    if dts.size:
        axis.semilogy(
            np.arange(1, dts.size + 1),
            dts,
            marker=markers[index % len(markers)],
            color=colors[index % len(colors)],
            markersize=3,
            linewidth=1.5,
        )
    axis.text(
        0.015,
        0.82,
        f"{row['strategy']}\nsteps={row['accepted_steps']}, jac={row['jacobians']}",
        transform=axis.transAxes,
        ha="left",
        va="top",
        fontsize=8,
        bbox={"boxstyle": "round,pad=0.25", "fc": "white", "ec": "0.85"},
    )
    axis.grid(True, which="both", alpha=0.3)
    if all_dts:
        axis.set_ylim(0.7 * min(all_dts), 1.35 * max(all_dts))
axes[len(axes) // 2].set_ylabel("dt [s]")
axes[-1].set_xlabel("accepted pseudo-step number")
title = selected[0]["scenario_title"] if selected else scenario_name
fig.suptitle(f"{title}: accepted timestep histories, growth factor {growth_factor:g}")
fig.savefig(output_path, dpi=200)
plt.close(fig)

def normalized_rows(rows: list[dict[str, object]], scenario_name: str, growth_factor: float) -> list[list[object]]:
selected = [
row for row in rows
if row["scenario"] == scenario_name and float(row["growth_factor"]) == growth_factor
]
fixed = next(row for row in selected if row["strategy"] == "fixed-growth")
table = []
for row in selected:
table.append(
[
row["strategy"],
row["jacobians"],
float(row["jacobians"]) / float(fixed["jacobians"]),
row["wall_time_s"],
float(row["wall_time_s"]) / float(fixed["wall_time_s"]),
row["accepted_steps"],
row["dt_max_s"],
]
)
return table

def write_report(
path: Path,
rows: list[dict[str, object]],
scenarios: list[Scenario],
plots_dir: Path,
) -> None:
lines = [
"# Flame timestep heuristic benchmark",
"",
f"Cantera: {ct.__version__} from {ct.__file__}",
"",
"The comparisons below are normalized to fixed-growth within each scenario. "
"The target is not the largest pseudo-time step; the target is the same final "
"steady solution with less nonlinear/Jacobian work and less fragile growth.",
"",
]

for scenario in scenarios:
    scenario_rows = [row for row in rows if row["scenario"] == scenario.name]
    if not scenario_rows:
        continue

    lines.extend(
        [
            f"## {scenario.title}",
            "",
            scenario.description,
            "",
        ]
    )

    growth_factors = sorted({float(row["growth_factor"]) for row in scenario_rows})
    if len(growth_factors) > 1:
        sweep_path = plots_dir / f"{scenario.name}_growth_sweep.png"
        lines.extend([f"Growth sweep: `{sweep_path}`", ""])

    for growth_factor in growth_factors:
        bars_path = plots_dir / f"{scenario.name}_{growth_tag(growth_factor)}_bars.png"
        histories_path = (
            plots_dir / f"{scenario.name}_{growth_tag(growth_factor)}_timestep_histories.png"
        )
        lines.extend(
            [
                f"### Growth factor {growth_factor:g}",
                "",
                f"Bars: `{bars_path}`",
                "",
                f"Timestep histories: `{histories_path}`",
                "",
            ]
        )
        lines.append(
            markdown_table(
                [
                    "Strategy",
                    "Jac",
                    "Jac / fixed",
                    "Wall [s]",
                    "Wall / fixed",
                    "Steps",
                    "max(dt) [s]",
                ],
                normalized_rows(rows, scenario.name, growth_factor),
            )
        )
        lines.append("")

lines.append("")
path.write_text("\n".join(lines), encoding="utf-8")

def parse_args() -> argparse.Namespace:
parser = argparse.ArgumentParser(description=doc)
parser.add_argument(
"--output-dir",
type=Path,
default=DEFAULT_OUTPUT_DIR,
help=f"output directory (default: {DEFAULT_OUTPUT_DIR})",
)
parser.add_argument("--loglevel", type=int, default=0)
parser.add_argument(
"--suite",
choices=("standard", "long", "all"),
default="standard",
help="scenario suite to run when --scenario is not specified",
)
parser.add_argument(
"--scenario",
action="append",
choices=[scenario.name for scenario in SCENARIOS],
help="specific scenario to run; may be provided more than once",
)
parser.add_argument(
"--skip-gri30",
action="store_true",
help="only run the quick H2/O2 growth-factor sweep",
)
return parser.parse_args()

def main() -> None:
args = parse_args()
output_dir = args.output_dir
if args.output_dir == DEFAULT_OUTPUT_DIR and (args.suite != "standard" or args.scenario):
tag = args.suite if not args.scenario else "custom"
output_dir = DEFAULT_OUTPUT_DIR.with_name(f"{DEFAULT_OUTPUT_DIR.name}_{tag}")
output_dir = output_dir.resolve()
plots_dir = output_dir / "plots"
histories_dir = output_dir / "timestep_histories"
profiles_dir = output_dir / "temperature_profiles"
for directory in (output_dir, plots_dir, histories_dir, profiles_dir):
directory.mkdir(parents=True, exist_ok=True)

if args.scenario:
    requested = set(args.scenario)
    scenarios = [scenario for scenario in SCENARIOS if scenario.name in requested]
elif args.suite == "all":
    scenarios = list(SCENARIOS)
else:
    scenarios = [
        scenario for scenario in SCENARIOS if args.suite in scenario.suites
    ]
if args.skip_gri30:
    scenarios = [
        scenario for scenario in scenarios if "gri30" not in scenario.name
    ]
rows: list[dict[str, object]] = []
arrays_by_key: dict[tuple[str, str, float], dict[str, np.ndarray]] = {}

print(f"Cantera {ct.__version__}: {ct.__file__}")
print(f"Output directory: {output_dir}")
print("Scenarios: " + ", ".join(scenario.name for scenario in scenarios))

for scenario in scenarios:
    print()
    print(f"Running {scenario.title}")
    for growth_factor in scenario.growth_factors:
        print(f"  growth factor {growth_factor:g}")
        for strategy in STRATEGIES:
            print(f"    {strategy.label}", flush=True)
            row, arrays = run_one(scenario, strategy, growth_factor, args.loglevel)
            rows.append(row)
            arrays_by_key[row_key(row)] = arrays

            suffix = f"{scenario.name}__growth{growth_factor:g}__{strategy.label}"
            write_timestep_history(histories_dir / f"{suffix}.csv", arrays["dts"])
            write_temperature_profile(
                profiles_dir / f"{suffix}.csv",
                arrays["grid"],
                arrays["temperature"],
            )
            print(
                f"      ok={row['success']} steps={row['accepted_steps']} "
                f"jac={row['jacobians']} wall={row['wall_time_s']:.2f}s",
                flush=True,
            )

write_csv(output_dir / "summary.csv", rows)
write_summary_md(output_dir / "summary.md", rows)
(output_dir / "metadata.json").write_text(
    json.dumps(
        {
            "cantera_version": ct.__version__,
            "cantera_module": ct.__file__,
            "scenarios": [
                {
                    "name": scenario.name,
                    "title": scenario.title,
                    "description": scenario.description,
                    "mechanism": scenario.mechanism,
                    "growth_factors": scenario.growth_factors,
                    "suites": scenario.suites,
                }
                for scenario in scenarios
            ],
            "strategies": [strategy.__dict__ for strategy in STRATEGIES],
        },
        indent=2,
    )
    + "\n",
    encoding="utf-8",
)

for scenario in scenarios:
    scenario_rows = [row for row in rows if row["scenario"] == scenario.name]
    if not scenario_rows:
        continue
    growth_factors = sorted({float(row["growth_factor"]) for row in scenario_rows})
    if len(growth_factors) > 1:
        plot_growth_sweep(
            rows,
            scenario.name,
            plots_dir / f"{scenario.name}_growth_sweep.png",
        )
    for growth_factor in growth_factors:
        plot_strategy_bars(
            rows,
            scenario.name,
            growth_factor,
            plots_dir / f"{scenario.name}_{growth_tag(growth_factor)}_bars.png",
        )
        plot_timestep_histories(
            rows,
            arrays_by_key,
            scenario.name,
            growth_factor,
            plots_dir / f"{scenario.name}_{growth_tag(growth_factor)}_timestep_histories.png",
        )
write_report(output_dir / "report.md", rows, scenarios, plots_dir)

print()
print(f"Wrote {output_dir / 'summary.csv'}")
print(f"Wrote {output_dir / 'summary.md'}")
if (output_dir / "report.md").exists():
    print(f"Wrote {output_dir / 'report.md'}")
print(f"Wrote plots in {plots_dir}")

if name == "main":
main()

Some of the plots for the cases that the script runs through.

This one is a methane case with a max growth factor of 20.
gri30_methane_burner_growth20_bars

This is the same case on a larger grid.

gri30_methane_burner_refined_growth20_bars

And this is a premixed flame case.
gri30_twin_premixed_growth20_bars

These are larger grids with larger mechanisms where some differences seem to be showing up.

@speth
Copy link
Copy Markdown
Member

speth commented May 19, 2026

I went back and ran some cases with larger mechanisms like the gri30 one, and it does seem like at least some subset of those heuristics seem to give a performance improvement.

Thanks for sharing that set of results and the script. I think what I was missing is that the "growth strategies" matter most when you substantially increase the time_step_growth_factor -- these results are all with a growth factor of 20, versus 2 in the example I was looking at (and where the default is 1.5). The combination of higher growth factor and one of these alternative strategies -- particularly steady-norm -- seems like it can reduce runtime by a meaningful amount. I think it makes sense to leave the defaults as-is for now and just have these options available. Once we get a little more comfortable with them and maybe do some more testing on different mechanisms / configurations, we can consider changing the default solver behavior.

@speth speth force-pushed the timestep_regridding branch from 7879319 to b797be0 Compare May 29, 2026 02:01
@speth
Copy link
Copy Markdown
Member

speth commented May 29, 2026

@wandadars - Thanks again for the additional exploration of the effect the different options have on convergence rate. I updated the conditions in the example based on this to show a case where a couple different time step growth strategies do better than the default, using the BurnerFlame configuration. Take a look and let me know if you have any other suggestions for that example. Otherwise, I'm ready to merge this.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants