Skip to content

UDF: pipe friction

Status: available. kind = "python" on any [pipe.friction] block.

Read the chapter overview first — the permission model, the PipeState object and the return-validation rules are shared.

# case.toml
[pipe.friction]
kind = "python"
path = "my_models:rough_wall"      # "importable.module:callable"

# my_models.py — must be importable by the interpreter running aeolus1d
import numpy as np

def rough_wall(state):
    f = 0.02
    s_mom = -(f / (2.0 * state.D_h)) * state.rho * np.abs(state.u) * state.u
    return s_mom, np.zeros_like(s_mom)

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

Return: (s_mom, s_E)

Two arrays of length N — the momentum and total-energy source terms per cell.

s_E is normally zeros. In the total-energy formulation the kinetic energy removed by s_mom becomes internal energy automatically by conservation. An s_E = u * s_mom sink double-debits it and silently sheds the frictional heat, driving the gas below the Fanno prediction — this was a real bug in the built-in closure once. Deliberate, non-adiabatic wall heat loss belongs in the heat closure, not here.

Sign convention: friction opposes the flow, so s_mom is negative for positive u. The built-in darcy_const is -(f / (2·D_h))·ρ·|u|·u, which is the model in the snippet above.

state.D_h is geometric — it is not [pipe.friction] D_h

Corrected 2026-08-17. The sentence above is true of the formula and misleading about the inputs: the two D_h are unrelated numbers that share a name.

  • [pipe.friction] D_h (FrictionSpec.D_h) is whatever you set. The built-in kernel reads this one.
  • state.D_h, handed to your callable, is derived from the pipe's cross-sectional area as sqrt(4A/π) — the diameter of a round pipe of that area. Nothing consults the spec field.

So the snippet above reproduces darcy_const only when the pipe's area corresponds to the D_h you configured. Measured on a 0.2 m duct at f = 0.02: with area = 1e-3 and D_h = 0.02 the built-in gives Δp = 14.098 Pa and the snippet 7.903 Pa — a factor of 1.784, which is exactly sqrt(4·1e-3/π) / 0.02. Set area = π·0.02²/4 = 3.141593e-4 and both give 144.318782 Pa, agreeing to 12 significant figures (the residual 5.8e-11 Pa is @njit-vs-NumPy summation order).

Update 2026-08-25 — you no longer have to hard-code it. The configured value is now handed to your callable as state.D_h_spec, a scalar, or None when the case set none. So the snippet above reproduces the built-in exactly by reading that instead:

D_h = state.D_h_spec if state.D_h_spec is not None else state.D_h
s_mom = -(f / (2.0 * D_h)) * state.rho * np.abs(state.u) * state.u

A different name and a different type, on purpose: an array against a scalar, and None against a silent geometric fallback, so the two cannot be swapped in by accident. Use state.D_h when you genuinely want the geometric diameter — a duct whose area varies along its length has no single configured one.

The [pipe.porous] seam has no equivalent: PorousSpec carries no configured hydraulic diameter, so there is nothing there to disagree with.

solver = "quasi_steady" is refused, not ignored

The quasi-steady network folds friction into a constant drag coefficient read from friction_kind, and never calls pipe.py_friction. Your closure would simply not run.

Rather than let that happen quietly, the build raises:

UserCodeError: pipe 'tube' sets friction kind='python' ('my_models:rough_wall'),
but the case runs solver='quasi_steady', whose network solve folds friction into
a constant drag read from `friction_kind` and never calls `pipe.py_friction`.
The closure would be silently ignored.

Use solver = "explicit" or "strang", or express the model as a constant via kind = "darcy_const".

This is stricter than the neighbouring colebrook kind, which warns and runs frictionless under QS. The difference is intent: a case that happens to use Colebrook probably was not written for QS, whereas writing kind = "python" and pointing it at your own module is an explicit request for that code to run.

Interaction with the compiled kernel

kind = "python" sets pipe.py_friction and records friction_kind = 4 as provenance only. Before each step, _py_closure_kinds substitutes none for any family carrying a Python override, so the @njit kernel does flux and area only and the Python source is overlaid afterwards. Friction is applied exactly once.

The integer code exists so a snapshot or GUI round-trip can still report which closure a pipe carries, rather than none.

Why friction went first

The seam was already consumed by the solver (solver/timeloop.py, _overlay_py_sources) and had zero producers — no built-in closure and no TOML path wrote pipe.py_friction. Adding one is therefore purely additive: it cannot change the result of any case that does not opt in.

Gotchas found while building it

Two bugs were caught by asserting numbers rather than reading code, and both would have produced plausible-looking wrong physics:

  • R_gas captured at build time gave the dataclass default 287.05 instead of the case's declared 287.0, so every T was wrong in the third digit — too small to look broken, plenty to matter inside an Arrhenius exponent. The build now passes the same R_eff it packs into every other closure.
  • D_h read from friction_params[1] is identically zero for this kind, because nothing packs that array when the closure is a Python callable. A model handed D_h = 0 either divides by zero or produces no friction at all. It is now derived per cell from the actual area.

Both are pinned by tests/unit/test_udf_friction_seam.py.