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, StatisticsMaking 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 runGrid 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 mDomain 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 = [[[2e4, 6e4], [4e4, 6.5e4], [4.5e4, 6.5e4], [4.5e4, 6e4], [2e4, 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, island2, 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}
∟4-element TopograpahyElement{Float64} listOcean 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 CAtmos 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 CFloe Creation
floe_settings = FloeSettings(
subfloe_point_generator = MonteCarloPointsGenerator(; npoints = 100, ntries = 10, err = 0.1),
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)77-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 (69737.26796, 57783.60867) m
⊢Height of 0.25 m
⊢Area of 1.568338356922e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (47265.80984, 23976.84359) m
⊢Height of 0.25 m
⊢Area of 1.0347634530574e8 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (65067.19204, 17742.94609) m
⊢Height of 0.25 m
⊢Area of 1.117016619884e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (66742.00503, 39061.04421) m
⊢Height of 0.25 m
⊢Area of 2.684660021698e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (68383.16346, 18448.59657) m
⊢Height of 0.25 m
⊢Area of 2.415875239623e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (31470.19764, 57177.11967) m
⊢Height of 0.25 m
⊢Area of 3.258321777488e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (29345.06611, 64673.74641) m
⊢Height of 0.25 m
⊢Area of 2.172916663444e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (70532.84407, 62356.96392) m
⊢Height of 0.25 m
⊢Area of 1.639235628319e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (30797.91042, 51695.28907) m
⊢Height of 0.25 m
⊢Area of 1.340857161312e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (38489.15398, 88371.05636) m
⊢Height of 0.25 m
⊢Area of 7.818403899649e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
⋮
Floe{Float64}
⊢Centroid of (46460.02913, 9299.60876) m
⊢Height of 0.25 m
⊢Area of 2.247414444534e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (40663.83274, 54231.56218) m
⊢Height of 0.25 m
⊢Area of 2.335732274489e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (32494.19391, 73891.10361) m
⊢Height of 0.25 m
⊢Area of 1.8180361744848e8 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (59282.94947, 34785.22907) m
⊢Height of 0.25 m
⊢Area of 1.5144841055465e8 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (31061.83997, 8164.60547) m
⊢Height of 0.25 m
⊢Area of 9.55393048531e6 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (48792.71312, 91358.78515) m
⊢Height of 0.25 m
⊢Area of 4.339968138094e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (29363.70946, 4664.58539) m
⊢Height of 0.25 m
⊢Area of 2.983913459886e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (50565.58158, 98180.34667) m
⊢Height of 0.25 m
⊢Area of 2.022809821757e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (47946.1723, 42658.36046) m
⊢Height of 0.25 m
⊢Area of 9.421356557589e7 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}
∟4-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: 77
⊢Total floe area: 3.44509642781778e9
∟Average floe height: 0.25Constants 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.1973333357660856e7
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 = 10
)RidgeRaftSettings{Float64}
⊢ ridge_raft_on = true
⊢ Δt = 10
⊢ 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
nout = 50
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(nout, 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
@time 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!
27.784282 seconds (195.41 M allocations: 15.294 GiB, 6.16% gc time, 16.96% compilation time)Plotting the Simulation
plot_sim(joinpath(dir, floe_fn), joinpath(dir, init_fn), Δt, joinpath(dir, "simple_strait.mp4"))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.