Wind- and convection-driven mixing in an ocean surface boundary layer

This example simulates mixing by three-dimensional turbulence in an ocean surface boundary layer driven by atmospheric winds and convection. It demonstrates:

• How to set-up a grid with varying spacing in the vertical direction
• How to use the SeawaterBuoyancy model for buoyancy with TEOS10EquationOfState.
• How to use a turbulence closure for large eddy simulation.
• How to use a function to impose a boundary condition.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, CairoMakie, SeawaterPolynomials, CUDA"

We start by importing all of the packages and functions that we'll need for this example.

using Oceananigans
using Oceananigans.Units

using CUDA
using Random
using Printf
using CairoMakie
using SeawaterPolynomials.TEOS10: TEOS10EquationOfState

The grid

We use 128²×64 grid points with 1 m grid spacing in the horizontal and varying spacing in the vertical, with higher resolution closer to the surface. Here we use a stretching function for the vertical nodes that maintains relatively constant vertical spacing in the mixed layer, which is desirable from a numerical standpoint:

Nx = Ny = 128    # number of points in each of horizontal directions
Nz = 64          # number of points in the vertical direction

Lx = Ly = 128    # (m) domain horizontal extents
Lz = 64          # (m) domain depth

refinement = 1.2 # controls spacing near surface (higher means finer spaced)
stretching = 12  # controls rate of stretching at bottom

# Normalized height ranging from 0 to 1
h(k) = (k - 1) / Nz

# Linear near-surface generator
ζ₀(k) = 1 + (h(k) - 1) / refinement

# Bottom-intensified stretching function
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))

# Generating function
z_interfaces(k) = Lz * (ζ₀(k) * Σ(k) - 1)

grid = RectilinearGrid(GPU(),
                       size = (Nx, Nx, Nz),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = z_interfaces)

128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CUDAGPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 128.0) regularly spaced with Δx=1.0
├── Periodic y ∈ [0.0, 128.0) regularly spaced with Δy=1.0
└── Bounded  z ∈ [-64.0, 0.0] variably spaced with min(Δz)=0.833413, max(Δz)=1.96618

We plot vertical spacing versus depth to inspect the prescribed grid stretching:

fig = Figure(size=(1200, 800))
ax = Axis(fig[1, 1], ylabel = "z (m)", xlabel = "Vertical spacing (m)")

lines!(ax, zspacings(grid, Center()))
scatter!(ax, zspacings(grid, Center()))

fig

Buoyancy that depends on temperature and salinity

We use the SeawaterBuoyancy model with the TEOS10 equation of state,

ρₒ = 1026 # kg m⁻³, average density at the surface of the world ocean
equation_of_state = TEOS10EquationOfState(reference_density=ρₒ)
buoyancy = SeawaterBuoyancy(; equation_of_state)

SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: BoussinesqEquationOfState{Float64}

Boundary conditions

We calculate the surface temperature flux associated with surface cooling of 200 W m⁻², reference density ρₒ, and heat capacity cᴾ,

Q = 200   # W m⁻², surface _heat_ flux
cᴾ = 3991 # J K⁻¹ kg⁻¹, typical heat capacity for seawater

Jᵀ = Q / (ρₒ * cᴾ) # K m s⁻¹, surface _temperature_ flux

4.884283985946938e-5

Finally, we impose a temperature gradient dTdz both initially (see "Initial conditions" section below) and at the bottom of the domain, culminating in the boundary conditions on temperature,

dTdz = 0.01 # K m⁻¹

T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᵀ),
                                bottom = GradientBoundaryCondition(dTdz))

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: GradientBoundaryCondition: 0.01
├── top: FluxBoundaryCondition: 4.88428e-5
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

Note that a positive temperature flux at the surface of the ocean implies cooling. This is because a positive temperature flux implies that temperature is fluxed upwards, out of the ocean.

For the velocity field, we imagine a wind blowing over the ocean surface with an average velocity at 10 meters u₁₀, and use a drag coefficient cᴰ to estimate the kinematic stress (that is, stress divided by density) exerted by the wind on the ocean:

u₁₀ = 10  # m s⁻¹, average wind velocity 10 meters above the ocean
cᴰ = 2e-3 # dimensionless drag coefficient
ρₐ = 1.2  # kg m⁻³, approximate average density of air at sea-level
τx = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²

-0.00023391812865497074

The boundary conditions on u are thus

u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τx))

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: -0.000233918
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

For salinity, S, we impose an evaporative flux of the form

@inline Jˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S # [salinity unit] m s⁻¹

where S is salinity. We use an evaporation rate of 1 millimeter per hour,

evaporation_rate = 1e-3 / hour # m s⁻¹

2.7777777777777776e-7

We build the Flux evaporation BoundaryCondition with the function Jˢ, indicating that Jˢ depends on salinity S and passing the parameter evaporation_rate,

evaporation_bc = FluxBoundaryCondition(Jˢ, field_dependencies=:S, parameters=evaporation_rate)

FluxBoundaryCondition: ContinuousBoundaryFunction Jˢ at (Nothing, Nothing, Nothing)

The full salinity boundary conditions are

S_bcs = FieldBoundaryConditions(top=evaporation_bc)

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: ContinuousBoundaryFunction Jˢ at (Nothing, Nothing, Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

Model instantiation

We fill in the final details of the model here, i.e., Coriolis forces, and the AnisotropicMinimumDissipation closure for large eddy simulation to model the effect of turbulent motions at scales smaller than the grid scale that are not explicitly resolved.

model = NonhydrostaticModel(; grid, buoyancy,
                            tracers = (:T, :S),
                            coriolis = FPlane(f=1e-4),
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs))

NonhydrostaticModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CUDAGPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: (T, S)
├── closure: AnisotropicMinimumDissipation{ExplicitTimeDiscretization, @NamedTuple{T::Float64, S::Float64}, Float64, Nothing}
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
└── coriolis: FPlane{Float64}(f=0.0001)

Note: To use the Smagorinsky-Lilly turbulence closure (with a constant model coefficient) rather than AnisotropicMinimumDissipation, use closure = SmagorinskyLilly() in the model constructor.

Initial conditions

Our initial condition for temperature consists of a linear stratification superposed with random noise damped at the walls, while our initial condition for velocity consists only of random noise.

# Random noise damped at top and bottom
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

# Temperature initial condition: a stable density gradient with random noise superposed.
Tᵢ(x, y, z) = 20 + dTdz * z + dTdz * model.grid.Lz * 1e-6 * Ξ(z)

# Velocity initial condition: random noise scaled by the friction velocity.
uᵢ(x, y, z) = sqrt(abs(τx)) * 1e-3 * Ξ(z)

# `set!` the `model` fields using functions or constants:
set!(model, u=uᵢ, w=uᵢ, T=Tᵢ, S=35)

Setting up a simulation

We set-up a simulation with an initial time-step of 10 seconds that stops at 2 hours, with adaptive time-stepping and progress printing.

simulation = Simulation(model, Δt=10, stop_time=2hours)

Simulation of NonhydrostaticModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 10 seconds
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 2 hours
├── Stop iteration: Inf
├── Wall time limit: Inf
├── Minimum relative step: 0.0
├── Callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => 4
│   ├── stop_iteration_exceeded => -
│   ├── wall_time_limit_exceeded => e
│   └── nan_checker => }
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entries

The TimeStepWizard helps ensure stable time-stepping with a Courant-Freidrichs-Lewy (CFL) number of 1.0.

wizard = TimeStepWizard(cfl=1, max_change=1.1, max_Δt=1minute)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

Callback of TimeStepWizard(cfl=1.0, max_Δt=60.0, min_Δt=0.0) on IterationInterval(10)

Nice progress messaging is helpful:

# Print a progress message
progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, max(|w|) = %.1e ms⁻¹, wall time: %s\n",
                                iteration(sim), prettytime(sim), prettytime(sim.Δt),
                                maximum(abs, sim.model.velocities.w), prettytime(sim.run_wall_time))

add_callback!(simulation, progress_message, IterationInterval(40))

We then set up the simulation:

Output

We use the JLD2Writer to save $x, z$ slices of the velocity fields, tracer fields, and eddy diffusivities. The prefix keyword argument to JLD2Writer indicates that output will be saved in ocean_wind_mixing_and_convection.jld2.

# Create a NamedTuple with eddy viscosity
eddy_viscosity = (; νₑ = model.diffusivity_fields.νₑ)

filename = "ocean_wind_mixing_and_convection"

simulation.output_writers[:slices] =
    JLD2Writer(model, merge(model.velocities, model.tracers, eddy_viscosity),
               filename = filename * ".jld2",
               indices = (:, grid.Ny/2, :),
               schedule = TimeInterval(1minute),
               overwrite_existing = true)

JLD2Writer scheduled on TimeInterval(1 minute):
├── filepath: ocean_wind_mixing_and_convection.jld2
├── 6 outputs: (u, v, w, T, S, νₑ)
├── array_type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 63.4 KiB

We're ready:

run!(simulation)

[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|w|) = 1.4e-05 ms⁻¹, wall time: 0 seconds
[ Info:     ... simulation initialization complete (23.380 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (5.509 seconds).
Iteration: 0040, time: 7 minutes, Δt: 8.309 seconds, max(|w|) = 3.1e-05 ms⁻¹, wall time: 29.791 seconds
Iteration: 0080, time: 11 minutes, Δt: 5.003 seconds, max(|w|) = 6.5e-03 ms⁻¹, wall time: 30.280 seconds
Iteration: 0120, time: 13.886 minutes, Δt: 4.018 seconds, max(|w|) = 2.2e-02 ms⁻¹, wall time: 30.788 seconds
Iteration: 0160, time: 16.518 minutes, Δt: 4.244 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 31.408 seconds
Iteration: 0200, time: 19.200 minutes, Δt: 4.051 seconds, max(|w|) = 2.5e-02 ms⁻¹, wall time: 31.889 seconds
Iteration: 0240, time: 21.872 minutes, Δt: 3.977 seconds, max(|w|) = 2.5e-02 ms⁻¹, wall time: 32.384 seconds
Iteration: 0280, time: 24.305 minutes, Δt: 3.684 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 32.867 seconds
Iteration: 0320, time: 26.589 minutes, Δt: 3.598 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 33.352 seconds
Iteration: 0360, time: 28.904 minutes, Δt: 3.413 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 33.946 seconds
Iteration: 0400, time: 31.117 minutes, Δt: 3.543 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 34.468 seconds
Iteration: 0440, time: 33.342 minutes, Δt: 3.387 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 34.952 seconds
Iteration: 0480, time: 35.481 minutes, Δt: 3.209 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 35.447 seconds
Iteration: 0520, time: 37.593 minutes, Δt: 3.215 seconds, max(|w|) = 3.6e-02 ms⁻¹, wall time: 35.925 seconds
Iteration: 0560, time: 39.638 minutes, Δt: 3.267 seconds, max(|w|) = 3.6e-02 ms⁻¹, wall time: 36.484 seconds
Iteration: 0600, time: 41.711 minutes, Δt: 3.061 seconds, max(|w|) = 3.4e-02 ms⁻¹, wall time: 36.939 seconds
Iteration: 0640, time: 43.714 minutes, Δt: 3.074 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 37.428 seconds
Iteration: 0680, time: 45.647 minutes, Δt: 3.014 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 37.913 seconds
Iteration: 0720, time: 47.523 minutes, Δt: 2.927 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 38.413 seconds
Iteration: 0760, time: 49.375 minutes, Δt: 2.794 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 38.955 seconds
Iteration: 0800, time: 51.187 minutes, Δt: 2.787 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 39.423 seconds
Iteration: 0840, time: 53.048 minutes, Δt: 2.887 seconds, max(|w|) = 4.2e-02 ms⁻¹, wall time: 39.944 seconds
Iteration: 0880, time: 54.949 minutes, Δt: 2.800 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 40.445 seconds
Iteration: 0920, time: 56.689 minutes, Δt: 2.767 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 41.037 seconds
Iteration: 0960, time: 58.461 minutes, Δt: 2.739 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 41.523 seconds
Iteration: 1000, time: 1.004 hours, Δt: 2.756 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 42.019 seconds
Iteration: 1040, time: 1.033 hours, Δt: 2.620 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 42.497 seconds
Iteration: 1080, time: 1.062 hours, Δt: 2.656 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 42.995 seconds
Iteration: 1120, time: 1.092 hours, Δt: 2.590 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 43.575 seconds
Iteration: 1160, time: 1.120 hours, Δt: 2.672 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 44.057 seconds
Iteration: 1200, time: 1.149 hours, Δt: 2.440 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 44.533 seconds
Iteration: 1240, time: 1.175 hours, Δt: 2.679 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 45.008 seconds
Iteration: 1280, time: 1.205 hours, Δt: 2.688 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 45.501 seconds
Iteration: 1320, time: 1.233 hours, Δt: 2.645 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 46.040 seconds
Iteration: 1360, time: 1.262 hours, Δt: 2.535 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 46.545 seconds
Iteration: 1400, time: 1.291 hours, Δt: 2.677 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 47.017 seconds
Iteration: 1440, time: 1.318 hours, Δt: 2.595 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 47.514 seconds
Iteration: 1480, time: 1.346 hours, Δt: 2.657 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 47.992 seconds
Iteration: 1520, time: 1.375 hours, Δt: 2.585 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 48.582 seconds
Iteration: 1560, time: 1.403 hours, Δt: 2.539 seconds, max(|w|) = 5.9e-02 ms⁻¹, wall time: 49.060 seconds
Iteration: 1600, time: 1.430 hours, Δt: 2.491 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 49.548 seconds
Iteration: 1640, time: 1.458 hours, Δt: 2.497 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 50.064 seconds
Iteration: 1680, time: 1.485 hours, Δt: 2.492 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 50.547 seconds
Iteration: 1720, time: 1.513 hours, Δt: 2.595 seconds, max(|w|) = 7.3e-02 ms⁻¹, wall time: 51.165 seconds
Iteration: 1760, time: 1.540 hours, Δt: 2.618 seconds, max(|w|) = 7.9e-02 ms⁻¹, wall time: 51.659 seconds
Iteration: 1800, time: 1.567 hours, Δt: 2.487 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 52.184 seconds
Iteration: 1840, time: 1.595 hours, Δt: 2.591 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 52.655 seconds
Iteration: 1880, time: 1.622 hours, Δt: 2.505 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 53.180 seconds
Iteration: 1920, time: 1.649 hours, Δt: 2.486 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 53.737 seconds
Iteration: 1960, time: 1.676 hours, Δt: 2.531 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 54.224 seconds
Iteration: 2000, time: 1.704 hours, Δt: 2.648 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 54.717 seconds
Iteration: 2040, time: 1.732 hours, Δt: 2.478 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 55.176 seconds
Iteration: 2080, time: 1.759 hours, Δt: 2.485 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 55.750 seconds
Iteration: 2120, time: 1.786 hours, Δt: 2.577 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 56.240 seconds
Iteration: 2160, time: 1.814 hours, Δt: 2.567 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 56.772 seconds
Iteration: 2200, time: 1.841 hours, Δt: 2.518 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 57.310 seconds
Iteration: 2240, time: 1.867 hours, Δt: 2.440 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 57.862 seconds
Iteration: 2280, time: 1.894 hours, Δt: 2.428 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 58.522 seconds
Iteration: 2320, time: 1.921 hours, Δt: 2.476 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 59.066 seconds
Iteration: 2360, time: 1.948 hours, Δt: 2.303 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 59.608 seconds
Iteration: 2400, time: 1.973 hours, Δt: 2.363 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 1.002 minutes
Iteration: 2440, time: 1.999 hours, Δt: 2.451 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 1.010 minutes
[ Info: Simulation is stopping after running for 1.010 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.0530265, min=-0.0583849, mean=-1.73644e-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.5886, 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.0134, 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.0178005, min=0.0, mean=0.00045223)

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.