Skip to content

User-defined functions (UDFs)

Write a closure in Python, name it from a case file, and have the solver call it per cell, per step.

This chapter is one section per seam. Read this page first — the permission model, the state object and the return contract are shared by all of them; the per-seam pages only cover what differs.

Status of each seam

The solver carries several py_* seams. They are at different stages, and this table is the honest version — a seam existing in the solver is not the same as a TOML path reaching it.

Seam Section TOML path for your own code?
Pipe friction Friction Yes — kind = "python", dotted path
Wall heat Wall heat No — seam is live, but only built-in closures produce it
Porous drag Porous Yes — kind = "python", dotted path
Junction BC Junction BC Partly — a fixed selector, not an arbitrary path
Surface reactions Reactions No — internal to the mechanism closures

Friction and porous were first because they were the two seams the solver consumed that nothing else wrote — no built-in closure and no TOML path. Adding an arm to either is purely additive: it cannot change the result of any case that does not opt in. The remaining seams are not in that position, and each section says why.

The permission model

Resolving a dotted path imports and executes the module it names. Case files get emailed, committed, attached to bug reports and pasted into issues, so a format that silently runs whatever it names is one you cannot accept from a stranger.

aeolus1d run case.toml --allow-user-code

Two properties that matter more than the flag itself:

  • Off by default. No flag, no import — you get a refusal naming the path, not a stack trace about a missing module.
  • The case file cannot grant it. The permission is process-level (AEOLUS_ALLOW_USER_CODE), never a key in the TOML. A file that could set allow_user_code = true about itself would be deciding whether to trust itself, which is not a check at all.

The module must be importable from your interpreter — PYTHONPATH, or pip-installed. aeolus1d deliberately does not add the case file's own directory to sys.path: a case's meaning should not depend on where the file happens to sit.

PipeState — what your function receives

Every per-cell UDF receives a PipeState (aeolus1d.closures.user_code.PipeState) describing the pipe's interior cells; ghosts are excluded.

field type meaning
rho, u, p array, length N density, velocity, static pressure
T array static temperature — precomputed, see below
x, A, D_h array cell centre, area, hydraulic diameter
dx float cell width
R_gas, gamma float per-pipe gas constant; cell-mean γ
t float simulation time
N int number of interior cells

All arrays are read-only. Mutating solver state from a closure is not a supported way to influence the solution, and silently succeeding at it would make runs irreproducible.

Do not recompute T

The one rule worth stating loudly.

T looks like one line — p / (rho * R) — and the trap is the p. The solver recovers cell pressure through a three-way branch: real-gas EOS, per-cell γ, then the scalar (γ−1)ρe. The scalar form is wrong for the first two. A closure that reaches for state.gamma gets the cell mean, which is exactly the two-γ seam that once hid 26 % of a cylinder's heat release — while conserving energy perfectly, so nothing downstream complained.

state.T is computed by the owner from the same pressure the flux kernel used. Use it. ADR-0013 applies at this boundary: the owner hands over primitives, the far side does not re-derive them.

The return contract

Each seam documents its own return shape, but the validation is common and deliberate: returns are checked, not trusted.

R[:, 1] += s_mom broadcasts — a scalar return would be added to every cell and produce a plausible wrong answer instead of an error. A UDF is the one place in the solver where the array on the right of that += was written by someone who has never seen the residual layout. So:

  • wrong shape → UserCodeError naming the shape it got and the one it wanted;
  • any non-finite value → UserCodeError at the source. NaN loses every comparison, so a downstream dt guard would pass it straight through (see ADR-0011);
  • an exception inside your model → wrapped in UserCodeError naming the closure, so it is obvious whose code failed.

Performance

Your closure is called once per pipe per step with whole arrays. Keep it vectorised — write np.abs(state.u), not a Python loop over state.N.

Measured at 200k cells, the same model three ways:

form time
njit kernel calling an njit function 0.58 ms
vectorised numpy 1.90 ms
per-cell Python loop 190 ms

A compiled tier that accepts an @njit function directly is the natural next step and is not built yet. When it is, it must compute T before dispatch for the reason above — the three-way branch must not be duplicated on the far side of the boundary.