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.3e-05 ms⁻¹, wall time: 0 seconds
[ Info:     ... simulation initialization complete (29.536 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (6.110 seconds).
Iteration: 0040, time: 7.000 minutes, Δt: 8.310 seconds, max(|w|) = 3.2e-05 ms⁻¹, wall time: 36.601 seconds
Iteration: 0080, time: 11.093 minutes, Δt: 4.975 seconds, max(|w|) = 6.5e-03 ms⁻¹, wall time: 37.277 seconds
Iteration: 0120, time: 13.967 minutes, Δt: 4.091 seconds, max(|w|) = 2.2e-02 ms⁻¹, wall time: 37.812 seconds
Iteration: 0160, time: 16.623 minutes, Δt: 4.177 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 38.393 seconds
Iteration: 0200, time: 19.274 minutes, Δt: 4.107 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 38.944 seconds
Iteration: 0240, time: 21.916 minutes, Δt: 3.829 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 39.540 seconds
Iteration: 0280, time: 24.376 minutes, Δt: 3.846 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 40.090 seconds
Iteration: 0320, time: 26.883 minutes, Δt: 3.542 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 40.631 seconds
Iteration: 0360, time: 29.178 minutes, Δt: 3.581 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 41.240 seconds
Iteration: 0400, time: 31.464 minutes, Δt: 3.525 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 41.799 seconds
Iteration: 0440, time: 33.647 minutes, Δt: 3.160 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 42.407 seconds
Iteration: 0480, time: 35.687 minutes, Δt: 3.176 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 42.975 seconds
Iteration: 0520, time: 37.757 minutes, Δt: 3.216 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 43.554 seconds
Iteration: 0560, time: 39.840 minutes, Δt: 3.020 seconds, max(|w|) = 3.2e-02 ms⁻¹, wall time: 44.102 seconds
Iteration: 0600, time: 41.767 minutes, Δt: 3.009 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 44.689 seconds
Iteration: 0640, time: 43.747 minutes, Δt: 2.870 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 45.234 seconds
Iteration: 0680, time: 45.590 minutes, Δt: 2.901 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 45.755 seconds
Iteration: 0720, time: 47.483 minutes, Δt: 2.923 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 46.335 seconds
Iteration: 0760, time: 49.342 minutes, Δt: 2.958 seconds, max(|w|) = 3.7e-02 ms⁻¹, wall time: 46.857 seconds
Iteration: 0800, time: 51.245 minutes, Δt: 2.881 seconds, max(|w|) = 3.7e-02 ms⁻¹, wall time: 47.414 seconds
Iteration: 0840, time: 53.143 minutes, Δt: 2.721 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 47.964 seconds
Iteration: 0880, time: 54.898 minutes, Δt: 2.743 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 48.541 seconds
Iteration: 0920, time: 56.696 minutes, Δt: 2.803 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 49.067 seconds
Iteration: 0960, time: 58.499 minutes, Δt: 2.772 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 49.593 seconds
Iteration: 1000, time: 1.005 hours, Δt: 2.795 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 50.171 seconds
Iteration: 1040, time: 1.036 hours, Δt: 2.791 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 50.726 seconds
Iteration: 1080, time: 1.066 hours, Δt: 2.607 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 51.306 seconds
Iteration: 1120, time: 1.093 hours, Δt: 2.636 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 51.846 seconds
Iteration: 1160, time: 1.121 hours, Δt: 2.741 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 52.416 seconds
Iteration: 1200, time: 1.151 hours, Δt: 2.724 seconds, max(|w|) = 6.6e-02 ms⁻¹, wall time: 52.962 seconds
Iteration: 1240, time: 1.180 hours, Δt: 2.601 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 53.468 seconds
Iteration: 1280, time: 1.207 hours, Δt: 2.669 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 54.044 seconds
Iteration: 1320, time: 1.236 hours, Δt: 2.693 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 54.602 seconds
Iteration: 1360, time: 1.266 hours, Δt: 2.687 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 55.181 seconds
Iteration: 1400, time: 1.294 hours, Δt: 2.579 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 55.745 seconds
Iteration: 1440, time: 1.322 hours, Δt: 2.538 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 56.335 seconds
Iteration: 1480, time: 1.350 hours, Δt: 2.642 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 56.880 seconds
Iteration: 1520, time: 1.379 hours, Δt: 2.583 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 57.485 seconds
Iteration: 1560, time: 1.407 hours, Δt: 2.613 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 58.070 seconds
Iteration: 1600, time: 1.434 hours, Δt: 2.524 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 58.644 seconds
Iteration: 1640, time: 1.462 hours, Δt: 2.488 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 59.220 seconds
Iteration: 1680, time: 1.488 hours, Δt: 2.490 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 59.773 seconds
Iteration: 1720, time: 1.516 hours, Δt: 2.595 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 1.006 minutes
Iteration: 1760, time: 1.543 hours, Δt: 2.448 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 1.015 minutes
Iteration: 1800, time: 1.569 hours, Δt: 2.521 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 1.025 minutes
Iteration: 1840, time: 1.597 hours, Δt: 2.400 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 1.034 minutes
Iteration: 1880, time: 1.622 hours, Δt: 2.319 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 1.044 minutes
Iteration: 1920, time: 1.648 hours, Δt: 2.451 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 1.054 minutes
Iteration: 1960, time: 1.673 hours, Δt: 2.487 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 1.063 minutes
Iteration: 2000, time: 1.700 hours, Δt: 2.562 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 1.073 minutes
Iteration: 2040, time: 1.728 hours, Δt: 2.436 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 1.082 minutes
Iteration: 2080, time: 1.754 hours, Δt: 2.597 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 1.092 minutes
Iteration: 2120, time: 1.782 hours, Δt: 2.521 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 1.101 minutes
Iteration: 2160, time: 1.809 hours, Δt: 2.421 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 1.111 minutes
Iteration: 2200, time: 1.835 hours, Δt: 2.448 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 1.121 minutes
Iteration: 2240, time: 1.862 hours, Δt: 2.368 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 1.130 minutes
Iteration: 2280, time: 1.887 hours, Δt: 2.365 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 1.140 minutes
Iteration: 2320, time: 1.913 hours, Δt: 2.436 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 1.149 minutes
Iteration: 2360, time: 1.939 hours, Δt: 2.441 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 1.159 minutes
Iteration: 2400, time: 1.966 hours, Δt: 2.308 seconds, max(|w|) = 5.9e-02 ms⁻¹, wall time: 1.168 minutes
Iteration: 2440, time: 1.991 hours, Δt: 2.334 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 1.178 minutes
[ Info: Simulation is stopping after running for 1.181 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.0472023, min=-0.0431702, mean=-1.01747e-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.0133, 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.0179819, min=0.0, mean=0.000399671)

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.