Array wrappers

The PencilArrays module defines types for handling MPI-distributed data.

The most important types are:

PencilArray: array wrapper including MPI decomposition information. Takes local indices starting at 1, regardless of the location of each MPI process on the global topology.
GlobalPencilArray: PencilArray wrapper that takes global indices, which generally don't start at 1. See also Global views.

Construction

An uninitialised PencilArray can be constructed from a Pencil instance as

pencil = Pencil(#= ... =#)
A = PencilArray{Float64}(undef, pencil)
parent(A)  # returns the Array wrapped by `A`

This allocates a new Array with the local dimensions and data type associated to the Pencil.

One can also construct a PencilArray wrapper from an existing AbstractArray, whose dimensions must be compatible with the Pencil configuration. For instance, the following works:

dims = size_local(pencil, MemoryOrder())  # dimensions must be in memory order!
data = zeros(dims)
A = PencilArray(pencil, data)

Note that data does not need to be a Array, but can be any subtype of AbstractArray.

It is also possible to construct higher dimensional arrays, as in:

data = zeros(dims..., 3, 2)
A = PencilArray(pencil, data)

This will construct a PencilArray where the rightmost dimensions (called extra dimensions in the PencilArrays API) will never be split among MPI processes.

Unlike the wrapped AbstractArray, the PencilArray wrapper takes indices in logical order. For instance, if the underlying permutation of the Pencil is (2, 3, 1), then A[i, j, k] points to the same value as parent(A)[j, k, i].

Global views

PencilArrays are accessed using local indices that start at 1, regardless of the location of the subdomain associated to the local process on the global grid. Sometimes it may be more convenient to use global indices describing the position of the local process in the domain. For this, the global_view function is provided that generates an OffsetArray wrapper taking global indices.

For more details, see for instance the gradient example in the PencilFFTs docs.

Global views

Regular PencilArrays have more functionality than global view wrappers. This includes broadcasting, which is currently not supported for global views. In general it should be preferred to work with PencilArrays.

Types

PencilArrays.PencilArray — Type

PencilArray(pencil::Pencil, data::AbstractArray{T,N})

Create array wrapper with pencil decomposition information.

The array dimensions and element type must be consistent with those of the given pencil.

Index permutations

If the Pencil has an associated index permutation, then data must have its dimensions permuted accordingly (in memory order).

Unlike data, the resulting PencilArray should be accessed with unpermuted indices (in logical order).

Example

Suppose pencil has local dimensions (10, 20, 30) before permutation, and has an asociated permutation (2, 3, 1). Then:

data = zeros(20, 30, 10)       # parent array (with dimensions in memory order)

u = PencilArray(pencil, data)  # wrapper with dimensions (10, 20, 30)
@assert size_local(u) === (10, 20, 30)

u[15, 25, 5]          # BoundsError (15 > 10 and 25 > 20)
u[5, 15, 25]          # correct
parent(u)[15, 25, 5]  # correct

Extra dimensions

The data array can have one or more extra dimensions to the right (slow indices), which are not affected by index permutations.

Example

dims = (20, 30, 10)
PencilArray(pencil, zeros(dims...))        # works (scalar)
PencilArray(pencil, zeros(dims..., 3))     # works (3-component vector)
PencilArray(pencil, zeros(dims..., 4, 3))  # works (4×3 tensor)
PencilArray(pencil, zeros(3, dims...))     # fails

PencilArray{T}(undef, pencil::Pencil, [extra_dims...])

Allocate an uninitialised PencilArray that can hold data in the local pencil.

Extra dimensions, for instance representing vector components, can be specified. These dimensions are added to the rightmost (slowest) indices of the resulting array.

Example

Suppose pencil has local dimensions (20, 10, 30). Then:

PencilArray{Float64}(undef, pencil)        # array dimensions are (20, 10, 30)
PencilArray{Float64}(undef, pencil, 4, 3)  # array dimensions are (20, 10, 30, 4, 3)

More examples:

julia> pen = Pencil((20, 10, 12), MPI.COMM_WORLD);

julia> u = PencilArray{Float64}(undef, pen);

julia> summary(u)
"20×10×12 PencilArray{Float64, 3}(::Pencil{3, 2, NoPermutation, Array})"

julia> PencilArray{Float64}(undef, pen, 4, 3) |> summary
"20×10×12×4×3 PencilArray{Float64, 5}(::Pencil{3, 2, NoPermutation, Array})"