Floes Bounded by Ocean Currents
This simulation has four open boundaries. The ocean is created such that there is a current in the middle of the domain that pushes the floes from left to right, and there are also currents at each of boundaries that push the floes back into the middle of the domain. The main point of this simulation is to highlight that the user can create a bounding box to control initial floe placement and that the initial set of floes does not need to span the entire domain.
using Subzero, CairoMakie, GeoInterfaceMakie
using JLD2, Random, StatisticsWARNING: Method definition deepcopy(Observables.Observable{T} where T) in module Makie at /home/runner/.julia/packages/Makie/frKXt/src/attributes.jl:51 overwritten in module MakieCore at /home/runner/.julia/packages/MakieCore/FHomX/src/attributes.jl:48.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
1 dependency had output during precompilation:
┌ Subzero → SubzeroMakieExt
│ [Output was shown above]
└User Inputs
const FT = Float64
const Lx = 1e5
const Ly = 1e5
const Δgrid = 2e3
const hmean = 0.25
const Δh = 0.0
const nfloes = 300
const concentrations = [0.4]
const Δt = 20
const nΔt = 10000;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 mDomain Creation
nboundary = OpenBoundary(North; grid)
sboundary = OpenBoundary(South; grid)
eboundary = OpenBoundary(East; grid)
wboundary = OpenBoundary(West; grid)
domain = Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary)Domain
⊢Northern boundary of type OpenBoundary{North, Float64}
⊢Southern boundary of type OpenBoundary{South, Float64}
⊢Eastern boundary of type OpenBoundary{East, Float64}
⊢Western boundary of type OpenBoundary{West, Float64}
∟0-element TopograpahyElement{Float64} listOcean Creation
Set ocean u and v velocities
ocean_uvels = zeros(FT, (grid.Nx + 1, grid.Ny + 1))
ocean_vvels = zeros(FT, (grid.Nx + 1, grid.Ny + 1))
for i in CartesianIndices(ocean_vvels)
r, c = Tuple(i)
if r ≤ 5 # set u vals
ocean_uvels[i] = 0.6
elseif r ≥ grid.Nx - 4
ocean_uvels[i] = -0.6
end
if c ≤ 5 # set v vals
ocean_vvels[i] = 0.6
elseif c ≥ grid.Ny - 4
ocean_vvels[i] = -0.6
end
end
ocean_uvels[10:40, 20:30] .= 0.3;Instatiate Ocean
ocean = Ocean(;
grid,
u = ocean_uvels,
v = ocean_vvels,
temp = 0,
)Ocean{Float64}
⊢Vector fields of dimension (51, 51)
⊢Tracer fields of dimension (51, 51)
⊢Average u-velocity of: 0.02757 m/s
⊢Average v-velocity of: -0.01176 m/s
∟Average temperature of: 0.0 CWe can then plot the ocean for a better understanding of the above setup.
fig = Figure();
ax1 = Axis(fig[1, 1]; title = "Ocean U-Velocities [m/s]", xticklabelrotation = pi/4)
ax2 = Axis(fig[2, 1]; title = "Ocean V-Velocities [m/s]", xticklabelrotation = pi/4)
xs = grid.x0:grid.Δx:grid.xf
ys = grid.y0:grid.Δy:grid.yf
u_hm = heatmap!(ax1, xs, ys, ocean.u)
Colorbar(fig[1, end+1], u_hm)
v_hm = heatmap!(ax2, xs, ys, ocean.v)
Colorbar(fig[2, end+1], v_hm)
resize_to_layout!(fig)
fig
Atmosphere Creation
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 CFloe Creation - bound floes within smaller part of the domain
floe_settings = FloeSettings(; subfloe_point_generator = SubGridPointsGenerator(; grid, npoint_per_cell = 2))
floe_generator = VoronoiTesselationFieldGenerator(; nfloes, concentrations, hmean, Δh)
floe_bounds = Subzero.make_polygon([[[0.1Lx, 0.1Ly], [0.1Lx, 0.9Ly], [0.9Lx, 0.9Ly], [0.9Lx, 0.1Ly], [0.1Lx, 0.1Ly]]], FT)
floe_arr = initialize_floe_field(FT; generator = floe_generator, domain, floe_bounds, rng = Xoshiro(1), floe_settings)307-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 (75325.88495, 40206.81959) m
⊢Height of 0.25 m
⊢Area of 1.045970665057e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (62072.25518, 41520.62709) m
⊢Height of 0.25 m
⊢Area of 1.455639393235e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (45066.363, 66914.48432) m
⊢Height of 0.25 m
⊢Area of 1.263261704127e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (20374.99757, 89306.77762) m
⊢Height of 0.25 m
⊢Area of 5.4339968306e6 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (29651.48144, 35958.23254) m
⊢Height of 0.25 m
⊢Area of 7.51792214909e6 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (84521.78057, 78134.70751) m
⊢Height of 0.25 m
⊢Area of 1.271087350074e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (15097.97839, 35744.96156) m
⊢Height of 0.25 m
⊢Area of 2.58324577994e6 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (57337.40069, 17177.33082) m
⊢Height of 0.25 m
⊢Area of 1.255752441478e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (87506.85285, 36330.00609) m
⊢Height of 0.25 m
⊢Area of 1.82656136452e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (39966.59921, 70573.02835) m
⊢Height of 0.25 m
⊢Area of 1.651632563835e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
⋮
Floe{Float64}
⊢Centroid of (23886.47948, 38272.67642) m
⊢Height of 0.25 m
⊢Area of 1.436282763712e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (66674.68703, 24956.90824) m
⊢Height of 0.25 m
⊢Area of 1.827558829176e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (83013.87468, 55155.07816) m
⊢Height of 0.25 m
⊢Area of 9.32990758421e6 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (68002.39549, 71805.82773) m
⊢Height of 0.25 m
⊢Area of 1.3343739539e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (49732.86916, 72562.30366) m
⊢Height of 0.25 m
⊢Area of 1.11735759329e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (57401.52829, 35518.89254) m
⊢Height of 0.25 m
⊢Area of 1.017558195603e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (53563.85422, 15312.54291) m
⊢Height of 0.25 m
⊢Area of 1.371207139475e7 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (56580.01033, 23321.68599) m
⊢Height of 0.25 m
⊢Area of 9.26794143818e6 m^2
∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)
Floe{Float64}
⊢Centroid of (13487.15251, 25025.90748) m
⊢Height of 0.25 m
⊢Area of 4.99539925511e6 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 OpenBoundary{North, Float64}
⊢Southern boundary of type OpenBoundary{South, Float64}
⊢Eastern boundary of type OpenBoundary{East, Float64}
⊢Western boundary of type OpenBoundary{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.02757 m/s
⊢Average v-velocity of: -0.01176 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: 307
⊢Total floe area: 2.56240695909283e9
∟Average floe height: 0.25Output Writer Creation
dir = "forcing_contained_floes"
init_fn, floe_fn = "contained_floes_init_state.jld2", "contained_floes.jld2"
initwriter = InitialStateOutputWriter(filename = init_fn, dir = dir, overwrite = true)
floewriter = FloeOutputWriter(50, filename = floe_fn, dir = dir, overwrite = true)
writers = OutputWriters(initwriter, floewriter)OutputWriters
⊢1 InitialStateOuputWriter(s)
⊢0 CheckpointOutputWriter(s)
⊢1 FloeOutputWriter(s)
∟0 GridOutputWriter(s)Simulation Creation
modulus = 1.5e3*(mean(sqrt.(floe_arr.area)) + minimum(sqrt.(floe_arr.area)))
consts = Constants(E = modulus)
simulation = Simulation(;
model = model,
consts = consts,
Δt = Δt,
nΔt = nΔt,
verbose = true,
writers = writers,
rng = Xoshiro(1),
floe_settings = floe_settings,
)Simulation
⊢Timestep: 20 seconds
⊢Runtime: 10000 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
5050 timesteps
5100 timesteps
5150 timesteps
5200 timesteps
5250 timesteps
5300 timesteps
5350 timesteps
5400 timesteps
5450 timesteps
5500 timesteps
5550 timesteps
5600 timesteps
5650 timesteps
5700 timesteps
5750 timesteps
5800 timesteps
5850 timesteps
5900 timesteps
5950 timesteps
6000 timesteps
6050 timesteps
6100 timesteps
6150 timesteps
6200 timesteps
6250 timesteps
6300 timesteps
6350 timesteps
6400 timesteps
6450 timesteps
6500 timesteps
6550 timesteps
6600 timesteps
6650 timesteps
6700 timesteps
6750 timesteps
6800 timesteps
6850 timesteps
6900 timesteps
6950 timesteps
7000 timesteps
7050 timesteps
7100 timesteps
7150 timesteps
7200 timesteps
7250 timesteps
7300 timesteps
7350 timesteps
7400 timesteps
7450 timesteps
7500 timesteps
7550 timesteps
7600 timesteps
7650 timesteps
7700 timesteps
7750 timesteps
7800 timesteps
7850 timesteps
7900 timesteps
7950 timesteps
8000 timesteps
8050 timesteps
8100 timesteps
8150 timesteps
8200 timesteps
8250 timesteps
8300 timesteps
8350 timesteps
8400 timesteps
8450 timesteps
8500 timesteps
8550 timesteps
8600 timesteps
8650 timesteps
8700 timesteps
8750 timesteps
8800 timesteps
8850 timesteps
8900 timesteps
8950 timesteps
9000 timesteps
9050 timesteps
9100 timesteps
9150 timesteps
9200 timesteps
9250 timesteps
9300 timesteps
9350 timesteps
9400 timesteps
9450 timesteps
9500 timesteps
9550 timesteps
9600 timesteps
9650 timesteps
9700 timesteps
9750 timesteps
9800 timesteps
9850 timesteps
9900 timesteps
9950 timesteps
10000 timesteps
sim done running!Plotting the Simulation
plot_sim(joinpath(dir, floe_fn), joinpath(dir, init_fn), Δt, joinpath(dir, "contained_floes.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.