GitHub

Working with grids

PencilArrays.jl includes functionality for conveniently working with the coordinates associated to a multidimensional domain. For this, the localgrid function can be used to construct an object describing the grid coordinates associated to the local MPI process. This object can be used to easily and efficiently perform operations that depend on the local coordinates.

Creating local grids

As an example, say we are performing a 3D simulation on a rectilinear grid, so that the coordinates of a grid point are given by $\bm{x}_{ijk} = (x_i, y_j, z_k)$ where x, y and z are separate one-dimensional arrays. For instance:

Nx, Ny, Nz = 65, 17, 21
xs = range(0, 1; length = Nx)
ys = range(-1, 1; length = Ny)
zs = range(0, 2; length = Nz)

Before continuing, let's create a domain decomposition configuration:

using MPI
using PencilArrays

MPI.Init()
comm = MPI.COMM_WORLD

pen = Pencil((Nx, Ny, Nz), comm)
Decomposition of 3D data
    Data dimensions: (65, 17, 21)
    Decomposed dimensions: (2, 3)
    Data permutation: NoPermutation()
    Array type: Array

Now, we can extract the local grid associated to the local MPI process:

grid = localgrid(pen, (xs, ys, zs))
LocalRectilinearGrid{3} with coordinates:
 (1) 0.0:0.015625:1.0
 (2) -1.0:0.125:1.0
 (3) 0.0:0.1:2.0

Note that this example was run on a single MPI process, which makes things somewhat less interesting, but the same applies to more processes. With more than one process, the local grid is a subset of the global grid defined by the coordinates (xs, ys, zs).

Using local grids

The grid object just created can be used to operate with PencilArrays. In particular, say we want to initialise a PencilArray to a function that depends on the domain coordinates, $u(x, y, z) = x + 2y + z^2$. This can be easily done using the broadcasting syntax (here we use the @. macro for convenience):

u = PencilArray{Float64}(undef, pen)  # construct PencilArray first
@. u = grid.x + 2 * grid.y + grid.z^2

Here, grid.x, grid.y and grid.z are a convenient way of extracting the three components of the grid. Alternatively, one can use the syntax grid[1], grid[2], etc..., which is in particularly useful when working in dimensions higher than 3.

Note that one could do the same as above using indexing instead of broadcasting:

for I ∈ eachindex(grid)
    x, y, z = grid[I]
    u[I] = x + 2y + z^2
end

Library

PencilArrays.LocalGrids.LocalRectilinearGrid — Type
LocalRectilinearGrid{N, Perm} <: AbstractLocalGrid{N, Perm}

Defines the local portion of a rectilinear grid in N dimensions.

A rectilinear grid is represented by a set of orthogonal coordinates (x, y, z, ...).

source
PencilArrays.LocalGrids.localgrid — Function
localgrid(x::PencilArray, args...)

Equivalent of localgrid(pencil(x), args...).

source
localgrid((xs, ys, ...), perm = NoPermutation()) -> LocalRectilinearGrid

Create a LocalRectilinearGrid from a set of orthogonal coordinates (xs, ys, ...), where each element is an AbstractVector.

Optionally, one can pass a static permutation (as in Permutation(2, 1, 3)) to change the order in which the coordinates are iterated.

source
localgrid(p::Pencil, (x_global, y_global, ...))      -> LocalRectilinearGrid
localgrid(u::PencilArray, (x_global, y_global, ...)) -> LocalRectilinearGrid

Create a LocalRectilinearGrid from a decomposition configuration and from a set of orthogonal global coordinates (x_global, y_global, ...).

In this case, each *_global is an AbstractVector describing the coordinates along one dimension of the global grid.

source

Index