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:

PartitionedArrays.jl/src/p_range.jl

Line 702 in 3349755

function local_range(p,np,n,ghost=false,periodic=false)

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.

[ZhangFengshun]

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:

PartitionedArrays.jl/src/p_range.jl

Line 702 in 3349755

function local_range(p,np,n,ghost=false,periodic=false)

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.

Hi~, I changed the "start-=1 and stop+=1" to "start-=2 and stop+=2" in the local_range function. This change works in the uniform_partition function but failed in the consistent! function.


The pvector of rank 2 should be [1, 1, 2, 2] as expected.

[ZhangFengshun]

This error is working as expected, and I find that the block_with_constant_size function should also be updated for consistency with the 2 ghost layers.

In this way, we got the right consistent! :

[fverdugo]

Hi @ZhangFengshun this looks promising, but we do not want to hard code 2 as the thickness of the ghost layer. Instead we want to use the value of variable 'ghost' already available in these functions. Now 'ghost' is assumed to be a bool, but we want to generalise it to an arbitrary integer. Feel free to propose the required changes.

[fverdugo]

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).

Done in #173