# Operations and averaging

`Field`s 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 `Field`s.
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,

```jldoctest operations
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 it only represents a computation.
Unlike `Field`s, `AbstractOperation`s do not have any data, and are associated only with minimal memory allocation.
`AbstractOperations` are generated by inflicting `Field`s with ordinary arithmetic expressions,

```jldoctest operations
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`,

```jldoctest operations
quadratic = c^2 + two_c + 1

# 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)
    ├── ^ 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 `Field`s, `AbstractOperations` have a location and a grid.
In addition to `BinaryOperation`s like the kind above, `UnaryOperation`s and `MultiaryOperation`s are also supported,

```jldoctest operations
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
```

```jldoctest operations
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,

```jldoctest operations
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 `Field`s, `AbstractOperations` are evaluated by indexing,

```jldoctest operations
@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 `Field`s, `AbstractOperation`s have a _location_. For example,

```jldoctest operations
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 `Field`s from one location to another is intrinsic to arithmetic on the staggered grid.
Consider the magnitude of the gradient of `c`:

```jldoctest operations
∇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

```jldoctest operations
∇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`,

```jldoctest operations
∇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`](@ref); some reductions include [`Average`](@ref), [`Integral`](@ref), and [`CumulativeIntegral`](@ref).
Integrals come with a second interface, introduced below, that spells them with Unicode integral signs and returns a computed `Field` instead of a `Reduction`.

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.


```jldoctest operations_avg_int
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
├── 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:

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

# output

true
```

The Unicode interface writes the same integral more compactly, in notation that mirrors the derivative operators `∂x`, `∂y`, and `∂z`.
[`∫dx`](@ref), [`∫dy`](@ref), and [`∫dz`](@ref) integrate over one dimension, [`∫∫dxdy`](@ref), [`∫∫dxdz`](@ref), and [`∫∫dydz`](@ref) integrate over two, and [`∫∫∫dxdydz`](@ref), aliased to [`∫dV`](@ref), integrates over all three.
Each of them builds the corresponding `Integral` and computes it into a `Field`:

```jldoctest operations_avg_int
∫∫dxdy(c) == Field(Integral(c, dims=(1, 2)))

# output

true
```

Because the result is a `Field` and not a `Reduction`, it can be used inside other abstract operations.
Here we divide the horizontal integral at each level by the integral over the whole domain, which for `c = 1` returns the inverse of the domain depth:

```jldoctest operations_avg_int
Field(∫∫dxdy(c) / ∫dV(c))

# output
1×1×5 Field{Nothing, Nothing, Center} reduced over dims = (1, 2) on LatitudeLongitudeGrid on CPU
├── grid: 60×10×5 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: BinaryOperation at (⋅, ⋅, Center)
├── 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=0.001, min=0.001, mean=0.001
```

The one-dimensional operators also take `cumulative=true`, which returns the running integral computed by [`CumulativeIntegral`](@ref), and `reverse=true`, which accumulates from the far end of the dimension:

```jldoctest operations_avg_int
∫ᶻc = ∫dz(c, cumulative=true)

# output
60×10×6 Field{Center, Center, Face} on LatitudeLongitudeGrid on CPU
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (60, 10, 6)
├── grid: 60×10×5 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── operand: CumulativeIntegral of BinaryOperation at (Center, Center, Center) over dims 3
├── status: time=0.0
└── data: 66×16×12 OffsetArray(::Array{Float64, 3}, -2:63, -2:13, -2:9) with eltype Float64 with indices -2:63×-2:13×-2:9
    └── max=1000.0, min=0.0, mean=500.0
```

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.

```jldoctest operations_avg_int
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
├── 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

```jldoctest operations_avg_int
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.

```jldoctest operations_avg_int
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
```

The `condition` and `mask` keyword arguments behave the same way with the Unicode operators, as in `∫∫dxdy(c, condition=cond)`.
