Full Subzero API documentation

YourGridType{FT} <: AbstractRectilinearGrid{FT}

Each simulation run with Subzero.jl must be run on a grid. A grid defines the points at which the ocean and atmosphere hold velocity data, as well as various other tracers, at grid points. The user must choose an existing subtype of AbstractRectilinearGrid or implement a new subtype to create a simulation.

Each grid implementation must define the number and dimensions of each grid cell. Right now, this assumes that the ocean and the atmosphere are on the same grid and that the grid is rectangular. We might not want this to be true in the future. Furthermore, each concrete implementation of AbstractRectilinearGrid that will be used to two-way couple with an ocean or an atmosphere must have a field called floe_locations that is a matrix of CellFloes, with one element for each grid point. The user should not worry about CellFloes, this is up to the developer to make sure that their grid contains and populates this field. Again, in the future, this might be related to the ocean and/or atmosphere rather than the grid. For more information, or if you're interested in working on this, see issue #107.

API

The following methods must be implemented for all subtypes:

• _get_grid_extent(grid::AbstractRectilinearGrid):

The _get_grid_extent function takes in a concrete subtype of AbstractRectilinearGrid and returns the grid's minimum x-value, maximum x-value, minimum y-value, and maximum y-value.

Given the original code was written for RegRectilinearGrid objects, the dispatch isn't fully implemented and other functions should be added to this list. If you were interested in adding a new subtype of AbstractRectilinearGrid, see issue #108.

RegRectilinearGrid{FT} <: AbstractRectilinearGrid{FT}

A concrete implementation of an AbstractRectilinearGrid that represents a tessellation of 2-dimensional Euclidean space into n-by-m congruent rectangles. Fields that hold float data are of type FT, a concrete subtype of AbstractFloat.

• Nx::Int: number of grid cells in the x-direction
• Ny::Int: number of grid cells in the y-direction
• Δx::FT: grid cell width
• Δy::FT: grid cell height
• x0::FT: value of first x grid line
• xf::FT: value of final x grid line
• y0::FT: value of first y grid line
• yf::FT: value of final y grid line
• floe_locations::Matrix{CellFloes}: Nx + 1 by Ny + 1 matrix of CellFloes

Here is how to construct a RegRectilinearGrid:

RegRectilinearGrid([FT = Float64]; x0, xf, y0, yf, Nx = nothing, Ny = nothing, Δx = nothing, Δy = nothing)

Positional arguments

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• x0::FT: value of first x grid line
• xf::FT: value of final x grid line
• y0::FT: value of first y grid line
• yf::FT: value of final y grid line
• Nx::Int: number of grid cells in the x-direction
• Ny::Int: number of grid cells in the y-direction
• Δx::FT: grid cell width
• Δy::FT: grid cell height

Examples

• Defining a RegRectilinearGrid using Nx and Ny.

julia> grid = RegRectilinearGrid(; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 10)
RegRectilinearGrid{Float64}
  ⊢x extent (0.0 to 500000.0) with 20 grid cells of size 25000.0 m
  ∟y extent (0.0 to 500000.0) with 10 grid cells of size 50000.0 m

• Defining a RegRectilinearGrid using Δx and Δy.

julia> grid = RegRectilinearGrid(Float32; x0 = -5e5, xf = 5e5, y0 = 0.0, yf = 5e5, Δx = 5e4, Δy = 5e4)
RegRectilinearGrid{Float32}
  ⊢x extent (-500000.0 to 500000.0) with 20 grid cells of size 50000.0 m
  ∟y extent (0.0 to 500000.0) with 10 grid cells of size 50000.0 m

• Error due to attemping to define a RegRectilinearGrid using Δx and Ny.

julia> grid = RegRectilinearGrid(; x0 = -5e5, xf = 5e5, y0 = 0.0, yf = 5e5, Δx = 5e4, Ny = 10)
ERROR: ArgumentError: To create a RegRectilinearGrid, either provide Δx and Δy OR Nx and Ny.
[...]

YourDirection <: AbstractDirection

Each domain within a Subzero.jl Model must have four (4) boundaries (subtypes of AbstractBoundary) where each of these boundaries is parametrically typed by the direction of the boundary. The user will first choose one of the four cardinal directions, the subtypes of AbstractDirection:

• North
• South
• East
• West

This abstract type is not meant to be extended by the user, unless the user wants to move away from a rectangular domain with assigned cardinal directions for each wall. This would a more major redesign and the user should check out the developer documentation.

_API

• _boundary_info_from_extent(D::Type{AbstractDirection}, FT::Type{AbstractFloat}, x0, xf, y0, yf)

The _boundary_info_from_extent function finds the values for the poly and val fields of a subtype of AbstractBoundary of direction D given a region of the extent defined by x0, xf, y0, and yf given that these are the minimum and maximum x and y values of the region desired. The returned extent will have values of type FT.

North<:AbstractDirection

A simple subtype of AbstractDirection used for parametrically typing a subtype of AbstractBoundary if that boundary is the northern boundary in a rectangular domain.

South<:AbstractDirection

A simple subtype of AbstractDirection used for parametrically typing a subtype of AbstractBoundary if that boundary is the southern boundary in a rectangular domain.

East<:AbstractDirection

A simple subtype of AbstractDirection used for parametrically typing a subtype of AbstractBoundary if that boundary is the eastern boundary in a rectangular domain.

West<:AbstractDirection

A simple subtype of AbstractDirection used for parametrically typing a subtype of AbstractBoundary if that boundary is the western boundary in a rectangular domain.

AbstractBoundary{D, FT} <: AbstractDomainElement{FT}

When running a Subzero simulation, each simulation has a rectangular domain that contains the ice pack. This domain is made up for four boundaries, one in each of the cardinal directions: North, South, East, and West. The user will pick four boundary types, each concrete subtypes of AbstractBoundary and assign each of them a direction, D and a float type FT, by typing them parametrically.

The user can pick from the currently implemented concrete subtypes of AbstractBoundary (OpenBoundary, CollisionBoundary, PeriodicBoundary, and MovingBoundary) or implement their own.

For now, each boundary wall is assumed to be a rectangle. For correct behavior, the corners of each of the boundary rectangles should overlap to make sure that no floes reach outside of the domain without touching a boundary. These rectangles should be represented as polygons with the poly field within each boundary. Furthermore, each boundary should have a val field that represents either the line y = val (for North and South walls) or x = val (for East and West walls) that denotes the line that marks the "inner-most" edge of the domain. If an ice floe passes over the boundary's val line, the floe interacts with the boundary.

 ________________
|__|____val___|__| <- North coordinates include corners
|  |          |  |
|  |          |  | <- East and west coordinates ALSO include corners
|  |          |  |

For ease of use, each boundary should have a constructor option where the user can provide an AbstractRectilinearGrid concrete type to define a boundary where val is at the edge of the grid so that the boundaries form a border right around the grid.

API

In addition to the below function, any new AbstractBoundary subtype must also implement AbstractDomainElement API functions.

• _update_boundary!(boundary::AbstractBoundary, Δt::Int)
• _periodic_compat(boundary1::AbstractBoundary, boundary2::AbstractBoundary)

The _update_boundary! function updates a boundary's position (by changing the val and poly fields) at every timestep. Currently, only MovingBoundary elements are updated, and their update depends on the length of the simulation's timestep, Δt.

The _periodic_compat function determines if a pair of functions are "periodically compatible", that is, could they form a pair on opposite edges of a domain to create a set of functioning periodic boundaries. It is called when creating a Domain.

OpenBoundary{D, FT} <: AbstractBoundary{D, FT}

A concrete subtype of AbstractBoundary that removes a floe from the simulation if any of the floe's vertices overlap with the OpenBoundary. This is meant to simulate the floe floating out of the simulation. The floe is immeditaly removed as there might not be ocean u and v values outside of the domian and thus we wouldn't know how to evolve the floe's trajectory. This boundary is direction D (AbstractDirection) of the domain and has field data of type FT <: AbstractFloat.

Fields

• poly::StaticQuadrilateral{FT}: rectangular polygon representing
  

the shape and location of the boundary on a given timestep with points of float type FT.

• val::FT: number of float type FT representing either the line y = val
  

(for North and South walls) or x = val (for East and West walls) that denotes the
line marking the inner-most edge of the domain.

Here is how to construct an OpenBoundary:

OpenBoundary(D, [FT = Float64]; grid  = nothing, x0 = nothing, xf = nothing, y0 = nothing, yf = nothing)

The user must specify which AbstractDirection the boundary is. The user also has an opportunity to specify which float type Float64 or Float32 will be used for field data. The user then must either provide an grid object (AbstractRectilinearGrid) for the domain to align with the grid edges, or provide four values (x0, xf, y0, and yf) that define the x and y-extents the user wants for the domain.

Positional arguments

• D::Type{<:AbstractDirection}: subtype of AbstractDirection used to represent
  

if a boundary is the North, South, East, or West boundary of a simulation

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• grid::AbstractRectilinearGrid: subtype of AbstractRectilinearGrid representing the grid used within a simulation
• x0::FT: minimum x-value of domain extent (region floes will be contained within)
• xf::FT: maximum x-value of domain extent (region floes will be contained within)
• y0::FT: minimum y-value of domain extent (region floes will be contained within)
• yf::FT: maximum y-value of domain extent (region floes will be contained within)

Examples

• Defining a Northern OpenBoundary with Float64 (default type) data using the grid keyword.

julia> g = RegRectilinearGrid(x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 20);

julia> OpenBoundary(North; grid = g)
OpenBoundary{North, Float64}
  ⊢polygon points are defined by the following set: (-250000.0, 750000.0), (750000.0, 500000.0), (750000.0, 750000.0), (-250000.0, 500000.0)
  ∟val is 500000.0

• Defining a Southern OpenBoundary with Float32 data using the x0, xf, y0 and yf keywords.

julia> OpenBoundary(South, Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5)
OpenBoundary{South, Float32}
  ⊢polygon points are defined by the following set: (750000.0f0, 0.0f0), (750000.0f0, -250000.0f0), (-250000.0f0, 0.0f0), (-250000.0f0, -250000.0f0)
  ∟val is 0.0

PeriodicBoundary <: AbstractBoundary

A concrete subtype of AbstractBoundary that moves a floe from one side of the domain to the opposite side of the domain (North to South and East to West and visa versa) when its centroid crosses the PeriodicBoundary, bringing the floe back into the domain. Due to this behavior, PeriodicBoundary pairs are required to form a valid domain. This boundary is direction D (AbstractDirection) of the domain and has field data of type FT <: AbstractFloat.

Fields

• poly::StaticQuadrilateral{FT}: rectangular polygon representing
  

the shape and location of the boundary on a given timestep with points of float type FT.

• val::FT: number of float type FT representing either the line y = val
  

(for North and South walls) or x = val (for East and West walls) that denotes the
line marking the inner-most edge of the domain.

Here is how to construct an PeriodicBoundary:

PeriodicBoundary(D, [FT = Float64]; grid  = nothing, x0 = nothing, xf = nothing, y0 = nothing, yf = nothing)

The user must specify which AbstractDirection the boundary is. The user also has an opportunity to specify which float type Float64 or Float32 will be used for field data. The user then must either provide an grid object (AbstractRectilinearGrid) for the domain to align with the grid edges, or provide four values (x0, xf, y0, and yf) that define the x and y-extents the user wants for the domain.

Positional arguments

• D::Type{<:AbstractDirection}: subtype of AbstractDirection used to represent
  

if a boundary is the North, South, East, or West boundary of a simulation

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• grid::AbstractRectilinearGrid: subtype of AbstractRectilinearGrid representing the grid used within a simulation
• x0::FT: minimum x-value of domain extent (region floes will be contained within)
• xf::FT: maximum x-value of domain extent (region floes will be contained within)
• y0::FT: minimum y-value of domain extent (region floes will be contained within)
• yf::FT: maximum y-value of domain extent (region floes will be contained within)

Examples

• Defining an Eastern PeriodicBoundary with Float32 data using the grid keyword.

julia> g = RegRectilinearGrid(Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 20);

julia> PeriodicBoundary(East; grid = g)
PeriodicBoundary{East, Float64}
  ⊢polygon points are defined by the following set: (750000.0, -250000.0), (500000.0, -250000.0), (750000.0, 750000.0), (500000.0, 750000.0)
  ∟val is 500000.0

• Defining a Western PeriodicBoundary with Float64 data using the x0, xf, y0 and yf keywords.

julia> PeriodicBoundary(West, Float64; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5)
PeriodicBoundary{West, Float64}
  ⊢polygon points are defined by the following set: (0.0, -250000.0), (-250000.0, 750000.0), (0.0, 750000.0), (-250000.0, -250000.0)
  ∟val is 0.0

CollisionBoundary <: AbstractBoundary

A concrete subtype of AbstractBoundary that calculates collision forces of a floe against the boundary if any of the floe's vertices overlap with the CollisionBoundary. This is meant to simulate any barrier that might stop a floe from flowing into or out of a given region, like the edges of a cove. With this type of wall, a CollisionBoundary is treated as an immovable, unbreakable floe for the purposes of calculations. This boundary is direction D (AbstractDirection) of the domain and has field data of type FT <: AbstractFloat.

Fields

• poly::StaticQuadrilateral{FT}: rectangular polygon representing
  

the shape and location of the boundary on a given timestep with points of float type FT.

• val::FT: number of float type FT representing either the line y = val
  

(for North and South walls) or x = val (for East and West walls) that denotes the
line marking the inner-most edge of the domain.

Here is how to construct an CollisionBoundary:

CollisionBoundary(D, [FT = Float64]; grid  = nothing, x0 = nothing, xf = nothing, y0 = nothing, yf = nothing)

The user must specify which AbstractDirection the boundary is. The user also has an opportunity to specify which float type Float64 or Float32 will be used for field data. The user then must either provide an grid object (AbstractRectilinearGrid) for the domain to align with the grid edges, or provide four values (x0, xf, y0, and yf) that define the x and y-extents the user wants for the domain.

Positional arguments

• D::Type{<:AbstractDirection}: subtype of AbstractDirection used to represent
  

if a boundary is the North, South, East, or West boundary of a simulation

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• grid::AbstractRectilinearGrid: subtype of AbstractRectilinearGrid representing the grid used within a simulation
• x0::FT: minimum x-value of domain extent (region floes will be contained within)
• xf::FT: maximum x-value of domain extent (region floes will be contained within)
• y0::FT: minimum y-value of domain extent (region floes will be contained within)
• yf::FT: maximum y-value of domain extent (region floes will be contained within)

Examples

• Defining an Northern CollisionBoundary with Float64 data using the grid keyword.

julia> g = RegRectilinearGrid(Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 20);

julia> CollisionBoundary(North; grid = g)
CollisionBoundary{North, Float64}
  ⊢polygon points are defined by the following set: (-250000.0, 750000.0), (750000.0, 500000.0), (750000.0, 750000.0), (-250000.0, 500000.0)
  ∟val is 500000.0

• Defining a Western CollisionBoundary with Float64 data using the x0, xf, y0 and yf keywords.

julia> CollisionBoundary(West, Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5)
CollisionBoundary{West, Float32}
  ⊢polygon points are defined by the following set: (0.0f0, -250000.0f0), (-250000.0f0, 750000.0f0), (0.0f0, 750000.0f0), (-250000.0f0, -250000.0f0)
  ∟val is 0.0

MovingBoundary <: AbstractBoundary

A concrete subtype of AbstractBoundary that can provide a compressive or shear force on the floes within the domain by moving parallel or perpendicular to the domain. The MovingBoundary calcuates the forces on the floe exactly the same as a CollisionBoundary, by acting as a immovable, unbreakable floe. This boundary is direction D (AbstractDirection) of the domain and has field data of type FT <: AbstractFloat.

If a North or South wall has a non-zero v velocity this will provide a compressive stress as the floe moves towards the center of the domain along the vector x = cx where (cx, cy) is the centroid of the domain (or decompresses with opposite sign velocities). The poly and val fields are updated to represent the movement of the domain at the user-provided velocities, u and v m/s.

Alternatively, if a North or South wall has a non-zero u velocity this will provide a shear stress as the domain "moves" along the line y = y0 (for North) or y = yf (for South). In this case, this only changes the frictional forces and we do not need up actually change the poly or val fields.

Fields

• poly::StaticQuadrilateral{FT}: rectangular polygon representing
  

the shape and location of the boundary on a given timestep with points of float type FT.

• val::FT: number of float type FT representing either the line y = val
  

(for North and South walls) or x = val (for East and West walls) that denotes the
line marking the inner-most edge of the domain.

• u::FT: boundary's u-velocity
• v::FT: boundary's v-velocity

Here is how to construct an MovingBoundary:

MovingBoundary(D, [FT = Float64]; u = 0.0, v = 0.0, grid  = nothing, x0 = nothing, xf = nothing, y0 = nothing, yf = nothing)

The user must specify which AbstractDirection the boundary is. The user also has an opportunity to specify which float type Float64 or Float32 will be used for field data. The user then must either provide an grid object (AbstractRectilinearGrid) for the domain to align with the grid edges, or provide four values (x0, xf, y0, and yf) that define the x and y-extents the user wants for the domain.

Positional arguments

• D::Type{<:AbstractDirection}: subtype of AbstractDirection used to represent
  

if a boundary is the North, South, East, or West boundary of a simulation

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• u::FT: boundary's u-velocity
• v::FT: boundary's v-velocity
• grid::AbstractRectilinearGrid: subtype of AbstractRectilinearGrid representing the grid used within a simulation
• x0::FT: minimum x-value of domain extent (region floes will be contained within)
• xf::FT: maximum x-value of domain extent (region floes will be contained within)
• y0::FT: minimum y-value of domain extent (region floes will be contained within)
• yf::FT: maximum y-value of domain extent (region floes will be contained within)

Examples

• Defining an Northern MovingBoundary with Float64 data using the grid keyword. Assigning u-velocity of 0.5m/s.

julia> g = RegRectilinearGrid(Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 20);

julia> MovingBoundary(North; u = 0.5, grid = g)
MovingBoundary{North, Float64}
  ⊢polygon points are defined by the following set: (-250000.0, 750000.0), (750000.0, 500000.0), (750000.0, 750000.0), (-250000.0, 500000.0)
  ⊢val is 500000.0
  ⊢u-velocity of 0.5 m/s
  ∟v-velocity of 0.0 m/s

• Defining a Southern MovingBoundary with Float32 data using the x0, xf, y0 and yf keywords.

Assigning u-velocity of 0.3 m/s and v-velocity of 0.25 m/s

julia> MovingBoundary(South, Float32; u = 0.3, v = 0.25, x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5)
MovingBoundary{South, Float32}
  ⊢polygon points are defined by the following set: (750000.0f0, 0.0f0), (750000.0f0, -250000.0f0), (-250000.0f0, 0.0f0), (-250000.0f0, -250000.0f0)
  ⊢val is 0.0
  ⊢u-velocity of 0.3 m/s
  ∟v-velocity of 0.25 m/s

TopographyElement{FT}<:AbstractDomainElement{FT}

A concrete subtype of AbstractDomainElement that represent topography in a given simulation. These TopographyElement objects will be treated as be treated as islands or coastline within the model in that they will not move or break due to floe interactions, but they will affect floes that collide with them.

Similar to Floe objects, we store the shape of each TopographyElement with a poly field. We also store the centroid and the shape's maximum radius (rmax) in order to calculate simple interaction checks with Floes.

Fields

• poly::Polys{FT}: Polygon used to represent the shape of a floe or topography
• centroid::Vector{FT}: Two-element vector meant to represent the (x, y) point that is the centroid of either a floe or topography
• rmax::FT: Float length representing the maximum radius of a floe or topography from the centroid to any given vertex

Here is how to construct an TopographyElement:

TopographyElement([FT = Float64]; poly::Polys)

The user must provide a GeoInterface polygon of the specific type defined by Polys. The user also has an opportunity to specify which float type Float64 or Float32 will be used for field data.

Positional arguments

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• poly::Polys{FT}: Polygon used to represent the shape of a floe or topography

Note: The user should NOT be using this constructor. The user should create a topography

field using initialize_topography_field rather than creating individual topography elements. This allows the abstraction of the use of the StructArrays package away from the user.

Examples

julia> import GeoInterface as GI;

julia> poly = GI.Polygon([[(0.0, 0.0), (0.0, 1e3), (1e3, 1e3), (1e3, 0.0), (0.0, 0.0)]]);

julia> TopographyElement(Float64; poly)
TopographyElement{Float64}
  ⊢centroid is (500.0, 500.0) in meters
  ∟maximum radius is 707.1067811865476 meters

julia> TopographyElement(Float32; poly)
TopographyElement{Float32}
  ⊢centroid is (500.0, 500.0) in meters
  ∟maximum radius is 707.1068 meters

initialize_topography_field([::Type{FT} = Float64]; polys = nothing, coords = nothing)

This function allows for easy initialization of a field of [TopographyElement(@ref)s and collects them into aStructArrayso that they can be passed into a [Model`](@ref).

This is the suggested way to create a topography field for a simulation. Do NOT construct individual TopographyElement objects as that method does not correct the field as a whole to ensure no topography polygons overlap and does not create a struct array that can be passed to a Model.

The user can create a topography field by either passing a list of polygons or by passing a list of coordiantes, which will then be made into polygons.

Positional arguments

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• polys::Vector{<:Polygon}: list of polygons meant to represent a field of floes or topography elements. Polygons can be any polygon type that supports GeoInterface.
• coords::Vector{<:PolyVec}: list of polygon coordinates meant to represent a field of floes or topography elements. PolyVec refers to a Vector{Vector{<:Points}} where the points can be tuples, vectors, or static vectors and the innermost vector refers to each ring of the polygon.

Examples

• Defining a topography field with coordinates

julia> coords = [[[(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (0.0, 0.0)]], [[(10.0, 0.0), (10.0, 1.0), (11.0, 1.0), (10.0, 0.0)]]];

julia> initialize_topography_field(Float64; coords)
2-element TopographyField{Float64} list:
 TopographyElement{Float64}
  ⊢centroid is (0.3333333333333333, 0.6666666666666666) in meters
  ∟maximum radius is 0.74535599249993 meters
 TopographyElement{Float64}
  ⊢centroid is (10.333333333333334, 0.6666666666666666) in meters
  ∟maximum radius is 0.7453559924999301 meters

• Defining a topography field with polygons

julia> import GeoInterface as GI;

julia> polys = [GI.Polygon(c) for c in coords];

julia> initialize_topography_field(Float32; polys)
2-element TopographyField{Float32} list:
 TopographyElement{Float32}
  ⊢centroid is (0.33333334, 0.6666667) in meters
  ∟maximum radius is 0.745356 meters
 TopographyElement{Float32}
  ⊢centroid is (10.333333, 0.6666667) in meters
  ∟maximum radius is 0.74535626 meters

• Creating an empty topography field without polys or coords

julia> initialize_topography_field(Float64)
0-element TopographyField{Float64} list

Domain{FT, NB, SB, EB, WB, TT}

A simulation Domain holds four (4) boundary elements (concrete subtypes of AbstractBoundary) and a list (potentially empty) of TopographyElement objects. There must be a wall typed by each cardinal direction, and each wall must hold data of the same float type FT. Additionally, the North boundary should be "higher" (on the y-axis) than the South boundary and the East should be futher to the "right" (on the x-axis) than the West boundary. Aditionally, the boundary walls should be overlapping as detailed in AbstractBoundary.

Additionally, the set of walls must be periodically compatible. This means that pairs of opposite boundaries (North and South AND East and West) both need to be periodic if one of them is periodic. This is because if a floe exits a periodic boundary, it must be able to re-enter the opposite boundary to fulfill its definition of periodic.

Fields

• north::NB: Northern boundary where NB <: AbstractBoundary{North, FT}
• south::SB: Southern boundary where SB <: AbstractBoundary{South, FT}
• east::EB: Eastern boundary where EB <: AbstractBoundary{East, FT}
• west::WB: Western boundary where WB <: AbstractBoundary{West, FT}
• topography::tT: Field of topography elements where TT <: TopographyField{FT}

Here is how to construct an MovingBoundary:

Domain(; north::NB, south::SB, east::EB, west::WB, topography::TT = nothing)

The user must provide the four (4) boundaries. They then have an option to provide topography. If no topography is provided, an empty TopographyField will be created with the same FT as the first of the boundaries (which should be shared amoung boundaries).

Keyword arguments

• north::NB: domain's northern boundary
• south::SB: domain's southern boundary
• east::EB: domain's eastern boundary
• west::WB: domain's western boundary
• topography::TT: domain's topography field, if there is one, else an empty field is created

Examples

• Creating a Domain with NO topography

julia> grid = RegRectilinearGrid(Float32; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 20);

julia> north = OpenBoundary(North, Float64; grid);

julia> south = OpenBoundary(South, Float64; grid);

julia> east = PeriodicBoundary(East, Float64; grid);

julia> west = PeriodicBoundary(West, Float64; grid);

julia> Domain(; north, south, east, west)
Domain
  ⊢Northern boundary of type OpenBoundary{North, Float64}
  ⊢Southern boundary of type OpenBoundary{South, Float64}
  ⊢Eastern boundary of type PeriodicBoundary{East, Float64}
  ⊢Western boundary of type PeriodicBoundary{West, Float64}
  ∟0-element TopograpahyElement{Float64} list

• Creating a Domain with topography

julia> import GeoInterface as GI;

julia> topo_polys = [GI.Polygon([[(1e4, 1e4), (1e4, 3e4), (3e4, 1e4), (1e4, 1e4)]]), GI.Polygon([[(8e4, 8e4), (8e4, 9e4), (9e4, 9e4), (9e4, 8e4), (8e4, 8e4)]])];

julia> topography = initialize_topography_field(; polys = topo_polys);

julia> Domain(; north, south, east, west, topography)
Domain
  ⊢Northern boundary of type OpenBoundary{North, Float64}
  ⊢Southern boundary of type OpenBoundary{South, Float64}
  ⊢Eastern boundary of type PeriodicBoundary{East, Float64}
  ⊢Western boundary of type PeriodicBoundary{West, Float64}
  ∟2-element TopograpahyElement{Float64} list

Ocean{FT}

Each simulation needs ocean values on the grid to drive the ice floes. The Ocean struct holds the needed fields to perform coupling, including vector fields (like u-velocity, v-velocity, x-stress, and y-stress) and tracer fields (like temperature, heatflux, and dissolved ice mass). We also hold and calculate the sea ice fraction in any given grid cell.

Right now, all values are stored on grid points as defined by an AbstractRectilinearGrid. This should eventually change to a c-grid. If interested in this change, see issue #30. This will require a big change in behavior, but is an important change for consistency and to match with Oceananigans for two-way coupling.

Coupling beyond velocity (for example thermodynamically) is not very well developed. If you are interested in working on this see issue #9. The hflx_factor field is related to thermodynamic coupling.

Fields

• u::Matrix{FT}: ocean u-velocities in the x-direction on each grid point
• v::Matrix{FT}: ocean v-velocities in the y-direction on each grid point
• temp::Matrix{FT}: ocean temperature on each grid point
• hflx_factor::Matrix{FT}: factor to calculate the ocean-atmosphere heat flux for a grid cell
• scells::Matrix{CellStresses}: used to accumulate stresses on the ocean per grid cell from floes in cell (see CellStresses)
• τx::Matrix{FT}: stress on the ocean in the x-direction on each grid point
• τy::Matrix{FT}: stress on the ocean in the y-direction on each grid point
• si_frac::Matrix{FT}: fraction of area in each grid cell (centered on grid points) that is covered in sea-ice (between 0-1)
• dissolved::Matrix{FT}: total mass from ice floes dissolved in each grid cell (centered on grid points)

Here is how to construct an Ocean:

Ocean([FT = Float64]; u, v, temp, grid = nothing)

Positional arguments

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• u::Union{Real, Matrix{Real}}: user can either provide a single u-velocity value (a constant field the size of the grid will be created) or provide a matrix field of u-velocity values
• v::Union{Real, Matrix{Real}}: user can either provide a single v-velocity value (a constant field the size of the grid will be created) or provide a matrix field of v-velocity values
• temp::Union{Real, Matrix{Real}}: user can either provide a single temperature value (a constant field the size of the grid will be created) or provide a matrix field of temperature values
• grid::AbstractRectilinearGrid: subtype of AbstractRectilinearGrid representing the grid used within a simulation

Examples

• Creating Ocean with constant u-velocity, v-velocity, and temp

julia> grid = RegRectilinearGrid(; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 10);

julia> Ocean(; u = 0.5, v = 0.25, temp = 0.0, grid)
Ocean{Float64}
  ⊢Vector fields of dimension (21, 11)
  ⊢Tracer fields of dimension (21, 11)
  ⊢Average u-velocity of: 0.5 m/s
  ⊢Average v-velocity of: 0.25 m/s
  ∟Average temperature of: 0.0 C

• Creating Ocean with user-defined u and v and a constant temp with Float32 data

julia> import Random.Xoshiro;

julia> seeded_rng = Xoshiro(100);

julia> u_field = rand(seeded_rng, grid.Nx + 1, grid.Ny + 1);

julia> v_field = rand(seeded_rng, grid.Nx + 1, grid.Ny + 1);

julia> Ocean(Float32; u = u_field, v = v_field, temp = 0.0, grid)
Ocean{Float32}
  ⊢Vector fields of dimension (21, 11)
  ⊢Tracer fields of dimension (21, 11)
  ⊢Average u-velocity of: 0.4856711 m/s
  ⊢Average v-velocity of: 0.5319979 m/s
  ∟Average temperature of: 0.0 C

• Trying to create an Ocean with a constant field and NO grid

julia> Ocean(; u = 0.2, v = 0.1, temp = 0.0)
ERROR: ArgumentError: To create a matrix from the constant value provided, you must provide a grid.
[...]

Atmos{FT}

Each simulation needs atmosphere/wind values on the grid to drive the ice floes. The Atmos struct holds the needed fields to perform coupling of vector fields (u-velocity and v-velocity) and tracer fields (temperature).

Right now, all values are stored on grid points as defined by an AbstractRectilinearGrid. This should eventually change to a c-grid. If interested in this change, see issue #30.

Coupling beyond velocity (for example thermodynamically) is not very well developed. If you are interested in working on this see issue #9. The hflx_factor field is related to thermodynamic coupling.

Fields

• u::Matrix{FT}: atmosphere/wind u-velocities in the x-direction on each grid point
• v::Matrix{FT}: atmosphere/wind v-velocities in the y-direction on each grid point
• temp::Matrix{FT}: atmosphere/wind temperature on each grid point

Here is how to construct an Atmos:

Atmos([FT = Float64]; u, v, temp, grid = nothing)

Positional arguments

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• u::Union{Real, Matrix{Real}}: user can either provide a single u-velocity value (a constant field the size of the grid will be created) or provide a matrix field of u-velocity values
• v::Union{Real, Matrix{Real}}: user can either provide a single v-velocity value (a constant field the size of the grid will be created) or provide a matrix field of v-velocity values
• temp::Union{Real, Matrix{Real}}: user can either provide a single temperature value (a constant field the size of the grid will be created) or provide a matrix field of temperature values
• grid::AbstractRectilinearGrid: subtype of AbstractRectilinearGrid representing the grid used within a simulation

Examples

• Creating Atmos with constant u-velocity, v-velocity, and temp

julia> grid = RegRectilinearGrid(; x0 = 0.0, xf = 5e5, y0 = 0.0, yf = 5e5, Nx = 20, Ny = 10);

julia> Atmos(; u = 0.5, v = 0.25, temp = 0.0, grid)
Atmos{Float64}
  ⊢Vector fields of dimension (21, 11)
  ⊢Tracer fields of dimension (21, 11)
  ⊢Average u-velocity of: 0.5 m/s
  ⊢Average v-velocity of: 0.25 m/s
  ∟Average temperature of: 0.0 C

• Creating Atmos with user-defined u and v and a constant temp with Float32 data

julia> import Random.Xoshiro;

julia> seeded_rng = Xoshiro(200);

julia> u_field = rand(seeded_rng, grid.Nx + 1, grid.Ny + 1);

julia> v_field = rand(seeded_rng, grid.Nx + 1, grid.Ny + 1);

julia> Atmos(Float32; u = u_field, v = v_field, temp = 0.0, grid)
Atmos{Float32}
  ⊢Vector fields of dimension (21, 11)
  ⊢Tracer fields of dimension (21, 11)
  ⊢Average u-velocity of: 0.4965465 m/s
  ⊢Average v-velocity of: 0.5340565 m/s
  ∟Average temperature of: 0.0 C

• Trying to create an Atmos with a constant field and NO grid

julia> Atmos(; u = 0.2, v = 0.1, temp = 0.0)
ERROR: ArgumentError: To create a matrix from the constant value provided, you must provide a grid.
[...]

Model which holds grid, ocean, atmos structs, each with the same underlying float type (either Float32 of Float64) and size. It also holds the domain information, which includes the topography and the boundaries. It holds a StructArray of floe structs, again each relying on the same underlying float type. Finally, it also holds the maximum floe id used thus far in the simulation. This should be the length of the floes array at the beginning of the run.

CellFloes{FT}

A CellFloes struct represents a single grid cell and accumulates the indices of floes with area in that grid cell. This is used to accumulate the forces of floes in a grid cell on the ocean below it. Due to the prevalence of periodic boundaries, the grid cell represented by a CellFloes object are centered on grid points (translated by Δx/2 and Δy/2 in the x and y directions for a RegRectilinearGrid), rather than on the grid cells defined by the grid object itself. Floes are recorded with their index in the list of floes. Furthermore, in a model with periodic boundaries, some floes may be in multiple grid cells on different edges of the domain if they pass through a periodic boundary. In these cases, the floe "linked" to with its index must be translated by a vector to get its "ghost" on the other side of the domain. This non-integer translation data is of float type FT. CellFloes are used with CellStresses objects, which aggregate the stress from each of the floes recorded in a given CellFloes object.

Fields

• floeidx: see keyword arguments
• Δx: see keyword arguments
• Δy: see keyword arguments

Here is how to construct a CellFloes object:

CellFloes([FT = Float64]; floeidx = nothing, Δx = nothing, Δy = nothing)

Positional arguments

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• floeidx::Vector{Int}: vector of floe indicies in the list of model floes for floes with area in the grid
• Δx::Vector{FT}: vector of x-translations for a floe at the corresponding index of the

floeidx vector to be in the cell represented with the CellFloes object.

• Δy::Vector{FT}: vector of y-translations for a floe at the corresponding index of the

floeidx vector to be in the cell represented with the CellFloes object.

CellStresses{FT}

Struct to collect stress from ice floes on ocean grid cells. One CellStresses corresponds to one grid cell and is "linked" to a corresponding CellFloes object. The CellFloes object records which floes are in a given grid cell and the CellStresses aggregates the total stress felt on the ocean from the ice (on each of the floe's sub-floe points generated by a subtype of AbstractSubFloePointsGenerator) in a given grid cell. This is needed ONLY for two-way coupling. Due to the prevalence of periodic boundaries, the grid cell represented by a CellStresses object are centered on grid points (translated by Δx/2 and Δy/2 in the x and y directions for a RegRectilinearGrid), rather than on the grid cells defined by the grid object itself. A CellStresses object holds a list of stresses on the ocean where each element is from a single floe within grid cell. The ith element of the τx and τy fields are measures of the total stress from the ith floe in the cell and the npoints field is a measure of how many of the ith floe's sub-floe points contributed to the running total of the stress. Again, each element in the list corresponds to one floe, which is denoted in the corresponding CellFloes matrix within the grid.

Fields

• τx::Vector{FT}: list of total x-stress caused by each floe in a cell on the ocean
• τy::Vector{FT}: list of total y-stress caused by each floe in a cell on the ocean
• npoints::Vector{Int}: list of total number of sub-floe points contributing to each individial floe's stress

Here is how to construct a CellStresses object:

CellStresses([FT = Float64]; τx = nothing, τy = nothing, npoints = nothing)

Positional arguments

• FT::Type{<:AbstractFloat}: Float type used to run the simulation, either Float64 (default) or Float32.

Keyword arguments

• τx::Vector{FT}: list of total x-stress caused by each floe in a cell on the ocean
• τy::Vector{FT}: list of total y-stress caused by each floe in a cell on the ocean
• npoints::Vector{Int}: list of total number of sub-floe points contributing to each individial floe's stress

TopographyField{FT}

Alias for StructArray type with TopographyElement elements with data of type FT.