Shear Flow with Periodic Boundaries

This simulation creates a periodic simulation with an ocean shear flow current. It is very similar to the example created in the tutorial. The tutorial explains the setup in greater detail. It is a good starter simulation for understanding how the code works and playing around with a basic setup.

This simulation shows how to setup a basic domain where all four boundaries are periodic. It also creates non-constant, shear ocean field. The ocean u-velocities are zero at the minimum and maximum y-extents (constant across x-values). The u values increase from the top and bottom of the domain towards the center to 0.5m/s, again constant across x-values, varying with y-values.

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

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 run

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 = PeriodicBoundary(North; grid)
sboundary = PeriodicBoundary(South; grid)
eboundary = PeriodicBoundary(East; grid)
wboundary = PeriodicBoundary(West; grid)
domain = Domain(; north = nboundary, south = sboundary, east = eboundary, west = wboundary)

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

Ocean Creation

half_nx = fld(grid.Nx, 2)
u_vec = [range(0, 0.5, length = half_nx); 0.5; range(0.5, 0, length = half_nx)]
u_row = transpose(u_vec)
uvels = repeat(u_row, outer = (grid.Ny + 1, 1))
ocean = Ocean(; u = uvels, v = 0, temp = 0, grid)

Ocean{Float64}
  ⊢Vector fields of dimension (51, 51)
  ⊢Tracer fields of dimension (51, 51)
  ⊢Average u-velocity of: 0.2549 m/s
  ⊢Average v-velocity of: 0.0 m/s
  ∟Average temperature of: 0.0 C

Atmosphere Creation

atmos = Atmos(; u = 0.0, v = 0.0, temp = -1.0, grid)

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

floe_settings = FloeSettings(subfloe_point_generator = SubGridPointsGenerator(; grid, npoint_per_cell = 2))
floe_generator = VoronoiTesselationFieldGenerator(; nfloes = 50, concentrations = [0.75], hmean, Δh)
floe_arr = initialize_floe_field(FT; generator = floe_generator, domain, rng = Xoshiro(1), floe_settings)

52-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 (62995.82957, 43616.14049) m
  ⊢Height of 0.25 m
  ⊢Area of 1.1690494652847e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (51929.76984, 11885.98576) m
  ⊢Height of 0.25 m
  ⊢Area of 3.6363702964023e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

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

 Floe{Float64}
  ⊢Centroid of (75571.0654, 63500.51622) m
  ⊢Height of 0.25 m
  ⊢Area of 3.2777213000065e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (12375.00366, 3874.26379) m
  ⊢Height of 0.25 m
  ⊢Area of 1.3325555598103e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (3218.84526, 47216.58624) m
  ⊢Height of 0.25 m
  ⊢Area of 1.1919235193645e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (4202.05882, 19564.111) m
  ⊢Height of 0.25 m
  ⊢Area of 1.1771165057383e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (11656.4435, 93622.62301) m
  ⊢Height of 0.25 m
  ⊢Area of 2.163885674736e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (88178.69607, 13424.52278) m
  ⊢Height of 0.25 m
  ⊢Area of 1.6783420286495e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (57001.03872, 63226.2207) m
  ⊢Height of 0.25 m
  ⊢Area of 1.060794008597e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 ⋮
 Floe{Float64}
  ⊢Centroid of (32195.79981, 36824.15059) m
  ⊢Height of 0.25 m
  ⊢Area of 1.9236323809466e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (66409.08867, 50729.61018) m
  ⊢Height of 0.25 m
  ⊢Area of 1.4576742987025e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

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

 Floe{Float64}
  ⊢Centroid of (95773.67734, 46711.12264) m
  ⊢Height of 0.25 m
  ⊢Area of 1.7966044078233e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (94894.16106, 65292.29207) m
  ⊢Height of 0.25 m
  ⊢Area of 1.6148405105348e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

 Floe{Float64}
  ⊢Centroid of (82415.80025, 17595.45063) m
  ⊢Height of 0.25 m
  ⊢Area of 1.1417469701719e8 m^2
  ∟Velocity of (u, v, ξ) of (0.0, 0.0, 0.0) in (m/s, m/s, rad/s)

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

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

 Floe{Float64}
  ⊢Centroid of (80274.89211, 73445.83944) m
  ⊢Height of 0.25 m
  ⊢Area of 4.023340890352e7 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 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.2549 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: 52
  ⊢Total floe area: 7.50266794432454e9
  ∟Average floe height: 0.25

Constants 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 = 2.0173906502770815e7

Output Creation

dir = "shear_flow"
init_fn, floe_fn = "shear_flow_init_state.jld2", "shear_flow_floes.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 Creation

simulation = Simulation(; model, consts, writers, Δt, nΔt, floe_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

run!(simulation)

Plotting the Simulation

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