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: 63.4 KiB

We're ready:

run!(simulation)

[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|w|) = 1.4e-05 ms⁻¹, wall time: 0 seconds
[ Info:     ... simulation initialization complete (20.692 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (8.425 seconds).
Iteration: 0040, time: 7.000 minutes, Δt: 8.309 seconds, max(|w|) = 2.9e-05 ms⁻¹, wall time: 30.128 seconds
Iteration: 0080, time: 11.183 minutes, Δt: 4.891 seconds, max(|w|) = 8.9e-03 ms⁻¹, wall time: 30.987 seconds
Iteration: 0120, time: 14 minutes, Δt: 4.085 seconds, max(|w|) = 2.2e-02 ms⁻¹, wall time: 31.471 seconds
Iteration: 0160, time: 16.691 minutes, Δt: 4.071 seconds, max(|w|) = 2.4e-02 ms⁻¹, wall time: 31.955 seconds
Iteration: 0200, time: 19.346 minutes, Δt: 4.138 seconds, max(|w|) = 2.8e-02 ms⁻¹, wall time: 32.477 seconds
Iteration: 0240, time: 21.889 minutes, Δt: 3.997 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 33.098 seconds
Iteration: 0280, time: 24.306 minutes, Δt: 3.813 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 33.639 seconds
Iteration: 0320, time: 26.674 minutes, Δt: 3.418 seconds, max(|w|) = 3.6e-02 ms⁻¹, wall time: 34.144 seconds
Iteration: 0360, time: 28.968 minutes, Δt: 3.740 seconds, max(|w|) = 2.8e-02 ms⁻¹, wall time: 34.633 seconds
Iteration: 0400, time: 31.234 minutes, Δt: 3.499 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 35.227 seconds
Iteration: 0440, time: 33.432 minutes, Δt: 3.189 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 35.853 seconds
Iteration: 0480, time: 35.540 minutes, Δt: 3.227 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 36.430 seconds
Iteration: 0520, time: 37.557 minutes, Δt: 3.021 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 37.025 seconds
Iteration: 0560, time: 39.492 minutes, Δt: 2.949 seconds, max(|w|) = 3.2e-02 ms⁻¹, wall time: 37.648 seconds
Iteration: 0600, time: 41.382 minutes, Δt: 2.884 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 38.248 seconds
Iteration: 0640, time: 43.242 minutes, Δt: 2.827 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 38.900 seconds
Iteration: 0680, time: 45.099 minutes, Δt: 3.022 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 39.467 seconds
Iteration: 0720, time: 47.049 minutes, Δt: 3.010 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 39.973 seconds
Iteration: 0760, time: 48.946 minutes, Δt: 2.813 seconds, max(|w|) = 4.0e-02 ms⁻¹, wall time: 40.491 seconds
Iteration: 0800, time: 50.762 minutes, Δt: 2.884 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 41.117 seconds
Iteration: 0840, time: 52.662 minutes, Δt: 2.809 seconds, max(|w|) = 4.0e-02 ms⁻¹, wall time: 41.601 seconds
Iteration: 0880, time: 54.472 minutes, Δt: 2.877 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 42.125 seconds
Iteration: 0920, time: 56.225 minutes, Δt: 2.722 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 42.625 seconds
Iteration: 0960, time: 58.045 minutes, Δt: 2.732 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 43.126 seconds
Iteration: 1000, time: 59.861 minutes, Δt: 2.700 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 43.679 seconds
Iteration: 1040, time: 1.026 hours, Δt: 2.657 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 44.168 seconds
Iteration: 1080, time: 1.055 hours, Δt: 2.650 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 44.657 seconds
Iteration: 1120, time: 1.084 hours, Δt: 2.730 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 45.148 seconds
Iteration: 1160, time: 1.113 hours, Δt: 2.619 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 45.683 seconds
Iteration: 1200, time: 1.141 hours, Δt: 2.604 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 46.303 seconds
Iteration: 1240, time: 1.169 hours, Δt: 2.423 seconds, max(|w|) = 5.9e-02 ms⁻¹, wall time: 46.787 seconds
Iteration: 1280, time: 1.197 hours, Δt: 2.620 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 47.363 seconds
Iteration: 1320, time: 1.225 hours, Δt: 2.635 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 47.865 seconds
Iteration: 1360, time: 1.254 hours, Δt: 2.570 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 48.356 seconds
Iteration: 1400, time: 1.282 hours, Δt: 2.514 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 48.945 seconds
Iteration: 1440, time: 1.309 hours, Δt: 2.612 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 49.493 seconds
Iteration: 1480, time: 1.337 hours, Δt: 2.480 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 50.053 seconds
Iteration: 1520, time: 1.365 hours, Δt: 2.536 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 50.533 seconds
Iteration: 1560, time: 1.393 hours, Δt: 2.502 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 51.113 seconds
Iteration: 1600, time: 1.421 hours, Δt: 2.637 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 51.624 seconds
Iteration: 1640, time: 1.450 hours, Δt: 2.573 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 52.210 seconds
Iteration: 1680, time: 1.477 hours, Δt: 2.614 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 52.785 seconds
Iteration: 1720, time: 1.505 hours, Δt: 2.643 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 53.365 seconds
Iteration: 1760, time: 1.533 hours, Δt: 2.538 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 54.035 seconds
Iteration: 1800, time: 1.562 hours, Δt: 2.547 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 54.582 seconds
Iteration: 1840, time: 1.589 hours, Δt: 2.454 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 55.163 seconds
Iteration: 1880, time: 1.615 hours, Δt: 2.521 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 55.730 seconds
Iteration: 1920, time: 1.642 hours, Δt: 2.436 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 56.313 seconds
Iteration: 1960, time: 1.669 hours, Δt: 2.382 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 57.000 seconds
Iteration: 2000, time: 1.694 hours, Δt: 2.384 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 57.567 seconds
Iteration: 2040, time: 1.721 hours, Δt: 2.505 seconds, max(|w|) = 5.9e-02 ms⁻¹, wall time: 58.138 seconds
Iteration: 2080, time: 1.749 hours, Δt: 2.442 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 58.743 seconds
Iteration: 2120, time: 1.775 hours, Δt: 2.494 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 59.400 seconds
Iteration: 2160, time: 1.802 hours, Δt: 2.422 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 59.986 seconds
Iteration: 2200, time: 1.828 hours, Δt: 2.448 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 1.009 minutes
Iteration: 2240, time: 1.855 hours, Δt: 2.441 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 1.019 minutes
Iteration: 2280, time: 1.881 hours, Δt: 2.406 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 1.028 minutes
Iteration: 2320, time: 1.908 hours, Δt: 2.479 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 1.039 minutes
Iteration: 2360, time: 1.933 hours, Δt: 2.435 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 1.048 minutes
Iteration: 2400, time: 1.960 hours, Δt: 2.436 seconds, max(|w|) = 6.8e-02 ms⁻¹, wall time: 1.058 minutes
Iteration: 2440, time: 1.987 hours, Δt: 2.401 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 1.067 minutes
[ Info: Simulation is stopping after running for 1.072 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.0470492, min=-0.0543595, mean=-4.86806e-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.5883, 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.0127, 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.0178075, min=0.0, mean=0.000423202)

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.