Wind- and convection-driven mixing in an ocean surface boundary layer
This example simulates mixing by three-dimensional turbulence in an ocean surface boundary layer driven by atmospheric winds and convection. It demonstrates:
- How to set-up a grid with varying spacing in the vertical direction
- How to use the
SeawaterBuoyancymodel for buoyancy withTEOS10EquationOfState. - How to use a turbulence closure for large eddy simulation.
- How to use a function to impose a boundary condition.
Install dependencies
First let's make sure we have all required packages installed.
using Pkg
pkg"add Oceananigans, CairoMakie, SeawaterPolynomials"We start by importing all of the packages and functions that we'll need for this example.
using Oceananigans
using Oceananigans.Units
using Random
using Printf
using CairoMakie
using SeawaterPolynomials.TEOS10: TEOS10EquationOfStateThe grid
We use 128²×64 grid points with 1 m grid spacing in the horizontal and varying spacing in the vertical, with higher resolution closer to the surface. Here we use a stretching function for the vertical nodes that maintains relatively constant vertical spacing in the mixed layer, which is desirable from a numerical standpoint:
Nx = Ny = 128 # number of points in each of horizontal directions
Nz = 64 # number of points in the vertical direction
Lx = Ly = 128 # (m) domain horizontal extents
Lz = 64 # (m) domain depth
refinement = 1.2 # controls spacing near surface (higher means finer spaced)
stretching = 12 # controls rate of stretching at bottom
# Normalized height ranging from 0 to 1
h(k) = (k - 1) / Nz
# Linear near-surface generator
ζ₀(k) = 1 + (h(k) - 1) / refinement
# Bottom-intensified stretching function
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))
# Generating function
z_interfaces(k) = Lz * (ζ₀(k) * Σ(k) - 1)
grid = RectilinearGrid(GPU(),
size = (Nx, Nx, Nz),
x = (0, Lx),
y = (0, Ly),
z = z_interfaces)128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CUDAGPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 128.0) regularly spaced with Δx=1.0
├── Periodic y ∈ [0.0, 128.0) regularly spaced with Δy=1.0
└── Bounded z ∈ [-64.0, 0.0] variably spaced with min(Δz)=0.833413, max(Δz)=1.96618We plot vertical spacing versus depth to inspect the prescribed grid stretching:
fig = Figure(size=(1200, 800))
ax = Axis(fig[1, 1], ylabel = "z (m)", xlabel = "Vertical spacing (m)")
lines!(ax, zspacings(grid, Center()))
scatter!(ax, zspacings(grid, Center()))
fig
Buoyancy that depends on temperature and salinity
We use the SeawaterBuoyancy model with the TEOS10 equation of state,
ρₒ = 1026 # kg m⁻³, average density at the surface of the world ocean
equation_of_state = TEOS10EquationOfState(reference_density=ρₒ)
buoyancy = SeawaterBuoyancy(; equation_of_state)SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: BoussinesqEquationOfState{Float64}Boundary conditions
We calculate the surface temperature flux associated with surface cooling of 200 W m⁻², reference density ρₒ, and heat capacity cᴾ,
Q = 200 # W m⁻², surface _heat_ flux
cᴾ = 3991 # J K⁻¹ kg⁻¹, typical heat capacity for seawater
Jᵀ = Q / (ρₒ * cᴾ) # K m s⁻¹, surface _temperature_ flux4.884283985946938e-5Finally, we impose a temperature gradient dTdz both initially (see "Initial conditions" section below) and at the bottom of the domain, culminating in the boundary conditions on temperature,
dTdz = 0.01 # K m⁻¹
T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᵀ),
bottom = GradientBoundaryCondition(dTdz))Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: GradientBoundaryCondition: 0.01
├── top: FluxBoundaryCondition: 4.88428e-5
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)Note that a positive temperature flux at the surface of the ocean implies cooling. This is because a positive temperature flux implies that temperature is fluxed upwards, out of the ocean.
For the velocity field, we imagine a wind blowing over the ocean surface with an average velocity at 10 meters u₁₀, and use a drag coefficient cᴰ to estimate the kinematic stress (that is, stress divided by density) exerted by the wind on the ocean:
u₁₀ = 10 # m s⁻¹, average wind velocity 10 meters above the ocean
cᴰ = 2e-3 # dimensionless drag coefficient
ρₐ = 1.2 # kg m⁻³, approximate average density of air at sea-level
τx = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²-0.00023391812865497074The boundary conditions on u are thus
u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τx))Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: -0.000233918
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)For salinity, S, we impose an evaporative flux of the form
@inline Jˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S # [salinity unit] m s⁻¹where S is salinity. We use an evaporation rate of 1 millimeter per hour,
evaporation_rate = 1e-3 / hour # m s⁻¹2.7777777777777776e-7We build the Flux evaporation BoundaryCondition with the function Jˢ, indicating that Jˢ depends on salinity S and passing the parameter evaporation_rate,
evaporation_bc = FluxBoundaryCondition(Jˢ, field_dependencies=:S, parameters=evaporation_rate)FluxBoundaryCondition: ContinuousBoundaryFunction Jˢ at (Nothing, Nothing, Nothing)The full salinity boundary conditions are
S_bcs = FieldBoundaryConditions(top=evaporation_bc)Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: ContinuousBoundaryFunction Jˢ at (Nothing, Nothing, Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)Model instantiation
We fill in the final details of the model here, i.e., Coriolis forces, and the AnisotropicMinimumDissipation closure for large eddy simulation to model the effect of turbulent motions at scales smaller than the grid scale that are not explicitly resolved.
model = NonhydrostaticModel(; grid, buoyancy,
tracers = (:T, :S),
coriolis = FPlane(f=1e-4),
closure = AnisotropicMinimumDissipation(),
boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs))NonhydrostaticModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CUDAGPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: (T, S)
├── closure: AnisotropicMinimumDissipation{ExplicitTimeDiscretization, @NamedTuple{T::Float64, S::Float64}, Float64, Nothing}
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
└── coriolis: FPlane{Float64}(f=0.0001)Note: To use the Smagorinsky-Lilly turbulence closure (with a constant model coefficient) rather than AnisotropicMinimumDissipation, use closure = SmagorinskyLilly() in the model constructor.
Initial conditions
Our initial condition for temperature consists of a linear stratification superposed with random noise damped at the walls, while our initial condition for velocity consists only of random noise.
# Random noise damped at top and bottom
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise
# Temperature initial condition: a stable density gradient with random noise superposed.
Tᵢ(x, y, z) = 20 + dTdz * z + dTdz * model.grid.Lz * 1e-6 * Ξ(z)
# Velocity initial condition: random noise scaled by the friction velocity.
uᵢ(x, y, z) = sqrt(abs(τx)) * 1e-3 * Ξ(z)
# `set!` the `model` fields using functions or constants:
set!(model, u=uᵢ, w=uᵢ, T=Tᵢ, S=35)Setting up a simulation
We set-up a simulation with an initial time-step of 10 seconds that stops at 2 hours, with adaptive time-stepping and progress printing.
simulation = Simulation(model, Δt=10, stop_time=2hours)Simulation of NonhydrostaticModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 10 seconds
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 2 hours
├── Stop iteration: Inf
├── Wall time limit: Inf
├── Minimum relative step: 0.0
├── Callbacks: OrderedDict with 4 entries:
│ ├── stop_time_exceeded => 4
│ ├── stop_iteration_exceeded => -
│ ├── wall_time_limit_exceeded => e
│ └── nan_checker => }
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entriesThe TimeStepWizard helps ensure stable time-stepping with a Courant-Freidrichs-Lewy (CFL) number of 1.0.
wizard = TimeStepWizard(cfl=1, max_change=1.1, max_Δt=1minute)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))Callback of TimeStepWizard(cfl=1.0, max_Δt=60.0, min_Δt=0.0) on IterationInterval(10)Nice progress messaging is helpful:
# Print a progress message
progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, max(|w|) = %.1e ms⁻¹, wall time: %s\n",
iteration(sim), prettytime(sim), prettytime(sim.Δt),
maximum(abs, sim.model.velocities.w), prettytime(sim.run_wall_time))
add_callback!(simulation, progress_message, IterationInterval(40))We then set up the simulation:
Output
We use the JLD2Writer to save $x, z$ slices of the velocity fields, tracer fields, and eddy diffusivities. The prefix keyword argument to JLD2Writer indicates that output will be saved in ocean_wind_mixing_and_convection.jld2.
# Create a NamedTuple with eddy viscosity
eddy_viscosity = (; νₑ = model.diffusivity_fields.νₑ)
filename = "ocean_wind_mixing_and_convection"
simulation.output_writers[:slices] =
JLD2Writer(model, merge(model.velocities, model.tracers, eddy_viscosity),
filename = filename * ".jld2",
indices = (:, grid.Ny/2, :),
schedule = TimeInterval(1minute),
overwrite_existing = true)JLD2Writer scheduled on TimeInterval(1 minute):
├── filepath: ocean_wind_mixing_and_convection.jld2
├── 6 outputs: (u, v, w, T, S, νₑ)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 63.4 KiBWe're ready:
run!(simulation)[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|w|) = 1.3e-05 ms⁻¹, wall time: 0 seconds
[ Info: ... simulation initialization complete (28.374 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (6.032 seconds).
Iteration: 0040, time: 7 minutes, Δt: 8.309 seconds, max(|w|) = 3.1e-05 ms⁻¹, wall time: 35.271 seconds
Iteration: 0080, time: 11.181 minutes, Δt: 4.842 seconds, max(|w|) = 7.5e-03 ms⁻¹, wall time: 35.832 seconds
Iteration: 0120, time: 14 minutes, Δt: 4.152 seconds, max(|w|) = 2.5e-02 ms⁻¹, wall time: 36.292 seconds
Iteration: 0160, time: 16.704 minutes, Δt: 4.274 seconds, max(|w|) = 2.5e-02 ms⁻¹, wall time: 36.814 seconds
Iteration: 0200, time: 19.277 minutes, Δt: 4.062 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 37.283 seconds
Iteration: 0240, time: 21.842 minutes, Δt: 3.819 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 37.777 seconds
Iteration: 0280, time: 24.243 minutes, Δt: 3.620 seconds, max(|w|) = 2.5e-02 ms⁻¹, wall time: 38.256 seconds
Iteration: 0320, time: 26.616 minutes, Δt: 3.489 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 38.782 seconds
Iteration: 0360, time: 28.974 minutes, Δt: 3.418 seconds, max(|w|) = 2.8e-02 ms⁻¹, wall time: 39.239 seconds
Iteration: 0400, time: 31.168 minutes, Δt: 3.309 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 39.717 seconds
Iteration: 0440, time: 33.274 minutes, Δt: 3.320 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 40.247 seconds
Iteration: 0480, time: 35.364 minutes, Δt: 3.278 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 40.743 seconds
Iteration: 0520, time: 37.434 minutes, Δt: 3.184 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 41.236 seconds
Iteration: 0560, time: 39.454 minutes, Δt: 2.955 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 41.698 seconds
Iteration: 0600, time: 41.300 minutes, Δt: 2.991 seconds, max(|w|) = 2.8e-02 ms⁻¹, wall time: 42.189 seconds
Iteration: 0640, time: 43.150 minutes, Δt: 2.988 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 42.661 seconds
Iteration: 0680, time: 44.993 minutes, Δt: 2.918 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 43.150 seconds
Iteration: 0720, time: 46.860 minutes, Δt: 2.838 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 43.613 seconds
Iteration: 0760, time: 48.771 minutes, Δt: 2.927 seconds, max(|w|) = 4.2e-02 ms⁻¹, wall time: 44.105 seconds
Iteration: 0800, time: 50.598 minutes, Δt: 2.749 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 44.576 seconds
Iteration: 0840, time: 52.355 minutes, Δt: 2.728 seconds, max(|w|) = 3.7e-02 ms⁻¹, wall time: 45.074 seconds
Iteration: 0880, time: 54.141 minutes, Δt: 2.823 seconds, max(|w|) = 4.2e-02 ms⁻¹, wall time: 45.539 seconds
Iteration: 0920, time: 55.948 minutes, Δt: 2.672 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 46.019 seconds
Iteration: 0960, time: 57.650 minutes, Δt: 2.661 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 46.478 seconds
Iteration: 1000, time: 59.416 minutes, Δt: 2.724 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 46.970 seconds
Iteration: 1040, time: 1.020 hours, Δt: 2.757 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 47.439 seconds
Iteration: 1080, time: 1.050 hours, Δt: 2.759 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 47.920 seconds
Iteration: 1120, time: 1.079 hours, Δt: 2.631 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 48.380 seconds
Iteration: 1160, time: 1.108 hours, Δt: 2.763 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 48.874 seconds
Iteration: 1200, time: 1.138 hours, Δt: 2.678 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 49.336 seconds
Iteration: 1240, time: 1.167 hours, Δt: 2.644 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 49.832 seconds
Iteration: 1280, time: 1.194 hours, Δt: 2.607 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 50.305 seconds
Iteration: 1320, time: 1.222 hours, Δt: 2.615 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 50.809 seconds
Iteration: 1360, time: 1.250 hours, Δt: 2.622 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 51.264 seconds
Iteration: 1400, time: 1.278 hours, Δt: 2.551 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 51.763 seconds
Iteration: 1440, time: 1.306 hours, Δt: 2.513 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 52.262 seconds
Iteration: 1480, time: 1.333 hours, Δt: 2.560 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 52.719 seconds
Iteration: 1520, time: 1.360 hours, Δt: 2.521 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 53.220 seconds
Iteration: 1560, time: 1.387 hours, Δt: 2.461 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 53.689 seconds
Iteration: 1600, time: 1.414 hours, Δt: 2.468 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 54.171 seconds
Iteration: 1640, time: 1.438 hours, Δt: 2.343 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 54.631 seconds
Iteration: 1680, time: 1.464 hours, Δt: 2.246 seconds, max(|w|) = 5.9e-02 ms⁻¹, wall time: 55.112 seconds
Iteration: 1720, time: 1.489 hours, Δt: 2.475 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 55.584 seconds
Iteration: 1760, time: 1.516 hours, Δt: 2.415 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 56.068 seconds
Iteration: 1800, time: 1.541 hours, Δt: 2.425 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 56.543 seconds
Iteration: 1840, time: 1.568 hours, Δt: 2.520 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 57.043 seconds
Iteration: 1880, time: 1.596 hours, Δt: 2.404 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 57.503 seconds
Iteration: 1920, time: 1.621 hours, Δt: 2.411 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 58.003 seconds
Iteration: 1960, time: 1.647 hours, Δt: 2.351 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 58.450 seconds
Iteration: 2000, time: 1.673 hours, Δt: 2.404 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 58.953 seconds
Iteration: 2040, time: 1.700 hours, Δt: 2.468 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 59.400 seconds
Iteration: 2080, time: 1.727 hours, Δt: 2.417 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 59.932 seconds
Iteration: 2120, time: 1.753 hours, Δt: 2.327 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 1.007 minutes
Iteration: 2160, time: 1.779 hours, Δt: 2.384 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 1.015 minutes
Iteration: 2200, time: 1.805 hours, Δt: 2.386 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 1.023 minutes
Iteration: 2240, time: 1.831 hours, Δt: 2.385 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 1.031 minutes
Iteration: 2280, time: 1.856 hours, Δt: 2.360 seconds, max(|w|) = 6.4e-02 ms⁻¹, wall time: 1.039 minutes
Iteration: 2320, time: 1.882 hours, Δt: 2.367 seconds, max(|w|) = 7.6e-02 ms⁻¹, wall time: 1.047 minutes
Iteration: 2360, time: 1.907 hours, Δt: 2.344 seconds, max(|w|) = 7.1e-02 ms⁻¹, wall time: 1.055 minutes
Iteration: 2400, time: 1.933 hours, Δt: 2.414 seconds, max(|w|) = 6.4e-02 ms⁻¹, wall time: 1.063 minutes
Iteration: 2440, time: 1.958 hours, Δt: 2.300 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 1.071 minutes
Iteration: 2480, time: 1.983 hours, Δt: 2.292 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 1.079 minutes
[ Info: Simulation is stopping after running for 1.085 minutes.
[ Info: Simulation time 2 hours equals or exceeds stop time 2 hours.
Turbulence visualization
We animate the data saved in ocean_wind_mixing_and_convection.jld2. We prepare for animating the flow by loading the data into FieldTimeSeries and defining functions for computing colorbar limits.
filepath = filename * ".jld2"
time_series = (w = FieldTimeSeries(filepath, "w"),
T = FieldTimeSeries(filepath, "T"),
S = FieldTimeSeries(filepath, "S"),
νₑ = FieldTimeSeries(filepath, "νₑ"))(w = 128×1×65×121 FieldTimeSeries{InMemory} located at (Center, Center, Face) of w at ocean_wind_mixing_and_convection.jld2
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: w
└── data: 134×1×71×121 OffsetArray(::Array{Float64, 4}, -2:131, 64:64, -2:68, 1:121) with eltype Float64 with indices -2:131×64:64×-2:68×1:121
└── max=0.0407118, min=-0.0522771, mean=-4.45716e-5, T = 128×1×64×121 FieldTimeSeries{InMemory} located at (Center, Center, Center) of T at ocean_wind_mixing_and_convection.jld2
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: T
└── data: 134×1×70×121 OffsetArray(::Array{Float64, 4}, -2:131, 64:64, -2:67, 1:121) with eltype Float64 with indices -2:131×64:64×-2:67×1:121
└── max=19.9958, min=0.0, mean=18.5885, S = 128×1×64×121 FieldTimeSeries{InMemory} located at (Center, Center, Center) of S at ocean_wind_mixing_and_convection.jld2
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: S
└── data: 134×1×70×121 OffsetArray(::Array{Float64, 4}, -2:131, 64:64, -2:67, 1:121) with eltype Float64 with indices -2:131×64:64×-2:67×1:121
└── max=35.0128, min=0.0, mean=33.0007, νₑ = 128×1×64×121 FieldTimeSeries{InMemory} located at (Center, Center, Center) of νₑ at ocean_wind_mixing_and_convection.jld2
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: νₑ
└── data: 134×1×70×121 OffsetArray(::Array{Float64, 4}, -2:131, 64:64, -2:67, 1:121) with eltype Float64 with indices -2:131×64:64×-2:67×1:121
└── max=0.0190465, min=0.0, mean=0.000405101)We start the animation at $t = 10$ minutes since things are pretty boring till then:
times = time_series.w.times
intro = searchsortedfirst(times, 10minutes)11We are now ready to animate using Makie. We use Makie's Observable to animate the data. To dive into how Observables work we refer to Makie.jl's Documentation.
n = Observable(intro)
wₙ = @lift time_series.w[$n]
Tₙ = @lift time_series.T[$n]
Sₙ = @lift time_series.S[$n]
νₑₙ = @lift time_series.νₑ[$n]
fig = Figure(size = (1800, 900))
axis_kwargs = (xlabel="x (m)",
ylabel="z (m)",
aspect = AxisAspect(grid.Lx/grid.Lz),
limits = ((0, grid.Lx), (-grid.Lz, 0)))
ax_w = Axis(fig[2, 1]; title = "Vertical velocity", axis_kwargs...)
ax_T = Axis(fig[2, 3]; title = "Temperature", axis_kwargs...)
ax_S = Axis(fig[3, 1]; title = "Salinity", axis_kwargs...)
ax_νₑ = Axis(fig[3, 3]; title = "Eddy viscocity", axis_kwargs...)
title = @lift @sprintf("t = %s", prettytime(times[$n]))
wlims = (-0.05, 0.05)
Tlims = (19.7, 19.99)
Slims = (35, 35.005)
νₑlims = (1e-6, 5e-3)
hm_w = heatmap!(ax_w, wₙ; colormap = :balance, colorrange = wlims)
Colorbar(fig[2, 2], hm_w; label = "m s⁻¹")
hm_T = heatmap!(ax_T, Tₙ; colormap = :thermal, colorrange = Tlims)
Colorbar(fig[2, 4], hm_T; label = "ᵒC")
hm_S = heatmap!(ax_S, Sₙ; colormap = :haline, colorrange = Slims)
Colorbar(fig[3, 2], hm_S; label = "g / kg")
hm_νₑ = heatmap!(ax_νₑ, νₑₙ; colormap = :thermal, colorrange = νₑlims)
Colorbar(fig[3, 4], hm_νₑ; label = "m s⁻²")
fig[1, 1:4] = Label(fig, title, fontsize=24, tellwidth=false)
fig
And now record a movie.
frames = intro:length(times)
@info "Making a motion picture of ocean wind mixing and convection..."
record(fig, filename * ".mp4", frames, framerate=8) do i
n[] = i
end[ Info: Making a motion picture of ocean wind mixing and convection...
This page was generated using Literate.jl.