Collision and State Storage

Recursive regularised BGK

The scalar collision is the recursive regularised BGK (RRBGK) operator, matching the fluid collision family and following the regularisation projection introduced by Latt and Chopard[1]. On D3Q7 only the rank-1 Hermite tensor (the vector flux \(q_\alpha\)) is non-trivial; on larger velocity sets higher-rank tensors are also reconstructed.

The non-equilibrium scalar flux is

\[q_\alpha^{neq} \equiv \sum_i g_i^{neq}\, c_{i\alpha} = \sum_i g_i\, c_{i\alpha} - \phi\, u_\alpha.\]

The reconstructed non-equilibrium populations and the post-collision state are

\[g_i^{neq} = \frac{w_i\, c_{i\alpha}}{c_{s,\phi}^2}\, q_\alpha^{neq}, \qquad g_i^{*} = g_i^{eq} + (1 - \omega_\phi)\, g_i^{neq}.\]

The relaxation rate \(\omega_\phi = 1/\tau_g\) is set node-locally from the molecular diffusivity (and the LES subgrid contribution when active, see LES Subgrid Coupling).

Post-collision macroscopics update analytically:

  • \(\phi^{*} = \phi\) (conserved by collision)

  • \(q_\alpha^{*,neq} = (1 - \omega_\phi)\, q_\alpha^{neq}\).

Macroscopic state storage

The scalar runs in memory-efficient mode, mirroring the fluid kernel. The full state stored per node is

Stored

Symbol

Count

Scalar concentration

\(\phi\)

1

Non-equilibrium scalar flux

\(q_x^{neq},\, q_y^{neq},\, q_z^{neq}\)

3

instead of the seven \(g_i\) populations. The regularised collision projects every non-equilibrium population onto the rank-1 subspace each step, which is what makes the four-macroscopic state lossless.

Key idea: storing four macroscopics, not seven populations

Because the regularised collision discards everything outside the rank-1 flux \(q_\alpha^{neq}\) at every step, \((\phi, q_x^{neq}, q_y^{neq}, q_z^{neq})\) reconstruct the full population set exactly. Storing the macroscopics rather than the populations is therefore lossless, and the scalar reuses the same memory-efficient kernel the fluid uses.

The per-step procedure for a scalar node is:

  1. Read \(\phi\) and the fluid velocity \(u_\alpha\).

  2. Compute \(g_i^{eq}\) from \((\phi, u_\alpha)\).

  3. Compute \(g_i^{neq}\) from \(q_\alpha^{neq}\).

  4. Apply collision: \(g_i^{*} = g_i^{eq} + (1 - \omega_\phi)\, g_i^{neq}\).

  5. Stream \(g_i^{*}\).

  6. Recompute \(\phi^{\text{new}} = \sum_i g_i^{\text{new}}\) and \(q_\alpha^{\text{new},neq} = \sum_i g_i^{\text{new}} c_{i\alpha} - \phi^{\text{new}} u_\alpha\).

Because the scalar reads the fluid velocity, the scalar step runs after the fluid step in every iteration.

Source term

A user-supplied source \(S(x, y, z, t, \phi)\) is added to the collision with the same half-step weighting nassu uses for the fluid body force Guo et al.[2], which keeps the recovered source second-order accurate and independent of the relaxation rate:

\[g_i \leftarrow g_i - \omega_\phi\,(g_i - g_i^{eq}) + \Delta t\,\Bigl(1 - \tfrac{\omega_\phi}{2}\Bigr)\, w_i\, S.\]

In the macroscopic (moment) view this is a half-step increment to the stored zeroth moment, applied before collision,

\[\phi = \sum_i g_i + \tfrac{\Delta t}{2}\, S,\]

exactly paralleling the \(\tfrac{\Delta t}{2\rho}\,\Delta F_\alpha\) velocity correction on the momentum side.

Why the half-step weighting and not a plain population increment

The naive addition \(\Delta g_i = w_i\, S\, \Delta t\) recovers \(S\) only to first order and leaves an error that depends on the relaxation rate \(\omega_\phi\), so the same source produces a different result as the diffusivity changes. The half-step form removes that dependence and keeps the recovered source second-order accurate.

The source is evaluated either as an Eulerian field from a symbolic expression provided in the configuration, or as a Lagrangian release at a point or surface; in the Lagrangian case the source is spread to the neighbouring Eulerian nodes via the same Roma 2-point Dirac kernel used by IBM (see Coupling with Other Modules).