Fields basics
Field
s and its relatives are core Oceananigans data structures. Field
s are more or less arrays of data
located on a grid
, whose entries correspond to the average value of some quantity over some finite-sized volume. Field
s also may contain boundary_conditions
, may be computed from an operand
or expression involving other fields, and may cover only a portion of the total indices
spanned by the grid.
Staggered grids and field locations
Oceananigans ocean-flavored fluids simulations rely fundamentally on "staggered grid" numerical methods.
Recall that grids represent a physical domain divided into finite volumes. For example, let's consider a horizontally-periodic, vertically-bounded grid of cells that divide up a cube with dimensions $1 \times 1 \times 1$:
using Oceananigans
grid = RectilinearGrid(topology = (Periodic, Periodic, Bounded),
size = (4, 5, 4),
halo = (1, 1, 1),
x = (0, 1),
y = (0, 1),
z = [0, 0.1, 0.3, 0.6, 1])
# output
4×5×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── Periodic x ∈ [0.0, 1.0) regularly spaced with Δx=0.25
├── Periodic y ∈ [0.0, 1.0) regularly spaced with Δy=0.2
└── Bounded z ∈ [0.0, 1.0] variably spaced with min(Δz)=0.1, max(Δz)=0.4
The cubic domain is divided into a "primary mesh" of $4 \times 5 \times 4 = 80$ cells, which are evenly spaced in $x, y$ but variably spaced in $z$. Now, in addition to the primary mesh, the grid defines also a set of "staggered" grids whose cells are shifted by half a cell width relative to the primary mesh. In other words, the staggered grid cells have a "location" in each direction – either Center
, and therefore co-located with the primary mesh, or Face
and located over the interfaces of the primary mesh. For example, the primary or Center
cell spacings in $z$ are
zspacings(grid, Center())
# output
4-element view(OffsetArray(::Vector{Float64}, 0:5), 1:4) with eltype Float64:
0.1
0.19999999999999998
0.3
0.4
corresponding to cell interfaces located at z = [0, 0.1, 0.3, 0.6, 1]
. But then for the grid which is staggered in z
relative to the primary mesh,
zspacings(grid, Face())
# output
5-element view(OffsetArray(::Vector{Float64}, -1:5), 1:5) with eltype Float64:
0.1
0.15000000000000002
0.24999999999999994
0.3500000000000001
0.3999999999999999
The cells for the vertically staggered grid have different spacings than the primary mesh. That's because the edges of the vertically-staggered mesh coincide with the nodes (the cell centers) of the primary mesh. The nodes of the primary mesh are
znodes(grid, Center(), with_halos=true)
# output
6-element OffsetArray(::Vector{Float64}, 0:5) with eltype Float64 with indices 0:5:
-0.05
0.05
0.2
0.44999999999999996
0.8
1.2
The center of the leftmost "halo cell" is z = -0.05
, while the center of the first cell from the left is z = 0.05
. This means that the width of the first cell on the vertically-staggered grid is 0.05 - (-0.05) = 0.1
– and so on. Finally, note that the nodes of the staggered mesh coincide with the cell interfaces of the primary mesh, so:
znodes(grid, Center())
# output
4-element view(OffsetArray(::Vector{Float64}, 0:5), 1:4) with eltype Float64:
0.05
0.2
0.44999999999999996
0.8
In a three-dimensional domain, there are $2³ = 8$ meshes – 1 primary mesh, and 7 meshes that are staggered to varying degrees from the primary mesh. This system of staggered grids is commonly used in fluid dynamics and was invented specifically for simulations of the atmosphere and ocean.
Constructing Fields at specified locations
Every Field
is associated with either the primary mesh or one of the staggered meshes by a three-dimensional "location" associated with each field. To build a fully-centered Field
, for example, we write
c = Field{Center, Center, Center}(grid)
# output
4×5×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 4×5×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 6×7×6 OffsetArray(::Array{Float64, 3}, 0:5, 0:6, 0:5) with eltype Float64 with indices 0:5×0:6×0:5
└── max=0.0, min=0.0, mean=0.0
Fully-centered fields also go by the alias CenterField
,
c == CenterField(grid)
# output
true
Many fluid dynamical variables are located at cell centers – for example, tracers like temperature and salinity. Another common type of Field
we encounter have cells located over the x
-interfaces of the primary grid,
u = Field{Face, Center, Center}(grid)
# output
4×5×4 Field{Face, Center, Center} on RectilinearGrid on CPU
├── grid: 4×5×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 6×7×6 OffsetArray(::Array{Float64, 3}, 0:5, 0:6, 0:5) with eltype Float64 with indices 0:5×0:6×0:5
└── max=0.0, min=0.0, mean=0.0
which also goes by the alias u = XFaceField(grid)
. The name u
is suggestive: in the Arakawa type-C grid ('C-grid' for short) used by Oceananigans, the x
-component of the velocity field is stored at Face, Center, Center
location.
The centers of the u
cells are shifted to the left relative to the c
cells:
@show xnodes(c)
@show xnodes(u)
nothing
# output
xnodes(c) = [0.125, 0.375, 0.625, 0.875]
xnodes(u) = [0.0, 0.25, 0.5, 0.75]
Notice that the first u
-node is at x=0
, the left end of the grid, but the last u
-node is at x=0.75
. Because the x
-direction is Periodic
, the XFaceField
u
has 4 cells in x
– the cell just right of x=0.75
is the same as the cell at x=0
.
Because the vertical direction is Bounded
, however, vertically-staggered fields have more vertical cells than CenterField
s:
w = Field{Center, Center, Face}(grid)
@show znodes(c)
@show znodes(w)
nothing
# output
znodes(c) = [0.05, 0.2, 0.44999999999999996, 0.8]
znodes(w) = [0.0, 0.1, 0.3, 0.6, 1.0]
Field
s at Center, Center, Face
are also called ZFaceField
, and the vertical velocity is a ZFaceField
on the C-grid. Let's visualize the situation:
using CairoMakie
fig = Figure(size=(600, 180))
ax = Axis(fig[1, 1], xlabel="x")
# Visualize the domain
lines!(ax, [0, 1], [0, 0], color=:gray)
xc = xnodes(c)
xu = xnodes(u)
scatter!(ax, xc, 0 * xc, marker=:circle, markersize=10, label="Cell centers")
scatter!(ax, xu, 0 * xu, marker=:vline, markersize=20, label="Cell interfaces")
ylims!(ax, -1, 1)
xlims!(ax, -0.1, 1.1)
hideydecorations!(ax)
hidexdecorations!(ax, ticklabels=false, label=false)
hidespines!(ax)
Legend(fig[0, 1], ax, nbanks=2, framevisible=false)
current_figure()
Setting Field
s
Field
s are full of 0's when they are created, which is not very exciting. The situation can be improved using set!
to change the values of a field. For example,
set!(c, 42)
# output
4×5×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 4×5×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 6×7×6 OffsetArray(::Array{Float64, 3}, 0:5, 0:6, 0:5) with eltype Float64 with indices 0:5×0:6×0:5
└── max=42.0, min=42.0, mean=42.0
Now c
is filled with 42
s (for this simple case, we could also have used c .= 42
). Let's confirm that:
c[1, 1, 1]
# output
42.0
Looks good. And
c[1:4, 1:5, 1]
# output
4×5 Matrix{Float64}:
42.0 42.0 42.0 42.0 42.0
42.0 42.0 42.0 42.0 42.0
42.0 42.0 42.0 42.0 42.0
42.0 42.0 42.0 42.0 42.0
Note that indexing into c
is the same as indexing into c.data
.
c[:, :, :] == c.data
# output
true
We can also set!
with arrays,
random_stuff = rand(size(c)...)
set!(c, random_stuff)
heatmap(view(c, :, :, 1))
or even use functions to set,
fun_stuff(x, y, z) = 2x
set!(c, fun_stuff)
# output
4×5×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 4×5×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 6×7×6 OffsetArray(::Array{Float64, 3}, 0:5, 0:6, 0:5) with eltype Float64 with indices 0:5×0:6×0:5
└── max=1.75, min=0.25, mean=1.0
and plot it
heatmap(view(c, :, :, 1))
For Field
s on three-dimensional grids, set!
functions must have arguments x, y, z
for RectilinearGrid
, or λ, φ, z
for LatitudeLongitudeGrid
and OrthogonalSphericalShellGrid
. But for Field
s on one- and two-dimensional grids, only the arguments that correspond to the non-Flat
directions must be included. For example, to set!
on a one-dimensional grid we write
# Make a field on a one-dimensional grid
one_d_grid = RectilinearGrid(size=7, x=(0, 7), topology=(Periodic, Flat, Flat))
one_d_c = CenterField(one_d_grid)
# The one-dimensional grid varies only in `x`
still_pretty_fun(x) = 3x
set!(one_d_c, still_pretty_fun)
# output
7×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 7×1×1 RectilinearGrid{Float64, Periodic, Flat, Flat} on CPU with 3×0×0 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 13×1×1 OffsetArray(::Array{Float64, 3}, -2:10, 1:1, 1:1) with eltype Float64 with indices -2:10×1:1×1:1
└── max=19.5, min=1.5, mean=10.5
A bit more about setting with functions
Let's return to the three-dimensional fun_stuff
case to investigate in more detail how set!
works with functions. The xnodes
of c
– the coordinates of the center of c
's finite volumes – are:
xc = xnodes(c)
@show xc
nothing # hide
# output
xc = [0.125, 0.375, 0.625, 0.875]
To set!
the values of c
we evaluate fun_stuff
at c
's nodes, producing
c[1:4, 1, 1]
# output
4-element Vector{Float64}:
0.25
0.75
1.25
1.75
This function-setting method is a first-order method for computing the finite volume of c
to fun_stuff
. Higher-order algorithms could be implemented – have a crack if you're keen.
As a result set!
can evaluate differently on Field
s at different locations:
u = XFaceField(grid)
set!(u, fun_stuff)
u[1:4, 1, 1]
# output
4-element Vector{Float64}:
0.0
0.5
1.0
1.5
Halo regions and boundary conditions
We built grid
with halo = (1, 1, 1)
, which means that the "interior" cells of the grid are surrounded by a "halo region" of cells that's one cell thick. The number of halo cells in each direction are stored in the properties Hx, Hy, Hz
, so,
(grid.Hx, grid.Hy, grid.Hz)
# output
(1, 1, 1)
set!
doesn't touch halo cells. Check out one of the two-dimensional slices of c
showing both the interior and the halo regions:
c[:, :, 1]
# output
6×7 OffsetArray(::Matrix{Float64}, 0:5, 0:6) with eltype Float64 with indices 0:5×0:6:
0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.25 0.25 0.25 0.25 0.25 0.0
0.0 0.75 0.75 0.75 0.75 0.75 0.0
0.0 1.25 1.25 1.25 1.25 1.25 0.0
0.0 1.75 1.75 1.75 1.75 1.75 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
The interior region is populated, but the surrounding halo regions are all 0. To remedy this situation we need to fill_halo_regions!
:
using Oceananigans.BoundaryConditions: fill_halo_regions!
fill_halo_regions!(c)
c[:, :, 1]
# output
6×7 OffsetArray(::Matrix{Float64}, 0:5, 0:6) with eltype Float64 with indices 0:5×0:6:
1.75 1.75 1.75 1.75 1.75 1.75 1.75
0.25 0.25 0.25 0.25 0.25 0.25 0.25
0.75 0.75 0.75 0.75 0.75 0.75 0.75
1.25 1.25 1.25 1.25 1.25 1.25 1.25
1.75 1.75 1.75 1.75 1.75 1.75 1.75
0.25 0.25 0.25 0.25 0.25 0.25 0.25
The way the halo regions are filled depends on c.boundary_conditions
:
c.boundary_conditions
# output
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: PeriodicBoundaryCondition
├── east: PeriodicBoundaryCondition
├── south: PeriodicBoundaryCondition
├── north: PeriodicBoundaryCondition
├── bottom: FluxBoundaryCondition: Nothing
├── top: FluxBoundaryCondition: Nothing
└── immersed: FluxBoundaryCondition: Nothing
Specifically for c
above, x
and y
are Periodic
while z
has been assigned the default "no-flux" boundary conditions for a Field
with Center
location in a Bounded
direction. For no-flux boundary conditions, the halo regions of c
are filled so that derivatives evaluated on the boundary return 0. To view only the interior cells of c
we use the function interior
,
interior(c, :, :, 1)
# output
4×5 view(::Array{Float64, 3}, 2:5, 2:6, 2) with eltype Float64:
0.25 0.25 0.25 0.25 0.25
0.75 0.75 0.75 0.75 0.75
1.25 1.25 1.25 1.25 1.25
1.75 1.75 1.75 1.75 1.75
Note that the indices of c
(and the indices of c.data
) are "offset" so that index 1
corresponds to the first interior cell. As a result,
c[1:4, 1:5, 1] == interior(c, :, :, 1)
# output
true
and more generally
typeof(c.data)
# output
OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}
Thus, for example, the x
-indices of c.data
vary from 1 - Hx
to Nx + Hx
– in this case, from 0
to 5
. The underlying array can be accessed with parent(c)
. But note that the "parent" array does not have offset indices, so
@show parent(c)[1:2, 2, 2]
@show c.data[1:2, 1, 1]
nothing
# output
(parent(c))[1:2, 2, 2] = [1.75, 0.25]
c.data[1:2, 1, 1] = [0.25, 0.75]