Python assembly: closures on a pipe¶

Part of the pure-Python assembly ladder. Read that page first — it covers Pipe, junctions, boundary conditions and the two surfaces. This one covers the source terms a pipe carries: friction, wall heat and porous drag, plus the seam for supplying your own.

Every number below was produced by the code above it, on one duct: L = 0.2 m, N = 24, A = 1e-3 m², 8 g/s of 800 K air in at the left, 1 bar out at the right, advanced to steady state at t = 0.05 s.

The three built-in seams¶

Closures are configured through their spec objects and materialised by build_network. Each is independent — a pipe can carry all three.

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

GAMMA, R_GAS = 1.4, 287.05

def duct(**closures):
    return PipeSpec(id="duct", length=0.2, n_cells=24, area=1.0e-3,
                    init=UniformInit(p=1.0e5, T=800.0, u=20.0), **closures)

def case_of(pipe):
    return Case(
        case=CaseHeader(name="closures", t_end=0.05, cfl=0.4),
        gas=GasSpec(gamma=GAMMA, R_gas=R_GAS),
        pipes=[pipe], 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)])

Now add one closure at a time:

bare     = run_case(case_of(duct()))[0]["duct"]
friction = run_case(case_of(duct(
    friction=FrictionSpec(kind="darcy_const", f=0.02, D_h=0.02))))[0]["duct"]
heat     = run_case(case_of(duct(
    heat=HeatSpec(kind="constant_twall", h=300.0, T_wall=400.0,
                  D_h=0.02))))[0]["duct"]
porous   = run_case(case_of(duct(
    porous=PorousSpec(kind="darcy_forchheimer",
                      alpha=5.0e3, beta=2.0e2))))[0]["duct"]

Pipe	Δp (Pa)	outlet T (K)
bare	0.0006	799.832
`+ friction` (`darcy_const`, f = 0.02)	14.098	799.834
`+ wall heat` (`constant_twall`, 400 K wall)	−99.543	492.170
`+ porous` (`darcy_forchheimer`)	31350.464	802.691

Three things in that table are worth reading rather than skimming.

Friction raises the gas temperature, 799.832 → 799.834 K. That is not a bug and it is not a wall-heat term: in the total-energy formulation the kinetic energy that s_mom removes becomes internal energy by conservation. Porous drag does the same thing harder — 802.691 K, because it removes 31 kPa rather than 14 Pa. If you write your own friction closure, this is the trap: subtracting u * s_mom from the energy equation as well double-debits it and silently sheds the frictional heat. See UDF: friction.

The bare pipe's Δp is 0.0006 Pa, not 0. That is the discrete residual of a converged solve, and it is the scale against which the other rows should be read — a "no closure" row that came out at exactly zero would more likely mean the run never started than that the physics is clean.

Wall heat gives a negative Δp. The gas cools from 800 K to 492 K, so it densifies, so ρu² falls and pressure rises along the duct — a Rayleigh-type effect, not a mistake. tests/integration/test_driver_steady_parity.py leans on the same behaviour.

Supplying your own closure¶

Two seams take user code today, and the honest status table is in the UDF overview: friction and porous have a kind = "python" path; wall heat does not — that seam is live but only built-in closures write it, for reasons that page records.

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

def rough_wall(state):
    """Reimplements darcy_const: -(f / (2*D_h)) * rho * |u| * u."""
    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)

spec = FrictionSpec(kind="python", path="my_models:rough_wall", D_h=0.02)

That function is the built-in darcy_const model written out, so the two ought to agree. On the duct above they do not — 14.098 Pa against 7.903 Pa:

A = 1e-3 m²  (the duct above)
  built-in dp =   14.097987 Pa
  python   dp =    7.903302 Pa
  max|dp_cell| = 6.335471e+00 Pa   bit-identical=False

state.D_h and FrictionSpec.D_h are different numbers

They have the same name and no relationship. FrictionSpec.D_h is the value you set. state.D_h — what a UDF receives — is geometric, derived from the pipe's area as sqrt(4A/π):

FrictionSpec.D_h            = 0.02
geometric sqrt(4A/π), A=1e-3 = 0.035682
ratio                        = 1.784124

1.784 is exactly 14.098 / 7.903. The closure was right; it was handed a hydraulic diameter 78 % larger than the one the built-in used, because the duct's area does not correspond to a 20 mm round pipe.

Give the pipe an area that does — A = πD²/4 = 3.141593e-04 m² — and the two agree to twelve significant figures:

A = π D²/4  (geometric D_h == 0.02 exactly)
  built-in dp =  144.318782 Pa
  python   dp =  144.318782 Pa
  max|dp_cell| = 5.820766e-11 Pa   bit-identical=False

The 5.8e-11 Pa residual is floating-point ordering, not physics: the built-in is a fused @njit kernel and the UDF is NumPy, so the two sum the same terms in different orders. Agreement to 12 figures is the correct expectation here; bit-identity is not.

Two things to take from this. Use the number you mean — if your closure wants a specified hydraulic diameter, write it in the closure rather than reading state.D_h, which answers a different question. And start by reproducing a built-in you can compare against, confirm you get its answer, then change the model: a closure wrong by a factor of 1.78 is much harder to notice once it is the only thing computing that term.

The permission model is not decorative¶

Resolving "my_models:rough_wall" imports and executes that module, so it is off unless you turn it on:

import os
os.environ["AEOLUS_ALLOW_USER_CODE"] = "1"    # the API equivalent of
                                              # the CLI's --allow-user-code

Without it the run refuses, rather than quietly falling back to no friction:

UserCodeError: case names a user closure ('my_models:rough_wall') but user
code is not enabled.

The permission is deliberately an environment variable and not a CaseHeader field — a case file that can grant itself permission is not a permission, since the whole point is that opening someone else's TOML must not run their code.

Where to go next¶

Next	Page
Species columns and chemistry on a pipe	Species and chemistry
Catalysts, filters, and the quasi-steady driver	Aftertreatment
Every closure kind and its fields	Closures — reference table
The shared UDF contract	User-defined functions