Partitioned multidimensional arrays (generalize PVector into PArray)

[stevengj]

For parallel stencil computations (e.g. finite-difference methods), it would be nice to have partitioned multidimensional arrays, with support for arbitrary thickness "ghost" overlap regions (so that you could loop over just the interior of each partition and the ghost regions would handle the stencil boundaries).

Not sure whether that is in scope for this package, or if it is something that should be implemented in another package on top of PVector?

It would be nice to at least broaden PVector to PArray (with PVector as an alias for PArray{1}) — same machinery, just multidimensional arrays as storage and either arrays of linear indices or arrays of CartesianIndex for the index_partition. Matrix-free stencil support could then be implemented in an add-on package.

[fverdugo]

Yes, this is in the scope definitively and multidimensional stencil computations are already supported via PVector with at most a layer of ghosts. For implicit methods, I guess this is what you need since you would need a PVector to represent A*x = b anyway. But I agree that for explicit or matrix free, then a multi dimensional array would be more natural.

n_per_dir = (10,10,10)
np_per_dir = (2,3,2)
ghost_per_dir = (true,true,true)
periodic_per_dir = (false,true,true)
ranks = LinearIndices(n_per_dir)
index_partition = uniform_partition(ranks,np_per_dir,n_per_dir,ghost_per_dir,periodic_per_dir)
v = pzeros(index_partition)

The Pvector v is the 1d version of a 10x10x10 multidimensional array split in 2x3x2 parts with ghost in all directions and periodic boundaries (influences the definition of ghost) in the 2nd and 3rd axis.

[stevengj]

Note that you might need a ghost layer of thickness more than 1 for higher-order finite-difference stencils.

Yes, I was indeed thinking of matrix-free and explicit methods.

[fverdugo]

Generalizing the thickness of the ghost layer should be very easy. I think it is as easy as changing the hardcoded 1 by ghost in start-=1 and stop+=1 in this function:

https://github.com/fverdugo/PartitionedArrays.jl/blob/3349755f92f80ec4d04020fcafcb6bb03859fd38/src/p_range.jl#L702

If I am not wrong, this would allow using ghost_per_dir = (2,2,2) in addition to ghost_per_dir = (true,true,true) which would correspond to ghost_per_dir = (1,1,1) . Luckily, this is orthogonal to the PArray thing.