—
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.reactor_models import PFRHomogeneousShellNet, RefractoryReactor
from bloc.utils import get_mechanism_path
@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()