Library

Documenting the public user interface.

Oceananigans.jl

Oceananigans.OceananigansModule

Main module for Oceananigans.jl – a Julia software for fast, friendly, flexible, data-driven, ocean-flavored fluid dynamics on CPUs and GPUs.

source

Abstract operations

Oceananigans.AbstractOperations.volumeConstant
volume = VolumeMetric()

Instance of VolumeMetric that generates BinaryOperations between AbstractFields and their cell volumes. Summing this BinaryOperation yields an integral of AbstractField over the domain.

Example

julia> using Oceananigans

julia> using Oceananigans.AbstractOperations: volume

julia> c = CenterField(RectilinearGrid(size=(2, 2, 2), extent=(1, 2, 3)));

julia> c .= 1;

julia> c_dV = c * volume
BinaryOperation at (Center, Center, Center)
├── grid: 2×2×2 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
    * at (Center, Center, Center)
    ├── 2×2×2 Field{Center, Center, Center} on RectilinearGrid on CPU
    └── Vᶜᶜᶜ at (Center, Center, Center)

julia> c_dV[1, 1, 1]
0.75

julia> sum(c_dV)
6.0
source
Oceananigans.AbstractOperations.ΔzConstant
Δz = ZSpacingMetric()

Instance of ZSpacingMetric that generates BinaryOperations between AbstractFields and the vertical grid spacing evaluated at the same location as the AbstractField.

Δx and Δy play a similar role for horizontal grid spacings.

Example

julia> using Oceananigans

julia> using Oceananigans.AbstractOperations: Δz

julia> c = CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 2, 3)));

julia> c_dz = c * Δz # returns BinaryOperation between Field and GridMetricOperation
BinaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
    * at (Center, Center, Center)
    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
    └── Δzᶜᶜᶜ at (Center, Center, Center)

julia> c .= 1;

julia> c_dz[1, 1, 1]
3.0
source
Oceananigans.AbstractOperations.AverageMethod
Average(field::AbstractField; condition = nothing, mask = 0, dims=:)

Return Reduction representing a spatial average of field over dims.

Over regularly-spaced dimensions this is equivalent to a numerical mean!.

Over dimensions of variable spacing, field is multiplied by the appropriate grid length, area or volume, and divided by the total spatial extent of the interval.

source
Oceananigans.AbstractOperations.BinaryOperationMethod
BinaryOperation{LX, LY, LZ}(op, a, b, ▶a, ▶b, grid)

Return an abstract representation of the binary operation op(▶a(a), ▶b(b)) on grid, where ▶a and ▶b interpolate a and b to locations (LX, LY, LZ).

source
Oceananigans.AbstractOperations.ConditionalOperationMethod
ConditionalOperation(operand::AbstractField;
                     func = identity,
                     condition = nothing,
                     mask = 0)

Return an abstract representation of a masking procedure applied when condition is satisfied on a field described by func(operand).

Positional arguments

  • operand: The AbstractField to be masked (it must have a grid property!)

Keyword arguments

  • func: A unary transformation applied element-wise to the field operand at locations where condition == true. Default is identity.

  • condition: either a function of (i, j, k, grid, operand) returning a Boolean, or a 3-dimensional Boolean AbstractArray. At locations where condition == false, operand will be masked by mask

  • mask: the scalar mask

condition_operand is a convenience function used to construct a ConditionalOperation

condition_operand(func::Function, operand::AbstractField, condition, mask) = ConditionalOperation(operand; func, condition, mask)

Example

julia> using Oceananigans

julia> using Oceananigans.Fields: condition_operand

julia> c = CenterField(RectilinearGrid(size=(2, 1, 1), extent=(1, 1, 1)));

julia> f(i, j, k, grid, c) = i < 2; d = condition_operand(cos, c, f, 10)
ConditionalOperation at (Center, Center, Center)
├── operand: 2×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── func: cos (generic function with 31 methods)
├── condition: f (generic function with 1 method)
└── mask: 10

julia> d[1, 1, 1]
1.0

julia> d[2, 1, 1]
10
source
Oceananigans.AbstractOperations.IntegralMethod
Integral(field::AbstractField; condition = nothing, mask = 0, dims=:)

Return a Reduction representing a spatial integral of field over dims.

Example

Compute the integral of $f(x, y, z) = x y z$ over the domain $(x, y, z) ∈ [0, 1] × [0, 1] × [0, 1]$. The analytical answer is $∭ x y z \, dx \, dy \, dz = 1/8$.

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(8, 8, 8), x=(0, 1), y=(0, 1), z=(0, 1));

julia> f = CenterField(grid);

julia> set!(f, (x, y, z) -> x * y * z)
8×8×8 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 14×14×14 OffsetArray(::Array{Float64, 3}, -2:11, -2:11, -2:11) with eltype Float64 with indices -2:11×-2:11×-2:11
    └── max=0.823975, min=0.000244141, mean=0.125

julia> ∫f = Integral(f)
sum! over dims (1, 2, 3) of BinaryOperation at (Center, Center, Center)
└── operand: BinaryOperation at (Center, Center, Center)
    └── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo

julia> ∫f = Field(Integral(f));

julia> compute!(∫f);

julia> ∫f[1, 1, 1]
0.125
source
Oceananigans.AbstractOperations.KernelFunctionOperationMethod
KernelFunctionOperation{LX, LY, LZ}(kernel_function, grid, arguments...)

Construct a KernelFunctionOperation at location (LX, LY, LZ) on grid with arguments.

kernel_function is called with

kernel_function(i, j, k, grid, arguments...)

Note that compute!(kfo::KernelFunctionOperation) calls compute! on all kfo.arguments.

Examples

Construct a KernelFunctionOperation that returns random numbers:

using Oceananigans

grid = RectilinearGrid(size=(1, 8, 8), extent=(1, 1, 1));

random_kernel_function(i, j, k, grid) = rand(); # use CUDA.rand on the GPU

kernel_op = KernelFunctionOperation{Center, Center, Center}(random_kernel_function, grid)

# output

KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: random_kernel_function (generic function with 1 method)
└── arguments: ()

Construct a KernelFunctionOperation using the vertical vorticity operator used internally to compute vertical vorticity on all grids:

using Oceananigans.Operators: ζ₃ᶠᶠᶜ # called with signature ζ₃ᶠᶠᶜ(i, j, k, grid, u, v)

model = HydrostaticFreeSurfaceModel(; grid);

u, v, w = model.velocities;

ζ_op = KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᶜ, grid, u, v)

# output

KernelFunctionOperation at (Face, Face, Center)
├── grid: 1×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: ζ₃ᶠᶠᶜ (generic function with 1 method)
└── arguments: ("1×8×8 Field{Face, Center, Center} on RectilinearGrid on CPU", "1×8×8 Field{Center, Face, Center} on RectilinearGrid on CPU")
source
Oceananigans.AbstractOperations.∂xMethod
∂x(L::Tuple, arg::AbstractField)

Return an abstract representation of an $x$-derivative acting on field arg followed by interpolation to L, where L is a 3-tuple of Faces and Centers.

source
Oceananigans.AbstractOperations.∂yMethod
∂y(L::Tuple, arg::AbstractField)

Return an abstract representation of a $y$-derivative acting on field arg followed by interpolation to L, where L is a 3-tuple of Faces and Centers.

source
Oceananigans.AbstractOperations.∂zMethod
∂z(L::Tuple, arg::AbstractField)

Return an abstract representation of a $z$-derivative acting on field arg followed by interpolation to L, where L is a 3-tuple of Faces and Centers.

source
Oceananigans.AbstractOperations.@binaryMacro
@binary op1 op2 op3...

Turn each binary function in the list (op1, op2, op3...) into a binary operator on Oceananigans.Fields for use in AbstractOperations.

Note: a binary function is a function with two arguments: for example, +(x, y) is a binary function.

Also note: a binary function in Base must be imported to be extended: use import Base: op; @binary op.

Example

julia> using Oceananigans, Oceananigans.AbstractOperations

julia> using Oceananigans.AbstractOperations: BinaryOperation, AbstractGridMetric, choose_location

julia> plus_or_times(x, y) = x < 0 ? x + y : x * y
plus_or_times (generic function with 1 method)

julia> @binary plus_or_times
Set{Any} with 6 elements:
  :+
  :/
  :^
  :-
  :*
  :plus_or_times

julia> c, d = (CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) for i = 1:2);

julia> plus_or_times(c, d)
BinaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
    plus_or_times at (Center, Center, Center)
    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
    └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
source
Oceananigans.AbstractOperations.@multiaryMacro
@multiary op1 op2 op3...

Turn each multiary operator in the list (op1, op2, op3...) into a multiary operator on Oceananigans.Fields for use in AbstractOperations.

Note that a multiary operator:

  • is a function with two or more arguments: for example, +(x, y, z) is a multiary function;
  • must be imported to be extended if part of Base: use import Base: op; @multiary op;
  • can only be called on Oceananigans.Fields if the "location" is noted explicitly; see example.

Example

julia> using Oceananigans, Oceananigans.AbstractOperations

julia> harmonic_plus(a, b, c) = 1/3 * (1/a + 1/b + 1/c)
harmonic_plus (generic function with 1 method)

julia> c, d, e = Tuple(CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) for i = 1:3);

julia> harmonic_plus(c, d, e) # before magic @multiary transformation
BinaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
    * at (Center, Center, Center)
    ├── 0.3333333333333333
    └── + at (Center, Center, Center)
        ├── / at (Center, Center, Center)
        │   ├── 1
        │   └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
        ├── / at (Center, Center, Center)
        │   ├── 1
        │   └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
        └── / at (Center, Center, Center)
            ├── 1
            └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU

julia> @multiary harmonic_plus
Set{Any} with 3 elements:
  :+
  :harmonic_plus
  :*

julia> harmonic_plus(c, d, e)
MultiaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
    harmonic_plus at (Center, Center, Center)
    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
    └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
source
Oceananigans.AbstractOperations.@unaryMacro
@unary op1 op2 op3...

Turn each unary function in the list (op1, op2, op3...) into a unary operator on Oceananigans.Fields for use in AbstractOperations.

Note: a unary function is a function with one argument: for example, sin(x) is a unary function.

Also note: a unary function in Base must be imported to be extended: use import Base: op; @unary op.

Example

julia> using Oceananigans, Oceananigans.Grids, Oceananigans.AbstractOperations

julia> square_it(x) = x^2
square_it (generic function with 1 method)

julia> @unary square_it
Set{Any} with 10 elements:
  :+
  :sqrt
  :square_it
  :cos
  :exp
  :interpolate_identity
  :-
  :tanh
  :sin
  :abs

julia> c = CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)));

julia> square_it(c)
UnaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
    square_it at (Center, Center, Center) via identity
    └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
source

Advection

Oceananigans.Advection.VectorInvariantMethod
VectorInvariant(; vorticity_scheme               = EnstrophyConserving(), 
                  vorticity_stencil              = VelocityStencil(),
                  vertical_scheme                = EnergyConserving(),
                  kinetic_energy_gradient_scheme = vertical_scheme)

Return a vector invariant momentum advection scheme.

Keyword arguments

  • vorticity_scheme: Scheme used for Center reconstruction of vorticity. Default: EnstrophyConserving(). Options: - UpwindBiased() - WENO() - EnergyConserving() - EnstrophyConserving()

  • vorticity_stencil: Stencil used for smoothness indicators for WENO schemes. Default: VelocityStencil(). Options: - VelocityStencil() (smoothness based on horizontal velocities) - DefaultStencil() (smoothness based on variable being reconstructed)

  • vertical_scheme: Scheme used for vertical advection of horizontal momentum. Default: EnergyConserving().

  • kinetic_energy_gradient_scheme: Scheme used for kinetic energy gradient reconstruction. Default: vertical_scheme.

  • divergence_scheme: Scheme used for divergence flux. Only upwinding schemes are supported. Default: vorticity_scheme.

  • upwinding: Treatment of upwinded reconstruction of divergence and kinetic energy gradient. Default: OnlySelfUpwinding(). Options: - CrossAndSelfUpwinding() - OnlySelfUpwinding() - VelocityUpwinding()

  • multi_dimensional_stencil : whether or not to use a horizontal two-dimensional stencil for the reconstruction of vorticity, divergence and kinetic energy gradient. Currently the "tangential" direction uses 5th-order centered WENO reconstruction.

Examples

julia> using Oceananigans

julia> VectorInvariant()
Vector Invariant, Dimension-by-dimension reconstruction 
 Vorticity flux scheme: 
 └── EnstrophyConserving{Float64} 
 Vertical advection / Divergence flux scheme: 
 └── EnergyConserving{Float64}
julia> using Oceananigans

julia> VectorInvariant(vorticity_scheme = WENO(), vertical_scheme = WENO(order = 3))
Vector Invariant, Dimension-by-dimension reconstruction 
 Vorticity flux scheme: 
 ├── WENO reconstruction order 5 
 └── smoothness ζ: Oceananigans.Advection.VelocityStencil()
 Vertical advection / Divergence flux scheme: 
 ├── WENO reconstruction order 3
 └── upwinding treatment: OnlySelfUpwinding 
 KE gradient and Divergence flux cross terms reconstruction: 
 └── Centered reconstruction order 2
 Smoothness measures: 
 ├── smoothness δU: FunctionStencil f = divergence_smoothness
 ├── smoothness δV: FunctionStencil f = divergence_smoothness
 ├── smoothness δu²: FunctionStencil f = u_smoothness
 └── smoothness δv²: FunctionStencil f = v_smoothness
      
source
Oceananigans.Advection.WENOType
WENO([FT=Float64;] 
     order = 5,
     grid = nothing, 
     zweno = true, 
     bounds = nothing)

Construct a weigthed essentially non-oscillatory advection scheme of order order.

Keyword arguments

  • order: The order of the WENO advection scheme. Default: 5
  • grid: (defaults to nothing)
  • zweno: When true implement a Z-WENO formulation for the WENO weights calculation. (defaults to true)

Examples

julia> using Oceananigans

julia> WENO()
WENO reconstruction order 5
 Smoothness formulation: 
    └── Z-weno  
 Boundary scheme: 
    └── WENO reconstruction order 3
 Symmetric scheme: 
    └── Centered reconstruction order 4
 Directions:
    ├── X regular 
    ├── Y regular 
    └── Z regular
julia> using Oceananigans

julia> Nx, Nz = 16, 10;

julia> Lx, Lz = 1e4, 1e3;

julia> chebychev_spaced_z_faces(k) = - Lz/2 - Lz/2 * cos(π * (k - 1) / Nz);

julia> grid = RectilinearGrid(size = (Nx, Nz), halo = (4, 4), topology=(Periodic, Flat, Bounded),
                              x = (0, Lx), z = chebychev_spaced_z_faces);

julia> WENO(grid; order=7)
WENO reconstruction order 7
 Smoothness formulation: 
    └── Z-weno  
 Boundary scheme: 
    └── WENO reconstruction order 5
 Symmetric scheme: 
    └── Centered reconstruction order 6
 Directions:
    ├── X regular 
    ├── Y regular 
    └── Z stretched
source
Oceananigans.Advection.div_UcMethod
div_uc(i, j, k, grid, advection, U, c)

Calculate the divergence of the flux of a tracer quantity $c$ being advected by a velocity field, $𝛁⋅(𝐯 c)$,

1/V * [δxᶜᵃᵃ(Ax * u * ℑxᶠᵃᵃ(c)) + δyᵃᶜᵃ(Ay * v * ℑyᵃᶠᵃ(c)) + δzᵃᵃᶜ(Az * w * ℑzᵃᵃᶠ(c))]

which ends up at the location ccc.

source
Oceananigans.Advection.div_𝐯uMethod
div_𝐯u(i, j, k, grid, advection, U, u)

Calculate the advection of momentum in the $x$-direction using the conservative form, $𝛁⋅(𝐯 u)$,

1/Vᵘ * [δxᶠᵃᵃ(ℑxᶜᵃᵃ(Ax * u) * ℑxᶜᵃᵃ(u)) + δy_fca(ℑxᶠᵃᵃ(Ay * v) * ℑyᵃᶠᵃ(u)) + δz_fac(ℑxᶠᵃᵃ(Az * w) * ℑzᵃᵃᶠ(u))]

which ends up at the location fcc.

source
Oceananigans.Advection.div_𝐯vMethod
div_𝐯v(i, j, k, grid, advection, U, v)

Calculate the advection of momentum in the $y$-direction using the conservative form, $𝛁⋅(𝐯 v)$,

1/Vʸ * [δx_cfa(ℑyᵃᶠᵃ(Ax * u) * ℑxᶠᵃᵃ(v)) + δyᵃᶠᵃ(ℑyᵃᶜᵃ(Ay * v) * ℑyᵃᶜᵃ(v)) + δz_afc(ℑxᶠᵃᵃ(Az * w) * ℑzᵃᵃᶠ(w))]

which ends up at the location cfc.

source
Oceananigans.Advection.div_𝐯wMethod
div_𝐯w(i, j, k, grid, advection, U, w)

Calculate the advection of momentum in the $z$-direction using the conservative form, $𝛁⋅(𝐯 w)$,

1/Vʷ * [δx_caf(ℑzᵃᵃᶠ(Ax * u) * ℑxᶠᵃᵃ(w)) + δy_acf(ℑzᵃᵃᶠ(Ay * v) * ℑyᵃᶠᵃ(w)) + δzᵃᵃᶠ(ℑzᵃᵃᶜ(Az * w) * ℑzᵃᵃᶜ(w))]

which ends up at the location ccf.

source

Architectures

Oceananigans.Architectures.CPUType
CPU <: AbstractArchitecture

Run Oceananigans on one CPU node. Uses multiple threads if the environment variable JULIA_NUM_THREADS is set.

source

Boundary conditions

Oceananigans.BoundaryConditions.BoundaryConditionMethod
BoundaryCondition(Classification::DataType, condition)

Construct a boundary condition of type BC with a number or array as a condition.

Boundary condition types include Periodic, Flux, Value, Gradient, and Open.

source
Oceananigans.BoundaryConditions.BoundaryConditionMethod
BoundaryCondition(Classification::DataType, condition::Function;
                  parameters = nothing,
                  discrete_form = false,
                  field_dependencies=())

Construct a boundary condition of type Classification with a function boundary condition.

By default, the function boudnary condition is assumed to have the 'continuous form' condition(ξ, η, t), where t is time and ξ and η vary along the boundary. In particular:

  • On x-boundaries, condition(y, z, t).
  • On y-boundaries, condition(x, z, t).
  • On z-boundaries, condition(x, y, t).

If parameters is not nothing, then function boundary conditions have the form func(ξ, η, t, parameters), where ξ and η are spatial coordinates varying along the boundary as explained above.

If discrete_form = true, the function condition is assumed to have the "discrete form",

condition(i, j, grid, clock, model_fields)

where i, and j are indices that vary along the boundary. If discrete_form = true and parameters is not nothing, the function condition is called with

condition(i, j, grid, clock, model_fields, parameters)
source
Oceananigans.BoundaryConditions.FieldBoundaryConditionsType
FieldBoundaryConditions(grid, location, indices=(:, :, :);
                        west     = default_auxiliary_bc(topology(grid, 1)(), location[1]()),
                        east     = default_auxiliary_bc(topology(grid, 1)(), location[1]()),
                        south    = default_auxiliary_bc(topology(grid, 2)(), location[2]()),
                        north    = default_auxiliary_bc(topology(grid, 2)(), location[2]()),
                        bottom   = default_auxiliary_bc(topology(grid, 3)(), location[3]()),
                        top      = default_auxiliary_bc(topology(grid, 3)(), location[3]()),
                        immersed = NoFluxBoundaryCondition())

Return boundary conditions for auxiliary fields (fields whose values are derived from a model's prognostic fields) on grid and at location.

Keyword arguments

Keyword arguments specify boundary conditions on the 6 possible boundaries:

  • west, left end point in the x-direction where i = 1
  • east, right end point in the x-direction where i = grid.Nx
  • south, left end point in the y-direction where j = 1
  • north, right end point in the y-direction where j = grid.Ny
  • bottom, right end point in the z-direction where k = 1
  • top, right end point in the z-direction where k = grid.Nz
  • immersed: boundary between solid and fluid for immersed boundaries

If a boundary condition is unspecified, the default for auxiliary fields and the topology in the boundary-normal direction is used:

  • PeriodicBoundaryCondition for Periodic directions
  • GradientBoundaryCondition(0) for Bounded directions and Centered-located fields
  • nothing for Bounded directions and Face-located fields
  • nothing for Flat directions and/or Nothing-located fields
source
Oceananigans.BoundaryConditions.FieldBoundaryConditionsType
FieldBoundaryConditions(; kwargs...)

Return a template for boundary conditions on prognostic fields.

Keyword arguments

Keyword arguments specify boundary conditions on the 7 possible boundaries:

  • west: left end point in the x-direction where i = 1
  • east: right end point in the x-direction where i = grid.Nx
  • south: left end point in the y-direction where j = 1
  • north: right end point in the y-direction where j = grid.Ny
  • bottom: right end point in the z-direction where k = 1
  • top: right end point in the z-direction where k = grid.Nz
  • immersed: boundary between solid and fluid for immersed boundaries

If a boundary condition is unspecified, the default for prognostic fields and the topology in the boundary-normal direction is used:

  • PeriodicBoundaryCondition for Periodic directions
  • NoFluxBoundaryCondition for Bounded directions and Centered-located fields
  • ImpenetrableBoundaryCondition for Bounded directions and Face-located fields
  • nothing for Flat directions and/or Nothing-located fields
source
Oceananigans.BoundaryConditions.FluxType
struct Flux <: AbstractBoundaryConditionClassification

A classification specifying a boundary condition on the flux of a field.

The sign convention is such that a positive flux represents the flux of a quantity in the positive direction. For example, a positive vertical flux implies a quantity is fluxed upwards, in the $+z$ direction.

Due to this convention, a positive flux applied to the top boundary specifies that a quantity is fluxed upwards across the top boundary and thus out of the domain. As a result, a positive flux applied to a top boundary leads to a reduction of that quantity in the interior of the domain; for example, a positive, upwards flux of heat at the top of the domain acts to cool the interior of the domain. Conversely, a positive flux applied to the bottom boundary leads to an increase of the quantity in the interior of the domain. The same logic holds for east, west, north, and south boundaries.

source
Oceananigans.BoundaryConditions.GradientType
struct Gradient <: AbstractBoundaryConditionClassification

A classification specifying a boundary condition on the derivative or gradient of a field. Also called a Neumann boundary condition.

source
Oceananigans.BoundaryConditions.OpenType
struct Open <: AbstractBoundaryConditionClassification

A classification that specifies the halo regions of a field directly.

For fields located at Faces, Open also specifies field value on the boundary.

Open boundary conditions are used to specify the component of a velocity field normal to a boundary and can also be used to describe nested or linked simulation domains.

source
Oceananigans.BoundaryConditions.ValueType
struct Value <: AbstractBoundaryConditionClassification

A classification specifying a boundary condition on the value of a field. Also called a Dirchlet boundary condition.

source

Buoyancy models

Oceananigans.BuoyancyModels.BuoyancyMethod
Buoyancy(; model, gravity_unit_vector=NegativeZDirection())

Construct a buoyancy given a buoyancy model. Optional keyword argument gravity_unit_vector can be used to specify the direction of gravity (default NegativeZDirection()). The buoyancy acceleration acts in the direction opposite to gravity.

Example

using Oceananigans

grid = RectilinearGrid(size=(1, 8, 8), extent=(1, 1, 1))

θ = 45 # degrees
g̃ = (0, -sind(θ), -cosd(θ))

buoyancy = Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=g̃)

model = NonhydrostaticModel(; grid, buoyancy, tracers=:b)

# output

NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = Tuple{Float64, Float64, Float64}
└── coriolis: Nothing
source
Oceananigans.BuoyancyModels.LinearEquationOfStateType
LinearEquationOfState([FT=Float64;] thermal_expansion=1.67e-4, haline_contraction=7.80e-4)

Return LinearEquationOfState for SeawaterBuoyancy with thermal_expansion coefficient and haline_contraction coefficient. The buoyancy perturbation $b$ for LinearEquationOfState is

\[ b = g (α T - β S),\]

where $g$ is gravitational acceleration, $α$ is thermal_expansion, $β$ is haline_contraction, $T$ is temperature, and $S$ is practical salinity units.

Default constants in units inverse Kelvin and practical salinity units for thermal_expansion and haline_contraction, respectively, are taken from Table 1.2 (page 33) of Vallis, "Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation" (2nd ed, 2017).

source
Oceananigans.BuoyancyModels.SeawaterBuoyancyType
SeawaterBuoyancy([FT = Float64;]
                 gravitational_acceleration = g_Earth,
                 equation_of_state = LinearEquationOfState(FT),
                 constant_temperature = nothing,
                 constant_salinity = nothing)

Return parameters for a temperature- and salt-stratified seawater buoyancy model with a gravitational_acceleration constant (typically called $g$), and an equation_of_state that related temperature and salinity (or conservative temperature and absolute salinity) to density anomalies and buoyancy.

Setting constant_temperature to something that is not nothing indicates that buoyancy depends only on salinity. For a nonlinear equation of state, the value provided constant_temperature is used as the temperature of the system. Vice versa, setting constant_salinity indicates that buoyancy depends only on temperature.

For a linear equation of state, the values of constant_temperature or constant_salinity are irrelevant.

Examples

The "TEOS10" equation of state, see https://www.teos-10.org

julia> using SeawaterPolynomials.TEOS10: TEOS10EquationOfState

julia> teos10 = TEOS10EquationOfState()
BoussinesqEquationOfState{Float64}:
    ├── seawater_polynomial: TEOS10SeawaterPolynomial{Float64}
    └── reference_density: 1020.0

Buoyancy that depends on both temperature and salinity

julia> using Oceananigans

julia> buoyancy = SeawaterBuoyancy(equation_of_state=teos10)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: BoussinesqEquationOfState{Float64}

Buoyancy that depends only on salinity with temperature held at 20 degrees Celsius

julia> salinity_dependent_buoyancy = SeawaterBuoyancy(equation_of_state=teos10, constant_temperature=20) 
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
├── constant_temperature: 20
└── equation_of_state: BoussinesqEquationOfState{Float64}

Buoyancy that depends only on temperature with salinity held at 35 psu

julia> temperature_dependent_buoyancy = SeawaterBuoyancy(equation_of_state=teos10, constant_salinity=35)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
├── constant_salinity: 35
└── equation_of_state: BoussinesqEquationOfState{Float64}
source
Oceananigans.BuoyancyModels.SeawaterBuoyancyType
SeawaterBuoyancy{FT, EOS, T, S} <: AbstractBuoyancyModel{EOS}

BuoyancyModels model for seawater. T and S are either nothing if both temperature and salinity are active, or of type FT if temperature or salinity are constant, respectively.

source
Oceananigans.BuoyancyModels.∂x_bMethod
∂x_b(i, j, k, grid, b::SeawaterBuoyancy, C)

Returns the $x$-derivative of buoyancy for temperature and salt-stratified water,

\[∂_x b = g ( α ∂_x T - β ∂_x S ) ,\]

where $g$ is gravitational acceleration, $α$ is the thermal expansion coefficient, $β$ is the haline contraction coefficient, $T$ is conservative temperature, and $S$ is absolute salinity.

Note: In Oceananigans, model.tracers.T is conservative temperature and model.tracers.S is absolute salinity.

Note that $∂_x T$ (∂x_T), $∂_x S$ (∂x_S), $α$, and $β$ are all evaluated at cell interfaces in x and cell centers in y and z.

source
Oceananigans.BuoyancyModels.∂y_bMethod
∂y_b(i, j, k, grid, b::SeawaterBuoyancy, C)

Returns the $y$-derivative of buoyancy for temperature and salt-stratified water,

\[∂_y b = g ( α ∂_y T - β ∂_y S ) ,\]

where $g$ is gravitational acceleration, $α$ is the thermal expansion coefficient, $β$ is the haline contraction coefficient, $T$ is conservative temperature, and $S$ is absolute salinity.

Note: In Oceananigans, model.tracers.T is conservative temperature and model.tracers.S is absolute salinity.

Note that $∂_y T$ (∂y_T), $∂_y S$ (∂y_S), $α$, and $β$ are all evaluated at cell interfaces in y and cell centers in x and z.

source
Oceananigans.BuoyancyModels.∂z_bMethod
∂z_b(i, j, k, grid, b::SeawaterBuoyancy, C)

Returns the vertical derivative of buoyancy for temperature and salt-stratified water,

\[∂_z b = N^2 = g ( α ∂_z T - β ∂_z S ) ,\]

where $g$ is gravitational acceleration, $α$ is the thermal expansion coefficient, $β$ is the haline contraction coefficient, $T$ is conservative temperature, and $S$ is absolute salinity.

Note: In Oceananigans, model.tracers.T is conservative temperature and model.tracers.S is absolute salinity.

Note that $∂_z T$ (∂z_T), $∂_z S$ (∂z_S), $α$, and $β$ are all evaluated at cell interfaces in z and cell centers in x and y.

source

Coriolis

Oceananigans.Coriolis.BetaPlaneType
BetaPlane([T=Float64;] f₀=nothing, β=nothing,
                       rotation_rate=Ω_Earth, latitude=nothing, radius=R_Earth)

Return a $β$-plane Coriolis parameter, $f = f₀ + β y$.

The user may specify both f₀ and β, or the three parameters rotation_rate, latitude (in degrees), and radius that specify the rotation rate and radius of a planet, and the central latitude (where $y = 0$) at which the β-plane approximation is to be made.

If f₀ and β are not specified, they are calculated from rotation_rate, latitude, and radius according to the relations f₀ = 2 * rotation_rate * sind(latitude) and β = 2 * rotation_rate * cosd(latitude) / radius.

By default, the rotation_rate and planet radius are assumed to be Earth's.

source
Oceananigans.Coriolis.BetaPlaneType
struct BetaPlane{T} <: AbstractRotation

A parameter object for meridionally increasing Coriolis parameter (f = f₀ + β y) that accounts for the variation of the locally vertical component of the rotation vector with latitude.

source
Oceananigans.Coriolis.ConstantCartesianCoriolisType
ConstantCartesianCoriolis([FT=Float64;] fx=nothing, fy=nothing, fz=nothing,
                                        f=nothing, rotation_axis=ZDirection(), 
                                        rotation_rate=Ω_Earth, latitude=nothing)

Return a parameter object for a constant rotation decomposed into the x, y, and z directions. In oceanography the components x, y, z correspond to the directions east, north, and up. This constant rotation can be specified in three different ways:

  • Specifying all components fx, fy and fz directly.
  • Specifying the Coriolis parameter f and (optionally) a rotation_axis (which defaults to the z direction if not specified).
  • Specifying latitude (in degrees) and (optionally) a rotation_rate in radians per second (which defaults to Earth's rotation rate).
source
Oceananigans.Coriolis.ConstantCartesianCoriolisType
struct ConstantCartesianCoriolis{FT} <: AbstractRotation

A Coriolis implementation that accounts for the locally vertical and possibly both local horizontal components of a constant rotation vector. This is a more general implementation of FPlane, which only accounts for the locally vertical component.

source
Oceananigans.Coriolis.FPlaneType
FPlane([FT=Float64;] f=nothing, rotation_rate=Ω_Earth, latitude=nothing)

Return a parameter object for constant rotation at the angular frequency f/2, and therefore with background vorticity f, around a vertical axis. If f is not specified, it is calculated from rotation_rate and latitude (in degrees) according to the relation f = 2 * rotation_rate * sind(latitude).

By default, rotation_rate is assumed to be Earth's.

Also called FPlane, after the "f-plane" approximation for the local effect of a planet's rotation in a planar coordinate system tangent to the planet's surface.

source
Oceananigans.Coriolis.HydrostaticSphericalCoriolisMethod
HydrostaticSphericalCoriolis([FT=Float64;]
                             rotation_rate = Ω_Earth,
                             scheme = EnergyConserving())

Return a parameter object for Coriolis forces on a sphere rotating at rotation_rate. By default, rotation_rate is assumed to be Earth's.

Keyword arguments

  • scheme: Either EnergyConserving() (default), EnstrophyConserving(), or ActiveCellEnstrophyConserving().
source
Oceananigans.Coriolis.NonTraditionalBetaPlaneType
NonTraditionalBetaPlane(FT=Float64;
                        fz=nothing, fy=nothing, β=nothing, γ=nothing,
                        rotation_rate=Ω_Earth, latitude=nothing, radius=R_Earth)

The user may directly specify fz, fy, β, γ, and radius or the three parameters rotation_rate, latitude (in degrees), and radius that specify the rotation rate and radius of a planet, and the central latitude (where $y = 0$) at which the non-traditional β-plane approximation is to be made.

If fz, fy, β, and γ are not specified, they are calculated from rotation_rate, latitude, and radius according to the relations fz = 2 * rotation_rate * sind(latitude), fy = 2 * rotation_rate * cosd(latitude), β = 2 * rotation_rate * cosd(latitude) / radius, and γ = - 4 * rotation_rate * sind(latitude) / radius.

By default, the rotation_rate and planet radius is assumed to be Earth's.

source
Oceananigans.Coriolis.NonTraditionalBetaPlaneType
struct NonTraditionalBetaPlane{FT} <: AbstractRotation

A Coriolis implementation that accounts for the latitudinal variation of both the locally vertical and the locally horizontal components of the rotation vector. The "traditional" approximation in ocean models accounts for only the locally vertical component of the rotation vector (see BetaPlane).

This implementation is based off of section 5 of Dellar (2011) and it conserves energy, angular momentum, and potential vorticity.

References

Dellar, P. (2011). Variations on a beta-plane: Derivation of non-traditional beta-plane equations from Hamilton's principle on a sphere. Journal of Fluid Mechanics, 674, 174-195. doi:10.1017/S0022112010006464

source

Diagnostics

Oceananigans.Diagnostics.CFLMethod
CFL(Δt [, timescale = Oceananigans.Advection.cell_advection_timescale])

Return an object for computing the Courant-Freidrichs-Lewy (CFL) number associated with time step Δt or TimeStepWizard and timescale.

See also AdvectiveCFL and DiffusiveCFL.

source
Oceananigans.Diagnostics.StateCheckerMethod
StateChecker(; schedule, fields)

Returns a StateChecker that logs field information (minimum, maximum, mean) for each field in a named tuple of fields when schedule actuates.

source
Oceananigans.Diagnostics.AdvectiveCFLMethod
AdvectiveCFL(Δt)

Return an object for computing the Courant-Freidrichs-Lewy (CFL) number associated with time step Δt or TimeStepWizard and the time scale for advection across a cell. The advective CFL is, e.g., $U Δt / Δx$.

Example

julia> using Oceananigans

julia> model = NonhydrostaticModel(grid = RectilinearGrid(size=(16, 16, 16), extent=(8, 8, 8)));

julia> Δt = 1.0;

julia> cfl = AdvectiveCFL(Δt);

julia> model.velocities.u .= π;

julia> cfl(model)
6.283185307179586
source
Oceananigans.Diagnostics.DiffusiveCFLMethod
DiffusiveCFL(Δt)

Returns an object for computing the diffusive Courant-Freidrichs-Lewy (CFL) number associated with time step Δt or TimeStepWizard and the time scale for diffusion across a cell associated with model.closure. The diffusive CFL, e.g., for viscosity is $ν Δt / Δx²$.

The maximum diffusive CFL number among viscosity and all tracer diffusivities is returned.

Example

julia> using Oceananigans

julia> model = NonhydrostaticModel(grid = RectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1)),
                                   closure = ScalarDiffusivity(; ν = 1e-2));

julia> Δt = 0.1;

julia> dcfl = DiffusiveCFL(Δt);

julia> dcfl(model)
0.256
source

Distributed

Oceananigans.DistributedComputations.DistributedType
Distributed(child_architecture = CPU(); 
            topology, 
            partition,
            devices = nothing, 
            communicator = MPI.COMM_WORLD)

Constructor for a distributed architecture that uses MPI for communications

Positional arguments

  • child_architecture: Specifies whether the computation is performed on CPUs or GPUs. Default: child_architecture = CPU().

Keyword arguments

  • synchronized_communication: if true, always use synchronized communication through ranks

  • ranks (required): A 3-tuple (Rx, Ry, Rz) specifying the total processors in the x, y and z direction. NOTE: support for distributed z direction is limited, so Rz = 1 is strongly suggested.

  • devices: GPU device linked to local rank. The GPU will be assigned based on the local node rank as such devices[node_rank]. Make sure to run --ntasks-per-node <= --gres=gpu. If nothing, the devices will be assigned automatically based on the available resources

  • communicator: the MPI communicator, MPI.COMM_WORLD. This keyword argument should not be tampered with if not for testing or developing. Change at your own risk!

source
Oceananigans.DistributedComputations.DistributedFFTBasedPoissonSolverMethod
DistributedFFTBasedPoissonSolver(global_grid, local_grid)

Return a FFT-based solver for the Poisson equation,

\[∇²φ = b\]

for Distributeditectures.

Supported configurations

We support two "modes":

  1. Vertical pencil decompositions: two-dimensional decompositions in $(x, y)$ for three dimensional problems that satisfy either Nz > Rx or Nz > Ry.

  2. One-dimensional decompositions in either $x$ or $y$.

Above, Nz = size(global_grid, 3) and Rx, Ry, Rz = architecture(local_grid).ranks.

Other configurations that are decomposed in $(x, y)$ but have too few Nz, or any configuration decomposed in $z$, are not supported.

Algorithm for two-dimensional decompositions

For two-dimensional decompositions (useful for three-dimensional problems), there are three forward transforms, three backward transforms, and four transpositions requiring MPI communication. In the schematic below, the first dimension is always the local dimension. In our implementation of the PencilFFTs algorithm, we require either Nz >= Rx, or Nz >= Ry, where Nz is the number of vertical cells, Rx is the number of ranks in $x$, and Ry is the number of ranks in $y$. Below, we outline the algorithm for the case Nz >= Rx. If Nz < Rx, but Nz > Ry, a similar algorithm applies with $x$ and $y$ swapped:

  1. first(storage) is initialized with layout $(z, x, y)$.
  2. Transform along $z$.

3 Transpose + communicate to storage[2] in layout $(x, z, y)$, which is distributed into (Rx, Ry) processes in $(z, y)$.

  1. Transform along $x$.
  2. Transpose + communicate to last(storage) in layout $(y, x, z)$, which is distributed into (Rx, Ry) processes in $(x, z)$.
  3. Transform in $y$.

At this point the three in-place forward transforms are complete, and we solve the Poisson equation by updating last(storage). Then the process is reversed to obtain first(storage) in physical space and with the layout $(z, x, y)$.

Restrictions

The algorithm for two-dimensional decompositions requires that Nz = size(global_grid, 3) is larger than either Rx = ranks[1] or Ry = ranks[2], where ranks are configured when building Distributed. If Nz does not satisfy this condition, we can only support a one-dimensional decomposition.

Algorithm for one-dimensional decompositions

This algorithm requires a one-dimensional decomposition with either Rx = 1 or Ry = 1, and is important to support two-dimensional transforms.

For one-dimensional decompositions, we place the decomposed direction last. If the number of ranks is Rh = max(Rx, Ry), this algorithm requires that both Nx > Rh and Ny > Rh. The resulting flow of transposes and transforms is similar to the two-dimensional case. It remains somewhat of a mystery why this succeeds (i.e., why the last transform is correctly decomposed).

source
Oceananigans.DistributedComputations.EqualType
Equal()

Return a type that partitions a direction equally among remaining processes.

Equal() can be used for only one direction. Other directions must either be unspecified, or specifically defined by Int, Fractional, or Sizes.

source
Oceananigans.DistributedComputations.PartitionMethod
Partition(; x = 1, y = 1, z = 1)

Return Partition representing the division of a domain in the x (first), y (second) and z (third) dimension

Keyword arguments:

  • x: partitioning of the first dimension
  • y: partitioning of the second dimension
  • z: partitioning of the third dimension

if supplied as positional arguments x will be the first argument, y the second and z the third

x, y and z can be:

  • x::Int: allocate x processors to the first dimension
  • Equal(): divide the domain in x equally among the remaining processes (not supported for multiple directions)
  • Fractional(ϵ₁, ϵ₂, ..., ϵₙ): divide the domain unequally among N processes. The total work is W = sum(ϵᵢ), and each process is then allocated ϵᵢ / W of the domain.
  • Sizes(n₁, n₂, ..., nₙ): divide the domain unequally. The total work is W = sum(nᵢ), and each process is then allocated nᵢ.

Examples:

julia> using Oceananigans; using Oceananigans.DistributedComputations

julia> Partition(1, 4)
Domain partitioning with (1, 4, 1) ranks
└── y-partitioning: 4

julia> Partition(x = Fractional(1, 2, 3, 4))
Domain partitioning with (4, 1, 1) ranks
└── x-partitioning: domain fractions: (0.1, 0.2, 0.3, 0.4)
source

Fields

Oceananigans.Fields.AbstractFieldType
AbstractField{LX, LY, LZ, G, T, N}

Abstract supertype for fields located at (LX, LY, LZ) and defined on a grid G with eltype T and N dimensions.

Note: we need the parameter T to subtype AbstractArray.

source
Oceananigans.Fields.BackgroundFieldMethod
BackgroundField(func; parameters=nothing)

Returns a BackgroundField to be passed to NonhydrostaticModel for use as a background velocity or tracer field.

If parameters is not provided, func must be callable with the signature

func(x, y, z, t)

If parameters is provided, func must be callable with the signature

func(x, y, z, t, parameters)
source
Oceananigans.Fields.FieldMethod
Field{LX, LY, LZ}(grid::AbstractGrid,
                  T::DataType=eltype(grid); kw...) where {LX, LY, LZ}

Construct a Field on grid with data type T at the location (LX, LY, LZ). Each of (LX, LY, LZ) is either Center or Face and determines the field's location in (x, y, z) respectively.

Keyword arguments

  • data :: OffsetArray: An offset array with the fields data. If nothing is provided the field is filled with zeros.
  • boundary_conditions: If nothing is provided, then field is created using the default boundary conditions via FieldBoundaryConditions.
  • indices: Used to prescribe where a reduced field lives on. For example, at which k index does a two-dimensional $x$-$y$ field lives on. Default: (:, :, :).

Example

A field at location (Face, Face, Center).

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(2, 3, 4), extent=(1, 1, 1));

julia> ω = Field{Face, Face, Center}(grid)
2×3×4 Field{Face, Face, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 8×9×10 OffsetArray(::Array{Float64, 3}, -2:5, -2:6, -2:7) with eltype Float64 with indices -2:5×-2:6×-2:7
    └── max=0.0, min=0.0, mean=0.0

Now, using indices we can create a two dimensional $x$-$y$ field at location (Face, Face, Center) to compute, e.g., the vertical vorticity $∂v/∂x - ∂u/∂y$ at the fluid's surface $z = 0$, which for Center corresponds to k = Nz.

julia> u = XFaceField(grid); v = YFaceField(grid);

julia> ωₛ = Field(∂x(v) - ∂y(u), indices=(:, :, grid.Nz))
2×3×1 Field{Face, Face, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 4:4)
├── operand: BinaryOperation at (Face, Face, Center)
├── status: time=0.0
└── data: 8×9×1 OffsetArray(::Array{Float64, 3}, -2:5, -2:6, 4:4) with eltype Float64 with indices -2:5×-2:6×4:4
    └── max=0.0, min=0.0, mean=0.0

julia> compute!(ωₛ)
2×3×1 Field{Face, Face, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 4:4)
├── operand: BinaryOperation at (Face, Face, Center)
├── status: time=0.0
└── data: 8×9×1 OffsetArray(::Array{Float64, 3}, -2:5, -2:6, 4:4) with eltype Float64 with indices -2:5×-2:6×4:4
    └── max=0.0, min=0.0, mean=0.0
source
Oceananigans.Fields.ReductionMethod
Reduction(reduce!, operand; dims)

Return a Reduction of operand with reduce!, along dims. Note that Reduction expects reduce! to operate in-place.

Example

using Oceananigans

Nx, Ny, Nz = 3, 3, 3

grid = RectilinearGrid(size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z=(0, 1),
                       topology=(Periodic, Periodic, Periodic))

c = CenterField(grid)

set!(c, (x, y, z) -> x + y + z)

max_c² = Field(Reduction(maximum!, c^2, dims=3))

compute!(max_c²)

max_c²[1:Nx, 1:Ny]

# output
3×3 Matrix{Float64}:
 1.36111  2.25     3.36111
 2.25     3.36111  4.69444
 3.36111  4.69444  6.25
source
Oceananigans.Fields.CenterFieldFunction
CenterField(grid, T=eltype(grid); kw...)

Return a Field{Center, Center, Center} on grid. Additional keyword arguments are passed to the Field constructor.

source
Oceananigans.Fields.PressureFieldsFunction
PressureFields(grid, bcs::NamedTuple)

Return a NamedTuple with pressure fields pHY′ and pNHS initialized as CenterFields on grid. Boundary conditions bcs may be specified via a named tuple of FieldBoundaryConditions.

source
Oceananigans.Fields.PressureFieldsMethod
PressureFields(proposed_pressures::NamedTuple{(:pHY′, :pNHS)}, grid, bcs)

Return a NamedTuple of pressure fields with, overwriting boundary conditions in proposed_tracer_fields with corresponding fields in the NamedTuple bcs.

source
Oceananigans.Fields.TendencyFieldsMethod
TendencyFields(grid, tracer_names;
               u = XFaceField(grid),
               v = YFaceField(grid),
               w = ZFaceField(grid),
               kwargs...)

Return a NamedTuple with tendencies for all solution fields (velocity fields and tracer fields), initialized on grid. Optional kwargs can be specified to assign data arrays to each tendency field.

source
Oceananigans.Fields.TracerFieldsMethod
TracerFields(tracer_names, grid, user_bcs)

Return a NamedTuple with tracer fields specified by tracer_names initialized as CenterFields on grid. Boundary conditions user_bcs may be specified via a named tuple of FieldBoundaryConditions.

source
Oceananigans.Fields.TracerFieldsMethod
TracerFields(tracer_names, grid; kwargs...)

Return a NamedTuple with tracer fields specified by tracer_names initialized as CenterFields on grid. Fields may be passed via optional keyword arguments kwargs for each field.

This function is used by OutputWriters.Checkpointer and TendencyFields. ```

source
Oceananigans.Fields.TracerFieldsMethod
TracerFields(proposed_tracers::NamedTuple, grid, bcs)

Return a NamedTuple of tracers, overwriting boundary conditions in proposed_tracers with corresponding fields in the NamedTuple bcs.

source
Oceananigans.Fields.VelocityFieldsFunction
VelocityFields(grid, user_bcs = NamedTuple())

Return a NamedTuple with fields u, v, w initialized on grid. Boundary conditions bcs may be specified via a named tuple of FieldBoundaryConditions.

source
Oceananigans.Fields.VelocityFieldsMethod
VelocityFields(proposed_velocities::NamedTuple{(:u, :v, :w)}, grid, bcs)

Return a NamedTuple of velocity fields, overwriting boundary conditions in proposed_velocities with corresponding fields in the NamedTuple bcs.

source
Oceananigans.Fields.XFaceFieldFunction
XFaceField(grid, T=eltype(grid); kw...)

Return a Field{Face, Center, Center} on grid. Additional keyword arguments are passed to the Field constructor.

source
Oceananigans.Fields.YFaceFieldFunction
YFaceField(grid, T=eltype(grid); kw...)

Return a Field{Center, Face, Center} on grid. Additional keyword arguments are passed to the Field constructor.

source
Oceananigans.Fields.ZFaceFieldFunction
ZFaceField(grid, T=eltype(grid); kw...)

Return a Field{Center, Center, Face} on grid. Additional keyword arguments are passed to the Field constructor.

source
Oceananigans.Fields.interpolateMethod
interpolate(at_node, from_field, from_loc, from_grid)

Interpolate from_field, at_node, on from_grid and at from_location, where at_node is a tuple of coordinates and and from_loc = (ℓx, ℓy, ℓz).

Note that this is a lower-level interpolate method defined for use in CPU/GPU kernels.

source
Oceananigans.Fields.regrid!Method
regrid!(a, b)

Regrid field b onto the grid of field a.

Example

Generate a tracer field on a vertically stretched grid and regrid it on a regular grid.

using Oceananigans

Nz, Lz = 2, 1.0
topology = (Flat, Flat, Bounded)

input_grid = RectilinearGrid(size=Nz, z = [0, Lz/3, Lz], topology=topology, halo=1)
input_field = CenterField(input_grid)
input_field[1, 1, 1:Nz] = [2, 3]

output_grid = RectilinearGrid(size=Nz, z=(0, Lz), topology=topology, halo=1)
output_field = CenterField(output_grid)

regrid!(output_field, input_field)

output_field[1, 1, :]

# output
4-element OffsetArray(::Vector{Float64}, 0:3) with eltype Float64 with indices 0:3:
 0.0
 2.333333333333333
 3.0
 0.0
source

Forcings

Oceananigans.Forcings.AdvectiveForcingMethod
AdvectiveForcing(u=ZeroField(), v=ZeroField(), w=ZeroField())

Build a forcing term representing advection by the velocity field u, v, w with an advection scheme.

Example

Using a tracer field to model sinking particles

using Oceananigans

# Physical parameters
gravitational_acceleration          = 9.81     # m s⁻²
ocean_density                       = 1026     # kg m⁻³
mean_particle_density               = 2000     # kg m⁻³
mean_particle_radius                = 1e-3     # m
ocean_molecular_kinematic_viscosity = 1.05e-6  # m² s⁻¹

# Terminal velocity of a sphere in viscous flow
Δb = gravitational_acceleration * (mean_particle_density - ocean_density) / ocean_density
ν = ocean_molecular_kinematic_viscosity
R = mean_particle_radius

w_Stokes = - 2/9 * Δb / ν * R^2 # m s⁻¹

settling = AdvectiveForcing(w=w_Stokes)

# output
AdvectiveForcing:
├── u: ZeroField{Int64}
├── v: ZeroField{Int64}
└── w: ConstantField(-1.97096)
source
Oceananigans.Forcings.ContinuousForcingType
ContinuousForcing{LX, LY, LZ, P, F, D, I, ℑ}

A callable object that implements a "continuous form" forcing function on a field at the location LX, LY, LZ with optional parameters.

source
Oceananigans.Forcings.ContinuousForcingMethod
ContinuousForcing(func; parameters=nothing, field_dependencies=())

Construct a "continuous form" forcing with optional parameters and optional field_dependencies on other fields in a model.

If neither parameters nor field_dependencies are provided, then func must be callable with the signature

func(X..., t)

where, on a three-dimensional grid with no Flat directions, X = (x, y, z) is a 3-tuple containing the east-west, north-south, and vertical spatial coordinates, and t is time.

Dimensions with Flat topology are omitted from the coordinate tuple X. For example, on a grid with topology (Periodic, Periodic, Flat), and with no parameters or field_dependencies, then func must be callable

func(x, y, t)

where x and y are the east-west and north-south coordinates, respectively. For another example, on a grid with topology (Flat, Flat, Bounded) (e.g. a single column), and for a forcing with no parameters or field_dependencies, then func must be callable with

func(z, t)

where z is the vertical coordinate.

If field_dependencies are provided, the signature of func must include them. For example, if field_dependencies=(:u, :S) (and parameters are not provided), and on a three-dimensional grid with no Flat dimensions, then func must be callable with the signature

func(x, y, z, t, u, S)

where u is assumed to be the u-velocity component, and S is a tracer. Note that any field which does not have the name u, v, or w is assumed to be a tracer and must be present in model.tracers.

If parameters are provided, then the last argument to func must be parameters. For example, if func has no field_dependencies but does depend on parameters, then on a three-dimensional grid it must be callable with the signature

func(x, y, z, t, parameters)

With field_dependencies=(:u, :v, :w, :c) and parameters and on a three-dimensional grid, then func must be callable with the signature

func(x, y, z, t, u, v, w, c, parameters)
source
Oceananigans.Forcings.DiscreteForcingMethod
DiscreteForcing(func; parameters=nothing)

Construct a "discrete form" forcing function with optional parameters. The forcing function is applied at grid point i, j, k.

When parameters are not specified, func must be callable with the signature

func(i, j, k, grid, clock, model_fields)

where grid is model.grid, clock.time is the current simulation time and clock.iteration is the current model iteration, and model_fields is a NamedTuple with u, v, w and the fields in model.tracers.

Note that the index end does not access the final physical grid point of a model field in any direction. The final grid point must be explicitly specified, as in model_fields.u[i, j, grid.Nz].

When parameters is specified, func must be callable with the signature.

func(i, j, k, grid, clock, model_fields, parameters)

Above, parameters is, in principle, arbitrary. Note, however, that GPU compilation can place constraints on typeof(parameters).

source
Oceananigans.Forcings.GaussianMaskType
GaussianMask{D}(center, width)

Callable object that returns a Gaussian masking function centered on center, with width, and varying along direction D, i.e.,

exp(-(D - center)^2 / (2 * width^2))

Examples

  • Create a Gaussian mask centered on z=0 with width 1 meter.
julia> mask = GaussianMask{:z}(center=0, width=1)
source
Oceananigans.Forcings.LinearTargetType
LinearTarget{D}(intercept, gradient)

Callable object that returns a Linear target function with intercept and gradient, and varying along direction D, i.e.,

intercept + D * gradient

Examples

  • Create a linear target function varying in z, equal to 0 at z=0 and with gradient 10⁻⁶:

    julia> target = LinearTarget{:z}(intercept=0, gradient=1e-6)
source
Oceananigans.Forcings.RelaxationMethod
Relaxation(; rate, mask=onefunction, target=zerofunction)

Returns a Forcing that restores a field to target(X..., t) at the specified rate, in the region mask(X...).

The functions onefunction and zerofunction always return 1 and 0, respectively. Thus the default mask leaves the whole domain uncovered, and the default target is zero.

Example

  • Restore a field to zero on a timescale of "3600" (equal to one hour if the time units of the simulation are seconds).
using Oceananigans

damping = Relaxation(rate = 1/3600)

# output
Relaxation{Float64, typeof(Oceananigans.Forcings.onefunction), typeof(Oceananigans.Forcings.zerofunction)}
├── rate: 0.0002777777777777778
├── mask: 1
└── target: 0
  • Restore a field to a linear z-gradient within the bottom 1/4 of a domain on a timescale of "60" (equal to one minute if the time units of the simulation are seconds).
dTdz = 0.001 # ⁰C m⁻¹, temperature gradient

T₀ = 20 # ⁰C, surface temperature at z=0

Lz = 100 # m, depth of domain

bottom_sponge_layer = Relaxation(; rate = 1/60,
                                   target = LinearTarget{:z}(intercept=T₀, gradient=dTdz),
                                   mask = GaussianMask{:z}(center=-Lz, width=Lz/4))

# output
Relaxation{Float64, GaussianMask{:z, Float64}, LinearTarget{:z, Float64}}
├── rate: 0.016666666666666666
├── mask: exp(-(z + 100.0)^2 / (2 * 25.0^2))
└── target: 20.0 + 0.001 * z
source
Oceananigans.Forcings.ForcingMethod
Forcing(func; parameters=nothing, field_dependencies=(), discrete_form=false)

Returns a forcing function added to the tendency of an Oceananigans model field.

If discrete_form=false (the default), and neither parameters nor field_dependencies are provided, then func must be callable with the signature

func(x, y, z, t)

where x, y, z are the east-west, north-south, and vertical spatial coordinates, and t is time. Note that this form is also default in the constructor for NonhydrostaticModel, so that Forcing is not needed.

If discrete_form=false (the default), and field_dependencies are provided, the signature of func must include them. For example, if field_dependencies=(:u, :S) (and parameters are not provided), then func must be callable with the signature

func(x, y, z, t, u, S)`

where u is assumed to be the u-velocity component, and S is a tracer. Note that any field which does not have the name u, v, or w is assumed to be a tracer and must be present in model.tracers.

If discrete_form=false (the default) and parameters are provided, then the last argument to func must be parameters. For example, if func has no field_dependencies but does depend on parameters, then it must be callable with the signature

func(x, y, z, t, parameters)

The object parameters is arbitrary in principle, however GPU compilation can place constraints on typeof(parameters).

With field_dependencies=(:u, :v, :w, :c) and parameters, then func must be callable with the signature

func(x, y, z, t, u, v, w, c, parameters)

If discrete_form=true then func must be callable with the "discrete form"

func(i, j, k, grid, clock, model_fields)

where i, j, k is the grid point at which the forcing is applied, grid is model.grid, clock.time is the current simulation time and clock.iteration is the current model iteration, and model_fields is a NamedTuple with u, v, w, the fields in model.tracers, and the fields in model.diffusivity_fields, each of which is an OffsetArrays (or NamedTuples of OffsetArrays depending on the turbulence closure) of field data.

When discrete_form=true and parameters is specified, func must be callable with the signature

func(i, j, k, grid, clock, model_fields, parameters)

Examples

using Oceananigans

# Parameterized forcing
parameterized_func(x, y, z, t, p) = p.μ * exp(z / p.λ) * cos(p.ω * t)

v_forcing = Forcing(parameterized_func, parameters = (μ=42, λ=0.1, ω=π))

# output
ContinuousForcing{NamedTuple{(:μ, :λ, :ω), Tuple{Int64, Float64, Irrational{:π}}}}
├── func: parameterized_func (generic function with 1 method)
├── parameters: (μ = 42, λ = 0.1, ω = π)
└── field dependencies: ()

Note that because forcing locations are regularized within the NonhydrostaticModel constructor:

grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
model = NonhydrostaticModel(grid=grid, forcing=(v=v_forcing,))

model.forcing.v

# output
ContinuousForcing{NamedTuple{(:μ, :λ, :ω), Tuple{Int64, Float64, Irrational{:π}}}} at (Center, Face, Center)
├── func: parameterized_func (generic function with 1 method)
├── parameters: (μ = 42, λ = 0.1, ω = π)
└── field dependencies: ()

After passing through the constructor for NonhydrostaticModel, the v-forcing location information is available and set to Center, Face, Center.

# Field-dependent forcing
growth_in_sunlight(x, y, z, t, P) = exp(z) * P

plankton_forcing = Forcing(growth_in_sunlight, field_dependencies=:P)

# output
ContinuousForcing{Nothing}
├── func: growth_in_sunlight (generic function with 1 method)
├── parameters: nothing
└── field dependencies: (:P,)
# Parameterized, field-dependent forcing
tracer_relaxation(x, y, z, t, c, p) = p.μ * exp((z + p.H) / p.λ) * (p.dCdz * z - c) 

c_forcing = Forcing(tracer_relaxation,
                    field_dependencies = :c,
                            parameters = (μ=1/60, λ=10, H=1000, dCdz=1))

# output
ContinuousForcing{NamedTuple{(:μ, :λ, :H, :dCdz), Tuple{Float64, Int64, Int64, Int64}}}
├── func: tracer_relaxation (generic function with 1 method)
├── parameters: (μ = 0.016666666666666666, λ = 10, H = 1000, dCdz = 1)
└── field dependencies: (:c,)
# Unparameterized discrete-form forcing function
filtered_relaxation(i, j, k, grid, clock, model_fields) =
    @inbounds - (model_fields.c[i-1, j, k] + model_fields.c[i, j, k] + model_fields.c[i+1, j, k]) / 3

filtered_forcing = Forcing(filtered_relaxation, discrete_form=true)

# output
DiscreteForcing{Nothing}
├── func: filtered_relaxation (generic function with 1 method)
└── parameters: nothing
# Discrete-form forcing function with parameters
masked_damping(i, j, k, grid, clock, model_fields, parameters) = 
    @inbounds - parameters.μ * exp(grid.zᵃᵃᶜ[k] / parameters.λ) * model_fields.u[i, j, k]

masked_damping_forcing = Forcing(masked_damping, parameters=(μ=42, λ=π), discrete_form=true)

# output
DiscreteForcing{NamedTuple{(:μ, :λ), Tuple{Int64, Irrational{:π}}}}
├── func: masked_damping (generic function with 1 method)
└── parameters: (μ = 42, λ = π)
source

Grids

Oceananigans.Grids.FlatType
Flat

Grid topology for flat dimensions, generally with one grid point, along which the solution is uniform and does not vary.

source
Oceananigans.Grids.LatitudeLongitudeGridType
LatitudeLongitudeGrid([architecture = CPU(), FT = Float64];
                      size,
                      longitude,
                      latitude,
                      z = nothing,
                      radius = R_Earth,
                      topology = nothing,
                      precompute_metrics = true,
                      halo = nothing)

Creates a LatitudeLongitudeGrid with coordinates (λ, φ, z) denoting longitude, latitude, and vertical coordinate respectively.

Positional arguments

  • architecture: Specifies whether arrays of coordinates and spacings are stored on the CPU or GPU. Default: CPU().

  • FT : Floating point data type. Default: Float64.

Keyword arguments

  • size (required): A 3-tuple prescribing the number of grid points each direction.

  • longitude (required), latitude (required), z (default: nothing): Each is either a:

    1. 2-tuple that specify the end points of the domain,
    2. one-dimensional array specifying the cell interface locations, or
    3. a single-argument function that takes an index and returns cell interface location.

    Note: the latitude and longitude coordinates extents are expected in degrees.

  • radius: The radius of the sphere the grid lives on. By default is equal to the radius of Earth.

  • topology: Tuple of topologies (Flat, Bounded, Periodic) for each direction. The vertical topology[3] must be Bounded, while the latitude-longitude topologies can be Bounded, Periodic, or Flat. If no topology is provided then, by default, the topology is (Periodic, Bounded, Bounded) if the latitudinal extent is 360 degrees or (Bounded, Bounded, Bounded) otherwise.

  • precompute_metrics: Boolean specifying whether to precompute horizontal spacings and areas. Default: true. When false, horizontal spacings and areas are computed on-the-fly during a simulation.

  • halo: A 3-tuple of integers specifying the size of the halo region of cells surrounding the physical interior. The default is 3 halo cells in every direction.

Examples

  • A default grid with Float64 type:
julia> using Oceananigans

julia> grid = LatitudeLongitudeGrid(size=(36, 34, 25),
                                    longitude = (-180, 180),
                                    latitude = (-85, 85),
                                    z = (-1000, 0))
36×34×25 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── longitude: Periodic λ ∈ [-180.0, 180.0) regularly spaced with Δλ=10.0
├── latitude:  Bounded  φ ∈ [-85.0, 85.0]   regularly spaced with Δφ=5.0
└── z:         Bounded  z ∈ [-1000.0, 0.0]  regularly spaced with Δz=40.0
  • A bounded spherical sector with cell interfaces stretched hyperbolically near the top:
julia> using Oceananigans

julia> σ = 1.1; # stretching factor

julia> Nz = 24; # vertical resolution

julia> Lz = 1000; # depth (m)

julia> hyperbolically_spaced_faces(k) = - Lz * (1 - tanh(σ * (k - 1) / Nz) / tanh(σ));

julia> grid = LatitudeLongitudeGrid(size=(36, 34, Nz),
                                    longitude = (-180, 180),
                                    latitude = (-20, 20),
                                    z = hyperbolically_spaced_faces,
                                    topology = (Bounded, Bounded, Bounded))
36×34×24 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── longitude: Bounded  λ ∈ [-180.0, 180.0] regularly spaced with Δλ=10.0
├── latitude:  Bounded  φ ∈ [-20.0, 20.0]   regularly spaced with Δφ=1.17647
└── z:         Bounded  z ∈ [-1000.0, -0.0] variably spaced with min(Δz)=21.3342, max(Δz)=57.2159
source
Oceananigans.Grids.RectilinearGridType
RectilinearGrid([architecture = CPU(), FT = Float64];
                size,
                x = nothing,
                y = nothing,
                z = nothing,
                halo = nothing,
                extent = nothing,
                topology = (Periodic, Periodic, Bounded))

Create a RectilinearGrid with size = (Nx, Ny, Nz) grid points.

Positional arguments

  • architecture: Specifies whether arrays of coordinates and spacings are stored on the CPU or GPU. Default: CPU().

  • FT: Floating point data type. Default: Float64.

Keyword arguments

  • size (required): A tuple prescribing the number of grid points in non-Flat directions. size is a 3-tuple for 3D models, a 2-tuple for 2D models, and either a scalar or 1-tuple for 1D models.

  • topology: A 3-tuple (TX, TY, TZ) specifying the topology of the domain. TX, TY, and TZ specify whether the x-, y-, and z directions are Periodic, Bounded, or Flat. The topology Flat indicates that a model does not vary in those directions so that derivatives and interpolation are zero. The default is topology = (Periodic, Periodic, Bounded).

  • extent: A tuple prescribing the physical extent of the grid in non-Flat directions, e.g., (Lx, Ly, Lz). All directions are constructed with regular grid spacing and the domain (in the case that no direction is Flat) is $0 ≤ x ≤ L_x$, $0 ≤ y ≤ L_y$, and $-L_z ≤ z ≤ 0$, which is most appropriate for oceanic applications in which $z = 0$ usually is the ocean's surface.

  • x, y, and z: Each of x, y, z are either (i) 2-tuples that specify the end points of the domain in their respect directions (in which case scalar values may be used in Flat directions), or (ii) arrays or functions of the corresponding indices i, j, or k that specify the locations of cell faces in the x-, y-, or z-direction, respectively. For example, to prescribe the cell faces in z we need to provide a function that takes k as argument and returns the location of the faces for indices k = 1 through k = Nz + 1, where Nz is the size of the stretched z dimension.

Note: Either extent, or all of x, y, and z must be specified.

  • halo: A tuple of integers that specifies the size of the halo region of cells surrounding the physical interior for each non-Flat direction. The default is 3 halo cells in every direction.

The physical extent of the domain can be specified either via x, y, and z keyword arguments indicating the left and right endpoints of each dimensions, e.g., x = (-π, π) or via the extent argument, e.g., extent = (Lx, Ly, Lz), which specifies the extent of each dimension in which case $0 ≤ x ≤ L_x$, $0 ≤ y ≤ L_y$, and $-L_z ≤ z ≤ 0$.

A grid topology may be specified via a tuple assigning one of Periodic, Bounded, and, Flat to each dimension. By default, a horizontally periodic grid topology (Periodic, Periodic, Bounded) is assumed.

Constants are stored using floating point values of type FT. By default this is Float64. Make sure to specify the desired FT if not using Float64.

Grid properties

  • (Nx, Ny, Nz) :: Int: Number of physical points in the $(x, y, z)$-direction.

  • (Hx, Hy, Hz) :: Int: Number of halo points in the $(x, y, z)$-direction.

  • (Lx, Ly, Lz) :: FT: Physical extent of the grid in the $(x, y, z)$-direction.

  • (Δxᶜᵃᵃ, Δyᵃᶜᵃ, Δzᵃᵃᶜ): Spacings in the $(x, y, z)$-directions between the cell faces. These are the lengths in $x$, $y$, and $z$ of Center cells and are defined at Center locations.

  • (Δxᶠᵃᵃ, Δyᵃᶠᵃ, Δzᵃᵃᶠ): Spacings in the $(x, y, z)$-directions between the cell centers. These are the lengths in $x$, $y$, and $z$ of Face cells and are defined at Face locations.

  • (xᶜᵃᵃ, yᵃᶜᵃ, zᵃᵃᶜ): $(x, y, z)$ coordinates of cell Centers.

  • (xᶠᵃᵃ, yᵃᶠᵃ, zᵃᵃᶠ): $(x, y, z)$ coordinates of cell Faces.

Examples

  • A grid with the default Float64 type:
julia> using Oceananigans

julia> grid = RectilinearGrid(size=(32, 32, 32), extent=(1, 2, 3))
32×32×32 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0)  regularly spaced with Δx=0.03125
├── Periodic y ∈ [0.0, 2.0)  regularly spaced with Δy=0.0625
└── Bounded  z ∈ [-3.0, 0.0] regularly spaced with Δz=0.09375
  • A grid with Float32 type:
julia> using Oceananigans

julia> grid = RectilinearGrid(Float32; size=(32, 32, 16), x=(0, 8), y=(-10, 10), z=(-π, π))
32×32×16 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 8.0)          regularly spaced with Δx=0.25
├── Periodic y ∈ [-10.0, 10.0)       regularly spaced with Δy=0.625
└── Bounded  z ∈ [-3.14159, 3.14159] regularly spaced with Δz=0.392699
  • A two-dimenisional, horizontally-periodic grid:
julia> using Oceananigans

julia> grid = RectilinearGrid(size=(32, 32), extent=(2π, 4π), topology=(Periodic, Periodic, Flat))
32×32×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── Periodic x ∈ [3.60072e-17, 6.28319) regularly spaced with Δx=0.19635
├── Periodic y ∈ [7.20145e-17, 12.5664) regularly spaced with Δy=0.392699
└── Flat z
  • A one-dimensional "column" grid:
julia> using Oceananigans

julia> grid = RectilinearGrid(size=256, z=(-128, 0), topology=(Flat, Flat, Bounded))
1×1×256 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── Flat x
├── Flat y
└── Bounded  z ∈ [-128.0, 0.0]    regularly spaced with Δz=0.5
  • A horizontally-periodic regular grid with cell interfaces stretched hyperbolically near the top:
julia> using Oceananigans

julia> σ = 1.1; # stretching factor

julia> Nz = 24; # vertical resolution

julia> Lz = 32; # depth (m)

julia> hyperbolically_spaced_faces(k) = - Lz * (1 - tanh(σ * (k - 1) / Nz) / tanh(σ));

julia> grid = RectilinearGrid(size = (32, 32, Nz),
                              x = (0, 64), y = (0, 64),
                              z = hyperbolically_spaced_faces)
32×32×24 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 64.0)   regularly spaced with Δx=2.0
├── Periodic y ∈ [0.0, 64.0)   regularly spaced with Δy=2.0
└── Bounded  z ∈ [-32.0, -0.0] variably spaced with min(Δz)=0.682695, max(Δz)=1.83091
  • A three-dimensional grid with regular spacing in $x$, cell interfaces at Chebyshev nodes in $y$, and cell interfaces hyperbolically stretched in $z$ near the top:
julia> using Oceananigans

julia> Nx, Ny, Nz = 32, 30, 24;

julia> Lx, Ly, Lz = 200, 100, 32; # (m)

julia> chebychev_nodes(j) = - Ly/2 * cos(π * (j - 1) / Ny);

julia> σ = 1.1; # stretching factor

julia> hyperbolically_spaced_faces(k) = - Lz * (1 - tanh(σ * (k - 1) / Nz) / tanh(σ));

julia> grid = RectilinearGrid(size = (Nx, Ny, Nz),
                              topology = (Periodic, Bounded, Bounded),
                              x = (0, Lx),
                              y = chebychev_nodes,
                              z = hyperbolically_spaced_faces)
32×30×24 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 200.0)  regularly spaced with Δx=6.25
├── Bounded  y ∈ [-50.0, 50.0] variably spaced with min(Δy)=0.273905, max(Δy)=5.22642
└── Bounded  z ∈ [-32.0, -0.0] variably spaced with min(Δz)=0.682695, max(Δz)=1.83091
source
Oceananigans.Grids.conformal_cubed_sphere_panelFunction
conformal_cubed_sphere_panel(architecture::AbstractArchitecture = CPU(),
                             FT::DataType = Float64;
                             size,
                             z,
                             topology = (Bounded, Bounded, Bounded),
                             ξ = (-1, 1),
                             η = (-1, 1),
                             radius = R_Earth,
                             halo = (1, 1, 1),
                             rotation = nothing)

Create a OrthogonalSphericalShellGrid that represents a section of a sphere after it has been conformally mapped from the face of a cube. The cube's coordinates are ξ and η (which, by default, both take values in the range $[-1, 1]$.

The mapping from the face of the cube to the sphere is done via the CubedSphere.jl package.

Positional arguments

  • architecture: Specifies whether arrays of coordinates and spacings are stored on the CPU or GPU. Default: CPU().

  • FT : Floating point data type. Default: Float64.

Keyword arguments

  • size (required): A 3-tuple prescribing the number of grid points each direction.

  • z (required): Either a

    1. 2-tuple that specify the end points of the $z$-domain,
    2. one-dimensional array specifying the cell interface locations, or
    3. a single-argument function that takes an index and returns cell interface location.
  • radius: The radius of the sphere the grid lives on. By default this is equal to the radius of Earth.

  • halo: A 3-tuple of integers specifying the size of the halo region of cells surrounding the physical interior. The default is 1 halo cells in every direction.

  • rotation :: Rotation: Rotation of the conformal cubed sphere panel about some axis that passes through the center of the sphere. If nothing is provided (default), then the panel includes the North Pole of the sphere in its center. For example, to construct a grid that includes tha South Pole we can pass either rotation = RotX(π) or rotation = RotY(π).

Examples

  • The default conformal cubed sphere panel grid with Float64 type:
julia> using Oceananigans, Oceananigans.Grids

julia> grid = conformal_cubed_sphere_panel(size=(36, 34, 25), z=(-1000, 0))
36×34×25 OrthogonalSphericalShellGrid{Float64, Bounded, Bounded, Bounded} on CPU with 1×1×1 halo and with precomputed metrics
├── centered at: North Pole, (λ, φ) = (0.0, 90.0)
├── longitude: Bounded  extent 90.0 degrees variably spaced with min(Δλ)=0.616164, max(Δλ)=2.58892
├── latitude:  Bounded  extent 90.0 degrees variably spaced with min(Δφ)=0.664958, max(Δφ)=2.74119
└── z:         Bounded  z ∈ [-1000.0, 0.0]  regularly spaced with Δz=40.0
  • The conformal cubed sphere panel that includes the South Pole with Float32 type:
julia> using Oceananigans, Oceananigans.Grids, Rotations

julia> grid = conformal_cubed_sphere_panel(Float32, size=(36, 34, 25), z=(-1000, 0), rotation=RotY(π))
36×34×25 OrthogonalSphericalShellGrid{Float32, Bounded, Bounded, Bounded} on CPU with 1×1×1 halo and with precomputed metrics
├── centered at: South Pole, (λ, φ) = (0.0, -90.0)
├── longitude: Bounded  extent 90.0 degrees variably spaced with min(Δλ)=0.616167, max(Δλ)=2.58891
├── latitude:  Bounded  extent 90.0 degrees variably spaced with min(Δφ)=0.664956, max(Δφ)=2.7412
└── z:         Bounded  z ∈ [-1000.0, 0.0]  regularly spaced with Δz=40.0
source
Oceananigans.Grids.minimum_xspacingMethod
minimum_xspacing(grid, ℓx, ℓy, ℓz)
minimum_xspacing(grid) = minimum_xspacing(grid, Center(), Center(), Center())

Return the minimum spacing for grid in $x$ direction at location ℓx, ℓy, ℓz.

Examples

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(2, 4, 8), extent=(1, 1, 1));

julia> minimum_xspacing(grid, Center(), Center(), Center())
0.5
source
Oceananigans.Grids.minimum_yspacingMethod
minimum_yspacing(grid, ℓx, ℓy, ℓz)
minimum_yspacing(grid) = minimum_yspacing(grid, Center(), Center(), Center())

Return the minimum spacing for grid in $y$ direction at location ℓx, ℓy, ℓz.

Examples

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(2, 4, 8), extent=(1, 1, 1));

julia> minimum_yspacing(grid, Center(), Center(), Center())
0.25
source
Oceananigans.Grids.minimum_zspacingMethod
minimum_zspacing(grid, ℓx, ℓy, ℓz)
minimum_zspacing(grid) = minimum_zspacing(grid, Center(), Center(), Center())

Return the minimum spacing for grid in $z$ direction at location ℓx, ℓy, ℓz.

Examples

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(2, 4, 8), extent=(1, 1, 1));

julia> minimum_zspacing(grid, Center(), Center(), Center())
0.125
source
Oceananigans.Grids.new_dataFunction
new_data(FT, arch, loc, topo, sz, halo_sz, indices)

Return an OffsetArray of zeros of float type FT on architecture, with indices corresponding to a field on a grid of size(grid) and located at loc.

source
Oceananigans.Grids.nodesMethod
nodes(grid, (ℓx, ℓy, ℓz); reshape=false, with_halos=false)
nodes(grid, ℓx, ℓy, ℓz; reshape=false, with_halos=false)

Return a 3-tuple of views over the interior nodes of the grid's native coordinates at the locations in loc=(ℓx, ℓy, ℓz) in x, y, z.

If reshape=true, the views are reshaped to 3D arrays with non-singleton dimensions 1, 2, 3 for x, y, z, respectively. These reshaped arrays can then be used in broadcast operations with 3D fields or arrays.

For RectilinearGrids the native coordinates are x, y, z; for curvilinear grids, like LatitudeLongitudeGrid or OrthogonalSphericalShellGrid the native coordinates are λ, φ, z.

See xnodes, ynodes, znodes, λnodes, and φnodes.

source
Oceananigans.Grids.offset_dataFunction
offset_data(underlying_data, grid::AbstractGrid, loc, indices=default_indices(length(loc)))

Return an OffsetArray that maps to underlying_data in memory, with offset indices appropriate for the data of a field on a grid of size(grid) and located at loc.

source
Oceananigans.Grids.total_sizeFunction
total_size(grid, loc)

Return the "total" size of a grid at loc. This is a 3-tuple of integers corresponding to the number of grid points along x, y, z.

source
Oceananigans.Grids.xnodesMethod
xnodes(grid, ℓx, ℓy, ℓz, with_halos=false)

Return the positions over the interior nodes on grid in the $x$-direction for the location ℓx, ℓy, ℓz. For Bounded directions, Face nodes include the boundary points.

See znodes for examples.

source
Oceananigans.Grids.xspacingsMethod
xspacings(grid, ℓx, ℓy, ℓz; with_halos=true)

Return the spacings over the interior nodes on grid in the $x$-direction for the location ℓx, ℓy, ℓz. For Bounded directions, Face nodes include the boundary points.

julia> using Oceananigans

julia> grid = LatitudeLongitudeGrid(size=(8, 15, 10), longitude=(-20, 60), latitude=(-10, 50), z=(-100, 0));

julia> xspacings(grid, Center(), Face(), Center())
16-element view(OffsetArray(::Vector{Float64}, -2:18), 1:16) with eltype Float64:
      1.0950562585518518e6
      1.1058578920188267e6
      1.1112718969963323e6
      1.1112718969963323e6
      1.1058578920188267e6
      1.0950562585518518e6
      1.0789196210678827e6
      1.0575265956426917e6
      1.0309814069457315e6
 999413.38046802
 962976.3124613502
 921847.720658409
 876227.979424229
 826339.3435524226
 772424.8654621692
 714747.2110712599
source
Oceananigans.Grids.ynodesMethod
ynodes(grid, ℓx, ℓy, ℓz, with_halos=false)

Return the positions over the interior nodes on grid in the $y$-direction for the location ℓx, ℓy, ℓz. For Bounded directions, Face nodes include the boundary points.

See znodes for examples.

source
Oceananigans.Grids.yspacingsMethod
yspacings(grid, ℓx, ℓy, ℓz; with_halos=true)

Return the spacings over the interior nodes on grid in the $y$-direction for the location ℓx, ℓy, ℓz. For Bounded directions, Face nodes include the boundary points.

julia> using Oceananigans

julia> grid = LatitudeLongitudeGrid(size=(20, 15, 10), longitude=(0, 20), latitude=(-15, 15), z=(-100, 0));

julia> yspacings(grid, Center(), Center(), Center())
222389.85328911748
source
Oceananigans.Grids.znodesMethod
znodes(grid, ℓx, ℓy, ℓz; with_halos=false)

Return the positions over the interior nodes on grid in the $z$-direction for the location ℓx, ℓy, ℓz. For Bounded directions, Face nodes include the boundary points.

julia> using Oceananigans

julia> horz_periodic_grid = RectilinearGrid(size=(3, 3, 3), extent=(2π, 2π, 1), halo=(1, 1, 1),
                                            topology=(Periodic, Periodic, Bounded));

julia> zC = znodes(horz_periodic_grid, Center())
3-element view(OffsetArray(::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, 0:4), 1:3) with eltype Float64:
 -0.8333333333333334
 -0.5
 -0.16666666666666666

julia> zC = znodes(horz_periodic_grid, Center(), Center(), Center())
3-element view(OffsetArray(::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, 0:4), 1:3) with eltype Float64:
 -0.8333333333333334
 -0.5
 -0.16666666666666666

julia> zC = znodes(horz_periodic_grid, Center(), Center(), Center(), with_halos=true)
-1.1666666666666667:0.3333333333333333:0.16666666666666666 with indices 0:4
source
Oceananigans.Grids.zspacingsMethod
zspacings(grid, ℓx, ℓy, ℓz; with_halos=true)

Return the spacings over the interior nodes on grid in the $z$-direction for the location ℓx, ℓy, ℓz. For Bounded directions, Face nodes include the boundary points.

julia> using Oceananigans

julia> grid = LatitudeLongitudeGrid(size=(20, 15, 10), longitude=(0, 20), latitude=(-15, 15), z=(-100, 0));

julia> zspacings(grid, Center(), Center(), Center())
10.0
source
Oceananigans.Grids.λnodesMethod
λnodes(grid::AbstractCurvilinearGrid, ℓx, ℓy, ℓz, with_halos=false)

Return the positions over the interior nodes on a curvilinear grid in the $λ$-direction for the location ℓλ, ℓφ, ℓz. For Bounded directions, Face nodes include the boundary points.

See znodes for examples.

source
Oceananigans.Grids.φnodesMethod
φnodes(grid::AbstractCurvilinearGrid, ℓx, ℓy, ℓz, with_halos=false)

Return the positions over the interior nodes on a curvilinear grid in the $φ$-direction for the location ℓλ, ℓφ, ℓz. For Bounded directions, Face nodes include the boundary points.

See znodes for examples.

source

Immersed boundaries

Logger

Oceananigans.Logger.OceananigansLoggerType
OceananigansLogger(stream::IO=stdout, level=Logging.Info; show_info_source=false)

Based on Logging.SimpleLogger, it tries to log all messages in the following format:

[yyyy/mm/dd HH:MM:SS.sss] log_level message [-@-> source_file:line_number]

where the source of the message between the square brackets is included only if show_info_source=true or if the message is not an info level message.

source

Models

Oceananigans.Models.seawater_densityMethod
seawater_density(model; temperature, salinity, geopotential_height)

Return a KernelFunctionOperation that computes the in-situ density of seawater with (gridded) temperature, salinity, and at geopotential_height. To compute the in-situ density, the 55 term polynomial approximation to the equation of state from Roquet et al. (2015) is used. By default the seawater_density extracts the geopotential height from the model to compute the in-situ density. To compute a potential density at some user chosen reference geopotential height, set geopotential_height to a constant for the density computation,

geopotential_height = 0 # sea-surface height
σ₀ = seawater_density(model; geopotential_height)

Note: seawater_density must be passed a BoussinesqEquationOfState to compute the density. See the relevant documentation for how to set SeawaterBuoyancy using a BoussinesqEquationOfState.

Example

Compute a density Field using the KernelFunctionOperation returned from seawater_density

julia> using Oceananigans, SeawaterPolynomials.TEOS10

julia> using Oceananigans.Models: seawater_density

julia> grid = RectilinearGrid(size=100, z=(-1000, 0), topology=(Flat, Flat, Bounded))
1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── Flat x
├── Flat y
└── Bounded  z ∈ [-1000.0, 0.0]   regularly spaced with Δz=10.0

julia> tracers = (:T, :S)
(:T, :S)

julia> eos = TEOS10EquationOfState()
BoussinesqEquationOfState{Float64}:
    ├── seawater_polynomial: TEOS10SeawaterPolynomial{Float64}
    └── reference_density: 1020.0

julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: BoussinesqEquationOfState{Float64}

julia> model = NonhydrostaticModel(; grid, buoyancy, tracers)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
└── coriolis: Nothing

julia> set!(model, S = 34.7, T = 0.5)

julia> density_operation = seawater_density(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: ρ (generic function with 2 methods)
└── arguments: ("BoussinesqEquationOfState{Float64}", "1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU", "1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU", "KernelFunctionOperation at (Center, Center, Center)")

julia> density_field = Field(density_operation)
1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 1×1×106 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:103) with eltype Float64 with indices 1:1×1:1×-2:103
    └── max=0.0, min=0.0, mean=0.0

julia> compute!(density_field)
1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 1×1×106 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:103) with eltype Float64 with indices 1:1×1:1×-2:103
    └── max=1032.38, min=1027.73, mean=1030.06

Values for temperature, salinity and geopotential_height can be passed to seawater_density to override the defaults that are obtained from the model.

source

Non-hydrostatic models

Oceananigans.Models.NonhydrostaticModels.NonhydrostaticModelMethod
NonhydrostaticModel(;          grid,
                              clock = Clock{eltype(grid)}(0, 0, 1),
                          advection = CenteredSecondOrder(),
                           buoyancy = nothing,
                           coriolis = nothing,
                       stokes_drift = nothing,
                forcing::NamedTuple = NamedTuple(),
                            closure = nothing,
    boundary_conditions::NamedTuple = NamedTuple(),
                            tracers = (),
                        timestepper = :QuasiAdamsBashforth2,
      background_fields::NamedTuple = NamedTuple(),
      particles::ParticlesOrNothing = nothing,

biogeochemistry::AbstractBGCOrNothing = nothing, velocities = nothing, pressures = nothing, diffusivityfields = nothing, pressuresolver = nothing, immersedboundary = nothing, auxiliaryfields = NamedTuple())

Construct a model for a non-hydrostatic, incompressible fluid on grid, using the Boussinesq approximation when buoyancy != nothing. By default, all Bounded directions are rigid and impenetrable.

Keyword arguments

  • grid: (required) The resolution and discrete geometry on which the model is solved. The architecture (CPU/GPU) that the model is solved on is inferred from the architecture of the grid. Note that the grid needs to be regularly spaced in the horizontal dimensions, $x$ and $y$.
  • advection: The scheme that advects velocities and tracers. See Oceananigans.Advection.
  • buoyancy: The buoyancy model. See Oceananigans.BuoyancyModels.
  • coriolis: Parameters for the background rotation rate of the model.
  • stokes_drift: Parameters for Stokes drift fields associated with surface waves. Default: nothing.
  • forcing: NamedTuple of user-defined forcing functions that contribute to solution tendencies.
  • closure: The turbulence closure for model. See Oceananigans.TurbulenceClosures.
  • boundary_conditions: NamedTuple containing field boundary conditions.
  • tracers: A tuple of symbols defining the names of the modeled tracers, or a NamedTuple of preallocated CenterFields.
  • timestepper: A symbol that specifies the time-stepping method. Either :QuasiAdamsBashforth2 or :RungeKutta3.
  • background_fields: NamedTuple with background fields (e.g., background flow). Default: nothing.
  • particles: Lagrangian particles to be advected with the flow. Default: nothing.
  • biogeochemistry: Biogeochemical model for tracers.
  • velocities: The model velocities. Default: nothing.
  • pressures: Hydrostatic and non-hydrostatic pressure fields. Default: nothing.
  • diffusivity_fields: Diffusivity fields. Default: nothing.
  • pressure_solver: Pressure solver to be used in the model. If nothing (default), the model constructor chooses the default based on the grid provide.
  • immersed_boundary: The immersed boundary. Default: nothing.
  • auxiliary_fields: NamedTuple of auxiliary fields. Default: nothing
source

Hydrostatic free-surface models

Oceananigans.Models.HydrostaticFreeSurfaceModels.HydrostaticFreeSurfaceModelMethod
HydrostaticFreeSurfaceModel(; grid,
                                         clock = Clock{eltype(grid)}(0, 0, 1),
                            momentum_advection = CenteredSecondOrder(),
                              tracer_advection = CenteredSecondOrder(),
                                      buoyancy = SeawaterBuoyancy(eltype(grid)),
                                      coriolis = nothing,
                                  free_surface = ImplicitFreeSurface(gravitational_acceleration = g_Earth),
                           forcing::NamedTuple = NamedTuple(),
                                       closure = nothing,
               boundary_conditions::NamedTuple = NamedTuple(),
                                       tracers = (:T, :S),
                 particles::ParticlesOrNothing = nothing,
         biogeochemistry::AbstractBGCOrNothing = nothing,
                                    velocities = nothing,
                                      pressure = nothing,
                            diffusivity_fields = nothing,
                              auxiliary_fields = NamedTuple(),
)

Construct a hydrostatic model with a free surface on grid.

Keyword arguments

  • grid: (required) The resolution and discrete geometry on which model is solved. The architecture (CPU/GPU) that the model is solve is inferred from the architecture of the grid.
  • momentum_advection: The scheme that advects velocities. See Oceananigans.Advection.
  • tracer_advection: The scheme that advects tracers. See Oceananigans.Advection.
  • buoyancy: The buoyancy model. See Oceananigans.BuoyancyModels.
  • coriolis: Parameters for the background rotation rate of the model.
  • forcing: NamedTuple of user-defined forcing functions that contribute to solution tendencies.
  • free_surface: The free surface model.
  • closure: The turbulence closure for model. See Oceananigans.TurbulenceClosures.
  • boundary_conditions: NamedTuple containing field boundary conditions.
  • tracers: A tuple of symbols defining the names of the modeled tracers, or a NamedTuple of preallocated CenterFields.
  • particles: Lagrangian particles to be advected with the flow. Default: nothing.
  • biogeochemistry: Biogeochemical model for tracers.
  • velocities: The model velocities. Default: nothing.
  • pressure: Hydrostatic pressure field. Default: nothing.
  • diffusivity_fields: Diffusivity fields. Default: nothing.
  • auxiliary_fields: NamedTuple of auxiliary fields. Default: nothing.
source
Oceananigans.Models.HydrostaticFreeSurfaceModels.ImplicitFreeSurfaceMethod
ImplicitFreeSurface(; solver_method=:Default, gravitational_acceleration=g_Earth, solver_settings...)

Return an implicit free-surface solver. The implicit free-surface equation is

\[\left [ 𝛁_h ⋅ (H 𝛁_h) - \frac{1}{g Δt^2} \right ] η^{n+1} = \frac{𝛁_h ⋅ 𝐐_⋆}{g Δt} - \frac{η^{n}}{g Δt^2} ,\]

where $η^n$ is the free-surface elevation at the $n$-th time step, $H$ is depth, $g$ is the gravitational acceleration, $Δt$ is the time step, $𝛁_h$ is the horizontal gradient operator, and $𝐐_⋆$ is the barotropic volume flux associated with the predictor velocity field $𝐮_⋆$, i.e.,

\[𝐐_⋆ = \int_{-H}^0 𝐮_⋆ \, 𝖽 z ,\]

where

\[𝐮_⋆ = 𝐮^n + \int_{t_n}^{t_{n+1}} 𝐆ᵤ \, 𝖽t .\]

This equation can be solved, in general, using the PreconditionedConjugateGradientSolver but other solvers can be invoked in special cases.

If $H$ is constant, we divide through out to obtain

\[\left ( ∇^2_h - \frac{1}{g H Δt^2} \right ) η^{n+1} = \frac{1}{g H Δt} \left ( 𝛁_h ⋅ 𝐐_⋆ - \frac{η^{n}}{Δt} \right ) .\]

Thus, for constant $H$ and on grids with regular spacing in $x$ and $y$ directions, the free surface can be obtained using the FFTBasedPoissonSolver.

solver_method can be either of:

By default, if the grid has regular spacing in the horizontal directions then the :FastFourierTransform is chosen, otherwise the :HeptadiagonalIterativeSolver.

source
Oceananigans.Models.HydrostaticFreeSurfaceModels.PrescribedVelocityFieldsMethod
PrescribedVelocityFields(; u = ZeroField(),
                           v = ZeroField(),
                           w = ZeroField(),
                           parameters = nothing)

Builds PrescribedVelocityFields with prescribed functions u, v, and w.

If isnothing(parameters), then u, v, w are called with the signature

u(x, y, z, t) = # something interesting

If !isnothing(parameters), then u, v, w are called with the signature

u(x, y, z, t, parameters) = # something parameterized and interesting

In the constructor for HydrostaticFreeSurfaceModel, the functions u, v, w are wrapped in FunctionField and associated with the model's grid and clock.

source
Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitFreeSurfaceType
SplitExplicitFreeSurface(; gravitational_acceleration = g_Earth, kwargs...)

Return a SplitExplicitFreeSurface representing an explicit time discretization of oceanic free surface dynamics with gravitational_acceleration.

Keyword Arguments

  • substeps: The number of substeps that divide the range (t, t + 2Δt), where Δt is the baroclinic timestep. Note that some averaging functions do not require substepping until 2Δt. The number of substeps is reduced automatically to the last index of averaging_weights for which averaging_weights > 0.

  • cfl: If set then the number of substeps are computed based on the advective timescale imposed from the barotropic gravity-wave speed, computed with depth grid.Lz. If fixed_Δt is provided then the number of substeps will adapt to maintain an exact cfl. If not the effective cfl will be always lower than the specified cfl provided that the baroclinic time step Δt_baroclinic < fixed_Δt

Needed keyword arguments

Either substeps or cfl (with grid) need to be prescribed.

  • grid: Used to compute the corresponding barotropic surface wave speed.

  • fixed_Δt: The maximum baroclinic timestep allowed. If fixed_Δt is a nothing and a cfl is provided, then the number of substeps will be computed on the fly from the baroclinic time step to maintain a constant cfl.

  • gravitational_acceleration: the gravitational acceleration (default: g_Earth)

  • averaging_kernel: function of τ used to average the barotropic transport U and free surface η within the barotropic advancement. τ is the fractional substep going from 0 to 2 with the baroclinic time step t + Δt located at τ = 1. This function should be centered at τ = 1, that is, $∑ (aₘ m /M) = 1$. By default the averaging kernel described by Shchepetkin and McWilliams (2005): https://doi.org/10.1016/j.ocemod.2004.08.002 is chosen.

  • timestepper: Time stepping scheme used for the barotropic advancement. Choose one of:

    • ForwardBackwardScheme() (default): η = f(U) then U = f(η),
    • AdamsBashforth3Scheme(): η = f(U, Uᵐ⁻¹, Uᵐ⁻²) then U = f(η, ηᵐ, ηᵐ⁻¹, ηᵐ⁻²).
source
Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitFreeSurfaceType
struct SplitExplicitFreeSurface

The split-explicit free surface solver.

  • η: The instantaneous free surface (ReducedField)

  • state: The entire state for the split-explicit solver (SplitExplicitState)

  • auxiliary: Parameters for timestepping split-explicit solver (NamedTuple)

  • gravitational_acceleration: Gravitational acceleration

  • settings: Settings for the split-explicit scheme (NamedTuple)

source

Shallow-water models

Oceananigans.Models.ShallowWaterModels.ShallowWaterModelMethod
ShallowWaterModel(; grid,
                    gravitational_acceleration,
                          clock = Clock{eltype(grid)}(0, 0, 1),
             momentum_advection = UpwindBiasedFifthOrder(),
               tracer_advection = WENO(),
                 mass_advection = WENO(),
                       coriolis = nothing,
            forcing::NamedTuple = NamedTuple(),
                        closure = nothing,
                     bathymetry = nothing,
                        tracers = (),
             diffusivity_fields = nothing,
boundary_conditions::NamedTuple = NamedTuple(),
            timestepper::Symbol = :RungeKutta3,
                    formulation = ConservativeFormulation())

Construct a shallow water model on grid with gravitational_acceleration constant.

Keyword arguments

  • grid: (required) The resolution and discrete geometry on which model is solved. The architecture (CPU/GPU) that the model is solve is inferred from the architecture of the grid.
  • gravitational_acceleration: (required) The gravitational acceleration constant.
  • clock: The clock for the model.
  • momentum_advection: The scheme that advects velocities. See Oceananigans.Advection. Default: UpwindBiasedFifthOrder().
  • tracer_advection: The scheme that advects tracers. See Oceananigans.Advection. Default: WENO().
  • mass_advection: The scheme that advects the mass equation. See Oceananigans.Advection. Default: WENO().
  • coriolis: Parameters for the background rotation rate of the model.
  • forcing: NamedTuple of user-defined forcing functions that contribute to solution tendencies.
  • closure: The turbulence closure for model. See Oceananigans.TurbulenceClosures.
  • bathymetry: The bottom bathymetry.
  • tracers: A tuple of symbols defining the names of the modeled tracers, or a NamedTuple of preallocated CenterFields.
  • diffusivity_fields: Stores diffusivity fields when the closures require a diffusivity to be calculated at each timestep.
  • boundary_conditions: NamedTuple containing field boundary conditions.
  • timestepper: A symbol that specifies the time-stepping method. Either :QuasiAdamsBashforth2 or :RungeKutta3 (default).
  • formulation: Whether the dynamics are expressed in conservative form (ConservativeFormulation(); default) or in non-conservative form with a vector-invariant formulation for the non-linear terms (VectorInvariantFormulation()).
Formulation-grid compatibility requirements

The ConservativeFormulation() requires RectilinearGrid. Use VectorInvariantFormulation() with LatitudeLongitudeGrid.

source
Oceananigans.Models.ShallowWaterModels.ShallowWaterScalarDiffusivityType
ShallowWaterScalarDiffusivity([FT::DataType=Float64;]
                              ν=0, ξ=0, discrete_form=false)

Return a scalar diffusivity for the shallow water model.

The diffusivity for the shallow water model is calculated as h * ν so that we get a viscous term in the form $h^{-1} 𝛁 ⋅ (h ν t)$, where $t$ is the 2D stress tensor plus a trace, i.e., $t = 𝛁𝐮 + (𝛁𝐮)^T - ξ I ⋅ (𝛁 ⋅ 𝐮)$.

With the VectorInvariantFormulation() (that evolves $u$ and $v$) we compute $h^{-1} 𝛁(ν h 𝛁 t)$, while with the ConservativeFormulation() (that evolves $u h$ and $v h$) we compute $𝛁 (ν h 𝛁 t)$.

source

Lagrangian particle tracking

Oceananigans.Models.LagrangianParticleTracking.LagrangianParticlesMethod
LagrangianParticles(particles::StructArray; restitution=1.0, tracked_fields::NamedTuple=NamedTuple(), dynamics=no_dynamics)

Construct some LagrangianParticles that can be passed to a model. The particles should be a StructArray and can contain custom fields. The coefficient of restitution for particle-wall collisions is specified by restitution.

A number of tracked_fields may be passed in as a NamedTuple of fields. Each particle will track the value of each field. Each tracked field must have a corresponding particle property. So if T is a tracked field, then T must also be a custom particle property.

dynamics is a function of (lagrangian_particles, model, Δt) that is called prior to advecting particles. parameters can be accessed inside the dynamics function.

source
Oceananigans.Models.LagrangianParticleTracking.LagrangianParticlesMethod
LagrangianParticles(; x, y, z, restitution=1.0, dynamics=no_dynamics, parameters=nothing)

Construct some LagrangianParticles that can be passed to a model. The particles will have initial locations x, y, and z. The coefficient of restitution for particle-wall collisions is specified by restitution.

dynamics is a function of (lagrangian_particles, model, Δt) that is called prior to advecting particles. parameters can be accessed inside the dynamics function.

source

MultiRegion

Oceananigans.MultiRegion.ConformalCubedSphereGridType
ConformalCubedSphereGrid(arch=CPU(), FT=Float64;
                         panel_size,
                         z,
                         horizontal_direction_halo = 1,
                         z_halo = horizontal_direction_halo,
                         horizontal_topology = FullyConnected,
                         z_topology = Bounded,
                         radius = R_Earth,
                         partition = CubedSpherePartition(; R = 1),
                         devices = nothing)

Return a ConformalCubedSphereGrid that comprises of six conformal_cubed_sphere_panel grids; we refer to each of these grids as a "panel". Each panel corresponds to a face of the cube.

The keyword arguments prescribe the properties of each of the panels. Only the topology in the vertical direction can be prescribed and that's done via the z_topology keyword argumet (default: Bounded). Topologies in both horizontal directions for a ConformalCubedSphereGrid are always FullyConnected.

Halo size in both horizontal dimensions must be equal; this is prescribed via the horizontal_halo :: Integer keyword argument. The number of halo points in the $z$-direction is prescribed by the z_halo :: Integer keyword argument.

The connectivity between the ConformalCubedSphereGrid panels is depicted below.

                          +==========+==========+
                          ∥    ↑     ∥    ↑     ∥
                          ∥    1W    ∥    1S    ∥
                          ∥←3N P5 6W→∥←5E P6 2S→∥
                          ∥    4N    ∥    4E    ∥
                          ∥    ↓     ∥    ↓     ∥
               +==========+==========+==========+
               ∥    ↑     ∥    ↑     ∥
               ∥    5W    ∥    5S    ∥
               ∥←1N P3 4W→∥←3E P4 6S→∥
               ∥    2N    ∥    2E    ∥
               ∥    ↓     ∥    ↓     ∥
    +==========+==========+==========+
    ∥    ↑     ∥    ↑     ∥
    ∥    3W    ∥    3S    ∥
    ∥←5N P1 2W→∥←1E P2 4S→∥
    ∥    6N    ∥    6E    ∥
    ∥    ↓     ∥    ↓     ∥
    +==========+==========+

The North Pole of the sphere lies in the center of panel 3 (P3) and the South Pole in the center of panel 6 (P6).

The partition keyword argument prescribes the partitioning in regions within each panel; see CubedSpherePartition. For example, a CubedSpherePartition(; R=2) implies that each of the panels are partitioned into 2 regions in each dimension; this adds up, e.g., to 24 regions for the whole sphere. In the depiction below, the intra-panel x, y indices are depicted in the center of each region and the overall region index is shown at the bottom right of each region.

                                                +==========+==========+==========+==========+
                                                ∥    ↑     |    ↑     ∥    ↑     |    ↑     ∥
                                                ∥          |          ∥          |          ∥
                                                ∥← (1, 2) →|← (2, 2) →∥← (1, 2) →|← (2, 2) →∥
                                                ∥          |          ∥          |          ∥
                                                ∥    ↓  19 |    ↓  20 ∥    ↓  23 |    ↓  24 ∥
                                                +-------- P 5 --------+-------- P 6 --------+
                                                ∥    ↑     |    ↑     ∥    ↑     |    ↑     ∥
                                                ∥          |          ∥          |          ∥
                                                ∥← (1, 1) →|← (2, 1) →∥← (1, 1) →|← (2, 1) →∥
                                                ∥          |          ∥          |          ∥
                                                ∥    ↓  17 |    ↓  18 ∥    ↓  21 |    ↓  22 ∥
                          +==========+==========+==========+==========+==========+==========+
                          ∥    ↑     |    ↑     ∥    ↑     |    ↑     ∥
                          ∥          |          ∥          |          ∥
                          ∥← (1, 2) →|← (2, 2) →∥← (1, 2) →|← (2, 2) →∥
                          ∥          |          ∥          |          ∥
                          ∥    ↓ 11  |    ↓  12 ∥    ↓  15 |    ↓  16 ∥
                          +-------- P 3 --------+-------- P 4 --------+
                          ∥    ↑     |    ↑     ∥    ↑     |    ↑     ∥
                          ∥          |          ∥          |          ∥
                          ∥← (1, 1) →|← (2, 1) →∥← (1, 1) →|← (2, 1) →∥
                          ∥          |          ∥          |          ∥
                          ∥    ↓  9  |    ↓  10 ∥    ↓  13 |    ↓  14 ∥
    +==========+==========+==========+==========+==========+==========+
    ∥    ↑     |    ↑     ∥    ↑     |    ↑     ∥
    ∥          |          ∥          |          ∥
    ∥← (1, 2) →|← (2, 2) →∥← (1, 2) →|← (2, 2) →∥
    ∥          |          ∥          |          ∥
    ∥    ↓   3 |    ↓   4 ∥    ↓   7 |    ↓   8 ∥
    +-------- P 1 --------+-------- P 2 --------+
    ∥    ↑     |    ↑     ∥    ↑     |    ↑     ∥
    ∥          |          ∥          |          ∥
    ∥← (1, 1) →|← (2, 1) →∥← (1, 1) →|← (2, 1) →∥
    ∥          |          ∥          |          ∥
    ∥    ↓   1 |    ↓   2 ∥    ↓   5 |    ↓   6 ∥
    +==========+==========+==========+==========+

Below, we show in detail panels 1 and 2 and the connectivity of each panel.

+===============+==============+==============+===============+
∥       ↑       |      ↑       ∥      ↑       |      ↑        ∥
∥      11W      |      9W      ∥      9S      |     10S       ∥
∥←19N (2, 1) 4W→|←3E (2, 2) 7W→∥←4E (2, 1) 8W→|←7E (2, 2) 13S→∥
∥       1N      |      2N      ∥      5N      |      6N       ∥
∥       ↓     3 |      ↓     4 ∥      ↓     7 |      ↓      8 ∥
+------------- P 1 ------------+------------ P 2 -------------+
∥       ↑       |      ↑       ∥      ↑       |      ↑        ∥
∥       3S      |      4S      ∥      7S      |      8S       ∥
∥←20N (1, 1) 2W→|←1E (2, 1) 5W→∥←2E (1, 1) 6W→|←5E (2, 1) 14S→∥
∥      23N      |     24N      ∥     24N      |     22N       ∥
∥       ↓     1 |      ↓     2 ∥      ↓     5 |      ↓      6 ∥
+===============+==============+==============+===============+

Example

using Oceananigans

grid = ConformalCubedSphereGrid(panel_size=(12, 12, 1), z=(-1, 0), radius=1)

We can find out all connectivities of the regions of our grid. For example, to determine the connectivites on the South boundary of each region we can call

using Oceananigans.MultiRegion: CubedSphereRegionalConnectivity, East, West, South, North, getregion

for region in 1:length(grid); println("panel ", region, ": ", getregion(grid.connectivity.connections, 3).south); end
source
Oceananigans.MultiRegion.ConformalCubedSphereGridType
ConformalCubedSphereGrid(filepath::AbstractString, arch::AbstractArchitecture=CPU(), FT=Float64;
                         Nz,
                         z,
                         panel_halo = (1, 1, 1),
                         panel_topology = (FullyConnected, FullyConnected, Bounded),
                         radius = R_Earth,
                         devices = nothing)

Load a ConformalCubedSphereGrid from filepath.

source
Oceananigans.MultiRegion.MultiRegionGridMethod
MultiRegionGrid(global_grid; partition = XPartition(2),
                             devices = nothing,
                             validate = true)

Split a global_grid into different regions handled by devices.

Positional Arguments

  • global_grid: the grid to be divided into regions.

Keyword Arguments

  • partition: the partitioning required. The implemented partitioning are XPartition (division along the $x$ direction) and YPartition (division along the $y$ direction).

  • devices: the devices to allocate memory on. If nothing is provided (default) then memorey is allocated on the the CPU. For GPU computation it is possible to specify the total number of GPUs or the specific GPUs to allocate memory on. The number of devices does not need to match the number of regions.

  • validate :: Boolean: Whether to validate devices; defautl: true.

Example

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(10, 12), extent=(1, 1), topology=(Bounded, Bounded, Flat))
10×12×1 RectilinearGrid{Float64, Bounded, Bounded, Flat} on CPU with 3×3×0 halo
├── Bounded  x ∈ [0.0, 1.0]       regularly spaced with Δx=0.1
├── Bounded  y ∈ [0.0, 1.0]       regularly spaced with Δy=0.0833333
└── Flat z

julia> multi_region_grid = MultiRegionGrid(grid, partition = XPartition(5))
┌ Warning: MultiRegion functionalities are experimental: help the development by reporting bugs or non-implemented features!
└ @ Oceananigans.MultiRegion ~/Research/OC.jl/src/MultiRegion/multi_region_grid.jl:53
MultiRegionGrid{Float64, Bounded, Bounded, Flat} partitioned on CPU(): 
├── grids: 2×12×1 RectilinearGrid{Float64, RightConnected, Bounded, Flat} on CPU with 3×3×0 halo
├── partitioning: Equal partitioning in X with (5 regions) 
└── devices: (CPU(), CPU(), CPU(), CPU(), CPU())
source

Operators

Oceananigans.Operators.divᶜᶜᶜMethod
divᶜᶜᶜ(i, j, k, grid, u, v, w)

Calculate the divergence $𝛁·𝐕$ of a vector field $𝐕 = (u, v, w)$,

1/V * [δxᶜᵃᵃ(Ax * u) + δxᵃᶜᵃ(Ay * v) + δzᵃᵃᶜ(Az * w)]

which ends up at the cell centers ccc.

source
Oceananigans.Operators.∇²ᶜᶜᶜMethod
∇²ᶜᶜᶜ(i, j, k, grid, c)

Calculate the Laplacian of $c$ via

1/V * [δxᶜᵃᵃ(Ax * ∂xᶠᵃᵃ(c)) + δyᵃᶜᵃ(Ay * ∂yᵃᶠᵃ(c)) + δzᵃᵃᶜ(Az * ∂zᵃᵃᶠ(c))]

which ends up at the location ccc.

source

Output readers

Oceananigans.OutputReaders.FieldDatasetMethod
FieldDataset(filepath;
             architecture=CPU(), grid=nothing, backend=InMemory(), metadata_paths=["metadata"])

Returns a Dict containing a FieldTimeSeries for each field in the JLD2 file located at filepath. Note that model output must have been saved with halos.

Keyword arguments

  • backend: Either InMemory() (default) or OnDisk(). The InMemory backend will

load the data fully in memory as a 4D multi-dimensional array while the OnDisk() backend will lazily load field time snapshots when the FieldTimeSeries is indexed linearly.

  • metadata_paths: A list of JLD2 paths to look for metadata. By default it looks in file["metadata"].

  • grid: May be specified to override the grid used in the JLD2 file.

source
Oceananigans.OutputReaders.FieldTimeSeriesMethod
FieldTimeSeries(path, name;
                backend = InMemory(),
                grid = nothing,
                iterations = nothing,
                times = nothing)

Return a FieldTimeSeries containing a time-series of the field name load from JLD2 output located at path.

Keyword arguments

  • backend: InMemory() to load data into a 4D array, OnDisk() to lazily load data from disk when indexing into FieldTimeSeries.

  • grid: A grid to associate with the data, in the case that the native grid was not serialized properly.

  • iterations: Iterations to load. Defaults to all iterations found in the file.

  • times: Save times to load, as determined through an approximate floating point comparison to recorded save times. Defaults to times associated with iterations. Takes precedence over iterations if times is specified.

source
Oceananigans.OutputReaders.FieldTimeSeriesMethod
FieldTimeSeries{LX, LY, LZ}(grid::AbstractGrid, times; kwargs...) where {LX, LY, LZ} =

Construct a FieldTimeSeries on grid and at times.

Keyword arguments

  • indices: spatial indices
  • backend: backend, InMemory(indices=Colon()) or OnDisk()
  • path: path to data for backend = OnDisk()
  • name: name of field for backend = OnDisk()
source

Output writers

Oceananigans.OutputWriters.AveragedTimeIntervalMethod
AveragedTimeInterval(interval; window=interval, stride=1)

Returns a schedule that specifies periodic time-averaging of output. The time window specifies the extent of the time-average, which reoccurs every interval.

output is computed and accumulated into the average every stride iterations during the averaging window. For example, stride=1 computs output every iteration, whereas stride=2 computes output every other iteration. Time-averages with longer strides are faster to compute, but less accurate.

The time-average of $a$ is a left Riemann sum corresponding to

\[⟨a⟩ = T⁻¹ \int_{tᵢ-T}^{tᵢ} a \mathrm{d} t \, ,\]

where $⟨a⟩$ is the time-average of $a$, $T$ is the time-window for averaging, and the $tᵢ$ are discrete times separated by the time interval. The $tᵢ$ specify both the end of the averaging window and the time at which output is written.

Example

using Oceananigans.OutputWriters: AveragedTimeInterval
using Oceananigans.Utils: days

schedule = AveragedTimeInterval(4days, window=2days)

# output
AveragedTimeInterval(window=2 days, stride=1, interval=4 days)

An AveragedTimeInterval schedule directs an output writer to time-average its outputs before writing them to disk:

using Oceananigans
using Oceananigans.OutputWriters: JLD2OutputWriter
using Oceananigans.Utils: minutes

model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)))

simulation = Simulation(model, Δt=10minutes, stop_time=30days)

simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities,
                                                          filename= "averaged_velocity_data.jld2",
                                                          schedule = AveragedTimeInterval(4days, window=2days, stride=2))

# output
JLD2OutputWriter scheduled on TimeInterval(4 days):
├── filepath: ./averaged_velocity_data.jld2
├── 3 outputs: (u, v, w) averaged on AveragedTimeInterval(window=2 days, stride=2, interval=4 days)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB
source
Oceananigans.OutputWriters.CheckpointerMethod
Checkpointer(model; schedule,
             dir = ".",
             prefix = "checkpoint",
             overwrite_existing = false,
             verbose = false,
             cleanup = false,
             properties = [:architecture, :grid, :clock, :coriolis,
                           :buoyancy, :closure, :timestepper, :particles])

Construct a Checkpointer that checkpoints the model to a JLD2 file on schedule. The model.clock.iteration is included in the filename to distinguish between multiple checkpoint files.

To restart or "pickup" a model from a checkpoint, specify pickup = true when calling run!, ensuring that the checkpoint file is in directory dir. See run! for more details.

Note that extra model properties can be safely specified, but removing crucial properties such as :velocities will make restoring from the checkpoint impossible.

The checkpointer attempts to serialize as much of the model to disk as possible, but functions or objects containing functions cannot be serialized at this time.

Keyword arguments

  • schedule (required): Schedule that determines when to checkpoint.

  • dir: Directory to save output to. Default: "." (current working directory).

  • prefix: Descriptive filename prefixed to all output files. Default: "checkpoint".

  • overwrite_existing: Remove existing files if their filenames conflict. Default: false.

  • verbose: Log what the output writer is doing with statistics on compute/write times and file sizes. Default: false.

  • cleanup: Previous checkpoint files will be deleted once a new checkpoint file is written. Default: false.

  • properties: List of model properties to checkpoint. This list must contain [:grid, :architecture, :timestepper, :particles]. Default: [:architecture, :grid, :clock, :coriolis, :buoyancy, :closure, :velocities, :tracers, :timestepper, :particles]

source
Oceananigans.OutputWriters.JLD2OutputWriterMethod
JLD2OutputWriter(model, outputs; filename, schedule,
                          dir = ".",
                      indices = (:, :, :),
                   with_halos = false,
                   array_type = Array{Float64},
                 max_filesize = Inf,
           overwrite_existing = false,
                         init = noinit,
                    including = [:grid, :coriolis, :buoyancy, :closure],
                      verbose = false,
                         part = 1,
                      jld2_kw = Dict{Symbol, Any}())

Construct a JLD2OutputWriter for an Oceananigans model that writes label, output pairs in outputs to a JLD2 file.

The argument outputs may be a Dict or NamedTuple. The keys of outputs are symbols or strings that "name" output data. The values of outputs are either AbstractFields, objects that are called with the signature output(model), or WindowedTimeAverages of AbstractFieldss, functions, or callable objects.

Keyword arguments

Filenaming

  • filename (required): Descriptive filename. ".jld2" is appended to filename in the file path if filename does not end in ".jld2".

  • dir: Directory to save output to. Default: "." (current working directory).

Output frequency and time-averaging

  • schedule (required): AbstractSchedule that determines when output is saved.

Slicing and type conversion prior to output

  • indices: Specifies the indices to write to disk with a Tuple of Colon, UnitRange, or Int elements. Indices must be Colon, Int, or contiguous UnitRange. Defaults to (:, :, :) or "all indices". If !with_halos, halo regions are removed from indices. For example, indices = (:, :, 1) will save xy-slices of the bottom-most index.

  • with_halos (Bool): Whether or not to slice halo regions from fields before writing output. Note, that to postprocess saved output (e.g., compute derivatives, etc) information about the boundary conditions is often crucial. In that case you might need to set with_halos = true.

  • array_type: The array type to which output arrays are converted to prior to saving. Default: Array{Float64}.

File management

  • max_filesize: The writer will stop writing to the output file once the file size exceeds max_filesize, and write to a new one with a consistent naming scheme ending in part1, part2, etc. Defaults to Inf.

  • overwrite_existing: Remove existing files if their filenames conflict. Default: false.

Output file metadata management

  • init: A function of the form init(file, model) that runs when a JLD2 output file is initialized. Default: noinit(args...) = nothing.

  • including: List of model properties to save with every file. Default: [:grid, :coriolis, :buoyancy, :closure]

Miscellaneous keywords

  • verbose: Log what the output writer is doing with statistics on compute/write times and file sizes. Default: false.

  • part: The starting part number used if max_filesize is finite. Default: 1.

  • jld2_kw: Dict of kwargs to be passed to jldopen when data is written.

Example

Write out 3D fields for $u$, $v$, $w$, and a tracer $c$, along with a horizontal average:

using Oceananigans
using Oceananigans.Utils: hour, minute

model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)), tracers=(:c,))
simulation = Simulation(model, Δt=12, stop_time=1hour)

function init_save_some_metadata!(file, model)
    file["author"] = "Chim Riggles"
    file["parameters/coriolis_parameter"] = 1e-4
    file["parameters/density"] = 1027
    return nothing
end

c_avg =  Field(Average(model.tracers.c, dims=(1, 2)))

# Note that model.velocities is NamedTuple
simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities,
                                                          filename = "some_data.jld2",
                                                          schedule = TimeInterval(20minute),
                                                          init = init_save_some_metadata!)

# output
JLD2OutputWriter scheduled on TimeInterval(20 minutes):
├── filepath: ./some_data.jld2
├── 3 outputs: (u, v, w)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB

and a time- and horizontal-average of tracer $c$ every 20 minutes of simulation time to a file called some_averaged_data.jld2

simulation.output_writers[:avg_c] = JLD2OutputWriter(model, (; c=c_avg),
                                                     filename = "some_averaged_data.jld2",
                                                     schedule = AveragedTimeInterval(20minute, window=5minute))

# output
JLD2OutputWriter scheduled on TimeInterval(20 minutes):
├── filepath: ./some_averaged_data.jld2
├── 1 outputs: c averaged on AveragedTimeInterval(window=5 minutes, stride=1, interval=20 minutes)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB
source
Oceananigans.OutputWriters.NetCDFOutputWriterMethod
NetCDFOutputWriter(model, outputs; filename, schedule
                                      dir = ".",
                               array_type = Array{Float64},
                                  indices = nothing,
                               with_halos = false,
                        global_attributes = Dict(),
                        output_attributes = Dict(),
                               dimensions = Dict(),
                       overwrite_existing = false,
                             deflatelevel = 0,
                                  verbose = false)

Construct a NetCDFOutputWriter that writes (label, output) pairs in outputs (which should be a Dict) to a NetCDF file, where label is a string that labels the output and output is either a Field (e.g. model.velocities.u) or a function f(model) that returns something to be written to disk. Custom output requires the spatial dimensions (a Dict) to be manually specified (see examples).

Keyword arguments

Filenaming

  • filename (required): Descriptive filename. ".nc" is appended to filename if filename does not end in ".nc".

  • dir: Directory to save output to.

Output frequency and time-averaging

  • schedule (required): AbstractSchedule that determines when output is saved.

Slicing and type conversion prior to output

  • indices: Tuple of indices of the output variables to include. Default is (:, :, :), which includes the full fields.

  • with_halos: Boolean defining whether or not to include halos in the outputs. Default: false. Note, that to postprocess saved output (e.g., compute derivatives, etc) information about the boundary conditions is often crucial. In that case you might need to set with_halos = true.

  • array_type: The array type to which output arrays are converted to prior to saving. Default: Array{Float64}.

  • dimensions: A Dict of dimension tuples to apply to outputs (required for function outputs).

File management

  • overwrite_existing: If false, NetCDFOutputWriter will be set to append to filepath. If true, NetCDFOutputWriter will overwrite filepath if it exists or create it if not. Default: false. See NCDatasets.jl documentation for more information about its mode option.

  • deflatelevel: Determines the NetCDF compression level of data (integer 0-9; 0 (default) means no compression and 9 means maximum compression). See NCDatasets.jl documentation for more information.

Miscellaneous keywords

  • global_attributes: Dict of model properties to save with every file. Default: Dict().

  • output_attributes: Dict of attributes to be saved with each field variable (reasonable defaults are provided for velocities, buoyancy, temperature, and salinity; otherwise output_attributes must be user-provided).

Examples

Saving the $u$ velocity field and temperature fields, the full 3D fields and surface 2D slices to separate NetCDF files:

using Oceananigans

grid = RectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1))

model = NonhydrostaticModel(grid=grid, tracers=:c)

simulation = Simulation(model, Δt=12, stop_time=3600)

fields = Dict("u" => model.velocities.u, "c" => model.tracers.c)

simulation.output_writers[:field_writer] =
    NetCDFOutputWriter(model, fields, filename="fields.nc", schedule=TimeInterval(60))

# output
NetCDFOutputWriter scheduled on TimeInterval(1 minute):
├── filepath: ./fields.nc
├── dimensions: zC(16), zF(17), xC(16), yF(16), xF(16), yC(16), time(0)
├── 2 outputs: (c, u)
└── array type: Array{Float64}
simulation.output_writers[:surface_slice_writer] =
    NetCDFOutputWriter(model, fields, filename="surface_xy_slice.nc",
                       schedule=TimeInterval(60), indices=(:, :, grid.Nz))

# output
NetCDFOutputWriter scheduled on TimeInterval(1 minute):
├── filepath: ./surface_xy_slice.nc
├── dimensions: zC(1), zF(1), xC(16), yF(16), xF(16), yC(16), time(0)
├── 2 outputs: (c, u)
└── array type: Array{Float64}
simulation.output_writers[:averaged_profile_writer] =
    NetCDFOutputWriter(model, fields,
                       filename = "averaged_z_profile.nc",
                       schedule = AveragedTimeInterval(60, window=20),
                       indices = (1, 1, :))

# output
NetCDFOutputWriter scheduled on TimeInterval(1 minute):
├── filepath: ./averaged_z_profile.nc
├── dimensions: zC(16), zF(17), xC(1), yF(1), xF(1), yC(1), time(0)
├── 2 outputs: (c, u) averaged on AveragedTimeInterval(window=20 seconds, stride=1, interval=1 minute)
└── array type: Array{Float64}

NetCDFOutputWriter also accepts output functions that write scalars and arrays to disk, provided that their dimensions are provided:

using Oceananigans

Nx, Ny, Nz = 16, 16, 16

grid = RectilinearGrid(size=(Nx, Ny, Nz), extent=(1, 2, 3))

model = NonhydrostaticModel(grid=grid)

simulation = Simulation(model, Δt=1.25, stop_iteration=3)

f(model) = model.clock.time^2; # scalar output

g(model) = model.clock.time .* exp.(znodes(Center, grid)) # vector/profile output

xC, yF = xnodes(grid, Center()), ynodes(grid, Face())

XC = [xC[i] for i in 1:Nx, j in 1:Ny]
YF = [yF[j] for i in 1:Nx, j in 1:Ny]

h(model) = @. model.clock.time * sin(XC) * cos(YF) # xy slice output

outputs = Dict("scalar" => f, "profile" => g, "slice" => h)

dims = Dict("scalar" => (), "profile" => ("zC",), "slice" => ("xC", "yC"))

output_attributes = Dict(
    "scalar"  => Dict("long_name" => "Some scalar", "units" => "bananas"),
    "profile" => Dict("long_name" => "Some vertical profile", "units" => "watermelons"),
    "slice"   => Dict("long_name" => "Some slice", "units" => "mushrooms")
)

global_attributes = Dict("location" => "Bay of Fundy", "onions" => 7)

simulation.output_writers[:things] =
    NetCDFOutputWriter(model, outputs,
                       schedule=IterationInterval(1), filename="things.nc", dimensions=dims, verbose=true,
                       global_attributes=global_attributes, output_attributes=output_attributes)

# output
NetCDFOutputWriter scheduled on IterationInterval(1):
├── filepath: ./things.nc
├── dimensions: zC(16), zF(17), xC(16), yF(16), xF(16), yC(16), time(0)
├── 3 outputs: (profile, slice, scalar)
└── array type: Array{Float64}
source
Oceananigans.OutputWriters.WindowedTimeAverageType
WindowedTimeAverage(operand, model=nothing; schedule)

Returns an object for computing running averages of operand over schedule.window and recurring on schedule.interval, where schedule is an AveragedTimeInterval. During the collection period, averages are computed every schedule.stride iteration.

operand may be a Oceananigans.Field or a function that returns an array or scalar.

Calling wta(model) for wta::WindowedTimeAverage object returns wta.result.

source

Simulations

Oceananigans.Simulations.CallbackType
Callback(func, schedule=IterationInterval(1);
         parameters=nothing, callsite=TimeStepCallsite())

Return Callback that executes func on schedule at the callsite with optional parameters. By default, schedule = IterationInterval(1) and callsite = TimeStepCallsite().

If isnothing(parameters), func(sim::Simulation) is called. Otherwise, func is called via func(sim::Simulation, parameters).

The callsite determines where Callback is executed. The possible values for callsite are

  • TimeStepCallsite(): after a time-step.

  • TendencyCallsite(): after tendencies are calculated, but before taking a time-step (useful for modifying tendency calculations).

  • UpdateStateCallsite(): within update_state!, after auxiliary variables have been computed (for multi-stage time-steppers, update_state! may be called multiple times per time-step).

source
Oceananigans.Simulations.SimulationMethod
Simulation(model; Δt,
           verbose = true,
           stop_iteration = Inf,
           stop_time = Inf,
           wall_time_limit = Inf)

Construct a Simulation for a model with time step Δt.

Keyword arguments

  • Δt: Required keyword argument specifying the simulation time step. Can be a Number for constant time steps or a TimeStepWizard for adaptive time-stepping.

  • stop_iteration: Stop the simulation after this many iterations.

  • stop_time: Stop the simulation once this much model clock time has passed.

  • wall_time_limit: Stop the simulation if it's been running for longer than this many seconds of wall clock time.

source
Oceananigans.Simulations.TimeStepWizardType
TimeStepWizard([FT=Float64;]
               cfl = 0.2,
               diffusive_cfl = Inf,
               max_change = 1.1,
               min_change = 0.5,
               max_Δt = Inf,
               min_Δt = 0.0,
               cell_advection_timescale = cell_advection_timescale,
               cell_diffusion_timescale = infinite_diffusion_timescale)

Callback function that adjusts the simulation time step to meet specified target values for advective and diffusive Courant-Friedrichs-Lewy (CFL) numbers (cfl and diffusive_cfl), subject to the limits

max(min_Δt, min_change * previous_Δt) ≤ new_Δt ≤ min(max_Δt, max_change * previous_Δt)

where new_Δt is the new time step calculated by the TimeStepWizard.

For more information on the CFL number, see its wikipedia entry.

Example

To use TimeStepWizard, insert it into a Callback and then add the Callback to a Simulation:

julia> simulation = Simulation(model, Δt=0.9, stop_iteration=100)

julia> wizard = TimeStepWizard(cfl=0.2)

julia> simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(4))

Then when run!(simulation) is invoked, the time-step simulation.Δt will be updated every 4 iterations.

(Note that the name :wizard is unimportant.)

source
Oceananigans.Simulations.add_callback!Method
add_callback!(simulation, callback::Callback; name = GenericName(), callback_kw...)

add_callback!(simulation, func, schedule=IterationInterval(1); name = GenericName(), callback_kw...)

Add Callback(func, schedule) to simulation.callbacks under name. The default GenericName() generates a name of the form :callbackN, where N is big enough for the name to be unique.

If name is supplied, it may be modified if simulation.callbacks[name] already exists.

callback_kw are passed to the constructor for Callback.

The callback (which contains a schedule) can also be supplied directly.

source
Oceananigans.Simulations.conjure_time_step_wizard!Function
conjure_time_step_wizard!(simulation, schedule=IterationInterval(5), wizard_kw...)

Add a TimeStepWizard built with wizard_kw as a Callback to simulation, called on schedule which is IterationInterval(5) by default.

source
Oceananigans.Simulations.run!Method
run!(simulation; pickup=false)

Run a simulation until one of simulation.stop_criteria evaluates true. The simulation will then stop.

Picking simulations up from a checkpoint

Simulations are "picked up" from a checkpoint if pickup is either true, a String, or an Integer greater than 0.

Picking up a simulation sets field and tendency data to the specified checkpoint, leaving all other model properties unchanged.

Possible values for pickup are:

  • pickup=true picks a simulation up from the latest checkpoint associated with the Checkpointer in simulation.output_writers.

  • pickup=iteration::Int picks a simulation up from the checkpointed file associated with iteration and the Checkpointer in simulation.output_writers.

  • pickup=filepath::String picks a simulation up from checkpointer data in filepath.

Note that pickup=true and pickup=iteration fails if simulation.output_writers contains more than one checkpointer.

source

Solvers

Oceananigans.Solvers.BatchedTridiagonalSolverMethod
BatchedTridiagonalSolver(grid;
                         lower_diagonal,
                         diagonal,
                         upper_diagonal,
                         scratch = arch_array(architecture(grid), zeros(eltype(grid), size(grid)...)),
                         tridiagonal_direction = ZDirection()
                         parameters = nothing)

Construct a solver for batched tridiagonal systems on grid of the form

                    bⁱʲ¹ ϕⁱʲ¹ + cⁱʲ¹ ϕⁱʲ²   = fⁱʲ¹,
    aⁱʲᵏ⁻¹ ϕⁱʲᵏ⁻¹ + bⁱʲᵏ ϕⁱʲᵏ + cⁱʲᵏ ϕⁱʲᵏ⁺¹ = fⁱʲᵏ,  k = 2, ..., N-1
    aⁱʲᴺ⁻¹ ϕⁱʲᴺ⁻¹ + bⁱʲᴺ ϕⁱʲᴺ               = fⁱʲᴺ,

or in matrix form

    ⎡ bⁱʲ¹   cⁱʲ¹     0       ⋯         0   ⎤ ⎡ ϕⁱʲ¹ ⎤   ⎡ fⁱʲ¹ ⎤
    ⎢ aⁱʲ¹   bⁱʲ²   cⁱʲ²      0    ⋯    ⋮   ⎥ ⎢ ϕⁱʲ² ⎥   ⎢ fⁱʲ² ⎥
    ⎢  0      ⋱      ⋱       ⋱              ⎥ ⎢   .  ⎥   ⎢   .  ⎥
    ⎢  ⋮                                0   ⎥ ⎢ ϕⁱʲᵏ ⎥   ⎢ fⁱʲᵏ ⎥
    ⎢  ⋮           aⁱʲᴺ⁻²   bⁱʲᴺ⁻¹   cⁱʲᴺ⁻¹ ⎥ ⎢      ⎥   ⎢   .  ⎥
    ⎣  0      ⋯      0      aⁱʲᴺ⁻¹    bⁱʲᴺ  ⎦ ⎣ ϕⁱʲᴺ ⎦   ⎣ fⁱʲᴺ ⎦

where a is the lower_diagonal, b is the diagonal, and c is the upper_diagonal.

Note the convention used here for indexing the upper and lower diagonals; this can be different from other implementations where, e.g., aⁱʲ² may appear at the second row, instead of aⁱʲ¹ as above.

ϕ is the solution and f is the right hand side source term passed to solve!(ϕ, tridiagonal_solver, f).

a, b, c, and f can be specified in three ways:

  1. A 1D array means, e.g., that aⁱʲᵏ = a[k].

  2. A 3D array means, e.g., that aⁱʲᵏ = a[i, j, k].

Other coefficient types can be implemented by extending get_coefficient.

source
Oceananigans.Solvers.FFTBasedPoissonSolverType
FFTBasedPoissonSolver(grid, planner_flag=FFTW.PATIENT)

Return an FFTBasedPoissonSolver that solves the "generalized" Poisson equation,

\[(∇² + m) ϕ = b,\]

where $m$ is a number, using a eigenfunction expansion of the discrete Poisson operator on a staggered grid and for periodic or Neumann boundary conditions.

In-place transforms are applied to $b$, which means $b$ must have complex-valued elements (typically the same type as solver.storage).

See solve! for more information about the FFT-based Poisson solver algorithm.

source
Oceananigans.Solvers.HeptadiagonalIterativeSolverMethod
HeptadiagonalIterativeSolver(coeffs;
                             grid,
                             iterative_solver = cg!,
                             maximum_iterations = prod(size(grid)),
                             tolerance = 1e-13,
                             reduced_dim = (false, false, false), 
                             placeholder_timestep = -1.0, 
                             preconditioner_method = :Default, 
                             preconditioner_settings = nothing,
                             template = arch_array(architecture(grid), zeros(prod(size(grid)))),
                             verbose = false)

Return a HeptadiagonalIterativeSolver to solve the problem A * x = b, provided that A is a symmetric matrix.

The solver relies on a sparse version of the matrix A that is stored in matrix_constructors.

In particular, given coefficients Ax, Ay, Az, C, D, the solved problem is

    Axᵢ₊₁ ηᵢ₊₁ + Axᵢ ηᵢ₋₁ + Ayⱼ₊₁ ηⱼ₊₁ + Ayⱼ ηⱼ₋₁ + Azₖ₊₁ ηₖ₊₁ + Azₖ ηₖ₋₁ 
    - 2 ( Axᵢ₊₁ + Axᵢ + Ayⱼ₊₁ + Ayⱼ + Azₖ₊₁ + Azₖ ) ηᵢⱼₖ 
    +   ( Cᵢⱼₖ + Dᵢⱼₖ/Δt^2 ) ηᵢⱼₖ  = b

To have the equation solved at location {Center, Center, Center}, the coefficients must be specified at:

  • Ax -> {Face, Center, Center}
  • Ay -> {Center, Face, Center}
  • Az -> {Center, Center, Face}
  • C -> {Center, Center, Center}
  • D -> {Center, Center, Center}

solver.matrix is precomputed with a placeholder timestep value of placeholder_timestep = -1.0.

The sparse matrix A can be constructed with:

  • SparseMatrixCSC(constructors...) for CPU
  • CuSparseMatrixCSC(constructors...) for GPU

The matrix constructors are calculated based on the pentadiagonal coeffients passed as an input to matrix_from_coefficients function.

To allow for variable time step, the diagonal term - Az / (g * Δt²) is only added later on and it is updated only when the previous time step changes (previous_Δt != Δt).

Preconditioning is done through the various methods implemented in Solvers/sparse_preconditioners.jl.

The iterative_solver used can is to be chosen from the IterativeSolvers.jl package. The default solver is a Conjugate Gradient (cg):

solver = HeptadiagonalIterativeSolver((Ax, Ay, Az, C, D); grid)
source
Oceananigans.Solvers.PreconditionedConjugateGradientSolverMethod
PreconditionedConjugateGradientSolver(linear_operation;
                                      template_field,
                                      maxiter = size(template_field.grid),
                                      reltol = sqrt(eps(eltype(template_field.grid))),
                                      abstol = 0,
                                      preconditioner = nothing)

Returns a PreconditionedConjugateGradientSolver that solves the linear equation $A x = b$ using a iterative conjugate gradient method with optional preconditioning.

The solver is used by calling

solve!(x, solver::PreconditionedConjugateGradientOperator, b, args...)

for solver, right-hand side b, solution x, and optional arguments args....

Arguments

  • linear_operation: Function with signature linear_operation!(p, y, args...) that calculates A * y and stores the result in p for a "candidate solution" y. args... are optional positional arguments passed from solve!(x, solver, b, args...).

  • template_field: Dummy field that is the same type and size as x and b, which is used to infer the architecture, grid, and to create work arrays that are used internally by the solver.

  • maxiter: Maximum number of iterations the solver may perform before exiting.

  • reltol, abstol: Relative and absolute tolerance for convergence of the algorithm. The iteration stops when norm(A * x - b) < tolerance.

  • preconditioner: Object for which precondition!(z, preconditioner, r, args...) computes z = P * r, where r is the residual. Typically P is approximately A⁻¹.

See solve! for more information about the preconditioned conjugate-gradient algorithm.

source
Oceananigans.Solvers.solve!Function
solve!(ϕ, solver::FFTBasedPoissonSolver, b, m=0)

Solves the "generalized" Poisson equation,

\[(∇² + m) ϕ = b,\]

where $m$ is a number, using a eigenfunction expansion of the discrete Poisson operator on a staggered grid and for periodic or Neumann boundary conditions.

In-place transforms are applied to $b$, which means $b$ must have complex-valued elements (typically the same type as solver.storage).

Alternative names for 'generalized' Poisson equation

Equation $(∇² + m) ϕ = b$ is sometimes referred to as the "screened Poisson" equation when $m < 0$, or the Helmholtz equation when $m > 0$.

source
Oceananigans.Solvers.solve!Method
solve!(ϕ, solver::BatchedTridiagonalSolver, rhs, args...)

Solve the batched tridiagonal system of linear equations with right hand side rhs and lower diagonal, diagonal, and upper diagonal coefficients described by the BatchedTridiagonalSolver solver. BatchedTridiagonalSolver uses a modified TriDiagonal Matrix Algorithm (TDMA).

The result is stored in ϕ which must have size (grid.Nx, grid.Ny, grid.Nz).

Implementation follows Press William et al. (1992); §2.4. Note that a slightly different notation from Press et al. is used for indexing the off-diagonal elements; see BatchedTridiagonalSolver.

Reference

Press William, H., Teukolsky Saul, A., Vetterling William, T., & Flannery Brian, P. (1992). Numerical recipes: the art of scientific computing. Cambridge University Press

source
Oceananigans.Solvers.solve!Method
solve!(x, solver::PreconditionedConjugateGradientSolver, b, args...)

Solve A * x = b using an iterative conjugate-gradient method, where A * x is determined by solver.linear_operation

See figure 2.5 in

The Preconditioned Conjugate Gradient Method in "Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods" Barrett et. al, 2nd Edition.

Given:

  • Linear Preconditioner operator M!(solution, x, other_args...) that computes M * x = solution
  • A matrix operator A as a function A();
  • A dot product function norm();
  • A right-hand side b;
  • An initial guess x; and
  • Local vectors: z, r, p, q

This function executes the psuedocode algorithm

β  = 0
r = b - A(x)
iteration  = 0

Loop:
     if iteration > maxiter
        break
     end

     ρ = r ⋅ z

     z = M(r)
     β = ρⁱ⁻¹ / ρ
     p = z + β * p
     q = A(p)

     α = ρ / (p ⋅ q)
     x = x + α * p
     r = r - α * q

     if |r| < tolerance
        break
     end

     iteration += 1
     ρⁱ⁻¹ = ρ
source

Stokes drift

Oceananigans.StokesDrifts.StokesDriftMethod
StokesDrift(; ∂z_uˢ=zerofunction, ∂y_uˢ=zerofunction, ∂t_uˢ=zerofunction, 
              ∂z_vˢ=zerofunction, ∂x_vˢ=zerofunction, ∂t_vˢ=zerofunction, 
              ∂x_wˢ=zerofunction, ∂y_wˢ=zerofunction, ∂t_wˢ=zerofunction, parameters=nothing)

Construct a set of functions of space and time for a Stokes drift velocity field corresponding to a surface gravity wave field with an envelope that (potentially) varies in the horizontal directions.

To resolve the evolution of the Lagrangian-mean momentum, we require all the components of the "psuedovorticity",

\[𝛁 × 𝐯ˢ = \hat{\boldsymbol{x}} (∂_y wˢ - ∂_z vˢ) + \hat{\boldsymbol{y}} (∂_z uˢ - ∂_x wˢ) + \hat{\boldsymbol{z}} (∂_x vˢ - ∂_y uˢ)\]

as well as the time-derivatives of $uˢ$, $vˢ$, and $wˢ$.

Note that each function (e.g., ∂z_uˢ) is generally a function of depth, horizontal coordinates, and time.Thus, the correct function signature depends on the grid, since Flat horizontal directions are omitted.

For example, on a grid with topology = (Periodic, Flat, Bounded) (and parameters=nothing), then, e.g., ∂z_uˢ is callable via ∂z_uˢ(x, z, t). When !isnothing(parameters), then ∂z_uˢ is callable via ∂z_uˢ(x, z, t, parameters). Similarly, on a grid with topology = (Periodic, Periodic, Bounded) and parameters=nothing, ∂z_uˢ is called via ∂z_uˢ(x, y, z, t).

Example

A wavepacket moving with the group velocity in the $x$-direction. We write the Stokes drift as:

\[uˢ(x, y, z, t) = A(x - cᵍ \, t, y) ûˢ(z)\]

with $A(ξ, η) = \exp{[-(ξ^2 + η^2) / 2δ^2]}$. We also assume $vˢ = 0$. If $𝐯ˢ$ represents the solenoidal component of the Stokes drift, then in this system from incompressibility requirement we have that $∂_z wˢ = - ∂_x uˢ = - (∂_ξ A) ûˢ$ and therefore, under the assumption that $wˢ$ tends to zero at large depths, we get $wˢ = - (∂_ξ A / 2k) ûˢ$.

using Oceananigans
using Oceananigans.Units

g = 9.81 # gravitational acceleration

ϵ = 0.1
λ = 100meters  # horizontal wavelength
const k = 2π / λ  # horizontal wavenumber
c = sqrt(g / k)  # phase speed
const δ = 400kilometers  # wavepacket spread
const cᵍ = c / 2  # group speed
const Uˢ = ϵ^2 * c

@inline A(ξ, η) = exp(- (ξ^2 + η^2) / 2δ^2)

@inline ∂ξ_A(ξ, η) = - ξ / δ^2 * A(ξ, η)
@inline ∂η_A(ξ, η) = - η / δ^2 * A(ξ, η)
@inline ∂η_∂ξ_A(ξ, η) = η * ξ / δ^4 * A(ξ, η)
@inline ∂²ξ_A(ξ, η) = (ξ^2 / δ^2 - 1) * A(ξ, η) / δ^2

@inline ûˢ(z) = Uˢ * exp(2k * z)
@inline uˢ(x, y, z, t) = A(x - cᵍ * t, y) * ûˢ(z)

@inline ∂z_uˢ(x, y, z, t) = 2k * A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂y_uˢ(x, y, z, t) = ∂η_A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂t_uˢ(x, y, z, t) = - cᵍ * ∂ξ_A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂x_wˢ(x, y, z, t) = - 1 / 2k * ∂²ξ_A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂y_wˢ(x, y, z, t) = - 1 / 2k * ∂η_∂ξ_A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂t_wˢ(x, y, z, t) = + cᵍ / 2k * ∂²ξ_A(x - cᵍ * t, y) * ûˢ(z)

stokes_drift = StokesDrift(; ∂z_uˢ, ∂t_uˢ, ∂y_uˢ, ∂t_wˢ, ∂x_wˢ, ∂y_wˢ)

# output

StokesDrift{Nothing}:
├── ∂x_vˢ: zerofunction
├── ∂x_wˢ: ∂x_wˢ
├── ∂y_uˢ: ∂y_uˢ
├── ∂y_wˢ: ∂y_wˢ
├── ∂z_uˢ: ∂z_uˢ
├── ∂z_vˢ: zerofunction
├── ∂t_uˢ: ∂t_uˢ
├── ∂t_vˢ: zerofunction
└── ∂t_wˢ: ∂t_wˢ
source
Oceananigans.StokesDrifts.UniformStokesDriftMethod
UniformStokesDrift(; ∂z_uˢ=zerofunction, ∂z_vˢ=zerofunction, ∂t_uˢ=zerofunction, ∂t_vˢ=zerofunction, parameters=nothing)

Construct a set of functions for a Stokes drift velocity field corresponding to a horizontally-uniform surface gravity wave field, with optional parameters.

If parameters=nothing, then the functions ∂z_uˢ, ∂z_vˢ, ∂t_uˢ, ∂t_vˢ must be callable with signature (z, t). If !isnothing(parameters), then functions must be callable with the signature (z, t, parameters).

To resolve the evolution of the Lagrangian-mean momentum, we require vertical-derivatives and time-derivatives of the horizontal components of the Stokes drift, and .

Examples

Exponentially decaying Stokes drift corresponding to a surface Stokes drift of uˢ(z=0) = 0.005 and decay scale h = 20:

using Oceananigans

@inline uniform_stokes_shear(z, t) = 0.005 * exp(z / 20)

stokes_drift = UniformStokesDrift(∂z_uˢ=uniform_stokes_shear)

# output

UniformStokesDrift{Nothing}:
├── ∂z_uˢ: uniform_stokes_shear
├── ∂z_vˢ: zerofunction
├── ∂t_uˢ: zerofunction
└── ∂t_vˢ: zerofunction

Exponentially-decaying Stokes drift corresponding to a surface Stokes drift of uˢ = 0.005 and decay scale h = 20, using parameters:

using Oceananigans

@inline uniform_stokes_shear(z, t, p) = p.uˢ * exp(z / p.h)

stokes_drift_parameters = (uˢ = 0.005, h = 20)
stokes_drift = UniformStokesDrift(∂z_uˢ=uniform_stokes_shear, parameters=stokes_drift_parameters)

# output

UniformStokesDrift with parameters (uˢ=0.005, h=20):
├── ∂z_uˢ: uniform_stokes_shear
├── ∂z_vˢ: zerofunction
├── ∂t_uˢ: zerofunction
└── ∂t_vˢ: zerofunction
source

Time steppers

Oceananigans.TimeSteppers.ClockType
mutable struct Clock{T<:Number}

Keeps track of the current time, iteration number, and time-stepping stage. The stage is updated only for multi-stage time-stepping methods. The time::T is either a number or a DateTime object.

source
Oceananigans.TimeSteppers.ClockMethod
Clock(; time, iteration=0, stage=1)

Returns a Clock object. By default, Clock is initialized to the zeroth iteration and first time step stage.

source
Oceananigans.TimeSteppers.QuasiAdamsBashforth2TimeStepperMethod
QuasiAdamsBashforth2TimeStepper(grid, tracers,
                                χ = 0.1;
                                implicit_solver = nothing,
                                Gⁿ = TendencyFields(grid, tracers),
                                G⁻ = TendencyFields(grid, tracers))

Return a 2nd-order quasi Adams-Bashforth (AB2) time stepper (QuasiAdamsBashforth2TimeStepper) on grid, with tracers, and AB2 parameter χ. The tendency fields Gⁿ and G⁻ can be specified via optional kwargs.

The 2nd-order quasi Adams-Bashforth timestepper steps forward the state Uⁿ by Δt via

Uⁿ⁺¹ = Uⁿ + Δt * [(3/2 + χ) * Gⁿ - (1/2 + χ) * Gⁿ⁻¹]

where Uⁿ is the state at the $n$-th timestep, Gⁿ is the tendency at the $n$-th timestep, and Gⁿ⁻¹ is the tendency at the previous timestep (G⁻).

First timestep

For the first timestep, since there are no saved tendencies from the previous timestep, the QuasiAdamsBashforth2TimeStepper performs an Euler timestep:

Uⁿ⁺¹ = Uⁿ + Δt * Gⁿ
source
Oceananigans.TimeSteppers.RungeKutta3TimeStepperMethod
RungeKutta3TimeStepper(grid, tracers;
                       implicit_solver = nothing,
                       Gⁿ = TendencyFields(grid, tracers),
                       G⁻ = TendencyFields(grid, tracers))

Return a 3rd-order Runge0Kutta timestepper (RungeKutta3TimeStepper) on grid and with tracers. The tendency fields Gⁿ and G⁻ can be specified via optional kwargs.

The scheme described by Le and Moin (1991). In a nutshel, the 3rd-order Runge Kutta timestepper steps forward the state Uⁿ by Δt via 3 substeps. A pressure correction step is applied after at each substep.

The state U after each substep m is

Uᵐ⁺¹ = Uᵐ + Δt * (γᵐ * Gᵐ + ζᵐ * Gᵐ⁻¹)

where Uᵐ is the state at the $m$-th substep, Gᵐ is the tendency at the $m$-th substep, Gᵐ⁻¹ is the tendency at the previous substep, and constants $γ¹ = 8/15$, $γ² = 5/12$, $γ³ = 3/4$, $ζ¹ = 0$, $ζ² = -17/60$, $ζ³ = -5/12$.

The state at the first substep is taken to be the one that corresponds to the $n$-th timestep, U¹ = Uⁿ, and the state after the third substep is then the state at the Uⁿ⁺¹ = U⁴.

source
Oceananigans.TimeSteppers.time_step!Method
time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt; euler=false, compute_tendencies=true)

Step forward model one time step Δt with a 2nd-order Adams-Bashforth method and pressure-correction substep. Setting euler=true will take a forward Euler time step. Setting compute_tendencies=false will not calculate new tendencies

source
Oceananigans.TimeSteppers.time_step!Method
time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt)

Step forward model one time step Δt with a 3rd-order Runge-Kutta method. The 3rd-order Runge-Kutta method takes three intermediate substep stages to achieve a single timestep. A pressure correction step is applied at each intermediate stage.

source

Turbulence closures

Oceananigans.TurbulenceClosures.AnisotropicMinimumDissipationMethod
AnisotropicMinimumDissipation([time_discretization = ExplicitTimeDiscretization, FT = Float64;]
                              C = 1/12, Cν = nothing, Cκ = nothing, Cb = nothing)

Return parameters of type FT for the AnisotropicMinimumDissipation turbulence closure.

Arguments

  • time_discretization: Either ExplicitTimeDiscretization() or VerticallyImplicitTimeDiscretization(), which integrates the terms involving only $z$-derivatives in the viscous and diffusive fluxes with an implicit time discretization. Default ExplicitTimeDiscretization().

  • FT: Float type; default Float64.

Keyword arguments

  • C: Poincaré constant for both eddy viscosity and eddy diffusivities. C is overridden for eddy viscosity or eddy diffusivity if or are set, respecitvely.

  • : Poincaré constant for momentum eddy viscosity.

  • : Poincaré constant for tracer eddy diffusivities. If one number or function, the same number or function is applied to all tracers. If a NamedTuple, it must possess a field specifying the Poncaré constant for every tracer.

  • Cb: Buoyancy modification multiplier (Cb = nothing turns it off, Cb = 1 was used by Abkar et al. (2016)). Note: that we do not subtract the horizontally-average component before computing this buoyancy modification term. This implementation differs from Abkar et al. (2016)'s proposal and the impact of this approximation has not been tested or validated.

By default: C = Cν = Cκ = 1/12, which is appropriate for a finite-volume method employing a second-order advection scheme, and Cb = nothing, which turns off the buoyancy modification term.

or may be numbers, or functions of x, y, z.

Examples

julia> using Oceananigans

julia> pretty_diffusive_closure = AnisotropicMinimumDissipation(C=1/2)
AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with:
           Poincaré constant for momentum eddy viscosity Cν: 0.5
    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: 0.5
                        Buoyancy modification multiplier Cb: nothing
julia> using Oceananigans

julia> const Δz = 0.5; # grid resolution at surface

julia> surface_enhanced_tracer_C(x, y, z) = 1/12 * (1 + exp((z + Δz/2) / 8Δz));

julia> fancy_closure = AnisotropicMinimumDissipation(Cκ=surface_enhanced_tracer_C)
AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with:
           Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: surface_enhanced_tracer_C
                        Buoyancy modification multiplier Cb: nothing
julia> using Oceananigans

julia> tracer_specific_closure = AnisotropicMinimumDissipation(Cκ=(c₁=1/12, c₂=1/6))
AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with:
           Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: (c₁ = 0.08333333333333333, c₂ = 0.16666666666666666)
                        Buoyancy modification multiplier Cb: nothing

References

Vreugdenhil C., and Taylor J. (2018), "Large-eddy simulations of stratified plane Couette flow using the anisotropic minimum-dissipation model", Physics of Fluids 30, 085104.

Verstappen, R. (2018), "How much eddy dissipation is needed to counterbalance the nonlinear production of small, unresolved scales in a large-eddy simulation of turbulence?", Computers & Fluids 176, pp. 276-284.

source
Oceananigans.TurbulenceClosures.ConvectiveAdjustmentVerticalDiffusivityType
ConvectiveAdjustmentVerticalDiffusivity([time_discretization = VerticallyImplicitTimeDiscretization(), FT=Float64;]
                                        convective_κz = 0,
                                        convective_νz = 0,
                                        background_κz = 0,
                                        background_νz = 0)

Return a convective adjustment vertical diffusivity closure that applies different values of diffusivity and/or viscosity depending whether the region is statically stable (positive or zero buoyancy gradient) or statically unstable (negative buoyancy gradient).

Arguments

  • time_discretization: Either ExplicitTimeDiscretization() or VerticallyImplicitTimeDiscretization(); default VerticallyImplicitTimeDiscretization().

  • FT: Float type; default Float64.

Keyword arguments

  • convective_κz: Vertical tracer diffusivity in regions with negative (unstable) buoyancy gradients. Either a single number, function, array, field, or tuple of diffusivities for each tracer.

  • background_κz: Vertical tracer diffusivity in regions with zero or positive (stable) buoyancy gradients.

  • convective_νz: Vertical viscosity in regions with negative (unstable) buoyancy gradients. Either a number, function, array, or field.

  • background_κz: Vertical viscosity in regions with zero or positive (stable) buoyancy gradients.

Example

julia> using Oceananigans

julia> cavd = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1)
ConvectiveAdjustmentVerticalDiffusivity{VerticallyImplicitTimeDiscretization}(background_κz=0.0 convective_κz=1 background_νz=0.0 convective_νz=0.0)
source
Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivityMethod
IsopycnalSkewSymmetricDiffusivity([time_disc=VerticallyImplicitTimeDiscretization(), FT=Float64;]
                                  κ_skew = 0,
                                  κ_symmetric = 0,
                                  isopycnal_tensor = SmallSlopeIsopycnalTensor(),
                                  slope_limiter = FluxTapering(1e-2))

Return parameters for an isopycnal skew-symmetric tracer diffusivity with skew diffusivity κ_skew and symmetric diffusivity κ_symmetric that uses an isopycnal_tensor model for for calculating the isopycnal slopes, and (optionally) applying a slope_limiter to the calculated isopycnal slope values.

Both κ_skew and κ_symmetric may be constants, arrays, fields, or functions of (x, y, z, t).

source
Oceananigans.TurbulenceClosures.RiBasedVerticalDiffusivityType
RiBasedVerticalDiffusivity([time_discretization = VerticallyImplicitTimeDiscretization(),
                           FT = Float64;]
                           Ri_dependent_tapering = HyperbolicTangentRiDependentTapering(),
                           ν₀  = 0.7,
                           κ₀  = 0.5,
                           κᶜᵃ = 1.7,
                           Cᵉⁿ = 0.1,
                           Cᵃᵛ = 0.6,
                           Ri₀ = 0.1,
                           Riᵟ = 0.4,
                           warning = true)

Return a closure that estimates the vertical viscosity and diffusivity from "convective adjustment" coefficients ν₀ and κ₀ multiplied by a decreasing function of the Richardson number, $Ri$.

Arguments

  • time_discretization: Either ExplicitTimeDiscretization() or VerticallyImplicitTimeDiscretization(), which integrates the terms involving only $z$-derivatives in the viscous and diffusive fluxes with an implicit time discretization. Default VerticallyImplicitTimeDiscretization().

  • FT: Float type; default Float64.

Keyword arguments

  • Ri_dependent_tapering: The $Ri$-dependent tapering. Options are: PiecewiseLinearRiDependentTapering(), HyperbolicTangentRiDependentTapering() (default), and ExponentialRiDependentTapering().

  • ν₀: Non-convective viscosity.

  • κ₀: Non-convective diffusivity for tracers.

  • κᶜᵃ: Convective adjustment diffusivity for tracers.

  • Cᵉⁿ: Entrainment coefficient for tracers.

  • Cᵃᵛ: Time-averaging coefficient for viscosity and diffusivity.

  • Ri₀: $Ri$ threshold for decreasing viscosity and diffusivity.

  • Riᵟ: $Ri$-width over which viscosity and diffusivity decreases to 0.

source
Oceananigans.TurbulenceClosures.ScalarBiharmonicDiffusivityType
ScalarBiharmonicDiffusivity(formulation = ThreeDimensionalFormulation(), FT = Float64;
                            ν = 0,
                            κ = 0,
                            discrete_form = false,
                            loc = (nothing, nothing, nothing),
                            parameters = nothing)

Return a scalar biharmonic diffusivity turbulence closure with viscosity coefficient ν and tracer diffusivities κ for each tracer field in tracers. If a single κ is provided, it is applied to all tracers. Otherwise κ must be a NamedTuple with values for every tracer individually.

Arguments

  • formulation:

    • HorizontalFormulation() for diffusivity applied in the horizontal direction(s)
    • VerticalFormulation() for diffusivity applied in the vertical direction,
    • ThreeDimensionalFormulation() (default) for diffusivity applied isotropically to all directions
  • FT: the float datatype (default: Float64)

Keyword arguments

  • ν: Viscosity. Number, AbstractArray, Field, or Function.

  • κ: Diffusivity. Number, three-dimensional AbstractArray, Field, Function, or NamedTuple of diffusivities with entries for each tracer.

  • discrete_form: Boolean; default: False.

When prescribing the viscosities or diffusivities as functions, depending on the value of keyword argument discrete_form, the constructor expects:

  • discrete_form = false (default): functions of the grid's native coordinates and time, e.g., (x, y, z, t) for a RectilinearGrid or (λ, φ, z, t) for a LatitudeLongitudeGrid.

  • discrete_form = true:

    • with loc = (nothing, nothing, nothing) (default): functions of (i, j, k, grid, ℓx, ℓy, ℓz) with ℓx, ℓy and ℓz either Face() or Center().
    • with loc = (ℓx, ℓy, ℓz) with ℓx, ℓy and ℓz either Face() or Center(): functions of (i, j, k, grid).
  • parameters: NamedTuple with parameters used by the functions that compute viscosity and/or diffusivity; default: nothing.

For examples see ScalarDiffusivity.

source
Oceananigans.TurbulenceClosures.ScalarDiffusivityType
ScalarDiffusivity(time_discretization = ExplicitTimeDiscretization(),
                  formulation = ThreeDimensionalFormulation(), FT = Float64;
                  ν = 0,
                  κ = 0,
                  discrete_form = false,
                  loc = (nothing, nothing, nothing),
                  parameters = nothing)

Return ScalarDiffusivity turbulence closure with viscosity ν and tracer diffusivities κ for each tracer field in tracers. If a single κ is provided, it is applied to all tracers. Otherwise κ must be a NamedTuple with values for every tracer individually.

Arguments

* `time_discretization`: either `ExplicitTimeDiscretization()` (default)
  or `VerticallyImplicitTimeDiscretization()`.

* `formulation`:
  - `HorizontalFormulation()` for diffusivity applied in the horizontal direction(s)
  - `VerticalFormulation()` for diffusivity applied in the vertical direction,
  - `ThreeDimensionalFormulation()` (default) for diffusivity applied isotropically to all directions

* `FT`: the float datatype (default: `Float64`)

Keyword arguments

* `ν`: Viscosity. `Number`, three-dimensional `AbstractArray`, `Field`, or `Function`.

* `κ`: Diffusivity. `Number`, `AbstractArray`, `Field`, `Function`, or
       `NamedTuple` of diffusivities with entries for each tracer.

* `discrete_form`: `Boolean`; default: `False`.

When prescribing the viscosities or diffusivities as functions, depending on the value of keyword argument discrete_form, the constructor expects:

* `discrete_form = false` (default): functions of the grid's native coordinates
  and time, e.g., `(x, y, z, t)` for a `RectilinearGrid` or
  `(λ, φ, z, t)` for a `LatitudeLongitudeGrid`.
  
* `discrete_form = true`: 
    - with `loc = (nothing, nothing, nothing)` (default):
      functions of `(i, j, k, grid, ℓx, ℓy, ℓz)` with `ℓx`, `ℓy`
      and `ℓz` either `Face()` or `Center()`.
    - with `loc = (ℓx, ℓy, ℓz)` with `ℓx`, `ℓy`
      and `ℓz` either `Face()` or `Center()`: functions of `(i, j, k, grid)`.

* `parameters`: `NamedTuple` with parameters used by the functions
  that compute viscosity and/or diffusivity; default: `nothing`.

Examples

julia> using Oceananigans

julia> ScalarDiffusivity(ν = 1000, κ=2000)
ScalarDiffusivity{ExplicitTimeDiscretization}(ν=1000.0, κ=2000.0)
julia> const depth_scale = 100;

julia> @inline ν(x, y, z) = 1000 * exp(z / depth_scale)
ν (generic function with 1 method)

julia> ScalarDiffusivity(ν = ν)
ScalarDiffusivity{ExplicitTimeDiscretization}(ν=ν (generic function with 1 method), κ=0.0)
julia> using Oceananigans.Grids: znode

julia> @inline function κ(i, j, k, grid, ℓx, ℓy, ℓz)
           z = znode(i, j, k, grid, ℓx, ℓy, ℓz)
           return 2000 * exp(z / depth_scale)
       end
κ (generic function with 1 method)

julia> ScalarDiffusivity(κ = κ, discrete_form = true)
ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0, κ=Oceananigans.TurbulenceClosures.DiscreteDiffusionFunction{Nothing, Nothing, Nothing, Nothing, typeof(κ)})
julia> @inline function another_κ(i, j, k, grid, p)
           z = znode(i, j, k, grid)
           return 2000 * exp(z / p.depth_scale)
       end
another_κ (generic function with 1 method)

julia> ScalarDiffusivity(κ = another_κ, discrete_form = true, loc = (Center, Center, Face), parameters = (; depth_scale = 120.0))
ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0, κ=Oceananigans.TurbulenceClosures.DiscreteDiffusionFunction{Center, Center, Face, NamedTuple{(:depth_scale,), Tuple{Float64}}, typeof(another_κ)})
source
Oceananigans.TurbulenceClosures.SmagorinskyLillyMethod
SmagorinskyLilly([time_discretization::TD = ExplicitTimeDiscretization(), FT=Float64;] C=0.16, Cb=1.0, Pr=1.0)

Return a SmagorinskyLilly type associated with the turbulence closure proposed by Lilly (1962), Smagorinsky (1958), Smagorinsky (1963), and Lilly (1966), which has an eddy viscosity of the form

νₑ = (C * Δᶠ)² * √(2Σ²) * √(1 - Cb * N² / Σ²)

and an eddy diffusivity of the form

κₑ = νₑ / Pr

where Δᶠ is the filter width, Σ² = ΣᵢⱼΣᵢⱼ is the double dot product of the strain tensor Σᵢⱼ, Pr is the turbulent Prandtl number, is the total buoyancy gradient, and Cb is a constant the multiplies the Richardson number modification to the eddy viscosity.

Arguments

  • time_discretization: Either ExplicitTimeDiscretization() or VerticallyImplicitTimeDiscretization(), which integrates the terms involving only $z$-derivatives in the viscous and diffusive fluxes with an implicit time discretization. Default ExplicitTimeDiscretization().

  • FT: Float type; default Float64.

Keyword arguments

  • C: Smagorinsky constant. Default value is 0.16 as obtained by Lilly (1966).

  • Cb: Buoyancy term multipler based on Lilly (1962) (Cb = 0 turns it off, Cb ≠ 0 turns it on. Typically, and according to the original work by Lilly (1962), Cb = 1 / Pr.)

  • Pr: Turbulent Prandtl numbers for each tracer. Either a constant applied to every tracer, or a NamedTuple with fields for each tracer individually.

References

Smagorinsky, J. "On the numerical integration of the primitive equations of motion for baroclinic flow in a closed region." Monthly Weather Review (1958)

Lilly, D. K. "On the numerical simulation of buoyant convection." Tellus (1962)

Smagorinsky, J. "General circulation experiments with the primitive equations: I. The basic experiment." Monthly Weather Review (1963)

Lilly, D. K. "The representation of small-scale turbulence in numerical simulation experiments." NCAR Manuscript No. 281, 0, (1966)

source
Oceananigans.TurbulenceClosures.TwoDimensionalLeithType
TwoDimensionalLeith(FT=Float64;
                    C=0.3, C_Redi=1, C_GM=1,
                    isopycnal_model=SmallSlopeIsopycnalTensor())

Return a TwoDimensionalLeith type associated with the turbulence closure proposed by Leith (1968) and Fox-Kemper and Menemenlis (2008) which has an eddy viscosity of the form

νₑ = (C * Δᶠ)³ * √(|∇ₕ ζ|² + |∇ₕ ∂w/∂z|²)

and an eddy diffusivity of the form...

where Δᶠ is the filter width, ζ = ∂v/∂x - ∂u/∂y is the vertical vorticity, and C is a model constant.

Keyword arguments

  • C: Model constant
  • C_Redi: Coefficient for down-gradient tracer diffusivity for each tracer. Either a constant applied to every tracer, or a NamedTuple with fields for each tracer individually.
  • C_GM: Coefficient for down-gradient tracer diffusivity for each tracer. Either a constant applied to every tracer, or a NamedTuple with fields for each tracer individually.

References

Leith, C. E. (1968). "Diffusion Approximation for Two‐Dimensional Turbulence", The Physics of Fluids 11, 671. doi: 10.1063/1.1691968

Fox‐Kemper, B., & D. Menemenlis (2008), "Can large eddy simulation techniques improve mesoscale rich ocean models?", in Ocean Modeling in an Eddying Regime, Geophys. Monogr. Ser., 177, pp. 319–337. doi: 10.1029/177GM19

source
Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretizationType
struct VerticallyImplicitTimeDiscretization <: AbstractTimeDiscretization

A vertically-implicit time-discretization of a TurbulenceClosure.

This implies that a flux divergence such as $𝛁 ⋅ 𝐪$ at the $n$-th timestep is time-discretized as

[∇ ⋅ q]ⁿ = [explicit_flux_divergence]ⁿ + [∂z (κ ∂z c)]ⁿ⁺¹
source

Utilities

Oceananigans.Utils.IterationIntervalMethod
IterationInterval(interval; offset=0)

Return a callable IterationInterval that "actuates" (schedules output or callback execution) whenever the model iteration (modified by offset) is a multiple of interval.

For example,

  • IterationInterval(100) actuates at iterations [100, 200, 300, ...].
  • IterationInterval(100, offset=-1) actuates at iterations [99, 199, 299, ...].
source
Oceananigans.Utils.KernelParametersMethod
KernelParameters(size, offsets)

Return parameters for kernel launching and execution that define (i) a tuple that defines the size of the kernel being launched and (ii) a tuple of offsets that offset loop indices. For example, offsets = (0, 0, 0) with size = (N, N, N) means all indices loop from 1:N. If offsets = (1, 1, 1), then all indices loop from 2:N+1. And so on.

Example

size = (8, 6, 4)
offsets = (0, 1, 2)
kp = KernelParameters(size, offsets)

# Launch a kernel with indices that range from i=1:8, j=2:7, k=3:6,
# where i, j, k are the first, second, and third index, respectively:
launch!(arch, grid, kp, kernel!; kernel_args...)

See the documentation for launch!.

source
Oceananigans.Utils.SpecifiedTimesMethod
SpecifiedTimes(times)

Return a callable TimeInterval that "actuates" (schedules output or callback execution) whenever the model's clock equals the specified values in times. For example,

  • SpecifiedTimes([1, 15.3]) actuates when model.clock.time is 1 and 15.3.
Sorting specified times

The specified times need not be ordered as the SpecifiedTimes constructor will check and order them in ascending order if needed.

source
Oceananigans.Utils.TimeIntervalType
struct TimeInterval <: AbstractSchedule

Callable TimeInterval schedule for periodic output or diagnostic evaluation according to model.clock.time.

source
Oceananigans.Utils.TimeIntervalMethod
TimeInterval(interval)

Return a callable TimeInterval that schedules periodic output or diagnostic evaluation on a interval of simulation time, as kept by model.clock.

source
Oceananigans.Utils.WallTimeIntervalMethod
WallTimeInterval(interval; start_time = time_ns() * 1e-9)

Return a callable WallTimeInterval that schedules periodic output or diagnostic evaluation on a interval of "wall time" while a simulation runs, in units of seconds.

The "wall time" is the actual real world time in seconds, as kept by an actual or hypothetical clock hanging on your wall.

The keyword argument start_time can be used to specify a starting wall time other than the moment WallTimeInterval is constructed.

source
Oceananigans.Utils.launch!Method
launch!(arch, grid, layout, kernel!, args...; kwargs...)

Launches kernel!, with arguments args and keyword arguments kwargs, over the dims of grid on the architecture arch. kernels run on the default stream

source
Oceananigans.Utils.pretty_filesizeFunction
pretty_filesize(s, suffix="B")

Convert a floating point value s representing a file size to a more human-friendly formatted string with one decimal places with a suffix defaulting to "B". Depending on the value of s the string will be formatted to show s using an SI prefix from bytes, kiB (1024 bytes), MiB (1024² bytes), and so on up to YiB (1024⁸ bytes).

source
Oceananigans.Utils.prettytimeFunction
prettytime(t, longform=true)

Convert a floating point value t representing an amount of time in SI units of seconds to a human-friendly string with three decimal places. Depending on the value of t the string will be formatted to show t in nanoseconds (ns), microseconds (μs), milliseconds (ms), seconds, minutes, hours, or days.

With longform=false, we use s, m, hrs, and d in place of seconds, minutes, and hours.

source
Oceananigans.Utils.with_tracersMethod
with_tracers(tracer_names, initial_tuple, tracer_default)

Create a tuple corresponding to the solution variables u, v, w, and tracer_names. initial_tuple is a NamedTuple that at least has fields u, v, and w, and may have some fields corresponding to the names in tracer_names. tracer_default is a function that produces a default tuple value for each tracer if not included in initial_tuple.

source
Oceananigans.Utils.work_layoutMethod
work_layout(grid, dims; include_right_boundaries=false, location=nothing)

Returns the workgroup and worksize for launching a kernel over dims on grid. The workgroup is a tuple specifying the threads per block in each dimension. The worksize specifies the range of the loop in each dimension.

Specifying include_right_boundaries=true will ensure the work layout includes the right face end points along bounded dimensions. This requires the field location to be specified.

For more information, see: https://github.com/CliMA/Oceananigans.jl/pull/308

source

Units

Oceananigans.Units.GiBConstant
GiB

A Float64 constant equal to 1024MiB. Useful for increasing the clarity of scripts, e.g. max_filesize = 50GiB.

source
Oceananigans.Units.KiBConstant
KiB

A Float64 constant equal to 1024.0. Useful for increasing the clarity of scripts, e.g. max_filesize = 250KiB.

source
Oceananigans.Units.MiBConstant
MiB

A Float64 constant equal to 1024KiB. Useful for increasing the clarity of scripts, e.g. max_filesize = 100MiB.

source
Oceananigans.Units.TiBConstant
TiB

A Float64 constant equal to 1024GiB. Useful for increasing the clarity of scripts, e.g. max_filesize = 2TiB.

source
Oceananigans.Units.dayConstant
day

A Float64 constant equal to 24hours. Useful for increasing the clarity of scripts, e.g. stop_time = 1day.

source
Oceananigans.Units.daysConstant
days

A Float64 constant equal to 24hours. Useful for increasing the clarity of scripts, e.g. stop_time = 7days.

source
Oceananigans.Units.hourConstant
hour

A Float64 constant equal to 60minutes. Useful for increasing the clarity of scripts, e.g. Δt = 1hour.

source
Oceananigans.Units.hoursConstant
hours

A Float64 constant equal to 60minutes. Useful for increasing the clarity of scripts, e.g. Δt = 3hours.

source
Oceananigans.Units.TimeType
Time(t)

Return a time "selector" at the continuous time t for linearly interpolating FieldTimeSeries.

Examples

# Interpolate `field_time_series` to `t=0.1`, returning `interpolated::Field`
interpolated = field_time_series[Time(0.1)]

# Interpolate `field_time_series` at `i, j, k` and `t=0.1`
interpolated_ijk = field_time_series[i, j, k, Time(0.1)]
source