Simple Strait Simulation

This simulation creates a north to south strait that ice can flow through, pushed by the ocean. The north and south boundaries form a periodic pair, so that the ice can endlessly flow through the strait. The east and west boundaries are collision bounds, but they are completely covered with topography forming the edges of the domain. This is a good simulation to understand how to setup topography and how to turn on fractures using the fracture settings.

using Subzero, CairoMakie, GeoInterfaceMakie
using JLD2, Random, Statistics

Making the Simulation

User Inputs

const FT = Float64  # Float type used to run simulation
const Lx = 1e5      # grid x-length
const Ly = 1e5      # grid y-length
const Δgrid = 2e3   # grid cell edge-size
const hmean = 0.25  # mean floe height
const Δh = 0.0      # difference in floe heights - here all floes are the same height
const Δt = 20       # timestep
const nΔt = 5000;    # number of timesteps to run

Grid Creation

grid = RegRectilinearGrid(; x0 = 0.0, xf = Lx, y0 = 0.0, yf = Ly, Δx = Δgrid, Δy = Δgrid)

RegRectilinearGrid{Float64}
  ⊢x extent (0.0 to 100000.0) with 50 grid cells of size 2000.0 m
  ∟y extent (0.0 to 100000.0) with 50 grid cells of size 2000.0 m

Domain Creation

nboundary = PeriodicBoundary(North; grid)
sboundary = PeriodicBoundary(South; grid)
eboundary = CollisionBoundary(East; grid)
wboundary = CollisionBoundary(West; grid)

island1 = [[[6e4, 4e4], [6e4, 4.5e4], [6.5e4, 4.5e4], [6.5e4, 4e4], [6e4, 4e4]]]
island2 = [[[4e4, 6e4], [4e4, 6.5e4], [4.5e4, 6.5e4], [4.5e4, 6e4], [4e4, 6e4]]]
topo1 = [[[0, 0.0], [0, 1e5], [2e4, 1e5], [3e4, 5e4], [2e4, 0], [0.0, 0.0]]]
topo2 = [[[8e4, 0], [7e4, 5e4], [8e4, 1e5], [1e5, 1e5], [1e5, 0], [8e4, 0]]]
topo_arr = initialize_topography_field(FT; coords = [island1, topo1, topo2])

domain = Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary, topography = topo_arr)

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

Ocean Creation

ocean = Ocean(; grid, u = 0.0, v = -0.3, temp = 0.0)

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

Atmos Creation

atmos = Atmos(; grid, u = 0.0, v = 0.0, temp = 0.0)

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

Floe Creation

floe_settings = FloeSettings(
    subfloe_point_generator = SubGridPointsGenerator(; grid, npoint_per_cell = 2),
    stress_calculator = DecayAreaScaledCalculator(),
)
floe_generator = VoronoiTesselationFieldGenerator(; nfloes = 75, concentrations = [0.7], hmean, Δh)
floe_arr = initialize_floe_field(FT; generator = floe_generator, domain, rng = Xoshiro(3), floe_settings = floe_settings)

80-element StructArray(::Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}, ::Vector{Vector{Float64}}, ::Vector{Vector{Vector{Vector{Float64}}}}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Vector{Float64}}, ::Vector{Vector{Float64}}, ::Vector{Vector{Float64}}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Subzero.Status}, ::Vector{Int64}, ::Vector{Int64}, ::Vector{Vector{Int64}}, ::Vector{Vector{Int64}}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Matrix{Float64}}, ::Vector{Float64}, ::Vector{Matrix{Float64}}, ::Vector{Int64}, ::Vector{Matrix{Float64}}, ::Vector{Matrix{Float64}}, ::Vector{Matrix{Float64}}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}) with eltype Floe{Float64}:
 Floe{Float64}
  ⊢Centroid of (65396.35599, 4713.68728) m
  ⊢Height of 0.25 m
  ⊢Area of 4.687847838451e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (64505.92977, 43741.21282) m
  ⊢Height of 0.25 m
  ⊢Area of 3.182988604443e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (23971.21364, 8073.47522) m
  ⊢Height of 0.25 m
  ⊢Area of 4.342911293146e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (63903.57815, 51357.35459) m
  ⊢Height of 0.25 m
  ⊢Area of 3.411767837845e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (58425.55193, 28344.95238) m
  ⊢Height of 0.25 m
  ⊢Area of 3.304893702652e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (50493.75949, 11431.05595) m
  ⊢Height of 0.25 m
  ⊢Area of 4.381029256633e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (40104.00673, 19175.56906) m
  ⊢Height of 0.25 m
  ⊢Area of 4.285299073719e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (28106.78589, 27547.34516) m
  ⊢Height of 0.25 m
  ⊢Area of 6.615948736958e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (60641.83658, 66671.61639) m
  ⊢Height of 0.25 m
  ⊢Area of 5.418484390923e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (31853.93175, 44759.71555) m
  ⊢Height of 0.25 m
  ⊢Area of 5.567470362437e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 ⋮
 Floe{Float64}
  ⊢Centroid of (58954.24081, 22044.11618) m
  ⊢Height of 0.25 m
  ⊢Area of 1.972455771574e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (36584.66568, 15754.80437) m
  ⊢Height of 0.25 m
  ⊢Area of 1.336973617971e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (44899.3221, 80441.18759) m
  ⊢Height of 0.25 m
  ⊢Area of 7.240710711909e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (33070.79801, 53957.47145) m
  ⊢Height of 0.25 m
  ⊢Area of 1.87496265058e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (60463.1156, 71466.66776) m
  ⊢Height of 0.25 m
  ⊢Area of 5.163696357609e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (70535.48964, 6364.46493) m
  ⊢Height of 0.25 m
  ⊢Area of 4.035277098086e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (72458.79837, 15292.72875) m
  ⊢Height of 0.25 m
  ⊢Area of 2.589619855022e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (29966.9007, 4099.96515) m
  ⊢Height of 0.25 m
  ⊢Area of 5.282326224194e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (64567.07101, 74715.22959) m
  ⊢Height of 0.25 m
  ⊢Area of 5.509833503247e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

Model Creation

model = Model(; grid, ocean, atmos, domain, floes = floe_arr)

Model{Float64, ...}
 ⊢RegRectilinearGrid{Float64}
  ⊢x extent (0.0 to 100000.0) with 50 grid cells of size 2000.0 m
  ∟y extent (0.0 to 100000.0) with 50 grid cells of size 2000.0 m
 ⊢Domain
  ⊢Northern boundary of type PeriodicBoundary{North, Float64}
  ⊢Southern boundary of type PeriodicBoundary{South, Float64}
  ⊢Eastern boundary of type CollisionBoundary{East, Float64}
  ⊢Western boundary of type CollisionBoundary{West, Float64}
  ∟3-element TopograpahyElement{Float64} list
 ⊢Ocean{Float64}
  ⊢Vector fields of dimension (51, 51)
  ⊢Tracer fields of dimension (51, 51)
  ⊢Average u-velocity of: 0.0 m/s
  ⊢Average v-velocity of: -0.3 m/s
  ∟Average temperature of: 0.0 C
 ⊢Atmos{Float64}
  ⊢Vector fields of dimension (51, 51)
  ⊢Tracer fields of dimension (51, 51)
  ⊢Average u-velocity of: 0.0 m/s
  ⊢Average v-velocity of: 0.0 m/s
  ∟Average temperature of: 0.0 C
 ∟Floe List:
  ⊢Number of floes: 80
  ⊢Total floe area: 3.52276779805866e9
  ∟Average floe height: 0.25

Constants Creation

modulus = 1.5e3*(mean(sqrt.(floe_arr.area)) + minimum(sqrt.(floe_arr.area)))
consts = Constants(; E = modulus)

Constants{Float64}
  ⊢ ρo = 1027.0
  ⊢ ρa = 1.2
  ⊢ Cd_io = 0.003
  ⊢ Cd_ia = 0.001
  ⊢ Cd_ao = 0.00125
  ⊢ f = 0.00014
  ⊢ turnθ = 0.2617993877991494
  ⊢ L = 293000.0
  ⊢ k = 2.14
  ⊢ ν = 0.3
  ⊢ μ = 0.2
  ⊢ E = 1.5124944922310926e7

Settings Creation

Fracture Settings

fracture_settings = FractureSettings(;
        fractures_on = true,
        criteria = HiblerYieldCurve(floe_arr),
        Δt = 75,
        npieces = 3,
        deform_on = false,
)

FractureSettings{HiblerYieldCurve{Float64}}
  ⊢ fractures_on = true
  ⊢ criteria = HiblerYieldCurve{Float64}
  ⊢ Δt = 75
  ⊢ deform_on = false
  ⊢ npieces = 3

Ridge Raft Settings

ridgeraft_settings = RidgeRaftSettings(;
    ridge_raft_on = true,
    Δt = 150
)

RidgeRaftSettings{Float64}
  ⊢ ridge_raft_on = true
  ⊢ Δt = 150
  ⊢ 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

Output Creation

dir = "simple_strait"
init_fn, floe_fn = "simple_strait_init_state.jld2", "simple_strait_floes.jld2"
initwriter = InitialStateOutputWriter(dir = dir, filename = init_fn, overwrite = true)
floewriter = FloeOutputWriter(50, dir = dir, filename = floe_fn, overwrite = true)
writers = OutputWriters(initwriter, floewriter)

OutputWriters
  ⊢1 InitialStateOuputWriter(s)
  ⊢0 CheckpointOutputWriter(s)
  ⊢1 FloeOutputWriter(s)
  ∟0 GridOutputWriter(s)

Simulation Creation

simulation = Simulation(; model, consts, writers, Δt, nΔt,
    floe_settings, fracture_settings, ridgeraft_settings,
    verbose = true, rng = Xoshiro(1))

Simulation
  ⊢Timestep: 20 seconds
  ⊢Runtime: 5000 timesteps
  ⊢RNG: Xoshiro(0xfff0241072ddab67, 0xc53bc12f4c3f0b4e, 0x56d451780b2dd4ba, 0x50a4aa153d208dd8, 0x3649a58b3b63d5db)
  ⊢verbose: true
  ⊢model
  ⊢consts
  ⊢floe_settings
  ⊢collision_settings
  ∟ ...

Running the Simulation

run!(simulation)

Plotting the Simulation

plot_sim(joinpath(dir, floe_fn), joinpath(dir, init_fn), Δt, joinpath(dir, "simple_strait.mp4"))

Note

Note that this is just using the built-in basic plotting. However, it is easy to write your own plotting code. See the source code for a basic outline.

This page was generated using Literate.jl.