LBMNode#
Module: maddening.nodes.lbm
Stability: experimental
Algorithm ID: MADD-NODE-011
Version: 1.1.0
From 0.4.0 this node’s algorithm ID is
MADD-NODE-011. Releases 0.1.0 to 0.3.1 gave itMADD-NODE-007, whichRigidBodyNodealso carried and keeps; the algorithm is unchanged by the renumbering.
Summary#
A lattice Boltzmann solver for weakly compressible, isothermal, laminar flow on a D2Q9 or D3Q19 lattice [@Kruger2017], with a BGK collision operator, Guo body forcing [@Guo2002], bounce-back walls on an arbitrary wall mask, and Zou-He pressure (density) boundaries on a chosen inlet and outlet face [@ZouHe1997; @HechtHarting2010].
Governing Equations#
The lattice Boltzmann equation with the BGK collision operator and a forcing term,
with the second-order equilibrium
and the moments
In the low-Mach limit this recovers the Navier-Stokes equations with kinematic viscosity \(\nu = c_s^2(\tau - \tfrac12)\) and the isothermal equation of state \(p = c_s^2\rho\). The node takes viscosity and sets \(\tau = \tfrac12 + \nu / c_s^2\).
Pressure boundary. On an inlet or outlet face the node imposes a density \(\rho_p = p / c_s^2\) and zero tangential velocity, and solves for the face-normal velocity \(u_n\) and the populations that stream into the domain through the face.
Discretization#
Lattice units, \(\Delta x = \Delta t = 1\); update ignores its dt argument. One call of update is, in order:
moments \(\rho\), \(\mathbf{u}\) from \(f\) (with the Guo half-force correction);
BGK collision plus the Guo source term \(S_q = (1 - \tfrac{1}{2\tau})\,w_q\left[\frac{\mathbf{e}_q - \mathbf{u}}{c_s^2} + \frac{\mathbf{e}_q\cdot\mathbf{u}}{c_s^4}\mathbf{e}_q\right]\cdot\mathbf{F}\);
streaming, periodic on every axis (
jnp.roll), or across a halo in the sharded path, whose halos at the edges of the global grid must be filled periodically (boundary="periodic") for it to be the same step (halo_boundary(), below);bounce-back: in every wall cell each population is replaced by its opposite;
the Zou-He closure on the inlet face, then on the outlet face (fluid cells only);
moments of the result, the velocity zeroed in wall cells, \(p = c_s^2 \rho\).
Declared order of accuracy: DiscretizationOrder(spatial=2.0, temporal=None), measured on a wall-free periodic domain by MADD-VER-007. There is no independent timestep, so no temporal order is declared. Walls are first order (see the limitations).
Zou-He pressure boundary#
Let \(n\) be the face axis and \(\sigma = +1\) on a *_min face, \(\sigma = -1\) on a *_max face, so that \(\sigma\hat{\mathbf{n}}\) points into the domain. After streaming, the directions on the face split into three sets:
Set |
Definition |
Meaning |
|---|---|---|
\(U\) (unknown) |
\(e_{qn} = \sigma\) |
arrive from outside the domain; rebuilt by the closure |
\(K\) (known) |
\(e_{qn} = -\sigma\) |
streamed out of the interior; the opposites of \(U\) |
\(T\) (tangential) |
\(e_{qn} = 0\) |
the rest population included |
With \(S_X = \sum_{q\in X} f_q\), the density and the normal momentum of the face are
and eliminating the unknown \(S_U\) gives the face-normal velocity
The known populations appear twice: once in the density and once, reflected, in the momentum. For each unknown \(q\) with opposite \(\bar q\), non-equilibrium bounce-back of the normal part plus a transverse-momentum correction gives
The rebuilt face has exactly density \(\rho_p\) and zero tangential momentum. This rests on four identities of the velocity set, which _zou_he_face_closure checks for each face before the closure runs, refusing a set that fails one: the opposites of \(U\) are exactly \(K\); \(2\sum_{U} w_q e_{qn}^2 / c_s^2 = 1\), so the bounce-back terms add up to \(\rho_p u_n\); \(\sum_U e_{qt} = \sum_U w_q e_{qn} e_{qt} = 0\), so the correction leaves density and normal momentum alone; and \(\sum_U e_{qt} e_{qt'} = 0\) for \(t \neq t'\), so each tangential axis is corrected on its own. On both supported lattices \(\sum_U e_{pt}^2 = 2\) on every face, so the coefficient is \(\tfrac12\).
D2Q9 (this module’s numbering: 1 \(+x\), 2 \(-x\), 3 \(+y\), 4 \(-y\), 5 \((+,+)\), 6 \((-,+)\), 7 \((+,-)\), 8 \((-,-)\)), x_min:
x_max rebuilds \(f_2\), \(f_6\), \(f_8\) from \(f_1\), \(f_7\), \(f_5\), with \(u = -1 + [f_0 + f_3 + f_4 + 2(f_1 + f_5 + f_7)]/\rho_p\) and the signs of the \(\rho_p u\) terms reversed. These are the equations of Zou and He [@ZouHe1997] with \(u_y = 0\). The \(y\) faces follow by exchanging the roles of the axes.
D3Q19, x_min: \(f_{(1,0,0)} = f_{(-1,0,0)} + \tfrac13\rho_p u\), and for the four edge directions
where \(N_y\) and \(N_z\) sum over the nine tangential populations, the \((0,\pm1,\pm1)\) edges included. This is the D3Q19 on-site closure of Hecht and Harting [@HechtHarting2010] with zero tangential velocity. The other five faces follow by symmetry.
Wall cells on a pressure face keep all of their populations. inlet_face and outlet_face must be different faces: the constructor refuses one face for both, because the outlet closure, applied second, would overwrite the inlet’s on every cell and inlet_pressure would be silently dropped. Perpendicular faces are allowed; they share the cells along their common edge (one corner cell in 2-D, a line of cells in 3-D), and because the inlet closure runs first and the outlet closure second, those shared cells carry the outlet pressure.
The closure is also used unchanged by the sharded step update_padded, which is correct only while the inlet/outlet axis is not sharded: it closes the edge plane of its own slab, so on a sharded face axis every seam between shards would be forced to the face pressure too. That is refused – see limitation 4.
Note
Before 0.4.0 the closure computed \(u_n = \sigma[1 - (S_T + S_K)/\rho_p]\), with the factor 2 on \(S_K\) missing, and had no transverse correction. The face density came out as \(\rho_p + S_K\): 15.4% high for \(p = 0.36\) on a unit-density lattice, on every face of both lattices. The pressure-driven channels measured carried 0.58 to 0.80 of the imposed pressure drop, and outlet_pressure_avg reported the wrong pressure. Recorded as MADD-ANO-020; MADD-VER-016 is the benchmark that would have caught it.
Implementation Mapping#
Equation Term |
Implementation |
Notes |
|---|---|---|
Velocity sets \(\mathbf{e}_q\), \(w_q\), opposites |
|
numpy arrays, concrete inside |
\(f_q^{\text{eq}}(\rho, \mathbf{u})\) |
|
Second-order Hermite equilibrium |
\(\rho\), \(\rho\mathbf{u} = \sum f_q\mathbf{e}_q + \tfrac12\mathbf{F}\) |
|
Guo half-force correction |
BGK collision, \(\tau = \tfrac12 + \nu/c_s^2\) |
|
|
Guo source term \(S_q\) |
|
[@Guo2002] |
Streaming, periodic |
|
|
Streaming across a halo (sharded) |
|
Slicing into halo-exchanged neighbours |
Periodic halo fill at the global edges (the sharded step is |
|
|
Pressure face on a sharded axis refused |
|
At trace time, from |
Bounce-back in wall cells |
|
|
Runtime wall mask |
|
|
Opposite-direction map \(\bar q\) |
|
\(\mathbf{e}_{\bar q} = -\mathbf{e}_q\) |
Sets \(U\), \(K\), \(T\) of a face |
|
By the sign of \(e_{qn}\) |
Lattice identities of the closure |
|
Refuses a velocity set that fails one; returns \(\sigma\) and the transverse coefficients |
\(u_n = \sigma[1 - (S_T + 2S_K)/\rho_p]\) |
|
[@ZouHe1997]; MADD-ANO-020 was the missing 2 |
Unknown populations \(f_q\) (bounce-back + transverse correction) |
|
[@ZouHe1997] (D2Q9), [@HechtHarting2010] (D3Q19) |
\(\rho_p = p/c_s^2\) from |
|
Inlet face first, then outlet |
\(p = c_s^2\rho\) |
|
|
Outlet pressure published to a coupled node |
|
Mean over the fluid cells of the outlet face, under the runtime wall mask |
Initial condition \(f = f^{\text{eq}}(1, \mathbf{0})\) |
|
The wall mask is carried as |
|
See below |
|
Parameter bounds (\(\nu > 0\)) |
|
Log transform |
Assumptions and Simplifications#
Weakly compressible, isothermal flow: \(p = c_s^2 \rho\), valid for Mach number \(\ll 1\)
BGK single-relaxation-time collision
Rigid, impermeable walls, modelled by bounce-back in wall cells
Pressure faces impose a density and zero tangential velocity; the face-normal velocity is whatever the interior delivers
Every axis without a pressure face or a wall is periodic
Validated Physical Regimes#
Parameter |
Verified Range |
Notes |
|---|---|---|
|
0.8 – 1.0 |
The values the benchmarks run at: 0.8 (MADD-VER-003, MADD-VER-007) and 1.0 (MADD-VER-016). \(\tau > 0.5\) is required; \(\tau \gg 1\) adds numerical diffusion. The 0.501 – 2.0 this row claimed until 0.4.0 was never run |
0 – 0.21 |
Stokes regime only. Measured on the benchmarks’ own runs: 0.21 and 0.19 for MADD-VER-016’s channel at \(H = 8\) and \(16\) (\(u_{mean} H / \nu\)), 0.21 for MADD-VER-003’s pipe (\(u_{mean} \cdot 2R / \nu\)). The 0 – 100 this row claimed until 0.4.0 was never run |
|
Inlet/outlet density difference |
up to 1% of the mean |
MADD-VER-016’s range; larger drops add compressibility error |
Known Limitations and Failure Modes#
Walls are first order, straight walls included. Wall cells collide as well as reflect, so the hydrodynamic wall does not sit on the half-way plane. In MADD-VER-016’s channel it sits about 0.1 lattice units from the wall node (0.097 at \(H = 8\), \(\tau = 1\); 0.12 at \(\tau = 0.8\)) instead of 0.5. The channel is effectively about 0.8 lattice units wider than nominal, and at \(H = 8\) the centreline velocity is 20% above Hagen-Poiseuille. The excess halves with each doubling of \(H\), and body-force driving shows nearly the same excess (+19.8% at \(H = 8\)).
Pressure faces reflect acoustic waves. A pressure-driven channel settles on a viscous time scale set by the slowest acoustic mode between the two faces, which is longer than the \(H^2/\nu\) of the velocity profile.
Entrance effect of a pressure face. Near a face the gradient departs from the linear profile. In the middle half of MADD-VER-016’s channel it is 0.91% steeper than \(\Delta p / L\) at \(H = 8\) and 0.18% at \(H = 16\).
Pressure faces and sharding.
update_paddedapplies the closure to the edge plane of its own slab, which is correct only while the inlet/outlet axis is not sharded. This is enforced: the first sharded step that imposes a pressure on a face whose axis is split across devices raises, naming the axis to shard instead (measured before the check on a 16x10 D2Q9 channel sharded along its streamwise axis on 2 devices: density 0.990 / 1.010 either side of the seam and a centreline velocity 2.10x the unsharded one, with no error). It fires at the first step rather than at construction because the faces act only when a pressure input arrives; a body-force-driven node may be sharded along any axis. Separately, the halos at the edges of the global grid must be filled periodically, because that is whatupdatedoes: the node declareshalo_boundary() == "periodic"andShardedStencilNoderefuses at construction any otherboundary, its own default"edge"included, saying to passboundary="periodic". The default used to run and moved that channel’s centreline velocity by 0.64% after 200 steps when sharded across its walls; the periodic fill matches the unsharded node to float32 rounding.Compressibility errors at Mach number above about 0.1; no turbulence model; BGK is less stable than MRT at high Reynolds number.
Stability Conditions#
\(\tau > \tfrac12\), which the constructor enforces through \(\nu > 0\). The equilibrium is a low-Mach expansion, so \(|\mathbf{u}| \lesssim 0.1\) in lattice units. \(\tau\) close to \(0.5\) with sharp gradients is prone to instability under BGK.
State Variables#
Field |
Shape |
Units |
Description |
|---|---|---|---|
|
|
lattice |
Distribution functions |
|
|
lattice |
\(\rho = \sum_q f_q\) |
|
|
lattice |
\(\mathbf{u}\), zero in wall cells |
|
|
lattice |
\(p = c_s^2\rho\) |
|
|
— |
Wall cells, |
Parameters#
Parameter |
Type |
Default |
Units |
Description |
|---|---|---|---|---|
|
tuple of int |
|
cells |
2 entries for D2Q9, 3 for D3Q19 |
|
float |
0.1 |
lattice |
\(\nu\); trainable; \(\tau = \tfrac12 + \nu/c_s^2\) |
|
str |
|
— |
|
|
bool array or None |
None |
— |
True in wall cells |
|
str |
|
— |
One of |
|
str |
|
— |
One of |
Boundary Inputs#
Field |
Shape |
Default |
Description |
|---|---|---|---|
|
scalar |
none (no BC) |
Pressure imposed on |
|
scalar |
none (no BC) |
Pressure imposed on |
|
|
zero |
Guo body force per unit volume |
|
|
the state’s mask |
Runtime wall mask override |
Boundary Fluxes#
Field |
Shape |
Units |
Description |
|---|---|---|---|
|
scalar |
lattice |
Mean pressure over the fluid cells of the outlet face (runtime wall mask); equals |
References#
[@Kruger2017] Krüger, T. et al. (2017). The Lattice Boltzmann Method: Principles and Practice. Springer. — BGK, equilibria, bounce-back.
[@Guo2002] Guo, Z., Zheng, C., Shi, B. (2002). Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys. Rev. E 65, 046308. — Forcing term and half-force velocity.
[@ZouHe1997] Zou, Q., He, X. (1997). On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Phys. Fluids 9, 1591. — The pressure closure (D2Q9).
[@HechtHarting2010] Hecht, M., Harting, J. (2010). Implementation of on-site velocity boundary conditions for D3Q19 lattice Boltzmann simulations. J. Stat. Mech. P01018. — The D3Q19 transverse-momentum correction.
Verification Evidence#
MADD-VER-007 — spatial order 2 by manufactured solution (body-forced Kolmogorov flow, periodic):
tests/verification/test_mms_order.py.MADD-VER-016 — pressure-driven Poiseuille through the Zou-He faces, in absolute terms, under refinement; also the pressure drop and body-force equivalence checks:
tests/verification/test_lbm_pressure_poiseuille.py.MADD-VER-003 — body-force Hagen-Poiseuille in a pipe on the sharded path:
tests/cloud/multigpu/test_lbm_poiseuille.py.Pressure face, cell by cell (density, tangential momentum, which populations change, fixed point, Zou and He’s D2Q9 equations written out, node pressure field, reported outlet pressure, runtime wall mask), the face rules (one face for both refused, shared cells of perpendicular faces at the outlet pressure) and the periodic halo fill:
tests/nodes/test_lbm_pressure_boundary.py.Sharded pressure channel equal to the unsharded one to float32 rounding under the default halo fill; a pressure face on a sharded axis refused, other halo fills refused, body-force driving unrestricted:
tests/cloud/multigpu/test_lbm_sharded_pressure_faces.py.Structural battery (finite, deterministic, jit-consistent, differentiable, every trainable parameter effective):
tests/verification/test_builtin_nodes_verified_lbm.py.
Changelog#
Version |
Date |
Change |
|---|---|---|
1.0.0 |
2026-03-16 |
Initial implementation |
1.1.0 |
2026-09-24 |
Zou-He closure corrected: the face now carries the prescribed density and zero tangential velocity (MADD-ANO-020). |