Building a model in pure Python¶

Every other tutorial starts from a TOML case. This one never touches one.

You need this page when aeolus1d is a library rather than an application — inside an optimiser that rebuilds the model each iteration, in a test harness, in a parameter sweep, in a notebook. There are two Python surfaces, and the difference between them is the first thing worth being clear about:

Surface	You write	Use it when
Declarative	`Case`, `PipeSpec`, `BCSpec` … then `run_case(case)`	You want the TOML model, built in Python. Every schema default, validation rule and closure kind still applies.
Component	`Pipe`, `Volume`, `Port`, `Junction`, BC callables … then `advance_network(...)`	You want the objects themselves — to inspect state mid-run, swap a boundary condition for your own callable, or drive the timeloop yourself.

They are not alternatives layered over different engines. run_case is the declarative surface calling the component one, and at the end of this page we build the same model both ways and check that they agree — they do, to 0.0 K, bit for bit.

Every bundled example already has a declarative builder

You do not have to transcribe a TOML case by hand. All 79 bundled examples ship a generated Python equivalent under examples/python/, guarded so it cannot drift from the TOML it mirrors. See Examples in Python. This page teaches the component surface, which is the one that cannot be generated.

Step 1 — one pipe¶

A Pipe is a finite-volume domain: L metres, N cells, Ng ghost cells either side, cross-sectional area A. set_uniform fills it with a constant state.

from aeolus1d._types import End
from aeolus1d.bc.base import transmissive_bc
from aeolus1d.components.pipe import Pipe
from aeolus1d.solver.network_timeloop import advance_network

GAMMA, R_GAS = 1.4, 287.05

tube = Pipe(L=1.0, N=200, gamma=GAMMA, Ng=2, A=1.0)
tube.set_uniform(rho=1.0, u=0.0, p=1.0e5)

Ng=2 because the default MUSCL reconstruction needs a two-cell stencil outside the domain; first-order runs are happy with Ng=1.

The conserved state is a plain (N, 3) array of (ρ, ρu, ρE), so a Riemann initial condition is an array write:

half = tube.N // 2
tube.U[:half, 0] = 1.0                     # ρ left
tube.U[half:, 0] = 0.125                   # ρ right
tube.U[:half, 2] = 1.0e5 / (GAMMA - 1.0)   # ρE left  (u = 0, so ρE = p/(γ-1))
tube.U[half:, 2] = 1.0e4 / (GAMMA - 1.0)   # ρE right

Read it back in primitive variables with primitives(), which returns (ρ, u, p):

rho, u, p = tube.primitives()

Boundary conditions are (pipe, end, callable) triples. The callable fills that end's ghost cells; transmissive_bc and wall_bc are ready-made, and anything with the same signature works — including your own function.

far = [(tube, End.LEFT, transmissive_bc), (tube, End.RIGHT, transmissive_bc)]
t = advance_network([tube], [], far, t_end=1.0e-3, cfl=0.8)

The three positional arguments are always the same: pipes, junctions, far-field BCs. A single pipe has no junctions, hence the empty list.

Running it:

after set_uniform: U.shape=(200, 3) dx=0.005 A_cell[0]=1
initial p: left=100000  right=10000
advanced to t=0.001 s
p range [29888.3, 99942] Pa
u max   296.168 m/s

Step 2 — two pipes and a plenum¶

A Junction is a Volume (the plenum gas) plus a list of Ports (where pipes meet it). Each port names the pipe, which end of it touches the plenum, and the port area.

import math

from aeolus1d.bc.base import wall_bc
from aeolus1d.components.junction import Junction, Port, Volume

D = 0.05
AREA = math.pi * 0.25 * D * D

a = Pipe(L=0.5, N=100, gamma=GAMMA, Ng=2, A=AREA)
a.set_uniform(rho=2.0e5 / (R_GAS * 400.0), u=0.0, p=2.0e5)
b = Pipe(L=0.5, N=100, gamma=GAMMA, Ng=2, A=AREA)
b.set_uniform(rho=1.0e5 / (R_GAS * 300.0), u=0.0, p=1.0e5)

plenum = Volume(V=5.0e-4, gamma=GAMMA, R_gas=R_GAS)
plenum.set_state(p=1.5e5, T=350.0)

junction = Junction(volume=plenum, ports=[
    Port(pipe=a, end=End.RIGHT, A_port=AREA),
    Port(pipe=b, end=End.LEFT,  A_port=AREA),
])
junction.__post_init__()

Call __post_init__() after constructing a Junction by hand

Junction is a dataclass whose __post_init__ wires the port bookkeeping (and builds neck pipes when neck_cells > 0). Constructing it with keyword arguments runs it automatically; the explicit call above is harmless and is written out because you must call it if you mutate ports afterwards.

The far-field list now covers only the two outer ends — the inner ends belong to the junction, and giving them a BC as well would be applying two boundary treatments to one face:

far = [(a, End.LEFT, wall_bc), (b, End.RIGHT, wall_bc)]
t = advance_network([a, b], [junction], far, t_end=2.0e-3, cfl=0.5)

plenum before: p=150000 Pa  m=0.00074651 kg
advanced to t=0.002 s
plenum after : p=147313 Pa  m=0.000715133 kg
total mass   = 0.00359661662135 kg

Both far ends are walls, so this network is closed — the total inventory is an invariant you can assert on, which is exactly what tests/integration/test_conservation_audit.py does:

import numpy as np

def inventory(pipes, junctions):
    m = sum(float(np.sum(p.U[:, 0] * p.A_cell) * p.dx) for p in pipes)
    m += sum(float(j.volume.m) for j in junctions)
    return m

Step 3 — driving it: real boundary conditions¶

Ready-made BCs come from factories named make_*_bc, each returning a callable you drop into the same (pipe, end, callable) triple. They live in aeolus1d.bc.* — mass-flow and velocity inlets, stagnation inlets, pressure outlets, non-reflecting, throttle, EGR valve, poppet valve, orifice, and transient (schedule-driven) versions of most.

from aeolus1d.bc.mass_flow_inlet import make_mass_flow_inlet_bc
from aeolus1d.bc.pressure_outlet import make_pressure_outlet_bc

duct = Pipe(L=0.2, N=24, gamma=GAMMA, Ng=2, A=1.0e-3)
duct.set_uniform(rho=1.0e5 / (R_GAS * 800.0), u=20.0, p=1.0e5)

inlet  = make_mass_flow_inlet_bc(mdot=8.0e-3, T0=800.0, R_gas=R_GAS, gamma=GAMMA)
outlet = make_pressure_outlet_bc(p_back=1.0e5, gamma=GAMMA)

far = [(duct, End.LEFT, inlet), (duct, End.RIGHT, outlet)]

advance_network takes an on_step hook that fires once per completed step — the seam for sampling, logging or a progress bar:

samples: list[float] = []
t = advance_network([duct], [], far, t_end=0.05, cfl=0.4,
                    on_step=lambda tt: samples.append(tt))

advanced to t=0.05 s in 8781 steps
outlet T = 799.832 K     mean u = 18.3673 m/s
mdot = rho*u*A = 0.008 kg/s

The recovered mass flow is the imposed 8 g/s, which is the cheapest sanity check available on a driven network: if the steady flux does not match what you prescribed, the boundary condition is not doing what you think.

Step 4 — the same model, declaratively¶

The identical network through the schema surface:

from aeolus1d.io.case import run_case
from aeolus1d.io.schema import (
    BCSpec, Case, CaseHeader, GasSpec, PipeSpec, UniformInit,
)

case = Case(
    case=CaseHeader(name="driven_duct", t_end=0.05, cfl=0.4),
    gas=GasSpec(gamma=GAMMA, R_gas=R_GAS),
    pipes=[PipeSpec(id="duct", length=0.2, n_cells=24, area=1.0e-3,
                    init=UniformInit(p=1.0e5, T=800.0, u=20.0))],
    junctions=[],
    bcs=[BCSpec(pipe="duct", end="left", kind="mass_flow_inlet",
                mdot=8.0e-3, T0=800.0),
         BCSpec(pipe="duct", end="right", kind="pressure_outlet",
                p_back=1.0e5)],
)
pipes, junctions, t_final = run_case(case)

advanced to t=0.05 s
outlet T = 799.832 K     mean u = 18.3673 m/s
max |dT| vs the hand-assembled network = 0.000000e+00 K

Bit-identical. That is the point of the page: run_case builds the same Pipe and calls the same advance_network, so choosing a surface is a question of what you need to reach, not of what the solver will do.

Which to pick:

Declarative gives you validation, defaults, closure dispatch by name, and every feature the TOML format covers — cylinders, turbochargers, controllers, aftertreatment. It is the right default, and examples/python/ has 79 worked instances of it.
Component gives you the objects. Take it when you need to hold a Pipe and read U mid-run, or supply a boundary condition that is a closure over your own state, or step the network yourself.

Nothing stops you mixing them: build with build_network(case) to get the objects the declarative way, then drive them by hand.

from aeolus1d.io.build import build_network

pipes, junctions, far = build_network(case)   # dict[str, Pipe], list, list
advance_network(list(pipes.values()), junctions, far, t_end=0.05, cfl=0.4)

The rest of the ladder¶

This page covered pipes, junctions, boundary conditions and the two surfaces. One page per family follows, each building on it:

Family	Page	What it adds
Closures	Closures on a pipe	friction, wall heat, porous drag, and writing your own
Species	Species and chemistry	primary + extra composition columns, and the layout invariant
Aftertreatment	Aftertreatment	monolith bricks, and when to change driver
Engines	Engines	cylinder, crank, valve ports, and the clock
Turbo + control	Turbo and control	compressors, turbines, shafts, controllers

Where to go next¶

Next	Page
A Python builder for every bundled example	Examples in Python
The functions and objects above, in reference form	API reference
The same models as TOML cases	Examples
What each component models, and its theory	Component index