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
├── run_wall_time: 0 seconds
├── run_wall_time / 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 => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: 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.closure_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.8 KiB

We're ready:

run!(simulation)

[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|w|) = 1.5e-05 ms⁻¹, wall time: 0 seconds
[ Info:     ... simulation initialization complete (13.020 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.375 seconds).
Iteration: 0040, time: 7.344 minutes, Δt: 7.912 seconds, max(|w|) = 4.2e-05 ms⁻¹, wall time: 17.446 seconds
Iteration: 0080, time: 11.439 minutes, Δt: 4.676 seconds, max(|w|) = 9.9e-03 ms⁻¹, wall time: 17.904 seconds
Iteration: 0120, time: 14.206 minutes, Δt: 4.134 seconds, max(|w|) = 2.4e-02 ms⁻¹, wall time: 18.404 seconds
Iteration: 0160, time: 16.922 minutes, Δt: 4.125 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 18.880 seconds
Iteration: 0200, time: 19.453 minutes, Δt: 3.929 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 19.371 seconds
Iteration: 0240, time: 21.995 minutes, Δt: 4.050 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 19.878 seconds
Iteration: 0280, time: 24.512 minutes, Δt: 3.682 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 20.278 seconds
Iteration: 0320, time: 26.917 minutes, Δt: 3.486 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 20.973 seconds
Iteration: 0360, time: 29.113 minutes, Δt: 3.340 seconds, max(|w|) = 3.4e-02 ms⁻¹, wall time: 21.750 seconds
Iteration: 0400, time: 31.329 minutes, Δt: 3.223 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 22.224 seconds
Iteration: 0440, time: 33.497 minutes, Δt: 3.100 seconds, max(|w|) = 3.2e-02 ms⁻¹, wall time: 22.806 seconds
Iteration: 0480, time: 35.610 minutes, Δt: 3.138 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 23.266 seconds
Iteration: 0520, time: 37.647 minutes, Δt: 3.289 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 23.696 seconds
Iteration: 0560, time: 39.741 minutes, Δt: 3.153 seconds, max(|w|) = 3.4e-02 ms⁻¹, wall time: 24.083 seconds
Iteration: 0600, time: 41.846 minutes, Δt: 3.144 seconds, max(|w|) = 3.2e-02 ms⁻¹, wall time: 24.481 seconds
Iteration: 0640, time: 43.868 minutes, Δt: 3.048 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 24.880 seconds
Iteration: 0680, time: 45.854 minutes, Δt: 3.046 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 25.383 seconds
Iteration: 0720, time: 47.824 minutes, Δt: 2.788 seconds, max(|w|) = 3.4e-02 ms⁻¹, wall time: 25.762 seconds
Iteration: 0760, time: 49.667 minutes, Δt: 2.853 seconds, max(|w|) = 3.6e-02 ms⁻¹, wall time: 26.317 seconds
Iteration: 0800, time: 51.464 minutes, Δt: 2.870 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 26.789 seconds
Iteration: 0840, time: 53.283 minutes, Δt: 2.772 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 27.210 seconds
Iteration: 0880, time: 55.043 minutes, Δt: 2.809 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 27.688 seconds
Iteration: 0920, time: 56.859 minutes, Δt: 2.636 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 28.154 seconds
Iteration: 0960, time: 58.604 minutes, Δt: 2.815 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 28.578 seconds
Iteration: 1000, time: 1.007 hours, Δt: 2.627 seconds, max(|w|) = 4.0e-02 ms⁻¹, wall time: 28.994 seconds
Iteration: 1040, time: 1.036 hours, Δt: 2.724 seconds, max(|w|) = 4.2e-02 ms⁻¹, wall time: 29.693 seconds
Iteration: 1080, time: 1.065 hours, Δt: 2.707 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 30.093 seconds
Iteration: 1120, time: 1.095 hours, Δt: 2.683 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 30.531 seconds
Iteration: 1160, time: 1.123 hours, Δt: 2.606 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 31.034 seconds
Iteration: 1200, time: 1.151 hours, Δt: 2.684 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 31.593 seconds
Iteration: 1240, time: 1.181 hours, Δt: 2.574 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 32.044 seconds
Iteration: 1280, time: 1.210 hours, Δt: 2.617 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 32.505 seconds
Iteration: 1320, time: 1.239 hours, Δt: 2.686 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 32.951 seconds
Iteration: 1360, time: 1.267 hours, Δt: 2.571 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 33.361 seconds
Iteration: 1400, time: 1.295 hours, Δt: 2.445 seconds, max(|w|) = 5.9e-02 ms⁻¹, wall time: 33.816 seconds
Iteration: 1440, time: 1.322 hours, Δt: 2.516 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 34.194 seconds
Iteration: 1480, time: 1.349 hours, Δt: 2.423 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 34.845 seconds
Iteration: 1520, time: 1.376 hours, Δt: 2.436 seconds, max(|w|) = 6.4e-02 ms⁻¹, wall time: 35.626 seconds
Iteration: 1560, time: 1.402 hours, Δt: 2.471 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 36.057 seconds
Iteration: 1600, time: 1.429 hours, Δt: 2.434 seconds, max(|w|) = 6.4e-02 ms⁻¹, wall time: 36.672 seconds
Iteration: 1640, time: 1.455 hours, Δt: 2.494 seconds, max(|w|) = 6.4e-02 ms⁻¹, wall time: 37.232 seconds
Iteration: 1680, time: 1.482 hours, Δt: 2.495 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 37.749 seconds
Iteration: 1720, time: 1.509 hours, Δt: 2.464 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 38.216 seconds
Iteration: 1760, time: 1.535 hours, Δt: 2.471 seconds, max(|w|) = 6.8e-02 ms⁻¹, wall time: 38.633 seconds
Iteration: 1800, time: 1.563 hours, Δt: 2.523 seconds, max(|w|) = 8.4e-02 ms⁻¹, wall time: 39.063 seconds
Iteration: 1840, time: 1.589 hours, Δt: 2.535 seconds, max(|w|) = 7.4e-02 ms⁻¹, wall time: 39.612 seconds
Iteration: 1880, time: 1.617 hours, Δt: 2.456 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 40.156 seconds
Iteration: 1920, time: 1.644 hours, Δt: 2.407 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 40.630 seconds
Iteration: 1960, time: 1.671 hours, Δt: 2.548 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 41.238 seconds
Iteration: 2000, time: 1.699 hours, Δt: 2.568 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 41.739 seconds
Iteration: 2040, time: 1.727 hours, Δt: 2.589 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 42.195 seconds
Iteration: 2080, time: 1.755 hours, Δt: 2.554 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 42.640 seconds
Iteration: 2120, time: 1.782 hours, Δt: 2.400 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 43.054 seconds
Iteration: 2160, time: 1.809 hours, Δt: 2.552 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 43.522 seconds
Iteration: 2200, time: 1.837 hours, Δt: 2.465 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 44.117 seconds
Iteration: 2240, time: 1.864 hours, Δt: 2.436 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 44.602 seconds
Iteration: 2280, time: 1.892 hours, Δt: 2.513 seconds, max(|w|) = 6.4e-02 ms⁻¹, wall time: 45.530 seconds
Iteration: 2320, time: 1.919 hours, Δt: 2.461 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 46.092 seconds
Iteration: 2360, time: 1.945 hours, Δt: 2.423 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 46.616 seconds
Iteration: 2400, time: 1.972 hours, Δt: 2.494 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 47.131 seconds
Iteration: 2440, time: 1.998 hours, Δt: 2.489 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 47.625 seconds
[ Info: Simulation is stopping after running for 47.653 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.0523118, min=-0.0490328, mean=-5.63376e-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.5887, 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.0124, 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.0166221, min=0.0, mean=0.000424559)

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.