<p>For anyone looking into this topic some hints.
Lattice-Boltzmann is generally <strong>bandwidth limited</strong>. This means its performance depends mainly on the amount of data that can be loaded from and written to memory.</p>

<ul>
<li><p>Use a highly efficient compiled programming language: <strong>C or C++</strong> are good choices for CPU-based implementations.</p></li>
<li><p>Choose an <strong>architecture with a high bandwidth</strong>. For a CPU this means high clock RAM and a lot of memory channels (quad-channel or more).</p></li>
<li><p>This makes it crucial to use an appropriate <strong>linear memory layout</strong> that makes effective use of <em>cache prefetching</em>: The data is arranged in memory in small portions, so called cache lines. Whenever a processor accesses an element the entire cache line (on modern architectures 64 Bytes) it lies in are loaded. This means 8 double or 16 float values are loaded at once! While I have not found this to be a problem for multi-core processors as they share the L3 cache this should lead to problems on many-core architectures as changes to the same cache line have to be kept in sync and problems arise when other processors are working on data that another processor is working on (<em>false sharing</em>). This can be avoided by introducing <strong>padding</strong>, meaning you add elements you won't use to fill the rest of the cache line. Assume you want to update a cell with a discretisation with 27 speeds for the D3Q27-lattice then in the case of doubles (8 Bytes) the data lies on 4 distinct cache lines. You should add 5 doubles of padding to match the 32 Bytes (4*8 Bytes).</p></li>
</ul>

<pre><code>unsigned int const PAD          = (64 - sizeof(double)*D3Q27.SPEEDS % 64); ///&lt; padding: number of doubles
size_t       const MEM_SIZE_POP = sizeof(double)*NZ*NY*NX*(SPEEDS+PAD);    ///&lt; amount of memory to be allocated

</code></pre>

<p>Most compilers naturally align the start of the array with a cache line so you don't have to take care of that.</p>

<ul>
<li>The linear indices are inconvenient for accessing. Therefore you should design the accessing as efficient as possible. You could write a wrapper class. In any case <strong>inline</strong> those <strong>functions</strong>, meaning every call is replaced by their definition in the code.</li>
</ul>

<pre><code>inline size_t const D3Q27_PopIndex(unsigned int const x, unsigned int const y, unsigned int const z, unsigned int const d)
{
    return (D3Q27.SPEEDS + D3Q27.PAD)*(NX*(NY*z + y) + x) + D3Q27.SPEEDS*p + d;
}
</code></pre>

<ul>
<li><p>Furthermore cache locality can be increased by maximising the ratio between computation and communication for example using three-dimensional <strong>spatial loop blocking</strong> (<a href="https://stackoverflow.com/questions/50847746/scaling-issues-with-openmp">Scaling issues with OpenMP</a>), meaning every code works on a cube of cells instead of a single cell.</p></li>
<li><p>Generally implementations make use of two distinct populations A and B and perform collision and streaming from one implementation into another. This means every value in memory exists twice, once pre- and once post-collision. There exist different strategies for <strong>recombining steps and storing in such a way that you only have to keep one population copy in memory</strong>. For instance see the <em>A-A pattern</em> as proposed by P. Bailey et al. - "Accelerating Lattice Boltzmann Fluid Flow Simulations Using Graphics Processors" (<a href="https://www2.cs.arizona.edu/people/pbailey/Accelerating_GPU_LBM.pdf" rel="nofollow noreferrer">https://www2.cs.arizona.edu/people/pbailey/Accelerating_GPU_LBM.pdf</a>) or the <em>Esoteric Twist</em> by M. Geier &amp; M. Schönherr - "Esoteric Twist: An Efficient in-Place Streaming Algorithmus for the Lattice Boltzmann Method on Massively Parallel Hardware" (<a href="https://pdfs.semanticscholar.org/ea64/3d63667900b60e6ff49f2746211700e63802.pdf" rel="nofollow noreferrer">https://pdfs.semanticscholar.org/ea64/3d63667900b60e6ff49f2746211700e63802.pdf</a>). I have implemented the first with the use of macros meaning every access of a population calls a macro of the form:</p></li>
</ul>

<pre><code>#define O_E(a,b)       a*odd + b*(!odd)

#define  READ_f_0      D3Q27_PopIndex(x,           y, z,           0, p)
#define  READ_f_1      D3Q27_PopIndex(O_E(x_m, x), y, z, O_E( 1,  2), p)
#define  READ_f_2      D3Q27_PopIndex(O_E(x_p, x), y, z, O_E( 2,  1), p)
...

#define  WRITE_f_0     D3Q27_PopIndex(x,           y, z,           0, p)
#define  WRITE_f_1     D3Q27_PopIndex(O_E(x_p, x), y, z, O_E( 2,  1), p)
#define  WRITE_f_2     D3Q27_PopIndex(O_E(x_m, x), y, z, O_E( 1,  2), p)
...
</code></pre>

<ul>
<li>If you have multiple interaction populations use <strong>grid merging</strong>. Lay the indices out linearly in memory and put two distinct populations side by side. The accessing of population p works then as follows:</li>
</ul>

<pre><code>inline size_t const D3Q27_PopIndex(unsigned int const x, unsigned int const y, unsigned int const z, unsigned int const d, unsigned int const p = 0)
{
    return (D3Q27.SPEEDS*D3Q27.NPOP + D3Q27.PAD)*(NX*(NY*z + y) + x) + D3Q27.SPEEDS*p + d;
}
</code></pre>

<ul>
<li><p>For a regular grid make the <em>algorithm as predictable as possible</em>. Let <strong>every cell perform collision and streaming</strong> and then do the <em>boundary conditions in reverse</em> afterwards. If you have many cells that do not contribute directly to the algorithm omit them with a <em>logical mask</em> that you can <em>store in the padding</em> as well!</p></li>
<li><p>Make everything <strong>know to the compiler at compilation time</strong>: Template for example boundary conditions with a function that takes care of index changes so you don't have to rewrite every boundary condition.</p></li>
<li><p>Modern architectures have registers that can perform <strong>SIMD</strong> operations, so the same instruction on multiple data. Some processors (AVX-512) can process up to 512 bits of data and thus 32 doubles almost as fast as a single number. This seems to be very attractive for LBM in particular ever since gathering and scattering have been introduced (<a href="https://en.wikipedia.org/wiki/Gather-scatter_(vector_addressing)" rel="nofollow noreferrer">https://en.wikipedia.org/wiki/Gather-scatter_(vector_addressing)</a>) but with the current bandwidth limitations (maybe it is worth it with DDR5 and processors with few cores) this is in my opinion not worth the hassle: The single core performance and parallel scaling is better (M. Wittmann et al. - "Lattice Boltzmann Benchmark Kernels as a Testbed for Performance Analysis" - <a href="https://arxiv.org/abs/1711.11468" rel="nofollow noreferrer">https://arxiv.org/abs/1711.11468</a>) but the overall algorithm performs not any better as it is bandwidth limited. So it only makes sense on architectures that are limited by the computing capacities rather than the bandwidth. On the Xeon Phi architecture the results seem to be remarkable Robertsen et al. - "High‐performance SIMD implementation of the lattice‐Boltzmann method on the Xeon Phi processor" (<a href="https://onlinelibrary.wiley.com/doi/abs/10.1002/cpe.5072" rel="nofollow noreferrer">https://onlinelibrary.wiley.com/doi/abs/10.1002/cpe.5072</a>).</p></li>
</ul>

<p>In my opinion most of this is not worth the effort for simple 2D implementations. Do the easy optimisations there, loop blocking, a linear memory layout but forget about the more complex access patterns. In 3D the effect can be enormous: I have achieved <em>up to 95% parallel scalability and an overall performance of over 150 Mlups with a D3Q19 on a modern 12-core processor</em>. For more performance switch to more adequate architectures like <em>GPUs</em> with <em>CUDA C</em> that are <em>optimised for bandwidth</em>.</p>
Improving memory layout for parallel computing
<p>I'm trying to optimize an algorithm (Lattice Boltzmann) for parallel computing using C++ AMP. And looking for some suggestions to optimize the memory layout, just found out that removing one parameter from the structure into another vector (the blocked vector) gave and increase of about 10%.</p>

<p>Anyone got any tips that can further improve this, or something i should take into consideration?
Below is the most time consuming function that is executed for each timestep, and the structure used for the layout.</p>

<pre><code>struct grid_cell {
//  int blocked;    // Define if blocked
    float n;        // North
    float ne;       // North-East
    float e;        // East
    float se;       // South-East
    float s;
    float sw;
    float w;
    float nw;
    float c;        // Center
};

int collision(const struct st_parameters param, vector&lt;struct grid_cell&gt; &amp;node, vector&lt;struct grid_cell&gt; &amp;tmp_node, vector&lt;int&gt; &amp;obstacle) {
    int x,y;
    int i = 0;
    float c_sq = 1.0f/3.0f;     // Square of speed of sound
    float w0 = 4.0f/9.0f;       // Weighting factors
    float w1 = 1.0f/9.0f;
    float w2 = 1.0f/36.0f;

    int chunk = param.ny/20;
    float total_density = 0;

    float u_x,u_y;              // Avrage velocities in x and y direction
    float u[9];                 // Directional velocities
    float d_equ[9];             // Equalibrium densities
    float u_sq;                 // Squared velocity
    float local_density;        // Sum of densities in a particular node

    for(y=0;y&lt;param.ny;y++) {
        for(x=0;x&lt;param.nx;x++) {
            i = y*param.nx + x; // Node index
            // Dont consider blocked cells
            if (obstacle[i] == 0) {
                // Calculate local density
                local_density = 0.0;
                local_density += tmp_node[i].n;
                local_density += tmp_node[i].e;
                local_density += tmp_node[i].s;
                local_density += tmp_node[i].w;
                local_density += tmp_node[i].ne;
                local_density += tmp_node[i].se;
                local_density += tmp_node[i].sw;
                local_density += tmp_node[i].nw;
                local_density += tmp_node[i].c;

                // Calculate x velocity component
                u_x = (tmp_node[i].e + tmp_node[i].ne + tmp_node[i].se - 
                      (tmp_node[i].w + tmp_node[i].nw + tmp_node[i].sw)) 
                      / local_density;
                // Calculate y velocity component
                u_y = (tmp_node[i].n + tmp_node[i].ne + tmp_node[i].nw - 
                      (tmp_node[i].s + tmp_node[i].sw + tmp_node[i].se)) 
                      / local_density;
                // Velocity squared
                u_sq = u_x*u_x +u_y*u_y;

                // Directional velocity components;
                u[1] =  u_x;        // East
                u[2] =        u_y;  // North
                u[3] = -u_x;        // West
                u[4] =      - u_y;  // South
                u[5] =  u_x + u_y;  // North-East
                u[6] = -u_x + u_y;  // North-West
                u[7] = -u_x - u_y;  // South-West
                u[8] =  u_x - u_y;  // South-East

                // Equalibrium densities
                // Zero velocity density: weight w0
                d_equ[0] = w0 * local_density * (1.0f - u_sq / (2.0f * c_sq));
                // Axis speeds: weight w1
                d_equ[1] = w1 * local_density * (1.0f + u[1] / c_sq
                                 + (u[1] * u[1]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[2] = w1 * local_density * (1.0f + u[2] / c_sq
                                 + (u[2] * u[2]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[3] = w1 * local_density * (1.0f + u[3] / c_sq
                                 + (u[3] * u[3]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[4] = w1 * local_density * (1.0f + u[4] / c_sq
                                 + (u[4] * u[4]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                // Diagonal speeds: weight w2
                d_equ[5] = w2 * local_density * (1.0f + u[5] / c_sq
                                 + (u[5] * u[5]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[6] = w2 * local_density * (1.0f + u[6] / c_sq
                                 + (u[6] * u[6]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[7] = w2 * local_density * (1.0f + u[7] / c_sq
                                 + (u[7] * u[7]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[8] = w2 * local_density * (1.0f + u[8] / c_sq
                                 + (u[8] * u[8]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));

                // Relaxation step
                node[i].c = (tmp_node[i].c + param.omega * (d_equ[0] - tmp_node[i].c));
                node[i].e = (tmp_node[i].e + param.omega * (d_equ[1] - tmp_node[i].e));
                node[i].n = (tmp_node[i].n + param.omega * (d_equ[2] - tmp_node[i].n));
                node[i].w = (tmp_node[i].w + param.omega * (d_equ[3] - tmp_node[i].w));
                node[i].s = (tmp_node[i].s + param.omega * (d_equ[4] - tmp_node[i].s));
                node[i].ne = (tmp_node[i].ne + param.omega * (d_equ[5] - tmp_node[i].ne));
                node[i].nw = (tmp_node[i].nw + param.omega * (d_equ[6] - tmp_node[i].nw));
                node[i].sw = (tmp_node[i].sw + param.omega * (d_equ[7] - tmp_node[i].sw));
                node[i].se = (tmp_node[i].se + param.omega * (d_equ[8] - tmp_node[i].se));
            }
        }
    }
    return 1;
}
</code></pre>