Moving Boundaries

This simulation has two moving boundaries, on the top and the bottom of the simulation. They push the floes inward towards the center of the domain. As the floes are pushed inward, they ridge and raft, gaining height and losing area. Users can also create shear moving boundaries, rather than the compression boundaries seen in this example.

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

User Inputs

const FT = Float64
const Lx = 1e5
const Ly = 1e5
const Δgrid = 2e3
const nfloes = 100
const concentrations = [1.0]
const hmean = 0.25
const Δh = 0.125
const Δt = 20
const nΔt = 1500;

Grid, Ocean, and Atmosphere instantiation

grid = RegRectilinearGrid(; x0 = 0.0, xf = Lx, y0 = 0.0, yf = Ly, Δx = Δgrid, Δy = Δgrid)
ocean = Ocean(; grid, u = 0.0, v = 0.0, temp = 0.0)
atmos = Atmos(; grid, u = 0.0, v = 0.0, temp = -1.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: -1.0 C

Domain creation

nboundary = MovingBoundary(North; grid, u = 0.0, v = -0.1)
sboundary = MovingBoundary(South; grid, u = 0.0, v = 0.1)
eboundary = PeriodicBoundary(East; grid)
wboundary = PeriodicBoundary(West; grid)

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

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

Floe creation

floe_settings = FloeSettings()
floe_generator = VoronoiTesselationFieldGenerator(; nfloes, concentrations, hmean, Δh)
floe_arr = initialize_floe_field(; generator = floe_generator, domain, floe_settings, rng = Xoshiro(1))

100-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 (95855.2678, 30617.93676) m
  ⊢Height of 0.32084 m
  ⊢Area of 5.340683388954e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (54279.50156, 5906.99736) m
  ⊢Height of 0.32212 m
  ⊢Area of 1.3309619199384e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (41873.72133, 56319.15856) m
  ⊢Height of 0.34801 m
  ⊢Area of 8.069710696096e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (6300.9056, 48641.41134) m
  ⊢Height of 0.20694 m
  ⊢Area of 8.464176301239e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (13679.72376, 85269.77911) m
  ⊢Height of 0.16617 m
  ⊢Area of 5.566189489785e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (15887.40277, 74903.53029) m
  ⊢Height of 0.31014 m
  ⊢Area of 8.001457075669e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (25677.36791, 97808.79568) m
  ⊢Height of 0.15723 m
  ⊢Area of 4.165318843636e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (69488.16244, 50360.74977) m
  ⊢Height of 0.18488 m
  ⊢Area of 1.8081555385474e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (30141.7917, 92544.59916) m
  ⊢Height of 0.17884 m
  ⊢Area of 7.823142319943e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (73310.94707, 32329.15148) m
  ⊢Height of 0.26112 m
  ⊢Area of 1.1671482155116e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 ⋮
 Floe{Float64}
  ⊢Centroid of (59787.08439, 80144.39138) m
  ⊢Height of 0.16783 m
  ⊢Area of 3.740521137656e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (27242.29102, 84916.37421) m
  ⊢Height of 0.22634 m
  ⊢Area of 7.780457000145e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (92132.08948, 23765.66271) m
  ⊢Height of 0.34643 m
  ⊢Area of 3.446492758606e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (90967.09841, 3535.09051) m
  ⊢Height of 0.18365 m
  ⊢Area of 1.283935295325e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (72663.51577, 70640.79198) m
  ⊢Height of 0.1695 m
  ⊢Area of 1.3380786033752e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (54676.46012, 51792.83996) m
  ⊢Height of 0.22213 m
  ⊢Area of 1.2427387696548e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (11473.61046, 57021.92083) m
  ⊢Height of 0.19646 m
  ⊢Area of 1.0521918803479e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (43352.04238, 79368.73698) m
  ⊢Height of 0.33844 m
  ⊢Area of 1.6911454198148e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (81700.77784, 23937.72996) m
  ⊢Height of 0.19744 m
  ⊢Area of 3.921164467272e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

You can also set the inital floe velocities manually like this:

floe_arr.u .= 0
floe_arr.v .= -0.01;

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 MovingBoundary{North, Float64}
  ⊢Southern boundary of type MovingBoundary{South, Float64}
  ⊢Eastern boundary of type PeriodicBoundary{East, Float64}
  ⊢Western boundary of type PeriodicBoundary{West, Float64}
  ∟0-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.0 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: -1.0 C
 ∟Floe List:
  ⊢Number of floes: 100
  ⊢Total floe area: 1.0e10
  ∟Average floe height: 0.24611

Output Writer Setup

dir = "moving_bounds"
init_fn, floe_fn = "moving_bounds_init_state.jld2", "moving_bounds.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 settings

modulus = 1.5e3*(mean(sqrt.(floe_arr.area)) + minimum(sqrt.(floe_arr.area)))
consts = Constants(; E = modulus, Cd_io = 0.0, f = 0.0, turnθ = 0.0)

ridgeraft_settings = RidgeRaftSettings(;
    ridge_raft_on = true,
    Δt = 150,
    domain_gain_probability = 0.5
)
weld_settings = WeldSettings(;
    weld_on = true,
    Δts = [150, 300, 600],  # weld at these specific timesteps
    Nxs = [2, 1, 1],  # split the domain into nx by ny sections and weld within each section
    Nys = [2, 2, 1],
)

WeldSettings{Float64}
  ⊢ weld_on = true
  ⊢ Δts = [600, 300, 150]
  ⊢ Nxs = [1, 1, 2]
  ⊢ Nys = [1, 2, 2]
  ⊢ min_weld_area = 1.0e6
  ⊢ max_weld_area = 2.0e9
  ⊢ welding_coeff = 150.0

Create Simulation

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

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, "moving_bounds.mp4"))

This page was generated using Literate.jl.