Lattice-Boltzmann Pointer-Shift In-Place Streaming - possible in AMRex?

[shreyasashok]

Hello all,

I am assessing whether AMReX is appropriate for a lattice-Boltzmann (LBM) code. In order to achieve efficient advection in LBM, there are many pointer-shift based algorithms that allow for the advection step to occur in-place without actually moving any memory around. An example of such an algorithm is here:

https://www.mdpi.com/2079-3197/5/2/19

I am struggling to understand how to implement this within AMReX. In order to achieve such an algorithm, I would need to allocate extra memory in my FArrayBox to allow room for the pointer shifting, and also have some way of altering the pointer to the zeroeth element of the data. Is there any way to achieve this without extensively modifying the code underneath?

[WeiqunZhang]

I probably will not have time to read the paper. Even if I try, I may not be able to understand it because I know almost nothing about LBM. Could you summarize what the algorithm does? What kind of memory operations you need? Are they going to be performed on regular data stored in MultiFab/FArrayBox or the special sparse data you mentioned in #3191.

[WeiqunZhang]

That probably means you need a 5x5 array. Let us assume the entire domain is NxN and it can be perfectly divided into a number of MxM Boxes. I suspect something like this might be what you need.

Box domain(IntVect(0), IntVect(N-1));
BoxArray ba(domain);
ba.maxSize(M);  // M is 3 in the example above
MultiFab mf(ba, DistributionMapping{ba}, ncomp, 1); // 1 ghost cell
for (MFIter mfi(mf); mfi.isValid(); ++mfi)
{
    Array4<Real> a = mf.array(mfi);
    // make changes to the meta-data inside a, which are all public
}

[WeiqunZhang]

Note that the operator() of Array4 is

        U& operator() (int i, int j, int k, int n) const noexcept {
#if defined(AMREX_DEBUG) || defined(AMREX_BOUND_CHECK)
            index_assert(i,j,k,n);
#endif
            return p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride];
        }

You can change p, begin, jstride, kstride and nstride. You probably also need to change end, because it's used in index_assert.

[WeiqunZhang]

I took your example literally and could not understand what it could possibly mean by shifting from (0,0) to (1,2). That's why I wrote out all the detail about which cells shifted to where. But I assume the shifting is to one of the 8 neighbors. So what's wrong with simply using expressions like a(i+1,j-1,k), which is equivalent to a(i,j,k) with the pointer inside a shifted?

[shreyasashok]

Hello Weiqun,

Thanks for the extensive help with this. I will play around with what you have suggested and see if it works for what I am trying to do. I am very new to AMReX, and my question may not be quite well-informed as a result.

The reason I might want to avoid using the shifted expression a(i+1,j-1,k) is because, won't the FillBoundary() routine mess up? Let's assume a 3x3 valid box with 1 ghost layer (5x5 total), as you suggested. On an even timestep, destination ghost cell (-1,-1)'s data should lie at pointer index (p+7). The source valid cell is at index (2,2), and its pointer index should be (p+7+18). But I assume the FillBoundary routine isn't going to know about my offset of 7. It would just fill destination pointer index (p) with source pointer index (p+18). That would not work.

An additional reason I would want to avoid the shifted expression is just for convenience. In this thread, I have only asked about pointer shifting of one floating-point quantity. In the actual LBM algorithm, I would have between 9 to 27 of these quantities, and all of the quantities would be shifted by a different amount. With a simple pointer shift, I can use the same compute kernel function for both my even timesteps and odd timesteps. If I have to do the shifting in the array expression, it might end up making the code unnecessarily complex-looking. However, I am entirely willing to drop this criteria if it makes everything else much easier.

[WeiqunZhang]

It seems to me that no matter which approach you take it's better to use two ghost cells. Then you only communicate once every two steps. This is more efficient in terms of communication. As for pointer shifting vs index shifting, I would argue the latter is better. I don't think you want to have 27 pointers in a single GPU kernel. It will use a lot of registers. (Launching 27 GPU kernel with one for each variable would be even worse.) As for different shifting patterns for different variables, if you create a function that returns a 3 ints struct (e.g., amrex::Dim3), you should be able to avoid duplication. You might also be able to use function template to reduce the runtime cost of index shifting. For example

if (step %2 == 0) {
    f<true>();
} else {
    f<false>();
}

template <bool> void f();

template<>
void f<true>()
{
    int comp = ...;
    amrex::Dim3 shift = compute_shift<IsOdd>(comp);
    a(i+shift.x,j+shift.y,...);  // The compiler sees shift is zero at compile time, so it will turn it into a(i,j,k) for IsOdd=true case.
}

template <> constexpr amrex:: Dim3 compute_shift<true>(int) { return {0,0,0}; }

You can probably template on the component index too, since the number of indices is fixed.

Take my opinion with a grain of salt, because I don't fully understand what you are trying to do.