## 3 Implementation  
The LBM implementation can directly follow Equation (3) but, in our version, we applied some of the techniques developed for higher performance: the _fused_ kernel and the _pull_ scheme [25], which are used on different hardware platforms, and a collection of low-level, GPU-only optimisations described below. The general idea of the presented implementation is shown in Algorithm 1. The algorithm contains three main operations. ```
forallnodesdo
1 Set initial values of \(\rho,\mathbf{v}\)
2 Set \(f_{i}^{pre}=f_{i}^{eq}\)
3 repeat
4 forallnodes at positions \(\mathbf{x}\)doin parallel
5 Kernelstream and collide
6 foralldirections \(i\)do
7 gather\(f_{i}^{\prime}\gets f_{i}^{pre}(\mathbf{x}-\mathbf{v}_{i})\)
8 compute\(\rho\) and \(\mathbf{v}\) from \(f_{i}^{\prime}\) using Equation (2)
9 foralldirections \(i\)do
10 compute\(f_{i}^{eq}\) from \(\rho\) and \(\mathbf{v}\) using Equation (6)
11 computenew\(f_{i}^{\prime\prime}\) from \(f_{i}^{\prime},f_{i}^{eq}\), RHS of Equations (3) and (4)
12 store\(f_{i}^{post}(\mathbf{x})\gets f_{i}^{\prime\prime}\)
13
14 forallnodes at positions \(\mathbf{x}\)doin parallel
15 swap\(f_{i}^{post}(\mathbf{x})\,\rightleftarrows\,f_{i}^{pre}(\mathbf{x})\)
16 untilend of simulation;
```  
**Algorithm 1** General idea of the presented LBM implementation. \(\rho,\mathbf{v},f_{i}^{\prime}\) and \(f_{i}^{\prime\prime}\) denote temporary values used only inside the kernel. Lines 1-3 are responsible for the initialisation of PDFs values \(f_{i}\) that are set to equilibrium state \(f_{i}^{eq}\) computed from initial values of velocity and density. Initialisation is done only once at the beginning of computations. Then, the simulation comes down to computing values of PDFs for the successive time steps. A single iteration of lines 5-15 corresponds to a single time step. During the time step computations, all nodes are processed in parallel. To avoid race conditions, we use two copies of \(f_{i}\) functions: \(f_{i}^{pre}=f_{i}(t)\) functions were computed during the previous time step and are only read, and \(f_{i}^{post}=f_{i}(t+\delta t)\) functions are computed during the current time step and are only written. There are also known parallel implementations that use a single copy of PDFs only [26] at the cost of increased code complexity. To minimise the memory bandwidth, we implemented the _fused_ kernel (lines 6-13), where collision and streaming are conducted in one step with a single read and write of \(f_{i}\) values. Additionally, we use the reversed order of collision and streaming, also known as the _pull_ scheme [25]. Direct implementation of Equation (3) computes the collision first and then _scatters_ the new, computed \(f_{i}\) values to neighbour nodes. In the pull approach, first, the \(f_{i}\) values from the previous time step are _gathered_ from neighbour nodes, then the collision is applied, and, eventually, the new \(f_{i}\) values are stored in the current node. This scheme keeps addresses of all writes to memory aligned what may additionally decrease memory traffic because unaligned memory writes often are more costly than unaligned reads (for example due to allocate-on-write policy). After processing of all nodes, the values of \(f_{i}^{pre}\) and \(f_{i}^{post}\) are exchanged in lines 14-15. However, since we have not found an efficient way to exchange fields, then we use two kernels with identical computations, but one of the kernels reads \(f_{i}^{pre}\) and writes \(f_{i}^{post}\) and the second kernel does the opposite. The operations shown in Algorithm 1 do not include support for boundary conditions, which is also present in the implemented kernel. To detect boundary nodes, we store in memory an additional field encoding each node type (fluid, solid, boundary type) along with a bitmask containing information about which neighbour nodes are present. We also use a separate field to store values for boundary nodes with fixed conditions, e.g., constant velocity.