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.7 KiB

We're ready:

run!(simulation)

[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|w|) = 1.3e-05 ms⁻¹, wall time: 0 seconds
[ Info:     ... simulation initialization complete (15.312 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (5.147 seconds).
Iteration: 0040, time: 7 minutes, Δt: 8.309 seconds, max(|w|) = 2.5e-05 ms⁻¹, wall time: 21.326 seconds
Iteration: 0080, time: 11.182 minutes, Δt: 4.833 seconds, max(|w|) = 7.9e-03 ms⁻¹, wall time: 21.860 seconds
Iteration: 0120, time: 13.957 minutes, Δt: 4.085 seconds, max(|w|) = 2.3e-02 ms⁻¹, wall time: 22.933 seconds
Iteration: 0160, time: 16.604 minutes, Δt: 3.877 seconds, max(|w|) = 2.3e-02 ms⁻¹, wall time: 23.496 seconds
Iteration: 0200, time: 19.069 minutes, Δt: 4.055 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 23.985 seconds
Iteration: 0240, time: 21.578 minutes, Δt: 3.980 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 24.474 seconds
Iteration: 0280, time: 24 minutes, Δt: 3.535 seconds, max(|w|) = 2.4e-02 ms⁻¹, wall time: 25.010 seconds
Iteration: 0320, time: 26.247 minutes, Δt: 3.634 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 25.518 seconds
Iteration: 0360, time: 28.513 minutes, Δt: 3.402 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 26.010 seconds
Iteration: 0400, time: 30.730 minutes, Δt: 3.376 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 26.500 seconds
Iteration: 0440, time: 32.821 minutes, Δt: 3.361 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 27.002 seconds
Iteration: 0480, time: 35 minutes, Δt: 3.326 seconds, max(|w|) = 3.6e-02 ms⁻¹, wall time: 27.517 seconds
Iteration: 0520, time: 37.167 minutes, Δt: 3.337 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 28.079 seconds
Iteration: 0560, time: 39.201 minutes, Δt: 3.000 seconds, max(|w|) = 3.4e-02 ms⁻¹, wall time: 28.575 seconds
Iteration: 0600, time: 41.195 minutes, Δt: 2.954 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 29.099 seconds
Iteration: 0640, time: 43.101 minutes, Δt: 3.056 seconds, max(|w|) = 4.0e-02 ms⁻¹, wall time: 29.594 seconds
Iteration: 0680, time: 45.099 minutes, Δt: 2.966 seconds, max(|w|) = 4.0e-02 ms⁻¹, wall time: 30.089 seconds
Iteration: 0720, time: 47 minutes, Δt: 2.829 seconds, max(|w|) = 3.7e-02 ms⁻¹, wall time: 30.542 seconds
Iteration: 0760, time: 48.911 minutes, Δt: 2.776 seconds, max(|w|) = 3.4e-02 ms⁻¹, wall time: 31.035 seconds
Iteration: 0800, time: 50.751 minutes, Δt: 2.878 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 31.531 seconds
Iteration: 0840, time: 52.560 minutes, Δt: 2.758 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 32.005 seconds
Iteration: 0880, time: 54.307 minutes, Δt: 2.721 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 32.510 seconds
Iteration: 0920, time: 56.094 minutes, Δt: 2.820 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 33.064 seconds
Iteration: 0960, time: 57.896 minutes, Δt: 2.734 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 33.572 seconds
Iteration: 1000, time: 59.676 minutes, Δt: 2.669 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 34.042 seconds
Iteration: 1040, time: 1.023 hours, Δt: 2.629 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 34.502 seconds
Iteration: 1080, time: 1.051 hours, Δt: 2.653 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 35.001 seconds
Iteration: 1120, time: 1.081 hours, Δt: 2.498 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 35.460 seconds
Iteration: 1160, time: 1.109 hours, Δt: 2.613 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 35.959 seconds
Iteration: 1200, time: 1.137 hours, Δt: 2.663 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 36.464 seconds
Iteration: 1240, time: 1.166 hours, Δt: 2.549 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 36.925 seconds
Iteration: 1280, time: 1.193 hours, Δt: 2.549 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 37.425 seconds
Iteration: 1320, time: 1.221 hours, Δt: 2.445 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 37.937 seconds
Iteration: 1360, time: 1.248 hours, Δt: 2.460 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 38.420 seconds
Iteration: 1400, time: 1.275 hours, Δt: 2.618 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 38.913 seconds
Iteration: 1440, time: 1.304 hours, Δt: 2.570 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 39.469 seconds
Iteration: 1480, time: 1.332 hours, Δt: 2.489 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 39.976 seconds
Iteration: 1520, time: 1.358 hours, Δt: 2.442 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 40.464 seconds
Iteration: 1560, time: 1.384 hours, Δt: 2.463 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 40.997 seconds
Iteration: 1600, time: 1.411 hours, Δt: 2.396 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 41.492 seconds
Iteration: 1640, time: 1.437 hours, Δt: 2.429 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 41.985 seconds
Iteration: 1680, time: 1.464 hours, Δt: 2.535 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 42.455 seconds
Iteration: 1720, time: 1.492 hours, Δt: 2.517 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 42.970 seconds
Iteration: 1760, time: 1.519 hours, Δt: 2.523 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 43.554 seconds
Iteration: 1800, time: 1.547 hours, Δt: 2.527 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 44.027 seconds
Iteration: 1840, time: 1.574 hours, Δt: 2.443 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 44.573 seconds
Iteration: 1880, time: 1.600 hours, Δt: 2.493 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 45.106 seconds
Iteration: 1920, time: 1.628 hours, Δt: 2.381 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 45.610 seconds
Iteration: 1960, time: 1.653 hours, Δt: 2.435 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 46.125 seconds
Iteration: 2000, time: 1.680 hours, Δt: 2.548 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 46.619 seconds
Iteration: 2040, time: 1.707 hours, Δt: 2.406 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 47.143 seconds
Iteration: 2080, time: 1.733 hours, Δt: 2.411 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 47.607 seconds
Iteration: 2120, time: 1.760 hours, Δt: 2.437 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 48.133 seconds
Iteration: 2160, time: 1.786 hours, Δt: 2.444 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 48.645 seconds
Iteration: 2200, time: 1.813 hours, Δt: 2.480 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 49.109 seconds
Iteration: 2240, time: 1.840 hours, Δt: 2.355 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 49.594 seconds
Iteration: 2280, time: 1.865 hours, Δt: 2.306 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 50.088 seconds
Iteration: 2320, time: 1.889 hours, Δt: 2.302 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 50.580 seconds
Iteration: 2360, time: 1.914 hours, Δt: 2.276 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 51.077 seconds
Iteration: 2400, time: 1.939 hours, Δt: 2.320 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 51.579 seconds
Iteration: 2440, time: 1.966 hours, Δt: 2.489 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 52.066 seconds
Iteration: 2480, time: 1.992 hours, Δt: 2.456 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 52.594 seconds
[ Info: Simulation is stopping after running for 52.761 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.0424928, min=-0.0529756, mean=4.27025e-6, T = 128×1×64×121 FieldTimeSeries{InMemory} located at (Center, Center, Center) of T at ocean_wind_mixing_and_convection.jld2
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: T
└── data: 134×1×70×121 OffsetArray(::Array{Float64, 4}, -2:131, 64:64, -2:67, 1:121) with eltype Float64 with indices -2:131×64:64×-2:67×1:121
    └── max=19.9958, min=0.0, mean=18.5885, S = 128×1×64×121 FieldTimeSeries{InMemory} located at (Center, Center, Center) of S at ocean_wind_mixing_and_convection.jld2
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: S
└── data: 134×1×70×121 OffsetArray(::Array{Float64, 4}, -2:131, 64:64, -2:67, 1:121) with eltype Float64 with indices -2:131×64:64×-2:67×1:121
    └── max=35.013, 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.0191207, min=0.0, mean=0.000414893)

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.