GitHub

Operations and averaging

Fields are more than mere vessels for data. They come equipped with a powerful infrastructure for assembling expression trees that represent arithmetic and discrete calculus operations. We also supply a system for computing reductions (averages, integrals, and cumulative integrals) of Fields. This infrastructure can be used to construct initial conditions, set up diagnostic calculations that are performed and saved to disk while a simulation runs, and also for post-processing.

We start by constructing a CenterField on a simple grid,

using Oceananigans

grid = RectilinearGrid(topology = (Periodic, Flat, Bounded),
                       size = (4, 4),
                       x = (0, 2π),
                       z = (-4, 0))

c = CenterField(grid)

periodic_but_decaying(x, z) = sin(x) * exp(z)
set!(c, periodic_but_decaying)

# output
4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 10×1×10 OffsetArray(::Array{Float64, 3}, -2:7, 1:1, -2:7) with eltype Float64 with indices -2:7×1:1×-2:7
    └── max=0.428882, min=-0.428882, mean=1.04083e-17

An AbstractOperation (or operation for short) differs from a Field in that only represents a computation. Unlike Fields, AbstractOperations do not have any data, and are associated only with minimal memory allocation. AbstractOperations are generated by inflicting Fields with ordinary arithmetic expressions,

two_c = 2 * c

# output
BinaryOperation at (Center, Center, Center)
├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
└── tree:
    * at (Center, Center, Center)
    ├── 2
    └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU

and even by chaining expressions together, which may themselves include AbstractOperations,

quadratic = c^2 + two_c + 1

# output
BinaryOperation at (Center, Center, Center)
├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
└── tree:
    + at (Center, Center, Center)
    ├── + at (Center, Center, Center)
    │   ├── ^ at (Center, Center, Center)
    │   │   ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
    │   │   └── 2
    │   └── * at (Center, Center, Center)
        │   ├── 2
        │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
    └── 1

Like Fields, AbstractOperations have a location and a grid. In addition to BinaryOperations like the kind above, UnaryOperations and MultiaryOperations are also supported,

cos_c = cos(c)

# output
UnaryOperation at (Center, Center, Center)
├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
└── tree:
    cos at (Center, Center, Center) via identity
    └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
four_c = c + c + c + c

# output
MultiaryOperation at (Center, Center, Center)
├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
└── tree:
    + at (Center, Center, Center)
    ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
    ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
    ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
    └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU

UnaryOperation, BinaryOperation and MultiaryOperation all have both an "operator", and between 1 and many. Last, and definitely not least, the fourth flavor of AbstractOperation represents a derivative,

dx_c = ∂x(c)

# output
Derivative at (Face, Center, Center)
├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
└── tree:
    ∂xᶠᶜᶜ at (Face, Center, Center) via identity
    └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
Note

The $x$-derivative of c is computed by invoking the function ∂xᶠᶜᶜ, where the superscript ᶠᶜᶜ indicates that the result of this function is located at (Face, Center, Center). This convention is used throughout Oceananigans source code. A third superscripted letter ᵃ stands for "any" location and applies to operators that are location-agnostic in the given direction.

Like Fields, AbstractOperations are evaluated by indexing,

@show c[1, 1, 1]
@show quadratic[1, 1, 1]
nothing

# output
c[1, 1, 1] = 0.02135277459201165
quadratic[1, 1, 1] = 1.0431614901668005

Also like Fields, AbstractOperations have a location. For example,

using Oceananigans: location

@show location(c)
@show location(dx_c)
nothing

# output
location(c) = (Center, Center, Center)
location(dx_c) = (Face, Center, Center)

Notice that the location of dx_c is shifted in x relative to c. Likewise, y-derivatives are shifted in y and z-derivatives are shifted in z.

Locations and interpolation

Reconstruction of Fields from one location to another is intrinsic to arithmetic on the staggered grid. Consider the magnitude of the gradient of c:

∇c² = ∂x(c)^2 + ∂z(c)^2

# output
BinaryOperation at (Face, Center, Center)
├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
└── tree:
    + at (Face, Center, Center)
    ├── ^ at (Face, Center, Center)
    │   ├── ∂xᶠᶜᶜ at (Face, Center, Center) via identity
    │   │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
    │   └── 2
    └── ^ at (Center, Center, Face)
        ├── ∂zᶜᶜᶠ at (Center, Center, Face) via identity
        │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
        └── 2

Because ∂x(c)^2 is located at (Face, Center, Center) and ∂z(c)^2 is located at (Center, Center, Face), a decision has to be made to compute ∇c². By default, AbstractOperations are reconstructed at the location of the first object in the expression. So

∇c²_ccf = ∂z(c)^2 + ∂x(c)^2

# output
BinaryOperation at (Center, Center, Face)
├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
└── tree:
    + at (Center, Center, Face)
    ├── ^ at (Center, Center, Face)
    │   ├── ∂zᶜᶜᶠ at (Center, Center, Face) via identity
    │   │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
    │   └── 2
    └── ^ at (Face, Center, Center)
        ├── ∂xᶠᶜᶜ at (Face, Center, Center) via identity
        │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
        └── 2

ends up at (Center, Center, Face). To control the location of an operation we use the macro @at,

∇c²_ccc = @at (Center, Center, Center) ∂x(c)^2 + ∂z(c)^2

# output
BinaryOperation at (Center, Center, Center)
├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
└── tree:
    + at (Center, Center, Center)
    ├── ^ at (Center, Center, Center)
    │   ├── ∂xᶠᶜᶜ at (Center, Center, Center) via ℑxᶜᵃᵃ
    │   │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
    │   └── 2
    └── ^ at (Center, Center, Center)
        ├── ∂zᶜᶜᶠ at (Center, Center, Center) via ℑzᵃᵃᶜ
        │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
        └── 2

Averages and integrals

An operation that reduces one or more dimensions of a field is called a Reduction; some reductions include Average, Integral, and CumulativeIntegral.

Let's demonstrate now how we can compute the average or the integral of a field over the grid or over some part of the grid. We start by creating a latitude-longitude grid that only goes up to 30 degrees latitude. Conveniently, with this latitude extent that grid covers half the total area of the sphere, i.e., 2π * grid.radius^2.

Let's try to estimate this area using Integral operation. We create a Field, we fill it with ones and we integrate it over the whole grid. We use a CenterField for the example below, but all location combinations work similarly.

using Oceananigans

grid = LatitudeLongitudeGrid(size=(60, 10, 5),
                             latitude = (-30, 30),
                             longitude = (0, 360),
                             z = (-1000, 0))

c = CenterField(grid)
set!(c, 1)

∫c = Field(Integral(c, dims=(1, 2)))

# output
1×1×5 Field{Nothing, Nothing, Center} reduced over dims = (1, 2) on LatitudeLongitudeGrid on CPU
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (1, 1, 5)
├── grid: 60×10×5 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── operand: Integral of BinaryOperation at (Center, Center, Center) over dims (1, 2)
├── status: time=0.0
└── data: 1×1×11 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:8) with eltype Float64 with indices 1:1×1:1×-2:8
    └── max=2.55032e14, min=2.55032e14, mean=2.55032e14

A few remarks: note that the ∫c has locations Nothing, Nothing, Center; this is because we have integrated in the first two dimensions and thus it's reduced over dims = (1, 2).

Above we see that the max, min and mean of the field are all the same. Let's check that these values are what we expect:

∫c[1, 1, 1] ≈ 2π * grid.radius^2 # area of spherical zone with |φ| ≤ 30ᵒ

# output

true

We can further have conditional reduced operations. Let's compute the above integral but only for North hemisphere.

First we need to define a condition which is a function with arguments (i, j, k, grid, field) that returns true or false. In this example we use Oceananigans.Grids.φnode to check whether the latitude is within the range we'd like.

using Oceananigans.Grids: φnode

cond(i, j, k, grid, c) = φnode(j, grid, Center()) ≥ 0

conditional_∫c = Field(Integral(c, dims=(1, 2), condition=cond)) # integrate only when condition is true

# output
1×1×5 Field{Nothing, Nothing, Center} reduced over dims = (1, 2) on LatitudeLongitudeGrid on CPU
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (1, 1, 5)
├── grid: 60×10×5 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── operand: Integral of ConditionalOperation of BinaryOperation at (Center, Center, Center) with condition cond (generic function with 1 method) over dims (1, 2)
├── status: time=0.0
└── data: 1×1×11 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:8) with eltype Float64 with indices 1:1×1:1×-2:8
    └── max=1.27516e14, min=1.27516e14, mean=1.27516e14

Above we have attached a condition to the operand. Now the operand is applied only when the condition is true. Let's see if that is 1/4 of the area of the sphere

conditional_∫c[1, 1, 1] ≈ π * grid.radius^2 # area of spherical zone with 0ᵒ ≤ φ ≤ 30ᵒ

# output
true

Another way to do the above is to provide the condition keyword argument with an array of booleans.

cond_array = trues(size(grid))
cond_array[:, 1:5, :] .= false # set the first half of the latitude range to false

conditional_∫c = Field(Integral(c, dims=(1, 2), condition=cond_array))

conditional_∫c[1, 1, 1] ≈ π * grid.radius^2 # area of spherical zone with 0ᵒ ≤ φ ≤ 30ᵒ

# output
true