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 SeawaterBuoyancy model for buoyancy with TEOS10EquationOfState.
• 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, CUDA"

We start by importing all of the packages and functions that we'll need for this example.

using Oceananigans
using Oceananigans.Units

using CUDA
using Random
using Printf
using CairoMakie
using SeawaterPolynomials.TEOS10: TEOS10EquationOfState

The 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.96618

We 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_ flux

4.884283985946938e-5

Finally, 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.00023391812865497074

The 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-7

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

The 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: 64.1 KiB

We'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 (18.349 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (7.340 seconds).
Iteration: 0040, time: 7 minutes, Δt: 8.309 seconds, max(|w|) = 2.7e-05 ms⁻¹, wall time: 26.575 seconds
Iteration: 0080, time: 11.183 minutes, Δt: 4.897 seconds, max(|w|) = 7.6e-03 ms⁻¹, wall time: 27.099 seconds
Iteration: 0120, time: 13.945 minutes, Δt: 3.730 seconds, max(|w|) = 2.8e-02 ms⁻¹, wall time: 27.579 seconds
Iteration: 0160, time: 16.534 minutes, Δt: 4.260 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 28.173 seconds
Iteration: 0200, time: 19.203 minutes, Δt: 3.996 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 28.675 seconds
Iteration: 0240, time: 21.771 minutes, Δt: 4.063 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 29.147 seconds
Iteration: 0280, time: 24.326 minutes, Δt: 3.826 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 29.644 seconds
Iteration: 0320, time: 26.803 minutes, Δt: 3.636 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 30.138 seconds
Iteration: 0360, time: 29.118 minutes, Δt: 3.550 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 30.777 seconds
Iteration: 0400, time: 31.333 minutes, Δt: 3.356 seconds, max(|w|) = 3.4e-02 ms⁻¹, wall time: 31.281 seconds
Iteration: 0440, time: 33.565 minutes, Δt: 3.322 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 31.782 seconds
Iteration: 0480, time: 35.672 minutes, Δt: 3.288 seconds, max(|w|) = 3.7e-02 ms⁻¹, wall time: 32.250 seconds
Iteration: 0520, time: 37.785 minutes, Δt: 3.184 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 32.735 seconds
Iteration: 0560, time: 39.829 minutes, Δt: 3.055 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 33.318 seconds
Iteration: 0600, time: 41.808 minutes, Δt: 2.999 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 33.786 seconds
Iteration: 0640, time: 43.806 minutes, Δt: 2.926 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 34.271 seconds
Iteration: 0680, time: 45.751 minutes, Δt: 3.103 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 34.755 seconds
Iteration: 0720, time: 47.696 minutes, Δt: 2.936 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 35.302 seconds
Iteration: 0760, time: 49.596 minutes, Δt: 2.897 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 35.890 seconds
Iteration: 0800, time: 51.408 minutes, Δt: 2.902 seconds, max(|w|) = 3.7e-02 ms⁻¹, wall time: 36.345 seconds
Iteration: 0840, time: 53.270 minutes, Δt: 2.795 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 36.804 seconds
Iteration: 0880, time: 55.092 minutes, Δt: 2.786 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 37.292 seconds
Iteration: 0920, time: 56.916 minutes, Δt: 2.704 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 37.745 seconds
Iteration: 0960, time: 58.721 minutes, Δt: 2.757 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 38.287 seconds
Iteration: 1000, time: 1.008 hours, Δt: 2.635 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 38.736 seconds
Iteration: 1040, time: 1.037 hours, Δt: 2.799 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 39.259 seconds
Iteration: 1080, time: 1.067 hours, Δt: 2.741 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 39.760 seconds
Iteration: 1120, time: 1.097 hours, Δt: 2.679 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 40.217 seconds
Iteration: 1160, time: 1.125 hours, Δt: 2.576 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 40.783 seconds
Iteration: 1200, time: 1.154 hours, Δt: 2.742 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 41.242 seconds
Iteration: 1240, time: 1.183 hours, Δt: 2.683 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 41.697 seconds
Iteration: 1280, time: 1.211 hours, Δt: 2.622 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 42.170 seconds
Iteration: 1320, time: 1.239 hours, Δt: 2.589 seconds, max(|w|) = 4.2e-02 ms⁻¹, wall time: 42.652 seconds
Iteration: 1360, time: 1.267 hours, Δt: 2.520 seconds, max(|w|) = 4.0e-02 ms⁻¹, wall time: 43.112 seconds
Iteration: 1400, time: 1.293 hours, Δt: 2.429 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 43.658 seconds
Iteration: 1440, time: 1.320 hours, Δt: 2.591 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 44.202 seconds
Iteration: 1480, time: 1.349 hours, Δt: 2.569 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 44.696 seconds
Iteration: 1520, time: 1.377 hours, Δt: 2.576 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 45.205 seconds
Iteration: 1560, time: 1.404 hours, Δt: 2.574 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 45.710 seconds
Iteration: 1600, time: 1.432 hours, Δt: 2.610 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 46.348 seconds
Iteration: 1640, time: 1.458 hours, Δt: 2.500 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 46.874 seconds
Iteration: 1680, time: 1.485 hours, Δt: 2.483 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 47.394 seconds
Iteration: 1720, time: 1.511 hours, Δt: 2.500 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 47.891 seconds
Iteration: 1760, time: 1.539 hours, Δt: 2.477 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 48.375 seconds
Iteration: 1800, time: 1.565 hours, Δt: 2.285 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 48.909 seconds
Iteration: 1840, time: 1.588 hours, Δt: 2.333 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 49.377 seconds
Iteration: 1880, time: 1.615 hours, Δt: 2.518 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 49.850 seconds
Iteration: 1920, time: 1.642 hours, Δt: 2.455 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 50.336 seconds
Iteration: 1960, time: 1.668 hours, Δt: 2.486 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 50.815 seconds
Iteration: 2000, time: 1.695 hours, Δt: 2.349 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 51.356 seconds
Iteration: 2040, time: 1.720 hours, Δt: 2.440 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 51.825 seconds
Iteration: 2080, time: 1.746 hours, Δt: 2.294 seconds, max(|w|) = 7.8e-02 ms⁻¹, wall time: 52.295 seconds
Iteration: 2120, time: 1.771 hours, Δt: 2.332 seconds, max(|w|) = 7.7e-02 ms⁻¹, wall time: 52.771 seconds
Iteration: 2160, time: 1.797 hours, Δt: 2.342 seconds, max(|w|) = 7.7e-02 ms⁻¹, wall time: 53.232 seconds
Iteration: 2200, time: 1.822 hours, Δt: 2.374 seconds, max(|w|) = 6.6e-02 ms⁻¹, wall time: 53.784 seconds
Iteration: 2240, time: 1.848 hours, Δt: 2.378 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 54.232 seconds
Iteration: 2280, time: 1.873 hours, Δt: 2.362 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 54.705 seconds
Iteration: 2320, time: 1.899 hours, Δt: 2.417 seconds, max(|w|) = 7.1e-02 ms⁻¹, wall time: 55.170 seconds
Iteration: 2360, time: 1.925 hours, Δt: 2.370 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 55.648 seconds
Iteration: 2400, time: 1.951 hours, Δt: 2.345 seconds, max(|w|) = 7.1e-02 ms⁻¹, wall time: 56.210 seconds
Iteration: 2440, time: 1.977 hours, Δt: 2.410 seconds, max(|w|) = 7.1e-02 ms⁻¹, wall time: 56.666 seconds
[ Info: Simulation is stopping after running for 57.100 seconds.
[ 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.0368549, min=-0.0599793, mean=-4.87989e-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.0131, 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.0204773, min=0.0, mean=0.00043399)

We start the animation at $t = 10$ minutes since things are pretty boring till then:

times = time_series.w.times
intro = searchsortedfirst(times, 10minutes)

11

We 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..."

CairoMakie.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.