# 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 {footcite:t}`latt2006lattice`.
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 {doc}`les_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.

```{admonition} Key idea: storing four macroscopics, not seven populations
---
class: note
---
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$.
1. Compute $g_i^{eq}$ from $(\phi, u_\alpha)$.
1. Compute $g_i^{neq}$ from $q_\alpha^{neq}$.
1. Apply collision: $g_i^{*} = g_i^{eq} + (1 - \omega_\phi)\, g_i^{neq}$.
1. Stream $g_i^{*}$.
1. 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 {footcite:t}`Guo2002-04`, 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.

```{admonition} Why the half-step weighting and not a plain population increment
---
class: important
---
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 {doc}`coupling`).

```{footbibliography}
```