Is there a possibility to have a Struct-of-Array layout for MultiFab?

[ToguyC]

As I saw in the code and changelog, the particle mode supports SoA and AoS, but it seems that the MultiFab does only support AoS.
Is there a possibility to have SoA using MultiFab? If not, is it planned to implement it in the future?

The way I understand it currently is as follows:

for (MFIter mfi(outputs); mfi.isValid(); ++mfi) {
    const Box &box = mfi.validbox();
    auto _mf = my_multifab.array(mfi);

    ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
        // do something with _mf(i, j, k, n) with n being in the range of the components
    });
}

When going into the type of _mf (Array4), the operator overload of () is in AoS mode.

template <class U=T, std::enable_if_t<!std::is_void_v<U>,int> = 0>
[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
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];
}

Thank you in advance for your responses and understanding!

[WeiqunZhang]

Could you clarify what you mean by AoS and SoA? For a MultiFab, T is Real, not a struct. So MultiFab with multiple components is already effectively a struct of array. Something like

struct multifab {
    3d_array a1;
    3d_array a2;
    ...
}

[ToguyC]

The way I understand the Array4 () operator overload, is that the i, j, and k indices are for x, y, and z respectively, and n for the component index.

So with this in mind, I see the p[...] access as getting the cell index first using i, j, and k, and then searching the n-th component with the +n*nstride. So an array of struct, struct being the n Real components.

I may be completely wrong about how it is implemented, and would be glad to be corrected if I am indeed wrong.

[WeiqunZhang]

I think you are wrong.

[AlexanderSinn]

nstride is larger than all the other strides, making it the outermost component (SoA). AoS would be if nstride is equal to 1 and istride > 1. 

[ToguyC]

oooh ok, with a large nstride it makes sens. I did not catch this in the code sorry. Thank you for giving this precision.

[AlexanderSinn]

Currently it is:

istride = 1
jstride = ilength
kstride = ilength * jlength
nstride = ilength * jlength * klength

[ToguyC]

yeah ok, I must have missed this initialization sorry for that. Helps me big time.

[WeiqunZhang]

(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)kstride+nnstride

Well, how the strides are set are just minor details. The expression above you quoted already shows i is the unit stride direction. So n cannot be the unit stride direction too.