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.10678 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.33333, 0.66667) in meters
  ∟maximum radius is 0.74536 meters
 TopographyElement{Float64}
  ⊢centroid is (10.33333, 0.66667) in meters
  ∟maximum radius is 0.74536 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.33333, 0.66667) in meters
  ∟maximum radius is 0.74536 meters
 TopographyElement{Float32}
  ⊢centroid is (10.33333, 0.66667) in meters
  ∟maximum radius is 0.74536 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.48567 m/s
  ⊢Average v-velocity of: 0.532 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.49655 m/s
  ⊢Average v-velocity of: 0.53406 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.
[...]

Floe{FT}

Each model has a list of floes that are advected and updated throughout the simulation. The Floe struct holds the fields needed to compute all floe calculations. However, to create a floe, only 2 arguments are neccessary and there are 2 more optional keyword arguments - all other fields will be set to good starting values for a simualation.

For those developers, here is how to construct a Floe:

Floe{FT}(shape, height; floe_settings = FloeSettings(), rng = Xoshiro(), kwargs...)

The FT type should be provided as shown in the examples below:

Positional arguments

• shape::Union{<:Polys, <:PolyVec}: either a polygon of type Polys or a vector of coordiantes of type PolyVec used to represent the shape of a floe.
• height:FT: height of the new floe to be created

Keyword arguments

• floe_settings::FloeSettings: floe settings for a simulation (see FloeSettings)
• rng:<AbstractRNG: random number generator to use to generate the subfloe points (if needed)
• kwargs...: additional field values can be passed in if desired

Examples

• Creating a Floe with coordiantes:

julia> coords = [[[0.0, 0.0], [0.0, 5.0], [10.0, 5.0], [10.0, 0.0], [0.0, 0.0]]];

julia> height = 1.0;

julia> Floe{Float64}(coords, height)
Floe{Float64}
  ⊢Centroid of (5.0, 2.5) m
  ⊢Height of 1.0 m
  ⊢Area of 50.0 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

• Creating a Floe with a polygon:

julia> poly = make_polygon(coords);

julia> Floe{Float32}(poly, height; u = 1.0, ξ = 0.02)
Floe{Float32}
  ⊢Centroid of (5.0, 2.5) m
  ⊢Height of 1.0 m
  ⊢Area of 50.0 m^2
  ∟Velocity of (u, v, ξ) of (1.0, 0.0, 0.02) in (m/s, m/s, rad/s)

Fields

The first catagory is physical properties. These have to do with the floe's physical shape.

Before listing the fields, one important thing to know is that a floe's coordinates are represented by a PolyVec, which is a shorthand for a vector of a vector of a vector of floats. This sounds complicated, but it is simply a way of representing a polygon's coordinates. A Polygon's coordinates are of the form below, where the xy-coordinates are the exterior border of the floe and the wz-coordinates, or any other following sets of coordinates, describe holes within the floe:

julia coords = [ [[x1, y1], [x2, y2], ..., [xn, yn], [x1, y1]], # Exterior vertices of the polygon represented as a list of cartesian points [[w1, z1], [w2, z2], ..., [wn, zn], [w1, z1]], # Interior holes of the polygon represented as a list of cartesian points ..., # Additional holes within the polygon represented as a list of cartesian points ] We will use the term PolyVec to describe this form of coordiantes and you will see it in the code if you take a look at the code base. It is also the form that floe coordinates are saved in output files.

The second catagory is sub-floe points. These are used for interpolation of the ocean and atmosphere onto the floe. They are a list of points within the floe. There are two ways they can be generated. One is for them are randomly generated, with the user providing an initial target number of points. These are generated with a monte carlo generator (explained below). The other way is for the points to be on a sub-grid within the floe. There are benefits and drawbacks to both strategies.

The third catagory is velocities and orientations. Floe's have both linear and angular velocity and keep track of the angle that they have rotated since the begining of the simulation. | Movement Fields| Meaning | Type | | ––––––– | –––––––––––––––- | ––––––––– | | u | floe's x-velocity in [m/s] | Float64 or Float32 | | v | floe's x-velocity in [m/s] | Float64 or Float32 | | ξ | floe's angular velocity in [rad/s]| Float64 or Float32 | | α | rotation from starting position<br> in [rad]| Float64 or Float32 |

The fourth catagory is status. These fields hold logistical information about each floe and through which process it originated. | Status Fields | Meaning | Type | | ––––––– | –––––––––––––––- | ––––––––– | | status | if the floe is still active in the simulation | Subzero.Status (see below)| | id | unique floe id for tracking the floe throughout the simulation | Int | | ghostid | if floe is not a ghost, `ghostid = 0, else it is in[1, 4]<br> as each floe can have up to 4 ghosts| Int | | parent_id | if floe is created from a fracture or the fusion of two floes,parent_idis a list of <br> the original floes'id`, else it is emtpy | Vector of Ints | | ghosts | indices of floe's ghost floes within the floe list| Vector of Ints |

A Status object has two fields: a tag and a fuse index list. There are currently three different tags: active, remove, and fuse. If a floe is active, it will continue in the simulation at the end of a timestep. If a floe's tag is remove, it will be removed at the end of the timestep. This ususally happens if a floe exits the domain, or becomes unstable for some reason. If a floe is marked at fuse, this means that is is overlapping with another floe by more than the user defined maximum overlap percent (see CollisionSettings for more information on this maximum overlap value. If a floe is marked for fusion, the index of the floe it is supposed to fuse with will be listed in the fuse_idx list.

The fifth catagory is forces and collisions. These fields hold information about the forces on each floe and the collisions it has been in. | Force Fields | Meaning | Type | | ––––––––- | ––––––––––––––––––––– | –––––––––––––- | | fxOA | x-force on floe from ocean and atmosphere in [N] | Float64 or Float32| | fyOA | y-force on floe from ocean and atmosphere in [N] | Float64 or Float32| | trqOA | torque on floe from ocean and atmosphere in [N m]| Float64 or Float32| | hflxfactor | coefficent of floe height to get heat flux directly <br> under floe in [W/m^3]| Float64 or Float32| | overarea | total overlap of floe from collisions in [m^2] | Float64 or Float32| | collisionforce | forces on floe from collisions in [N] | Float64 or Float32| | collisiontrq | torque on floe from collisions in [N m] | Float64 or Float32| | interactions | each row holds one collision's information, see below for more information | nx7 Matrix of Float64 or Float32 <br> where n is the number of collisions| | stressaccum | stress accumulated over the floe over past timesteps given StressCalculator, where it is of the form [xx yx; xy yy] | 2x2 Matrix{AbstractFloat}| | stress_instant | instantaneous stress on floe in current timestep | 2x2 Matrix{AbstractFloat} | strain | strain on floe where it is of the form [ux vx; uy vy] | 2x2 Matrix of Float64 or Float32| | damage | damage a floe assumulates through interactions (currently unused; see DamageStressCalculator)| Float64 or Float32|

The interactions field is a matrix where every row is a different collision that the floe has experienced in the given timestep. There are then seven columns, which are as follows:

• floeidx, which is the index of the floe that the current floe collided with
• xforce, which is the force in the x-direction caused by the collision
• yforce, which is the force in the y-direction caused by the collision
• xpoint, which is the x-coordinate of the collision point, which is the x-centroid of the overlap between the floes
• ypoint, which is the y-coordinate of the collision point, which is the y-centroid of the overlap between the floes
• torque, which is the torque caused by the collision
• overlap, which is the overlap area between the two floes in the collision.

You can use these column names to access columns of interactions. For example: floe.interactions[:, xforce] gives the list of x-force values for all collisions a given floe was involved in as that timestep.

The fifth catagory is timesteping values. These are values from the previous timestep that are needed for the current timestepping scheme. | Previous Value Fields | Meaning | Type | | ––––––––––- | ––––––––––––––––––––––- | –––––––––-| | pdxdt | previous timestep x-velocity (u) in [m/s] | Float64 or Float32| | pdydt | previous timestep y-velocity (v) in [m/s] | Float64 or Float32| | pdudt | previous timestep x-acceleration in [m/s^2] | Float64 or Float32| | pdvdt | previous timestep y-acceleration in [m/s^2] | Float64 or Float32| | pdαdt | previous timestep angular-velocity in [rad/s] | Float64 or Float32| | pdξdt | previous timestep time angular acceleration in [rad/s^2] | Float64 or Float32|

YourFloeFieldGenerator{FT} <: AbstractFloeFieldGenerator{FT}

Each simulation run with requires a field of ice floes, which is simply a list of Floe structs. However, there are various ways to create a floe field, and various configurations of floes that user might want to create.

To make it easier for the user to define a new floe creation functionality, this process uses multiple dispatch on subtypes of the AbstractFloeFieldGenerator type. The generators hold the information needed to generate a floe field - and the method is determined by the generator type. Right now, there are two subtypes implemented.

The first is the CoordinateListFieldGenerator. This generator takes in a list of floe coordinates (along with other arguments) to create a list of floes. The second is the VoronoiTesselationFieldGenerator. This generator tesselates the domain with floes and then randomly keeps various floes to match the user-provided floe concentration.

Users and developers can create new generators to easily implement new floe generators. For example, it would be good to have a generator that creates floes with a user-defined FSD.

API

The following methods must be implemented for all subtypes:

• _initialize_floe_field!(::Type{FT}, floe_arr, generator::CoordinateListFieldGenerator, domain; kwargs...):

CoordinateListFieldGenerator{FT} <: AbstractFloeFieldGenerator{FT}

A concrete implementatin of AbstractFloeFieldGenerator that generates a floe field from a list of floe coordinates. The floe coordinates are turned into polygons. The heights are a random distribution around hmean between hmean - Δh and hmean + Δh.

The fields are:

• polys::Vector{Polys{FT}}: list of polygons in the form of PolyVecs.
• hmean::FT: mean height of the generated floe field.
• Δh::FT: range of floe heights around the mean (i.e. random distribution around hmean between hmean - Δh and hmean + Δh)

Here is how to construct a CoordinateListFieldGenerator:

CoordinateListFieldGenerator([FT = Float64]; coords, hmean, Δh)

Positional arguments

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

Keyword arguments

• coords::Vector{PolyVec{FT}}: list of polygon coordinates in the form of PolyVecs
• hmean::FT: mean height of the generated floe field.
• Δh::FT: range of floe heights around the mean (i.e. random distribution around hmean between hmean - Δh and hmean + Δh)
• kwargs...

Examples

julia> coords = [[[(0.0, 0.0), (0.0, 100.0), (100.0, 100.0), (0.0, 0.0)]], [[(150.0, 0.0), (150.0, 100.0), (200.0, 100.0), (150.0, 0.0)]]];

julia> generator = CoordinateListFieldGenerator(Float64; coords, hmean = 1, Δh = 0.25)
CoordinateListFieldGenerator:
  ⊢Number of polygons: 2
  ⊢Mean height: 1.0
  ∟Height range: 0.25

VoronoiTesselationFieldGenerator{FT} <: AbstractFloeFieldGenerator{FT}

A concrete implementatin of AbstractFloeFieldGenerator that generates a floe field using voronoi tesselation according to user input number of floes and concentrations. The heights are a random distribution around hmean between hmean - Δh and hmean + Δh.

The fields are:

• nfloes::Int:number of floes to try to create - note you might not end up with this number of floes. Topography in domain and multiple concentrations can decrease number of floes created
• concentrations::Matrix{FT}: matrix of concentrations to fill domain. If size(concentrations) = N, M then split the domain into NxM cells, each to be filled with the corresponding concentration. If concentration is below 0, it will default to 0. If it is above 1, itwill default to 1
• hmean::FT: mean height of the generated floe field.
• Δh::FT: range of floe heights around the mean (i.e. random distribution around hmean between hmean - Δh and hmean + Δh)

Here is how to construct a VoronoiTesselationFieldGenerator:

VoronoiTesselationFieldGenerator([FT = Float64]; nfloes, concentrations, hmean, Δh)

Positional arguments

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

Keyword arguments

• nfloes::Int:number of floes to try to create (see above for more info)
• concentrations::Matrix{FT}: matrix of concentrations to fill domain (see above for more info)
• hmean::FT: mean height of the generated floe field.
• Δh::FT: range of floe heights around the mean (i.e. random distribution around hmean between hmean - Δh and hmean + Δh)
• floe_bounds
• kwargs...

Examples

julia> nfloes = 20;

julia> concentrations = [0.5];

julia> generator = VoronoiTesselationFieldGenerator(Float64; nfloes, concentrations, hmean = 1, Δh = 0.25)
VoronoiTesselationFieldGenerator:
  ⊢Requested number of floes: 20
  ⊢Requested concentrations [0.5;;]
  ⊢Mean height: 1.0
  ∟Height range: 0.25

initialize_floe_field(FT; generator, domain, supress_warnings, floe_settings, rng, kwargs...)

Create a field of floes using the provided generator. Regardless of generator, the floe array will be a vector of Floe structs, but the exact shape/location/number of the floes will be determined by the generator, which must be a subtype of the AbstractFloeFieldGenerator. After floe creation, the unique, numerical ids of the floes within the floe field will be set.

Positional arguments

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

Keyword arguments

• generator::AbstractFloeFieldGenerator: generator type that determines how the floe field is generated
• domain::Domain: simulation domain
• floe_settings::FloeSettings: floe settings that determine individual floe characteristics - also needed to create a Simulation object
• supress_warnings::Bool: boolean flag on if warnings regarding floe area and centroid position should checked
• rng::RandomNumberGenerator:: random number generator needed to randomly assign floe characteristics (e.g height), depending on the generator
• floe_bounds::Polys: bounding box for floes to be generated within - only used for VoronoiTesselationFieldGenerator
• kwargs...: other keywords night be needed if new generator types are implemented.

Model{FT, GT, DT, FLT}

A simulation Model holds a concrete subtype of AbstractRectilinearGrid, a Domain, an Ocean, an Atmos, and a list of Floe objects. These are all of the physical pieces of the simulation.

To sucessfully create a model, the domain must fit within the grid, the ocean and atmosphere must be the same dimensions (this could be changed later once the code can handle different sized grids), and all of the model elements must be of the same float type FT, either Float64 or Float32.

Fields

• grid::GT: model's grid where GT <: AbstractRectilinearGrid{FT}
• ocean::Ocean{FT}: model's ocean
• atmos::Atmos{FT}: model's atmosphere
• domain::DT: model's domain where DT <: Domain{FT, NB, SB, EB, WB, TT}
• floes::FLT: models' list of floes where FLT <: <:StructArray{<:Floe{FT}}

Here is how to construct a Model:

Model(; grid::GT, ocean::Ocean{FT}, atmos::Atmos{FT}, domain::DT, floes::FLT)

Keyword arguments

• grid::GT: model's grid
• ocean::Ocean{FT}: model's ocean
• atmos::Atmos{FT}: model's atmosphere
• domain::DT: model's domain
• floes::FLT: models' list of floes

Examples

• Creating a Model

```jldoctest model

julia> using Random

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

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

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

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

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

julia> domain = Domain(; north, south, east, west);

julia> ocean = Ocean(Float64; u = 0.5, v = 0.25, temp = 0.0, grid);

julia> atmos = Atmos(Float64; u = 0.0, v = 0.1, temp = 0.0, grid);

julia> floe_settings = FloeSettings(Float64);

julia> floes = initializefloefield(Float64, 3, [0.5], domain, 0.25, 0; floe_settings, rng = Xoshiro(1));

julia> Model(; grid, domain, ocean, atmos, floes) Model{Float64, ...}

⊢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 20 grid cells of size 25000.0 m

⊢Domain ⊢Northern boundary of type CollisionBoundary{North, Float64} ⊢Southern boundary of type CollisionBoundary{South, Float64} ⊢Eastern boundary of type CollisionBoundary{East, Float64} ⊢Western boundary of type CollisionBoundary{West, Float64} ∟0-element TopograpahyElement{Float64} list

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

⊢Atmos{Float64} ⊢Vector fields of dimension (21, 21) ⊢Tracer fields of dimension (21, 21) ⊢Average u-velocity of: 0.0 m/s ⊢Average v-velocity of: 0.1 m/s ∟Average temperature of: 0.0 C

⊢Floe List: ⊢Number of floes: 3 ⊢Total floe area: 1.375278018545777e11 ∟Average floe height: 0.25

```

Constants{FT}

Each simulation needs a list of physical constants to be used in various simulation calculations. All constants have default values that will be used if the user does not provide an alternative.

Fields

• ρo::FT: Ocean density (default = 1027.0)
• ρa::FT: Air density (default = 1.2)
• Cd_io::FT: Ice-ocean drag coefficent (default = 3e-3)
• Cd_ia::FT: Ice-atmosphere drag coefficent (default = 1e-3)
• Cd_ao::FT: Atmosphere-ocean momentum drag coefficient (default = 1.25e-3)
• f::FT: Ocean coriolis frequency (default = 1.4e-4)
• turnθ::FT: Ocean turn angle (default = 15π/180)
• L::FT: Latent heat of freezing Joules/kg
• k::FT: Thermal conductivity of surface iceW/(m*K)
• ν::FT: Poisson's ratio (default = 0.3)
• μ::FT: Coefficent of friction (default = 0.2)
• E::FT: Young's Modulus (default = 6e6)

Here is how to construct a Constants:

Constants([FT = Float64]; kwargs...)

Positional arguments

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

Keyword arguments

• Each of the above fields is an optional keyword argument.

Examples

• Creating Constants with a zero Coriolis frequency

julia> Constants(; f = 0)
Constants{Float64}
  ⊢ ρo = 1027.0
  ⊢ ρa = 1.2
  ⊢ Cd_io = 0.003
  ⊢ Cd_ia = 0.001
  ⊢ Cd_ao = 0.00125
  ⊢ f = 0.0
  ⊢ turnθ = 0.2617993877991494
  ⊢ L = 293000.0
  ⊢ k = 2.14
  ⊢ ν = 0.3
  ⊢ μ = 0.2
  ⊢ E = 6.0e6

• Creating Constants with all default constants, but Float32 float-type

julia> Constants(Float32)
Constants{Float32}
  ⊢ ρo = 1027.0
  ⊢ ρa = 1.2
  ⊢ Cd_io = 0.003
  ⊢ Cd_ia = 0.001
  ⊢ Cd_ao = 0.00125
  ⊢ f = 0.00014
  ⊢ turnθ = 0.2617994
  ⊢ L = 293000.0
  ⊢ k = 2.14
  ⊢ ν = 0.3
  ⊢ μ = 0.2
  ⊢ E = 6.0e6

FloeSettings{FT, GT, CT}

When you create a floe or a set of floes, you have the option to create a floe settings object. This set of settings controls certian floe fields and calculations.

Fields

• ρi::FT: floe's density (Default = 920.0 g/L)
• min_floe_area::FT: minimum floe area (Default = 1e6 m^2)
• min_floe_height::FT: minimum floe height (Default = 0.1 m)
• max_floe_height::FT: maximum floe height (Default = 10.0 m)
• min_aspect_ratio::FT: minimum ratio between floe x-length and y-length by maximum coordiante values (Default = 0.05)
• maximum_ξ::FT: the absolute maximum rotational velocity a floe can reach before it is capped at maximum_ξ (Default = 1e-5 rad/s)
• subfloe_point_generator::GT: subtype of AbstractSubFloePointsGenerator, which generates floe's subfloe points. These points which determines the method of subfloe point generation is used for each floe (Default = MonteCarloPointsGenerator)
• stress_calculator::CT: subtype of AbstractStressCalculator, which generates the calculator for stress, which in turn determines the method of calculating current stress of floes during the simulation (Default = DecayAreaScaledCalculator)

If any of the minimum values are exceeded, a floe is removed in the course of the simulation. If any of the maximum values are reached, the value is capped at the given value.

Here is how to construct a FloeSettings:

FloeSettings([FT = Float64]; kwargs..)

Positional arguments

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

Keyword arguments

• Each of the above fields is an optional keyword argument.

Examples

• Creating default FloeSettings

julia> FloeSettings()
FloeSettings{Float64, MonteCarloPointsGenerator{Float64}, DecayAreaScaledCalculator{Float64}}
  ⊢ ρi = 920.0
  ⊢ min_floe_area = 1.0e6
  ⊢ min_floe_height = 0.1
  ⊢ max_floe_height = 10.0
  ⊢ min_aspect_ratio = 0.05
  ⊢ maximum_ξ = 1.0e-5
  ⊢ subfloe_point_generator = MonteCarloPointsGenerator{Float64}(1000, 10, 0.1)
  ⊢ stress_calculator = DecayAreaScaledCalculator{Float64}(0.2, 0.0)

• Creating a Float32 FloeSettings with a zero minimum floe area and a SubGridPointsGenerator generator

julia> subfloe_point_generator =  SubGridPointsGenerator(Float32; Δg = 1000);

julia> FloeSettings(Float32; min_floe_area = 0, subfloe_point_generator)
FloeSettings{Float32, SubGridPointsGenerator{Float32}, DecayAreaScaledCalculator{Float32}}
  ⊢ ρi = 920.0
  ⊢ min_floe_area = 0.0
  ⊢ min_floe_height = 0.1
  ⊢ max_floe_height = 10.0
  ⊢ min_aspect_ratio = 0.05
  ⊢ maximum_ξ = 1.0e-5
  ⊢ subfloe_point_generator = SubGridPointsGenerator{Float32}(1000.0f0)
  ⊢ stress_calculator = DecayAreaScaledCalculator{Float32}(0.2f0, 0.0f0)

CouplingSettings

Settings needed for coupling within the model.

Fields

• coupling_on::Bool: if true, the model will be coupled with the simulation's ocean and atmosphere (Default = true)
• Δt::Int: determines how many simulation timesteps between calculating ocean and atmospheric forces on the floes (Default = 10)
• Δd::Int: Δd number of buffer grid cells on each side of floe for interpolation using sub-floe points - see AbstractSubFloePointsGenerator (Default = 1)
• two_way_coupling_on::Bool: if true, then the simulation calculates the stress the ice/atmosphere put on the ocean so that the user could couple to Oceananigans (Default = False)

Keyword arguments

• Each of the above fields is an optional keyword argument.

Examples

• Creating default CouplingSettings

julia> CouplingSettings()
CouplingSettings(true, 10, 1, false)

• Creating CouplingSettings with a two_way_coupling_on and a larger interpolation buffer.

julia> CouplingSettings(Δd = 2, two_way_coupling_on = true)
CouplingSettings(true, 10, 2, true)

CollisionSettings{FT}

Settings needed for collisions within the model.

Fields

• collisions_on::Bool: if true, collisions between floes and the domain elements will occur (Default = true)
• floe_floe_max_overlap::FT: defines the percentage of overlap allowed between floes before marking them for ridging/rafting (Default = 0.55)
• floe_domain_max_overlap::FT: defines the percentage of overlap allowed between floes and DomainElements before removing the floes from the simulation.

Both floe_floe_max_overlap and floe_domain_max_overlap should be between 0-1 and if a value < 0 is given or a value > 1 is given when collisions_on is true they will be set to 0 and 1 respectively.

Positional arguments

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

Keyword arguments

• Each of the above fields is an optional keyword argument.

Examples

• Creating default CollisionSettings

julia> CollisionSettings()
CollisionSettings{Float64}(true, 0.55, 0.75)

• Creating a Float32 CollisionSettings with a higher floe_floe_max_overlap

julia> CollisionSettings(Float32;  floe_floe_max_overlap = 0.65)
CollisionSettings{Float32}(true, 0.65f0, 0.75f0)

FractureSettings{CT}

Settings needed for fractures within the model.

Fields

• fractures_on::Bool: if true, floes will fracture, else they will not (Default = false)
• criteria::CT: defines which fracture criteria (subtype of AbstractFractureCriteria) are used to determine which floes to fracture
• Δt::Int: determines how many simulation timesteps between attempting floe fracture (Default = 0)
• deform_on::Bool: if true, then the floe will be deformed around floe primarily causing the fracture, identified by the largest overlap area on the most recent set of collisions (Default = False)
• npieces::Int: how many pieces to try to split a fractured floe into (Default = 3)

Keyword arguments

• Each of the above fields is an optional keyword argument.

Examples

• Creating default FractureSettings

julia> FractureSettings()
FractureSettings{NoFracture}
  ⊢ fractures_on = false
  ⊢ criteria = NoFracture
  ⊢ Δt = 0
  ⊢ deform_on = false
  ⊢ npieces = 3

• Creating a Float32 FractureSettings with fractures on

julia> mohrs_cone_criteria = MohrsCone(Float32);

julia> FractureSettings(; fractures_on = true, criteria = mohrs_cone_criteria)
FractureSettings{MohrsCone{Float32}}
  ⊢ fractures_on = true
  ⊢ criteria = MohrsCone{Float32}
  ⊢ Δt = 0
  ⊢ deform_on = false
  ⊢ npieces = 3

abstract type AbstractFractureCriteria

Abstract super type for fracture criteria, which determines if a floe fractures given its physical condition. Right now, fracture criteria define a set of vertices that define a shape in principal stress space. The minimum and maximum eigenvalues of a floe's stress field are then a point within the principal stress space. If that stress point falls outside of the criteria-verticies defined vertices, it is a stress big enough to fracture the floe. Otherwise the floe will not be fractured. If another form of fracture criteria is desired, some additional code may need to be pulled into helper fuctions within the fracture.jl file and then dispatched off of using subtypes of AbstractFractureCriteria. This is encuarged if It is helpful!

Right now, the subtypes serve as dispatch types for the following two methods.

API

The following methods must be implemented for all subtypes:

• __update_criteria!(criteria::AbstractFractureCriteria, floes::StructArray{<:Floe})
• _determine_fractures(criteria::AbstractFractureCriteria, floes::StructArray{<:Floe}, floe_settings::FloeSettings)

__update_criteria! is called in the fracture_floes! function and takes the fracture criteria and the current ice floe pack and updates the fracture criteria given the state of the ice floe pack if needed.

_determine_fractures is called in the fracture_floes! and takes in the fracture criteria, the current ice floe pack, and the floe settings and determines which of the ice floes fracture. It returns a Boolean vector equal in length to the floes list. If the ith element of the returned vector is true then the ith floe in floes will fracture.

Abstract type for fracture criteria. Each struct of this type must have a vertices field representing the criteria in principal stress space. For a given polygon, the minimum and maximum eigenvalues of its stress field will be its location in principal stress space. If that stress point falls outside of the criteria-verticies defined polygon it is a stress great enough to fracture the floe. Otherwise the floe will not be fractured. Each fracture criteria type must also have an updatecriteria! function defined that is used to update the criteria each timestep. If the criteria does not need to be updated, this function can be empty.

See the existing concrete subtypes: NoFracture, MohrsCone, and HiblerYieldCurve.

NoFracture<:AbstractFractureCriteria

Type of AbstractFractureCriteria representing when fracturing functionality is turned off. If this is the type provided to the simulation's FractureSettings, then fractures will not occur. As fracturing will not occur, this subtype depends on the default API functions.

HiblerYieldCurve{FT} <: AbstractFractureCriteria

Concrete subtype of AbstractFractureCriteria that calculates Hibler's Elliptical Yield curve using parameters pstar, c, and the current floe field. This is an elliptical yield curve that determines if a floe fractures based off if its stress in principal stress space is inside or outside of that elliptical yield curve.

Fields

• pstar::AbstractFloat: parameter used to tune ellipse for optimal fracturing (Default = 2.25e5)
• c::AbstractFloat: parameter used to tune ellipse for optimal fracturing (Default = 20)
• poly::Polys{FT}: polygon that defines the yield curve in principal stress space

Here is how to construct a HiblerYieldCurve object:

HiblerYieldCurve([FT = Float64]; floes, pstar, c)

Positional arguments

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

Keyword arguments

• floes::StructArray{Floes}: models's list of floes
• pstar::AbstractFloat: parameter used to tune ellipse for optimal fracturing (Default = 2.25e5)
• c::AbstractFloat: parameter used to tune ellipse for optimal fracturing (Default = 20)

Examples

• Creating default HiblerYieldCurve

julia> using Random

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

julia> north, south, east, west = CollisionBoundary(North, Float64; grid), CollisionBoundary(South, Float64; grid), CollisionBoundary(East, Float64; grid), CollisionBoundary(West, Float64; grid);

julia> domain = Domain(; north, south, east, west);

julia> floes = initialize_floe_field(Float64, 3, [0.5], domain, 0.25, 0; floe_settings = FloeSettings(Float64), rng = Xoshiro(1));

julia> hibler = HiblerYieldCurve(Float64; floes)
HiblerYieldCurve{Float64}
  ⊢ pstar: 225000.0
  ⊢ c: 20.0
  ⊢ yield curve area: 2.48338091663081e9
  ⊢ yield curve centroid: (-28125.0, -28125.0)

julia> hibler = HiblerYieldCurve(Float32; floes, c = 24)
HiblerYieldCurve{Float32}
  ⊢ pstar: 225000.0
  ⊢ c: 24.0
  ⊢ yield curve area: 2.4833810505e9
  ⊢ yield curve centroid: (-28124.99995, -28125.00086)

MohrsCone{FT} <: AbstractFractureCriteria

Concrete subtype of AbstractFractureCriteria that creates a conical yield curve in principal stress space. This conical yield curve determines if a floe fractures based off if its stress in principal stress space is inside or outside of that elliptical yield curve.

Fields

• poly::Polys{FT}: polygon that defines the yield curve in principal stress space

Here is how to construct a MohrsCone object:

MohrsCone([FT = Float64]; kwargs...)

Positional arguments

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

Keyword arguments

• q::AbstractFloat: based on the coefficient of internal friction ($µ_i$) by $(μ_i^2 + 1)^(1/2) + μ_i^2$ (Default = 5.2)
• σc::AbstractFloat: uniaxial compressive strength (Default = 2.5e5)
• σ11::AbstractFloat: negative of the x-coordinate of one vertex of cone (triangle in 2D) and negative of the y-coordinate of adjacend vertex in principal stress space (Default = -3.375e4)
• σ1::AbstractFloat: x-coordiante of first point in cone (Default = nothing)
• σ2::AbstractFloat: y-coordiante of first point in cone (Default = nothing)
• σ22::AbstractFloat: y-coordinate of one vertex of cone and negative of the x-coordinate of adjacend vertex in principal stress space (Default = nothing)

Examples

• Creating default MohrsCone

julia> MohrsCone()
MohrsCone{Float64}
  ⊢ Points: (59523.80952, 59523.80952), (33750.0, -74500.0), (-74500.0, 33750.0), (59523.80952, 59523.80952)

• Creating Float32 MohrsCone

julia> MohrsCone(Float32)
MohrsCone{Float32}
  ⊢ Points: (59523.81f0, 59523.81f0), (33750.0f0, -74500.0f0), (-74500.0f0, 33750.0f0), (59523.81f0, 59523.81f0)

SimplificationSettings{FT}

Settings needed for floe simplification within the simulation.

Fields

• smooth_vertices_on::Bool: if true then floe's with more vertices than max_vertices will be simplified every Δt_smooth timesteps (Default = True)
• max_vertices::Int: total number of vertices a floe can have without triggering smoothing if smooth_vertices_on = True (Default = 30)
• tol::FT: tolerance in meters of the the Douglas–Peucker algorithm, which is used to smooth edges of polygons (Default = 100)
• Δt_smooth::Int: determines how many simulation timesteps between attempting floe smoothing (Default = 20)

Positional arguments

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

Keyword arguments

• Each of the above fields is an optional keyword argument.

Examples

• Creating default SimplificationSettings

julia> SimplificationSettings()
SimplificationSettings{Float64}(true, 30, 100.0, 20)

• Creating a Float32 SimplificationSettings with smooth_vertices_on off

julia> SimplificationSettings(Float32; smooth_vertices_on = false)
SimplificationSettings{Float32}(false, 30, 100.0f0, 20)

RidgeRaftSettings{FT}

Settings needed for ridging and rafting within the simulation.

Fields

• ridge_raft_on::Bool: if true, ridging and rafting should be turned on in the simulation (Default = False)
• Δt: determines how many simulation timesteps between attempting ridging and rafting (Default = 0)
• ridge_probability::FT: the probability a floe ridges with another floe/domain if it meets all other criteria
• raft_probability::FT: the probability a floe rafts with another floe/domain if it meets all other criteria
• min_overlap_frac::FT: the minimum overlap area fraction between a floe and another floe/domain for that floe to ridge or raft
• min_ridge_height::FT: the minimum floe height to ridge with a floe/domain
• max_floe_ridge_height::FT: the maximum floe height to ridge with another floe
• max_domain_rdige_height::FT: maximum floe height to ridge with a domain element
• max_floe_raft_height::FT: maximum floe height to raft with another floe
• max_domain_raft_height::FT: maximum floe height to raft with a domain element
• domain_gain_probability::FT: the probalility that a floe that rafts with a domain element keeps all of its mass (0) or if that mass is removed and lost to the domain element (1).

Positional arguments

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

Keyword arguments

• Each of the above fields is an optional keyword argument.

Examples

• Creating default RidgeRaftSettings

julia> RidgeRaftSettings()
RidgeRaftSettings{Float64}
  ⊢ ridge_raft_on = false
  ⊢ Δt = 0
  ⊢ ridge_probability = 0.95
  ⊢ raft_probability = 0.95
  ⊢ min_overlap_frac = 0.01
  ⊢ min_ridge_height = 0.2
  ⊢ max_floe_ridge_height = 5.0
  ⊢ max_domain_ridge_height = 1.25
  ⊢ max_floe_raft_height = 0.25
  ⊢ max_domain_raft_height = 0.25
  ⊢ domain_gain_probability = 1.0

• Creating a Float32 RidgeRaftSettings with ridging and rafting turned on

julia> RidgeRaftSettings(Float32; ridge_raft_on = true)
RidgeRaftSettings{Float32}
  ⊢ ridge_raft_on = true
  ⊢ Δt = 0
  ⊢ ridge_probability = 0.95
  ⊢ raft_probability = 0.95
  ⊢ min_overlap_frac = 0.01
  ⊢ min_ridge_height = 0.2
  ⊢ max_floe_ridge_height = 5.0
  ⊢ max_domain_ridge_height = 1.25
  ⊢ max_floe_raft_height = 0.25
  ⊢ max_domain_raft_height = 0.25
  ⊢ domain_gain_probability = 1.0

WeldSettings{FT}

Settings needed for welding within the simulation.

Fields

• weld_on::Bool: a boolean flag for if welding should be turned on in the simulation (Default = False)
• Δts::Vector{Int}: a list of multiples of timesteps during which welding code will run, welding will be run at multiples of all elements, each with domain split (Default = Int[]) into corresponding Nx and Ny values
• Nxs::Vector{Int}: a list of number of x-directional bins to split the domain into at corresponding timesteps (Default = Int[])
• Nys::Vector{Int}: a list of number of x-directional bins to split the domain into at corresponding timesteps (Default = Int[])
• min_weld_area::FT: minimum area a weld can create for two floes to weld (Default = 1e6)
• max_weld_area::FT: maximum area a weld can create for two floes to weld (Default = 2e9)
• welding_coeff::FT: non-dimensional parameter, multiplied by ratio of overlap between two floes to original floe area to determin probability that a floe will merge. The larger this is, the more likely floes are to weld. Probability with 5% overlap is welding_coeff * (0.05) > rand() (Default = 150)

Positional arguments

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

Keyword arguments

• Each of the above fields is an optional keyword argument.

Examples

• Creating default WeldSettings

julia> WeldSettings()
WeldSettings{Float64}
  ⊢ weld_on = false
  ⊢ Δts = Int64[]
  ⊢ Nxs = Int64[]
  ⊢ Nys = Int64[]
  ⊢ min_weld_area = 1.0e6
  ⊢ max_weld_area = 2.0e9
  ⊢ welding_coeff = 150.0

• Creating a Float32 WeldSettings with welding turned on

julia> WeldSettings(Float32; weld_on = true, Δts = [100, 200], Nxs = [2, 1], Nys = [1, 2])
WeldSettings{Float32}
  ⊢ weld_on = true
  ⊢ Δts = [200, 100]
  ⊢ Nxs = [1, 2]
  ⊢ Nys = [2, 1]
  ⊢ min_weld_area = 1.0e6
  ⊢ max_weld_area = 2.0e9
  ⊢ welding_coeff = 150.0

AbstractOutputWriter

An abstract type for output writers that provide data from simulation runs.

Right now, there are four types of output writers: InitialStateOutputWriter, CheckpointOutputWriter, FloeOutputWriter, and GridOutputWriter.

They currently don't dispatch off of anything.

InitialStateOuputWriter <: AbstractOutputWriter

Concrete subtype of AbstractOutputWriter that records the intial state of the simulation so that you can re-load it later. Writes JLD2 file with initial simulation state to the filepath specified. If overwrite is true, and there is a file of the given name at the filepath, that file will be overwritten.

Fields

• filepath::String: filepath to the output files
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)

Keyword arguments

• dir::String: directory the files should be saved in (Default = ".")
• filename::String: name of output file (Default = "initial_state.jld2")
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)
• jld2_kw::Dict: dictionary of keywords for the jldopen function

Here is how to construct an InitialStateOutputWriter:

InitialStateOutputWriter(; kwargs...)

InitialStateOutputWriter <: AbstractOutputWriter

InitialStateOutputWriter can also be created using an existing InitialStateOutputWriter, copying all fields unless the new field values are explicity provided through keyword arguments.

Positional arguments

• writer::InitialStateOutputWriter: existing InitialStateOutputWriter

Keyword arguments

• dir::String: directory the files should be saved in (Default = ".")
• filename::String: name of output file (Default = "initial_state.jld2")
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)
• kwargs...

CheckpointOutputWriter <: AbstractOutputWriter

Checkpoint subtype of AbstractOutputWriter that holds information for outputting checkpoint information used for restarting the simulation from a point where the writer saved data. Checkpoint data is saved every Δtout timesteps to filepath. If the given file doesn't end in ".jld2", the extension will be appended. If overwrite is true then if there is already a file of the given name, it will be overwriten. Else it will thrown an error.

Fields

• Δtout::Int: number of timesteps between output
• filepath::String: filepath to the output files
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)

Positional arguments

• Δtout::Int: number of timesteps between output

Keyword arguments

• dir::String: directory the files should be saved in (Default = ".")
• filename::String: name of output file (Default = "checkpoint.jld2")
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)
• jld2_kw::Dict: dictionary of keywords for the jldopen function

Here is how to construct an CheckpointOutputWriter:

CheckpointOutputWriter(Δtout; kwargs...)

CheckpointOutputWriter <: AbstractOutputWriter

A CheckpointOutputWriter can also be created from an existing CheckpointOutputWriter, copying all fields unless new field values are explicity provided either as the optional argument Δtout or through keyword arguments.

Positional arguments

• writer::InitialStateOutputWriter: OPTIONAL argument - if present, the dir, filename, and overwrite of the provided writer will be used for the new writer, rather than the above argument
• Δtout::Int: number of timesteps between output

Keyword arguments

• dir::String: directory the files should be saved in (Default = ".")
• filename::String: name of output file (Default = "checkpoint.jld2")
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)
• kwargs...

Here is how to construct a new CheckpointOutputWriter using an existing CheckpointOutputWriter:

CheckpointOutputWriter(existing_checkpointer, Δtout; kwargs...)

GridOutputWriter{FT} <: AbstractOutputWriter

Grid subtype of AbstractOutputWriter that holds information for outputting floe data on the grid. This output does not need to be on the grid defined for the model. This grid can be coarser, or more fine, as defined by the xg and yg fields. Output on this scale will be saved to the file defined by filepath every Δtout timesteps. Data will be collected in the data field during calculation for easier writing to the NetCDF file. Only outputs within the outputs list will be saved. There is a limited number of floe outputs that can be calculated by the GridOutputWriter.

Fields

• outputs::Vector{Symbol}: list of field names (as symbols) to output - call get_known_grid_outputs() to see options
• Δtout::Int: number of timesteps between output
• filepath::String: filepath to the output files
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)
• xg::Vector{FT}: x-grid points on output grid
• yg::Vector{FT}: y-grid points on output grid
• data::Array{FT, 3}: 3D output grid of size nx, ny, and nt where nt is the number of timesteps saved
• average::Bool: if true, average gridded data over timesteps between outputs, else just calculate at output timestep

Positional arguments

• Δtout::Int: number of timesteps between output
• grid::AbstractRectilinearGrid: subtype of AbstractRectilinearGrid representing the grid used within a simulation
• dims::Tuple{Int, Int}: output new grid dimensions for these calculations (rows -> ny, cols -> nx)

Keyword arguments

• outputs::Vector{Symbol}: list of field names (as symbols) to output - call get_known_grid_outputs() to see options
• dir::String: directory the files should be saved in (Default = ".")
• filename::String: name of output file (Default = "gridded_data.nc")
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)
• average::Bool: if true, average gridded data over timesteps between outputs, else just calculate at output timestep

Here is how to construct an GridOutputWriter:

GridOutputWriter(Δtout, grid, dims; kwargs...)

GridOutputWriter <: AbstractOutputWriter

GridOutputWriter can also be created using existing GridOutputWriter, copying all fields unless the new field values are explicity provided through keyword arguments.

Positional arguments

• writer::InitialStateOutputWriter: OPTIONAL argument - if present, the dir, filename, and overwrite of the provided writer will be used for the new writer, rather than the above argument
• Δtout::Int: number of timesteps between output

Keyword arguments

• grid::AbstractRectilinearGrid: subtype of AbstractRectilinearGrid representing the grid used within a simulation
• dims::Tuple{Int, Int}: output new grid dimensions for these calculations (rows -> ny, cols -> nx)
• outputs::Vector{Symbol}: list of field names (as symbols) to output - call get_known_grid_outputs() to see options
• dir::String: directory the files should be saved in (Default = ".")
• filename::String: name of output file (Default = "gridded_data.nc")
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)
• average::Bool: if true, average gridded data over timesteps between outputs, else just calculate at output timestep

Here is how to construct a new GridOutputWriter using an existing GridOutputWriter:

GridOutputWriter(existing_gridwriter, Δtout; kwargs...)

FloeOutputWriter <: AbstractOutputWriter

Floe subtype of AbstractOutputWriter that holds information for outputting floe information from model throughout simulation. Output will be saved to the file defined by filename every Δtout timesteps. Only outputs within the outputs list will be saved. The outputs field takes in a list of symbols corresponding to floe fields. For example, if you want the floe output writer to output the floes centroid and coordinates then outputs = [:centroid, :coords]. If you want all floe fields then you can simply omit the outputs field all together and all floe fields will be output.

File will be saved as a JLD2 file to filepath. If the given file doesn't end in ".jld2", the extension will be appended. If overwrite is true then if there is already a file of the given name, it will be overwriten. Else it will thrown an error.

Fields

• Δtout::Int: number of timesteps between output
• outputs::Vector{Symbol}: list of Floe field names (as symbols) to output
• filepath::String: filepath to the output files
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)

Positional arguments

• Δtout::Int: number of timesteps between output

Keyword arguments

• outputs::Vector{Symbol}: list of Floe field names (as symbols) to output (Default outputs all Floe fields)
• dir::String: directory the files should be saved in (Default = ".")
• filename::String: name of output file (Default = "floes.jld2")
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)
• jld2_kw::Dict: dictionary of keywords for the jldopen function

Here is how to construct an FloeOutputWriter:

FloeOutputWriter(Δtout; kwargs...)

FloeOutputWriter <: AbstractOutputWriter

A FloeOutputWriter can also be created using an existing FloeOutputWriter, copying all fields unless the new field values are explicity provided through keyword arguments.

Positional arguments

• writer::InitialStateOutputWriter: OPTIONAL argument - if present, the dir, filename, and overwrite of the provided writer will be used for the new writer, rather than the above argument
• Δtout::Int: number of timesteps between output

Keyword arguments

• outputs::Vector{Symbol}: list of Floe field names (as symbols) to output
• dir::String: directory the files should be saved in (Default = ".")
• filename::String: name of output file (Default = "floes.jld2")
• overwrite::Bool: if true, existing files with the same filepath will be overwritten (Default = false)
• kwargs...

Here is how to construct a new FloeOutputWriter using an existing FloeOutputWriter:

FloeOutputWriter(existing_floewriter, Δtout; kwargs...)

OutputWriters{FT}

Structure to hold all types of output writers a user might want. Fields are vectors so that more than one of each type of output writer can be defined, and so that a default OutputWriter object doesn't create default output writer fields, which would create files, but rather empty lists of output writers. If any of the fields is not provided, the default is just an empty list.

Fields

• initialwriters::StructVector{<:InitialStateOutputWriter}: vector of InitialStateOutputWriters
• floewriters::StructVector{<:FloeOutputWriter}: vector of FloeOutputWriters
• gridwriters::StructVector{<:GridOutputWriter}: vector of GridOutputWriters
• checkpointwriters::StructVector{<:CheckpointOutputWriter}: vector of CheckpointOutputWriters

Positional arguments

• All output writers (one per argument) - they will be aggregated and put into lists by type

SubzeroLogger(; sim = nothing, filename = " ", messages_per_tstep = 1)

Logger for Subzero. Logs unique messages messages_per_tstep times per timestep to prevent overwhelming number of messages timesteps from multiple floes triggering the same log event.

Fields

• stream::IO: logs are written to this IO
• min_level::Logging.LogLevel: minimum log event level to write
• message_limits::Dict{Any, Int}: dictionary with message IDs for key whose values are the number of times that message can still be written in current timestep - current timestep is stored in same dictionary under key 'tstep'
• messages_per_tstep::Int: maximum number of times a given message should be written per timestep

Keyword arguments

• sim::Simulation: Subzero simulation - used to get simulation name (Default = nothing), which is then used to name log file
• filename::String: if simulation isn't provided, a filename (+ path) must be provided to save log file
• messages_per_tstep::Int: maximum number of times a given message should be written per timestep

Simulation{FT, MT, CT, PT, ST, RT, OT}

Simulation which holds a model and the parameters, settings, and output writers needed for running the simulation.

Only keyword arguments are used!

Fields / Keyword Arguments

General

• model::MT: Model to simulate
• consts::Constants{FT}: Constants used in Simulation (default = Constants())
• rng::RT: Random number generator (default = Xoshiro())
• verbose::Bool: String output printed during run (Default = false)
• name::String: Simulation name for printing/saving (Default = "sim")

Timesteping Information

• Δt::Int: Simulation timestep in seconds
• nΔt::Int: Total timesteps simulation runs for

Physical Processes

• floe_settings::FloeSettings{FT, PT, ST}: Settings that control floe size/mass/etc - no default!
• coupling_settings::CouplingSettings: Settings that control coupling between floes/ocean/atmosphere (Default = CouplingSettings())
• collision_settings::CollisionSettings{FT}: Settings that control floe collisions with other floes and the domain (Default = CollisionSettings())
• fracture_settings::FractureSettings{CT}: Settings that control floe fracturing (Default = FractureSettings())
• simp_settings::SimplificationSettings{FT}: Settings that control the simplification of floes (Default = SimplificationSettings())
• ridgeraft_settings::RidgeRaftSettings{FT}: Settings that control floe ridging and rafting (Default = RidgeRaftSettings())
• weld_settings::WeldSettings{FT}: Settings that control floe welding (Default = WeldSettings())

Output Writers

• writers::OT: Simulation output writers (Default = OutputWriters())

run!(sim; logger, messages_per_tstep, start_tstep)

Run given simulation and generate output for given output writers. Simulation calculations will be done with Floats of type FT (Float64 of Float32).

Positional arguments

• sim::Simulation: simulation to be run

Keyword arguments

• logger::AbstractLogger: logger for simulation (Default = Nothing, which triggers use of SubzeroLogger
• messages_per_tstep::Int"` number of messages to print per timestep if using default SubzeroLogger, else not needed (Default = 1)
• start_tstep::Int: which timestep to start the simulation on (Default = 0)

Returns

• None. The simulation will be run and outputs will be saved in the output folder.

restart!(initial_state_fn, checkpointer_fn, new_nΔt, new_output_writers; start_tstep)

Continue the simulation run started with the given initial state and floe file for an additional new_nΔt timesteps and with the new outputwriters provided. The simulation will restart with a recorded timestep of `starttstep`.

Note that this restart! function may not fit your needs and you may need to write your own. This function is meant to act as a simplest case and as a template for users to write their own restart functions.

Positional arguments

- `initial_state_fn::String`: file path to previously run simulation's initial state file
- `checkpointer_fn::String`: file path to previously run simulation's checkpointer file.
The simulation will restart right after the checkpoint captured in this file
- `new_output_writers::OutputWriters`: new output writers for the new simulation - new ones are required to the output writes to new, unique files.

Keyword arguments

- `start_tstep::Int`: which timestep to start the simulation on (Default = 0)

Returns

- None. The simulation will be run and outputs will be saved in the output folder.

timestep_sim!(sim, tstep, start_tstep)

Run one step of the simulation and write output.

Positional arguments

• sim::Simulation: simulation to be run
• tstep::Int: simulation's current timestep
• start_tstep::Int: timestep simulation started on (Default = 0)

Returns

• None. Simulation advances by one timestep.

write_data!(sim, tstep, start_tstep)

Writes data for the simulation's writers that are due to write at given tstep.

Positional arguments

• sim::Simulation: simulation to be run
• tstep::Int: current simulation timestep
• start_tstep::Int: starting timestep of the simulation

write_init_state_data!(sim, tstep)

Save initial simulation state.

Positional arguments

• sim::Simulation: simulation to be run
• tstep::Int: current simulation timestep

write_checkpoint_data!(sim, tstep)

Writes model's floe, ocean, and atmosphere data to JLD2 file. Data can be used to restart simulation run. Writes floes, ocean, and atmosphere to JLD2 file with name writer.fn for current timestep, which will be the group in the JLD2 file.

Positional arguments

• sim::Simulation: simulation to be run
• tstep::Int: current simulation timestep

write_floe_data!(sim, tstep)

Writes desired FloeOutputWriter data to JLD2 file. Writes desired fields writer.outputs to JLD2 file with name writer.fn for current timestep, which will be the group in the JLD2 file.

Positional arguments

• writers::StructArray{FloeWriter}: list of floe writers
• floes::StructArray{Floe}: simulation floes
• tstep::Int: current simulation timestep

write_grid_data!(sim, tstep)

Writes desired GridOutputWriter data to NetCDF file. Writes desired fields writer.outputs to file with name writer.fn for current timestep.

Positional arguments

• writers::StructArray{GridWriter}: list of grid writers
• floes::StructArray{Floe}: simulation floes
• topography::StructArray{TopographyElement}: simulation's list of topography elements
• tstep::Int: current simulation timestep

calc_eulerian_data!(floes, topography, writer, istep)

Calculate floe data averaged on grid defined by GridOutputWriter for current timestep (istep). Saved in writer.data field

Positional arguments

• floes::StructArray{Floe}: simulation floes
• topography::StructArray{TopographyElement}: simulation's list of topography elements
• writer::GridWriter: grid output writer
• istep::Int: current simulation timestep

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.