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{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 (89555.82733, 38616.18324) m
  ⊢Height of 0.25932 m
  ⊢Area of 1.2298298616685e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (52885.42069, 85826.75206) m
  ⊢Height of 0.18286 m
  ⊢Area of 1.7299318429344e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (35293.34871, 68756.55616) m
  ⊢Height of 0.29368 m
  ⊢Area of 3.537042552694e7 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.2723 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 (30195.73779, 79536.44682) m
  ⊢Height of 0.13275 m
  ⊢Area of 1.2303225783375e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (18397.89085, 7039.53287) m
  ⊢Height of 0.13982 m
  ⊢Area of 2.0558124326232e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (4941.54183, 93999.36068) m
  ⊢Height of 0.34323 m
  ⊢Area of 1.2058749330977e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (21160.30731, 48289.74816) m
  ⊢Height of 0.22165 m
  ⊢Area of 9.671165986542e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (35114.07558, 34585.78407) m
  ⊢Height of 0.34837 m
  ⊢Area of 1.9716468356947e8 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.16931 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 (43352.04238, 79368.73698) m
  ⊢Height of 0.3151 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 (18889.16779, 59086.92537) m
  ⊢Height of 0.35483 m
  ⊢Area of 7.062015356011e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (59363.03367, 15133.69564) m
  ⊢Height of 0.16173 m
  ⊢Area of 1.7714013447201e8 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.13614 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 (31617.05085, 72817.9764) m
  ⊢Height of 0.34364 m
  ⊢Area of 8.589904359324e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (35497.94221, 62229.04896) m
  ⊢Height of 0.22605 m
  ⊢Area of 6.780069623846e7 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (29858.25427, 50601.17202) m
  ⊢Height of 0.19664 m
  ⊢Area of 1.1881598836428e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (15471.61722, 38850.54989) m
  ⊢Height of 0.29884 m
  ⊢Area of 1.2210692342086e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (2783.48131, 58392.43519) m
  ⊢Height of 0.32945 m
  ⊢Area of 9.075220506041e7 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.25256

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)

sim is running!
0 timesteps
50 timesteps
100 timesteps
150 timesteps
200 timesteps
250 timesteps
300 timesteps
350 timesteps
400 timesteps
450 timesteps
500 timesteps
550 timesteps
600 timesteps
650 timesteps
700 timesteps
750 timesteps
800 timesteps
850 timesteps
900 timesteps
950 timesteps
1000 timesteps
1050 timesteps
1100 timesteps
1150 timesteps
1200 timesteps
1250 timesteps
1300 timesteps
1350 timesteps
1400 timesteps
1450 timesteps
1500 timesteps
1550 timesteps
1600 timesteps
1650 timesteps
1700 timesteps
1750 timesteps
1800 timesteps
1850 timesteps
1900 timesteps
1950 timesteps
2000 timesteps
2050 timesteps
2100 timesteps
2150 timesteps
2200 timesteps
2250 timesteps
2300 timesteps
2350 timesteps
2400 timesteps
2450 timesteps
2500 timesteps
2550 timesteps
2600 timesteps
2650 timesteps
2700 timesteps
2750 timesteps
2800 timesteps
2850 timesteps
2900 timesteps
2950 timesteps
3000 timesteps
3050 timesteps
3100 timesteps
3150 timesteps
3200 timesteps
3250 timesteps
3300 timesteps
3350 timesteps
3400 timesteps
3450 timesteps
3500 timesteps
3550 timesteps
3600 timesteps
3650 timesteps
3700 timesteps
3750 timesteps
3800 timesteps
3850 timesteps
3900 timesteps
3950 timesteps
4000 timesteps
4050 timesteps
4100 timesteps
4150 timesteps
4200 timesteps
4250 timesteps
4300 timesteps
4350 timesteps
4400 timesteps
4450 timesteps
4500 timesteps
4550 timesteps
4600 timesteps
4650 timesteps
4700 timesteps
4750 timesteps
4800 timesteps
4850 timesteps
4900 timesteps
4950 timesteps
5000 timesteps
sim done running!

Plotting the Simulation

plot_sim(joinpath(dir, floe_fn), joinpath(dir, init_fn), Δt, joinpath(dir, "moving_bounds.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.