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.2e-05 ms⁻¹, wall time: 0 seconds
[ Info:     ... simulation initialization complete (35.397 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (12.377 seconds).
Iteration: 0040, time: 7 minutes, Δt: 8.309 seconds, max(|w|) = 2.8e-05 ms⁻¹, wall time: 48.782 seconds
Iteration: 0080, time: 11.183 minutes, Δt: 4.810 seconds, max(|w|) = 7.5e-03 ms⁻¹, wall time: 49.415 seconds
Iteration: 0120, time: 13.949 minutes, Δt: 4.083 seconds, max(|w|) = 2.1e-02 ms⁻¹, wall time: 49.967 seconds
Iteration: 0160, time: 16.598 minutes, Δt: 4.037 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 50.597 seconds
Iteration: 0200, time: 19.261 minutes, Δt: 4.087 seconds, max(|w|) = 2.8e-02 ms⁻¹, wall time: 51.212 seconds
Iteration: 0240, time: 21.959 minutes, Δt: 3.924 seconds, max(|w|) = 3.9e-02 ms⁻¹, wall time: 51.836 seconds
Iteration: 0280, time: 24.459 minutes, Δt: 3.850 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 52.519 seconds
Iteration: 0320, time: 26.823 minutes, Δt: 3.579 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 53.172 seconds
Iteration: 0360, time: 29.177 minutes, Δt: 3.564 seconds, max(|w|) = 2.8e-02 ms⁻¹, wall time: 53.865 seconds
Iteration: 0400, time: 31.464 minutes, Δt: 3.473 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 54.487 seconds
Iteration: 0440, time: 33.604 minutes, Δt: 3.227 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 55.101 seconds
Iteration: 0480, time: 35.673 minutes, Δt: 3.180 seconds, max(|w|) = 3.2e-02 ms⁻¹, wall time: 55.634 seconds
Iteration: 0520, time: 37.754 minutes, Δt: 3.137 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 56.217 seconds
Iteration: 0560, time: 39.757 minutes, Δt: 3.167 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 56.730 seconds
Iteration: 0600, time: 41.721 minutes, Δt: 2.842 seconds, max(|w|) = 3.6e-02 ms⁻¹, wall time: 57.332 seconds
Iteration: 0640, time: 43.663 minutes, Δt: 2.978 seconds, max(|w|) = 3.2e-02 ms⁻¹, wall time: 57.988 seconds
Iteration: 0680, time: 45.662 minutes, Δt: 3.028 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 58.636 seconds
Iteration: 0720, time: 47.519 minutes, Δt: 2.879 seconds, max(|w|) = 3.6e-02 ms⁻¹, wall time: 59.276 seconds
Iteration: 0760, time: 49.394 minutes, Δt: 2.981 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 59.809 seconds
Iteration: 0800, time: 51.237 minutes, Δt: 2.858 seconds, max(|w|) = 4.0e-02 ms⁻¹, wall time: 1.007 minutes
Iteration: 0840, time: 53 minutes, Δt: 2.646 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 1.015 minutes
Iteration: 0880, time: 54.761 minutes, Δt: 2.691 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 1.024 minutes
Iteration: 0920, time: 56.555 minutes, Δt: 2.705 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 1.033 minutes
Iteration: 0960, time: 58.306 minutes, Δt: 2.648 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 1.044 minutes
Iteration: 1000, time: 1.001 hours, Δt: 2.737 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 1.054 minutes
Iteration: 1040, time: 1.031 hours, Δt: 2.708 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 1.065 minutes
Iteration: 1080, time: 1.060 hours, Δt: 2.661 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 1.076 minutes
Iteration: 1120, time: 1.088 hours, Δt: 2.676 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 1.087 minutes
Iteration: 1160, time: 1.117 hours, Δt: 2.482 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 1.098 minutes
Iteration: 1200, time: 1.145 hours, Δt: 2.609 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 1.109 minutes
Iteration: 1240, time: 1.173 hours, Δt: 2.438 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 1.119 minutes
Iteration: 1280, time: 1.200 hours, Δt: 2.505 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 1.128 minutes
Iteration: 1320, time: 1.227 hours, Δt: 2.483 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 1.140 minutes
Iteration: 1360, time: 1.254 hours, Δt: 2.461 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 1.151 minutes
Iteration: 1400, time: 1.281 hours, Δt: 2.527 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 1.161 minutes
Iteration: 1440, time: 1.308 hours, Δt: 2.350 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 1.173 minutes
Iteration: 1480, time: 1.334 hours, Δt: 2.461 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 1.183 minutes
Iteration: 1520, time: 1.361 hours, Δt: 2.392 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 1.194 minutes
Iteration: 1560, time: 1.387 hours, Δt: 2.387 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 1.205 minutes
Iteration: 1600, time: 1.412 hours, Δt: 2.388 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 1.218 minutes
Iteration: 1640, time: 1.438 hours, Δt: 2.368 seconds, max(|w|) = 5.9e-02 ms⁻¹, wall time: 1.228 minutes
Iteration: 1680, time: 1.465 hours, Δt: 2.556 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 1.241 minutes
Iteration: 1720, time: 1.493 hours, Δt: 2.466 seconds, max(|w|) = 6.6e-02 ms⁻¹, wall time: 1.251 minutes
Iteration: 1760, time: 1.519 hours, Δt: 2.533 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 1.263 minutes
Iteration: 1800, time: 1.547 hours, Δt: 2.544 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 1.275 minutes
Iteration: 1840, time: 1.574 hours, Δt: 2.576 seconds, max(|w|) = 7.8e-02 ms⁻¹, wall time: 1.285 minutes
Iteration: 1880, time: 1.602 hours, Δt: 2.505 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 1.297 minutes
Iteration: 1920, time: 1.629 hours, Δt: 2.412 seconds, max(|w|) = 5.9e-02 ms⁻¹, wall time: 1.307 minutes
Iteration: 1960, time: 1.656 hours, Δt: 2.467 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 1.319 minutes
Iteration: 2000, time: 1.683 hours, Δt: 2.580 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 1.328 minutes
Iteration: 2040, time: 1.710 hours, Δt: 2.441 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 1.340 minutes
Iteration: 2080, time: 1.736 hours, Δt: 2.405 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 1.350 minutes
Iteration: 2120, time: 1.762 hours, Δt: 2.233 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 1.361 minutes
Iteration: 2160, time: 1.786 hours, Δt: 2.302 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 1.373 minutes
Iteration: 2200, time: 1.812 hours, Δt: 2.414 seconds, max(|w|) = 7.0e-02 ms⁻¹, wall time: 1.384 minutes
Iteration: 2240, time: 1.838 hours, Δt: 2.418 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 1.395 minutes
Iteration: 2280, time: 1.865 hours, Δt: 2.450 seconds, max(|w|) = 7.5e-02 ms⁻¹, wall time: 1.406 minutes
Iteration: 2320, time: 1.892 hours, Δt: 2.488 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 1.418 minutes
Iteration: 2360, time: 1.918 hours, Δt: 2.484 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 1.429 minutes
Iteration: 2400, time: 1.945 hours, Δt: 2.376 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 1.442 minutes
Iteration: 2440, time: 1.970 hours, Δt: 2.264 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 1.453 minutes
Iteration: 2480, time: 1.995 hours, Δt: 2.429 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 1.465 minutes
[ Info: Simulation is stopping after running for 1.469 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.0432293, min=-0.0553175, mean=5.80752e-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.0137, 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.0212619, min=0.0, mean=0.000413142)

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.