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, Statistics

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 m

Domain 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} list

Ocean 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 C

We 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 C

Floe 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]]])
floe_arr = initialize_floe_field(FT; generator = floe_generator, domain, floe_bounds, rng = Xoshiro(1), floe_settings)

320-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 (57343.94252, 56153.45354) m
  ⊢Height of 0.25 m
  ⊢Area of 8.20740008326e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (52282.64229, 78921.63708) m
  ⊢Height of 0.25 m
  ⊢Area of 4.0771340159e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (56566.74808, 78374.7878) m
  ⊢Height of 0.25 m
  ⊢Area of 6.71268591182e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (76061.22462, 52722.59485) m
  ⊢Height of 0.25 m
  ⊢Area of 7.95333047605e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

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

 Floe{Float64}
  ⊢Centroid of (73037.57458, 67935.344) m
  ⊢Height of 0.25 m
  ⊢Area of 4.88912094424e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

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

 Floe{Float64}
  ⊢Centroid of (52182.3424, 69284.88788) m
  ⊢Height of 0.25 m
  ⊢Area of 5.70493841252e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (43325.68348, 67845.786) m
  ⊢Height of 0.25 m
  ⊢Area of 4.14183799699e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (26933.15605, 14593.5524) m
  ⊢Height of 0.25 m
  ⊢Area of 8.0584996622e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 ⋮
 Floe{Float64}
  ⊢Centroid of (83492.78126, 74177.95234) m
  ⊢Height of 0.25 m
  ⊢Area of 2.23634838744e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

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

 Floe{Float64}
  ⊢Centroid of (36930.80809, 29634.47731) m
  ⊢Height of 0.25 m
  ⊢Area of 8.58452100931e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (59076.23254, 33396.96246) m
  ⊢Height of 0.25 m
  ⊢Area of 4.33932101366e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

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

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

 Floe{Float64}
  ⊢Centroid of (29165.30001, 40237.9128) m
  ⊢Height of 0.25 m
  ⊢Area of 7.94405959064e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (25204.28209, 18098.60202) m
  ⊢Height of 0.25 m
  ⊢Area of 9.51716491854e6 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (31761.18771, 80197.95193) m
  ⊢Height of 0.25 m
  ⊢Area of 1.675365799582e7 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: 320
  ⊢Total floe area: 2.56031838868727e9
  ∟Average floe height: 0.25

Output 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)

Plotting the Simulation

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