MultiFab & Array4

[ax3l]

Reimplementation of the successful MFiter concepts from: https://github.com/sayerhs/pyamrex for MultiFab, avoiding copies. Reimplementation of yield continuations with an iterator.

Start new concepts for Array4 - first via __array_interface__: #17 / #9 (comment)

• we will make Array4(Real/Int) createable only by MultiFabs (amrex::Box indexing semantics by default for everything: index -1 is ghost)
• add explicit asarray / as_numpy / as_cupy functions for non-owning data from it: (numpy/Python Buffer Interface semantics for everything: index -1 is end of array)
• make data copies explicit or fail if coming from device (check also pytorch & tensorflow interfaces for good vibes)

Refs.:

• https://github.com/sayerhs/pyamrex/blob/master/amrex/amrex_base.pyx#L502-L511
• https://github.com/sayerhs/pyamrex/blob/master/amrex/cpp/multifab.pxi
• https://sayerhs.github.io/pyamrex/index.html#amrex.amrex_base.MultiFab

# assume a MultiFab mfab with 3 components
ngv = mfab.nGrowVect

# Looping
for mfi in mfab:
    bx = mfi.tilebox().grow(ngv)
    marr = mfab.array(mfi) # Array4: added __array_interface__ here

    # index by index assignment
    # - this is AMReX Array4, F-order indices
    # - even though we iterate by fastest varying index, such loops are naturally very slow in Python
    for i, j, k in bx:
        mar[i, j, k, 0] = 10.0 * i
        mar[i, j, k, 1] = 10.0 * j
        mar[i, j, k, 2] = 10.0 * k

    # numpy representation: non-copying view, including the
    # guard/ghost region
    #   note: in numpy, indices are in C-order!
    marr_np = np.array(marr, copy=False)

    # all components at all indices set at once to 42
    marr_np[:, :, :, :] = 42.
    # first component at all indices set at once to 12
    marr_np[0, :, :, :] = 12.

    # check the values at start/end are the same between MultiFab and numpy view in it: first component
    assert(marr_np[0, 0, 0, 0] == marr[bx.small_end])
    assert(marr_np[0, -1, -1, -1] == marr[bx.big_end])

[ax3l]

A comment by @WeiqunZhang

MultiFab contains a vector of ints. These ints are global box index that can be used to access Fab, MulitiFab::operator[](int K). There is a public function that returns the vector. Maybe this can be used to build an iterator in python without involving MFIter.
FabArrayBase::indexArray() if I recall it correctly.

[ax3l]

Windows memory violation: Cannot find any issue on Linux with address + undefined sanitizers so far... (besides #29)

update: I used a too-small integer type to store the pointer: numpy/numpy#21248