Skip to content

Tutorial 4 — Post-processing: from raw trace to engine numbers

Time: ~30 minutes. Prerequisites: Tutorial 3 (Engine end-to-end) — you should have an HDF5 trace from a multi-cylinder engine run on disk and the basic verbs of the tracer / examples/plots/ pattern down. Anchored on: any cylinder p–V trace; the worked example uses examples/52_knock_event.toml (single-cylinder SI with Wiebe + Livengood-Wu knock controller — produces a clean firing trace with deliberate knock events) so the numbers in this page are reproducible.

After this tutorial you'll have used every layer of the post-processing toolkit:

  • aeolus1d.engine_post — heat release inversion + CA10/50/90 + polytropic n + IMEP_gross / PMEP / volumetric efficiency / COV of IMEP
  • aeolus1d.analysis.acoustics — SPL / dB(A), 1/3-octave bands, engine-order analysis, muffler transmission loss
  • aeolus1d.analysis.system_reports — turbo operating-point overlay, energy-balance Sankey numbers, cumulative emissions

The full reference for each function — formulas, kwargs, references — lives in docs/post-processing.md. This page walks you through the typical workflow against a real run.

Step 1 — Run the case + capture the trace

The post-processing layer reads from the HDF5 Tracer output, so you need a run with --trace-dt set fine enough to resolve the cycle. 52_knock_event.toml turns the crank at 1800 rpm, so one 720° four-stroke cycle spans ~66.7 ms; a 50 µs trace step then gives ~1300 samples per cycle — plenty for IMEP and heat-release analysis.

Extend t_end so several cycles resolve

The bundled case ships with a short [case] t_end = 0.02 (a knock-seed smoke value — less than one full 720° cycle). The per-cycle / run-level metrics below need at least two complete cycles, so before running, bump t_end in the TOML to cover a handful of cycles — e.g. t_end = 0.3 (≈ 4.5 cycles at 1800 rpm). There is no --t-end CLI override; the duration lives in the [case] block (or the runner dock's t_end field in the GUI).

aeolus1d run examples/52_knock_event.toml \
    --out /tmp/knock.h5 \
    --trace-dt 5e-5

Open 52_knock_event.toml, tick Trace every Δt in the runner dock and set it to 5.0e-5, raise the t_end (s) field to 0.3, then run with F5. The GUI writes the HDF5 into a per-run temp directory (the run banner prints the full Output: path).

There is no field-selection flag: --trace-dt records every pipe cell's primitives (rho, u, p, plus any species Y_<name>) and every volume's scalar state (a cylinder's p / T / m / V) at each sample. You pick the fields you want at read time by indexing the on-disk arrays — p for the cylinder p–V work below, rho / u / p at an exhaust cell for the acoustics + system-reports paths.

Step 2 — Cycle metrics + heat release

Every post-processing function takes plain NumPy arrays — the on-disk HDF5 → array unpacking is the same one you'd do for a custom plot:

import h5py
import numpy as np
from aeolus1d.engine_post import (
    cycle_metrics, engine_metrics,
    heat_release_analysis, cov_imep,
)

# A cylinder is a zero-D *volume* in the network, so its trace lives
# in the time_series/junctions group — one indexed junction_NN group
# per volume, tagged with attrs kind="cylinder" and the user-facing
# id. Scan for the one you want. (The GUI's HDF5 reader does this
# partitioning for you; here we do it by hand.)
def cylinder_group(ts, cyl_id):
    jroot = ts["junctions"]
    for name in jroot:
        g = jroot[name]
        kind = g.attrs.get("kind")
        cid = g.attrs.get("id")
        # h5py may hand back bytes for string attrs — normalise.
        kind = kind.decode() if isinstance(kind, bytes) else kind
        cid = cid.decode() if isinstance(cid, bytes) else cid
        if kind == "cylinder" and cid == cyl_id:
            return g
    raise KeyError(f"no cylinder {cyl_id!r} in trace")

with h5py.File("/tmp/knock.h5") as h5:
    ts = h5["time_series"]
    t = ts["t"][:]
    # Cylinder p(t) and V(t) — scalar-per-sample trajectories on the
    # cylinder volume. V is the slider-crank V(θ) the cylinder tracks.
    cyl = cylinder_group(ts, "cyl1")
    p_cyl = cyl["p"][:]
    V_cyl = cyl["V"][:]

# omega_crank is a fixed engine parameter, not a traced field — read
# it from the case ([[cylinder]] omega_crank in 52_knock_event.toml =
# 188.49555921538757 rad/s = 1800 rpm).
omega_crank = 188.49555921538757

# Theta in the 720° four-stroke convention.
theta = (t * omega_crank * 180.0 / np.pi) % 720.0

# 52_knock_event geometry: 80 mm bore × 90 mm stroke → V_disp ≈ 452 cc.
V_disp = np.pi / 4.0 * 0.080**2 * 0.090

# Run-level rollup with the Phase 7g.2 fields populated.
em = engine_metrics(
    t, p_cyl, V_cyl, theta_deg=theta,
    V_disp=V_disp,
    omega_crank_rad_s=omega_crank,
    ivc_deg=230.0, evo_deg=490.0,            # intake close / exhaust open → IMEP_gross + PMEP
    m_fuel_per_cycle=5.0e-5, LHV=44.0e6,    # 50 mg/cycle → ISFC + η_th
    m_air_trapped_kg=V_disp * 1.6,          # ≈ ρ_intake(1.4 bar boost) × V_disp
    rho_intake_ref_kg_m3=1.6,
)

print(f"IMEP_net    : {em.imep_mean_pa/1e5:.2f} bar")
print(f"IMEP_gross  : {em.imep_gross_mean_pa/1e5:.2f} bar")
print(f"PMEP        : {em.pmep_mean_pa/1e5:.3f} bar")
print(f"η_vol       : {em.eta_volumetric:.3f}")
print(f"η_thermal   : {em.eta_thermal:.3f}")
print(f"ISFC        : {em.isfc_g_kwh:.0f} g/kWh")
print(f"COV(IMEP)   : {em.cov_imep*100:.2f} %")
print(f"Cycles      : {em.n_cycles_resolved}")

For a cycle-by-cycle look:

cycles = cycle_metrics(
    t, p_cyl, V_cyl, theta_deg=theta,
    V_disp=V_disp, omega_crank_rad_s=omega_crank,
    ivc_deg=230.0, evo_deg=490.0,
)
for k, c in enumerate(cycles):
    print(f"cycle {k:>2d} : IMEP={c.imep_pa/1e5:.2f} bar  "
          f"p_max={c.p_max_pa/1e5:.1f} bar @ {c.p_max_theta_deg:.1f}°")

You'll see one line per resolved cycle; on a knock-laden run the later cycles' p_max walks down as the closed-loop spark-retard controller pulls timing.

For the inverse heat-release problem — given p(θ) and V(θ), recover dQ/dθ + the MFB curve + the canonical CA10/50/90 combustion-phasing metrics:

# Pick the cycle you want to analyse — typically the last
# converged one (steady cycle-to-cycle behaviour).
last = cycles[-1]
i0 = next(i for i, ti in enumerate(t)
          if (ti * omega_crank * 180/np.pi) % 720.0 >= last.p_max_theta_deg - 50)
# … grab one cycle's slice of (theta_c, p_c, V_c) here …

hr = heat_release_analysis(
    p_c, V_c, theta_c,
    gamma=1.30,                  # mixed-burnt-unburnt at peak T
    ivc_deg=230.0, evo_deg=490.0,
    smooth_window=3,
)

print(f"CA10 / CA50 / CA90 : "
      f"{hr.ca10_deg:.1f}° / {hr.ca50_deg:.1f}° / {hr.ca90_deg:.1f}°")
print(f"Burn duration      : {hr.burn_duration_deg:.1f}°")
print(f"Q_total            : {hr.Q_total_j:.0f} J")
print(f"n_compression      : {hr.n_compression:.3f}")
print(f"n_expansion        : {hr.n_expansion:.3f}")

The polytropic exponents are the diagnostic dial here:

  • close to γ on both sides → almost no heat loss / blowby (the closed cycle is nearly isentropic, as a calibrated SI engine should be at part load)
  • compression n significantly less than γ → mid-stroke heat loss to the cylinder walls
  • expansion n less than compression n → late-cycle heat loss (or a slow burn that hasn't completed)

The bundled examples/plots/plot_post_engine_diagnostics.py runs all of the above plus renders a four-panel figure (p–V indicator, p(θ) with CA10/50/90 overlay, dQ/dθ + cumulative Q, MFB curve) — read it as a worked example.

Step 3 — Acoustics from a probe trace

"Probe" the exhaust at the tail-pipe outlet (or, for this single-cylinder demo, at the exhaust-pipe outlet) and you can run the same NVH post-processing chain GT-POST and AVL Cruise M ship for muffler tuning. In raw HDF5 a probe is just a cell of a pipe — there is no dedicated probe group; pipe fields are stored as (n_samples, n_cells) arrays, so index the cell you want. (In the GUI you drop a probe marker on the network and it maps to a pipe fractional position for you.)

from aeolus1d.analysis.acoustics import (
    spl_db, spl_db_a, power_spectrum,
    octave_band_levels, order_analysis,
)

with h5py.File("/tmp/knock.h5") as h5:
    # Tail-pipe pressure = the outlet cell of the exhaust pipe.
    p_tail = h5["time_series/pipes/exhaust/p"][:, -1]

fs = 1.0 / (t[1] - t[0])
# Fixed-speed run → the crank frequency is constant. (For an RPM
# sweep you would track ω(t) yourself; it is not a traced field.)
omega = np.full(t.shape, omega_crank)

print(f"SPL    : {spl_db(p_tail):.1f} dB re 20 µPa")
print(f"SPL(A) : {spl_db_a(p_tail, fs):.1f} dB(A)")

# Spectrum + 1/3-octave band reduction.
freqs, amp = power_spectrum(p_tail, fs, window="hann")
for b in octave_band_levels(freqs, amp, n_per_octave=3):
    print(f"  {b.centre_hz:>7.1f} Hz   {b.spl_db:5.1f} dB")

The dB(A) figure is the regulatory noise-pollution metric (IEC 61672 A-weighting attenuates the low-frequency content the human ear is insensitive to; the firing fundamental at 50 Hz on a 4-stroke 4-cylinder at 1500 rpm gets ~19 dB knocked off). dB linear is the engineering number to compare against an unweighted SPL meter.

For an RPM-sweep run (e.g. an open-loop ramp from 1500 to 4500 rpm), the canonical NVH plot is the order map:

track = order_analysis(
    p_tail, sample_rate_hz=fs,
    omega_crank_t_s=t, omega_crank_rad_s=omega,
    orders=np.array([0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0]),
    window_s=0.05, overlap=0.5,
)
# track.amplitude_pa[i, j] is the amplitude of orders[j] at
# rpm[i] — pcolormesh(track.times_s, track.orders, track.amplitude_pa.T)
# gives the GT-POST-style frequency × order × time map.

Engine "orders" are integer / half-integer multiples of crank rotation frequency. Order 2 is the dominant gas-pressure harmonic of a 4-stroke 4-cylinder; order 0.5 is the firing frequency of a single-cylinder 4-stroke; order 1 captures once-per-rev imbalances.

For muffler design — given probe pressures at the upstream and downstream sides of an acoustic element:

from aeolus1d.analysis.acoustics import transmission_loss

with h5py.File("/tmp/exhaust_with_muffler.h5") as h5:
    # Two probe cells straddling the acoustic element: the outlet cell
    # of the upstream pipe and the inlet cell of the downstream pipe.
    p_in = h5["time_series/pipes/pipe_upstream/p"][:, -1]
    p_out = h5["time_series/pipes/pipe_downstream/p"][:, 0]

freqs, tl_db = transmission_loss(p_in, p_out, fs)
# tl_db[k] > 0 → element attenuates at freqs[k]
# tl_db[k] < 0 → resonance amplifies

insertion_loss(p_baseline, p_with_element, fs) is the more end-user-relevant metric: the noise reduction the element actually delivers in situ, including the impedance-match interaction with the rest of the system.

The bundled examples/plots/plot_post_acoustics.py exercises the SPL / dB(A), 1/3-octave, order-analysis, and transmission-loss functions on synthetic traces — useful as a self-contained demo that runs in milliseconds without an aeolus1d simulation.

Step 4 — System-level reports

For end-of-run rollups:

from aeolus1d.analysis.system_reports import (
    energy_balance, cumulative_emissions, is_steady_state,
)

# Per-cycle first-law partition, ready to render as a Sankey
# diagram. Plug in the per-cycle measurements from your run:
eb = energy_balance(
    m_fuel_per_cycle_kg=5.0e-5,
    LHV_J_per_kg=44.0e6,
    indicated_work_j=em.work_mean_j,
    pumping_work_j=em.pmep_mean_pa * V_disp,
    exhaust_enthalpy_j=750.0,            # ∫ h_exh · ṁ from a tail-pipe cell
    wall_heat_j=320.0,                   # ∫ h_w · A_w · (T_w − T_g) from cylinder wall
)
print(eb.fractions)
# {'indicated_work': 0.34, 'pumping_work': 0.05, 'exhaust_enthalpy': 0.39,
#  'wall_heat': 0.15, 'residual': 0.07}

The residual collects everything that wasn't explicitly accounted for (blow-by, radiation, mechanical friction not modelled, integration error). A small residual (< 10 % of fuel-in) means the energy balance closes; a large one points at a missing modelling term.

For total emissions over a test cycle:

# Cumulative emissions integrate ṁ · Y past a tail-pipe cell, so this
# needs a *reacting* exhaust case (52_knock_event carries no species).
# ṁ is not a traced field — reconstruct it from ρ·u·A at the cell.
with h5py.File("/tmp/exhaust_reacting.h5") as h5:
    t_exh = h5["time_series/t"][:]
    pipe = h5["time_series/pipes/tailpipe"]
    A_cell = h5["pipes/tailpipe/A_cell"][-1]      # static geometry, outlet cell
    mdot_tail = pipe["rho"][:, -1] * pipe["u"][:, -1] * A_cell
    Y_co = pipe["Y_CO"][:, -1]
    Y_nox = pipe["Y_NO"][:, -1]

totals = cumulative_emissions(
    times_s=t_exh, mdot_kg_s=mdot_tail,
    Y_species={"CO": Y_co, "NOx": Y_nox},
)
print(f"Total CO emitted : {totals['CO']*1000:.3f} g")
print(f"Total NOx emitted: {totals['NOx']*1000:.3f} g")

For convergence detection — useful when sweeping through operating points and you want each one to settle before recording its steady-state numbers:

imep_history = [c.imep_pa for c in cycles]
report = is_steady_state(imep_history, cov_threshold=0.02, window=5)
if report.converged:
    print(f"Converged after {len(imep_history)} cycles "
          f"(COV(IMEP) = {report.cov_imep_window*100:.2f} %)")

Step 5 — Turbo operating-point overlay

For turbocharged cases, the standard sizing chart is the operating-point trajectory overlaid on the compressor map's surge / choke envelope. You record the compressor's per-step state during the run via an on_step recorder, then plot:

from dataclasses import dataclass
from aeolus1d.analysis.system_reports import (
    turbo_map_envelope, turbo_operating_points,
)

# Build a per-step recorder on advance_network(on_step=...)
@dataclass
class CompSnap:
    last_pr: float
    last_mdot_corr: float   # the CORRECTED flow the map is drawn on
    last_n_corr: float
    last_eta: float
    last_regime: str

times: list[float] = []
history: list[CompSnap] = []
def record(t):
    times.append(t)
    history.append(CompSnap(
        comp.last_pr, comp.last_mdot_corr,
        comp.last_n_corr, comp.last_eta, comp.last_regime,
    ))

# After the run:
op = turbo_operating_points(np.asarray(times), history)
env = turbo_map_envelope(comp.map)

# Plot:
#   env.{mdot,pr}_{surge,choke}_per_line — static surge / choke loci
#   (op.mdot_corr, op.pr) — trajectory through them
#   colour-code each scatter point by op.regime[i] to highlight
#   surge / choke / normal transitions

The Phase 7c-surge last_regime flag carries through every trajectory point so you can visually distinguish runs that hit surge from those that stayed in the normal envelope. The bundled examples/plots/plot_post_turbo_overlay.py is a worked example with synthetic data.

Step 6 — Drive-cycle integration

The fuel-economy / cumulative-emissions number you'd report for a regulatory test cycle (NEDC / WLTP / FTP) is ∫ ṁ_fuel · dt across the full schedule. Aeolus's drive_cycle_runner does the integration; you supply a function that maps (rpm, torque) to instantaneous rates:

from aeolus1d.analysis.system_reports import (
    DriveCyclePoint, drive_cycle_runner,
)

# 1. Build a steady-state engine map ONCE — one aeolus1d run per
#    operating point (or use a measured map / regression).
#    Returns at minimum {fuel_rate_kg_s, power_w, ...}.
steady_state_map = build_my_engine_op_map()    # your code

def engine_op(rpm, torque_nm):
    return steady_state_map.lookup(rpm, torque_nm)

# 2. Author the drive cycle (NEDC / WLTP / FTP / your own).
schedule = [
    DriveCyclePoint(t_s=0,     rpm=800,  torque_nm=0),
    DriveCyclePoint(t_s=10,    rpm=1500, torque_nm=60),
    # ... full cycle ...
]

# 3. Integrate.
res = drive_cycle_runner(schedule, operating_point_callable=engine_op)
print(f"Total fuel : {res.total_fuel_kg*1000:.1f} g")
print(f"CO / NOx   : {res.emissions_total_kg}")

The intent: build the steady-state engine map once and drive any cycle through it without re-running the 1-D solver every step. For a typical NEDC test cycle (~20 minutes of simulated time at 1-second resolution = ~1200 segments × ~50 ms per evaluation through a regression-fit map), this runs in seconds rather than the 20-minute real-time the actual aeolus1d run would need.

What you've done

  • Loaded a tracer HDF5 and pulled the cylinder p(t) / V(t) + probe pressure series into NumPy arrays.
  • Reduced a multi-cycle p-V trace to per-cycle and run-level engine metrics — IMEP, ISFC, η_thermal, η_volumetric, COV(IMEP), PMEP, IMEP_gross.
  • Inverted the first-law balance to recover the mass-fraction- burned curve + CA10/50/90 + polytropic compression / expansion exponents from the simulated p-V trace.
  • Computed dB / dB(A) sound pressure levels, 1/3-octave band reductions, and engine-order maps from a probe pressure trace.
  • Sized a muffler via transmission-loss / insertion-loss curves.
  • Built a per-cycle Sankey energy balance, integrated cumulative emissions past a tail-pipe probe, and run a regulatory drive- cycle through a steady-state engine map.
  • Plotted a turbo trajectory on the compressor map's surge / choke envelope.

That covers every layer of the post-processing toolkit. The reference page has the complete formula list + per-function kwargs + literature references; this tutorial is the workflow story.

Where to go from here

  • Custom probes. With --trace-dt the Tracer records every pipe cell's primitives (rho / u / p + any species Y_<name>) and every volume's scalar state — there is no field-selection flag. You "probe" a location by indexing the (n_samples, n_cells) pipe array at the cell you want; in the GUI you drop a probe marker on the network and it resolves to a pipe fractional position for you. The post-processing chain takes plain NumPy arrays, so every function is composable with arbitrary user data.
  • Optimise against measured data. Pair the cycle-resolved metrics with aeolus1d optimize (see the optimization example) to fit Wiebe parameters / SCR rate constants / DOC washcoat density to a measured trace.
  • Run the post-processing inside the GUI. Drop a probe in the GUI's network view, run, then open the 2-D plot dock — every aeolus1d.engine_post and aeolus1d.analysis.acoustics function works the same way against the live HDF5 the runner loads. The GUI's plot dock can also overlay multiple runs for case-sweep comparison.