—
Benchmark three PFR approaches on the same adiabatic case.
Methods compared: 1) Lagrangian particle march with one 0D constant-pressure reactor 2) Steady-state approximation with a chain of CSTRs 3) Bloc RefractoryReactor solved through PFRHomogeneousShellNet
- Run from repository root:
conda run -n bloc python examples/benchmark_pfr_three_approaches.py
Inspired from https://cantera.org/dev/examples/python/reactors/pfr.html “Plug flow reactor modeling approaches”
from __future__ import annotations
import math
import time
from dataclasses import dataclass
from pathlib import Path
import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
from bloc.paths import get_mechanism_path
from bloc.reactors import PFRHomogeneousShellNet, RefractoryReactor
@dataclass
class MethodResult:
"""Container for one method output and runtime."""
name: str
elapsed_s: float
x: np.ndarray
t: np.ndarray
T: np.ndarray
Y: dict[str, np.ndarray]
def _species_indices(gas: ct.Solution, names: list[str]) -> dict[str, int]:
"""Return available species indices for requested names."""
out: dict[str, int] = {}
for name in names:
if name in gas.species_names:
out[name] = gas.species_index(name)
return out
def run_lagrangian(
mechanism: str,
T_in_K: float,
P_in_Pa: float,
X_in: str,
length_m: float,
diameter_m: float,
mdot_kg_s: float,
species_to_track: list[str],
) -> MethodResult:
"""Run Method 1: one 0D reactor, time marched and mapped to x."""
t0 = time.perf_counter()
gas = ct.Solution(mechanism)
gas.TPX = T_in_K, P_in_Pa, X_in
reactor = ct.IdealGasConstPressureMoleReactor(gas, clone=False)
sim = ct.ReactorNet([reactor])
sim.preconditioner = ct.AdaptivePreconditioner()
area_m2 = math.pi * (diameter_m / 2) ** 2
states = ct.SolutionArray(gas, extra=["t", "x"])
states.append(gas.state, t=sim.time, x=0.0) # type: ignore[call-arg]
x = 0.0
dx_long = length_m / 60 # fewer Lagrangian cells for faster run
dt = 1e-9
while x < length_m:
uz = mdot_kg_s / max(gas.density * area_m2, 1e-20)
dx = min(uz * dt * 1.1, dx_long, length_m - x)
dt = dx / max(uz, 1e-20)
sim.advance(sim.time + dt)
x += dx
states.append(gas.state, t=sim.time, x=x) # type: ignore[call-arg]
idx = _species_indices(gas, species_to_track)
y_profiles = {sp: np.array(states.Y[:, i], dtype=float) for sp, i in idx.items()}
return MethodResult(
name="Lagrangian 0D",
elapsed_s=time.perf_counter() - t0,
x=np.array(states.x, dtype=float),
t=np.array(states.t, dtype=float),
T=np.array(states.T, dtype=float),
Y=y_profiles,
)
def run_cstr_chain_steady(
mechanism: str,
T_in_K: float,
P_in_Pa: float,
X_in: str,
length_m: float,
diameter_m: float,
mdot_kg_s: float,
n_cells: int,
species_to_track: list[str],
) -> MethodResult:
"""Run Method 2: closed-cell Euler march — N equal spatial cells, each advanced for τ.
Each cell is a **closed** constant-pressure reactor (no MFC/PressureController).
The cell receives the outlet state of the previous cell as its initial condition,
then advances for the local residence time τ = ρ·V_cell / ṁ.
This discretises the PFR ODE identically to the Lagrangian method, but on a
uniform spatial grid. It converges to the Lagrangian for any N and avoids the
Damköhler-number limitation of the open-CSTR/solve_steady() approach (which
requires Da_cell << 1, i.e. chemistry slower than τ_cell, to avoid premature
equilibration per cell).
"""
t0 = time.perf_counter()
area_m2 = math.pi * (diameter_m / 2) ** 2
dz = length_m / n_cells
cell_volume = area_m2 * dz
gas = ct.Solution(mechanism)
gas.TPX = T_in_K, P_in_Pa, X_in
idx = _species_indices(gas, species_to_track)
x = np.zeros(n_cells + 1)
t = np.zeros(n_cells + 1)
T = np.zeros(n_cells + 1)
y_profiles = {sp: np.zeros(n_cells + 1) for sp in idx}
x[0] = 0.0
t[0] = 0.0
T[0] = float(gas.T)
for sp, j in idx.items():
y_profiles[sp][0] = float(gas.Y[j])
tau = 0.0
for i in range(1, n_cells + 1):
# Closed constant-pressure reactor: one fresh ReactorNet per cell so
# there is no state leaking between cells.
reactor = ct.IdealGasConstPressureReactor(gas, clone=True)
sim = ct.ReactorNet([reactor])
# Residence time: mass of this fluid element / mass-flow rate.
tau_cell = gas.density * cell_volume / max(mdot_kg_s, 1e-30)
sim.advance(tau_cell)
# Propagate outlet state into `gas` for the next cell.
gas.TDY = reactor.phase.TDY
x[i] = i * dz
tau += tau_cell
t[i] = tau
T[i] = float(gas.T)
for sp, j in idx.items():
y_profiles[sp][i] = float(gas.Y[j])
return MethodResult(
name=f"CSTR chain (N={n_cells})",
elapsed_s=time.perf_counter() - t0,
x=x,
t=t,
T=T,
Y=y_profiles,
)
def run_design_pfr_net(
mechanism: str,
T_in_K: float,
P_in_Pa: float,
X_in: str,
length_m: float,
diameter_m: float,
mdot_kg_s: float,
species_to_track: list[str],
) -> MethodResult:
"""Run Method 3: Bloc RefractoryReactor via PFRHomogeneousShellNet."""
t0 = time.perf_counter()
gas = ct.Solution(mechanism)
gas.TPX = T_in_K, P_in_Pa, X_in
reactor = RefractoryReactor(gas, clone=False)
reactor._meta = {
"design_type": "RefractoryReactor",
"mechanism": mechanism,
"length": length_m,
"diameter": diameter_m,
"mass_flow_rate": mdot_kg_s,
"adiabatic": True,
}
net = PFRHomogeneousShellNet([reactor], meta=dict(reactor._meta))
net.advance(time=1.0)
states = reactor._states
idx = _species_indices(gas, species_to_track)
y_profiles = {sp: np.array(states.Y[:, i], dtype=float) for sp, i in idx.items()}
return MethodResult(
name="PFRHomogeneousShellNet",
elapsed_s=time.perf_counter() - t0,
x=np.array(states.x, dtype=float),
t=np.array(states.t, dtype=float),
T=np.array(states.T, dtype=float),
Y=y_profiles,
)
def _compare_against_reference(
ref: MethodResult, other: MethodResult
) -> dict[str, float]:
"""Compute profile error metrics with interpolation on reference x-grid."""
out: dict[str, float] = {}
T_other = np.interp(ref.x, other.x, other.T)
dT = T_other - ref.T
out["T_max_abs_K"] = float(np.max(np.abs(dT)))
out["T_rmse_K"] = float(np.sqrt(np.mean(dT**2)))
common_species = sorted(set(ref.Y).intersection(other.Y))
for sp in common_species:
y_other = np.interp(ref.x, other.x, other.Y[sp])
dy = y_other - ref.Y[sp]
out[f"Y_{sp}_max_abs"] = float(np.max(np.abs(dy)))
out[f"Y_{sp}_rmse"] = float(np.sqrt(np.mean(dy**2)))
return out
def _plot_results(
results: list[MethodResult], output_dir: Path, species: list[str]
) -> Path:
"""Create overlay plots for temperature and species profiles."""
nrows = 1 + len(species)
fig, axes = plt.subplots(nrows, 1, figsize=(10, 3 * nrows), sharex=True)
if nrows == 1:
axes = [axes]
for res in results:
axes[0].plot(res.x, res.T, label=res.name, linewidth=2)
axes[0].set_ylabel("Temperature [K]")
axes[0].grid(True, alpha=0.3)
axes[0].legend()
for row, sp in enumerate(species, start=1):
for res in results:
if sp in res.Y:
axes[row].plot(res.x, res.Y[sp], label=res.name, linewidth=2)
axes[row].set_ylabel(f"Y({sp}) [-]")
axes[row].grid(True, alpha=0.3)
axes[-1].set_xlabel("Axial position x [m]")
fig.suptitle(
"PFR benchmark: Lagrangian vs CSTR chain vs PFRHomogeneousShellNet (adiabatic)"
)
fig.tight_layout()
output_dir.mkdir(parents=True, exist_ok=True)
fig_path = output_dir / "benchmark_pfr_three_approaches.png"
fig.savefig(fig_path, dpi=180)
plt.close(fig)
return fig_path
def main() -> None:
"""Run benchmark and print performance + agreement summary."""
mechanism = get_mechanism_path("Fincke_GRC.yaml")
T_in_K = 1800.0
P_in_Pa = 1e5
X_in = "CH4:0.5,H2:0.5"
length_m = 1.0
diameter_m = 0.1
mdot_kg_s = 1e-3
n_cells_cstr = 1000
species_to_track = ["CH4", "H2", "C(s)"]
lag = run_lagrangian(
mechanism,
T_in_K,
P_in_Pa,
X_in,
length_m,
diameter_m,
mdot_kg_s,
species_to_track,
)
cstr = run_cstr_chain_steady(
mechanism,
T_in_K,
P_in_Pa,
X_in,
length_m,
diameter_m,
mdot_kg_s,
n_cells_cstr,
species_to_track,
)
design = run_design_pfr_net(
mechanism,
T_in_K,
P_in_Pa,
X_in,
length_m,
diameter_m,
mdot_kg_s,
species_to_track,
)
methods = [lag, cstr, design]
ref = lag
cmp_cstr = _compare_against_reference(ref, cstr)
cmp_design = _compare_against_reference(ref, design)
print("=== Runtime ===")
for m in methods:
print(f"{m.name:24s}: {m.elapsed_s:8.4f} s")
print("\n=== Outlet temperature [K] ===")
for m in methods:
print(f"{m.name:24s}: {m.T[-1]:10.6f}")
print("\n=== Agreement vs Lagrangian (profile errors) ===")
print("CSTR chain:")
for k, v in cmp_cstr.items():
print(f" {k:22s}: {v:.6e}")
print("PFRHomogeneousShellNet:")
for k, v in cmp_design.items():
print(f" {k:22s}: {v:.6e}")
print(
"\nNote: CSTR-chain accuracy is controlled by N. "
"Higher N reduces numerical diffusion at higher runtime."
)
fig_path = _plot_results(
methods,
output_dir=Path("examples") / "outputs",
species=[sp for sp in species_to_track if sp in lag.Y],
)
print(f"\nSaved plot: {fig_path}")
if __name__ == "__main__":
main()