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.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.7 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 (13.617 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (5.016 seconds).
Iteration: 0040, time: 7.344 minutes, Δt: 7.912 seconds, max(|w|) = 3.9e-05 ms⁻¹, wall time: 19.527 seconds
Iteration: 0080, time: 11.441 minutes, Δt: 4.706 seconds, max(|w|) = 8.2e-03 ms⁻¹, wall time: 19.958 seconds
Iteration: 0120, time: 14.205 minutes, Δt: 3.909 seconds, max(|w|) = 2.2e-02 ms⁻¹, wall time: 20.414 seconds
Iteration: 0160, time: 16.821 minutes, Δt: 3.841 seconds, max(|w|) = 2.4e-02 ms⁻¹, wall time: 20.834 seconds
Iteration: 0200, time: 19.346 minutes, Δt: 3.998 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 21.329 seconds
Iteration: 0240, time: 21.886 minutes, Δt: 3.793 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 21.851 seconds
Iteration: 0280, time: 24.239 minutes, Δt: 3.763 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 22.355 seconds
Iteration: 0320, time: 26.624 minutes, Δt: 3.561 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 22.892 seconds
Iteration: 0360, time: 28.924 minutes, Δt: 3.473 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 23.304 seconds
Iteration: 0400, time: 31.054 minutes, Δt: 3.327 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 23.764 seconds
Iteration: 0440, time: 33.221 minutes, Δt: 3.387 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 24.179 seconds
Iteration: 0480, time: 35.368 minutes, Δt: 3.244 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 24.631 seconds
Iteration: 0520, time: 37.481 minutes, Δt: 3.232 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 25.055 seconds
Iteration: 0560, time: 39.506 minutes, Δt: 3.019 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 25.554 seconds
Iteration: 0600, time: 41.521 minutes, Δt: 3.131 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 26.040 seconds
Iteration: 0640, time: 43.520 minutes, Δt: 3.109 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 26.571 seconds
Iteration: 0680, time: 45.479 minutes, Δt: 2.987 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 27.095 seconds
Iteration: 0720, time: 47.451 minutes, Δt: 2.892 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 27.503 seconds
Iteration: 0760, time: 49.334 minutes, Δt: 2.861 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 27.979 seconds
Iteration: 0800, time: 51.141 minutes, Δt: 2.800 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 28.393 seconds
Iteration: 0840, time: 52.956 minutes, Δt: 2.809 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 28.843 seconds
Iteration: 0880, time: 54.778 minutes, Δt: 2.755 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 29.265 seconds
Iteration: 0920, time: 56.613 minutes, Δt: 2.802 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 29.718 seconds
Iteration: 0960, time: 58.414 minutes, Δt: 2.732 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 30.127 seconds
Iteration: 1000, time: 1.004 hours, Δt: 2.785 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 30.569 seconds
Iteration: 1040, time: 1.034 hours, Δt: 2.758 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 30.996 seconds
Iteration: 1080, time: 1.064 hours, Δt: 2.723 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 31.453 seconds
Iteration: 1120, time: 1.094 hours, Δt: 2.613 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 31.878 seconds
Iteration: 1160, time: 1.123 hours, Δt: 2.718 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 32.333 seconds
Iteration: 1200, time: 1.151 hours, Δt: 2.582 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 32.784 seconds
Iteration: 1240, time: 1.180 hours, Δt: 2.647 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 33.196 seconds
Iteration: 1280, time: 1.208 hours, Δt: 2.544 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 33.643 seconds
Iteration: 1320, time: 1.236 hours, Δt: 2.544 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 34.056 seconds
Iteration: 1360, time: 1.264 hours, Δt: 2.555 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 34.495 seconds
Iteration: 1400, time: 1.290 hours, Δt: 2.404 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 34.906 seconds
Iteration: 1440, time: 1.316 hours, Δt: 2.459 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 35.340 seconds
Iteration: 1480, time: 1.342 hours, Δt: 2.327 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 35.755 seconds
Iteration: 1520, time: 1.367 hours, Δt: 2.218 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 36.208 seconds
Iteration: 1560, time: 1.392 hours, Δt: 2.349 seconds, max(|w|) = 6.8e-02 ms⁻¹, wall time: 36.623 seconds
Iteration: 1600, time: 1.418 hours, Δt: 2.372 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 37.065 seconds
Iteration: 1640, time: 1.444 hours, Δt: 2.358 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 37.476 seconds
Iteration: 1680, time: 1.471 hours, Δt: 2.366 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 37.921 seconds
Iteration: 1720, time: 1.496 hours, Δt: 2.362 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 38.386 seconds
Iteration: 1760, time: 1.523 hours, Δt: 2.397 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 38.795 seconds
Iteration: 1800, time: 1.550 hours, Δt: 2.414 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 39.235 seconds
Iteration: 1840, time: 1.576 hours, Δt: 2.423 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 39.658 seconds
Iteration: 1880, time: 1.603 hours, Δt: 2.428 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 40.109 seconds
Iteration: 1920, time: 1.630 hours, Δt: 2.565 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 40.518 seconds
Iteration: 1960, time: 1.658 hours, Δt: 2.583 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 40.963 seconds
Iteration: 2000, time: 1.685 hours, Δt: 2.484 seconds, max(|w|) = 7.0e-02 ms⁻¹, wall time: 41.384 seconds
Iteration: 2040, time: 1.712 hours, Δt: 2.407 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 41.826 seconds
Iteration: 2080, time: 1.738 hours, Δt: 2.391 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 42.251 seconds
Iteration: 2120, time: 1.765 hours, Δt: 2.478 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 42.677 seconds
Iteration: 2160, time: 1.792 hours, Δt: 2.500 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 43.101 seconds
Iteration: 2200, time: 1.819 hours, Δt: 2.373 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 43.561 seconds
Iteration: 2240, time: 1.846 hours, Δt: 2.521 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 44.008 seconds
Iteration: 2280, time: 1.873 hours, Δt: 2.577 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 44.432 seconds
Iteration: 2320, time: 1.901 hours, Δt: 2.564 seconds, max(|w|) = 7.4e-02 ms⁻¹, wall time: 44.892 seconds
Iteration: 2360, time: 1.929 hours, Δt: 2.530 seconds, max(|w|) = 7.2e-02 ms⁻¹, wall time: 45.299 seconds
Iteration: 2400, time: 1.955 hours, Δt: 2.330 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 45.745 seconds
Iteration: 2440, time: 1.981 hours, Δt: 2.444 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 46.156 seconds
[ Info: Simulation is stopping after running for 46.499 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.0450738, min=-0.0561239, mean=1.4673e-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.5884, 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.0195889, min=0.0, mean=0.000422781)

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.