Eady turbulence example

In this example, we initialize a random velocity field and observe its viscous, turbulent decay in a two-dimensional domain. This example demonstrates:

  • How to use a tuple of turbulence closures
  • How to use hyperdiffusivity
  • How to implement background velocity and tracer distributions
  • How to build Fields that compute output

Install dependencies

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

using Pkg
pkg"add Oceananigans, Plots"

The Eady problem

The "Eady problem" simulates the baroclinic instability problem proposed by Eric Eady in the classic paper "Long waves and cyclone waves," Tellus (1949). The Eady problem is a simple, canonical model for the generation of mid-latitude atmospheric storms and the ocean eddies that permeate the world sea.

In the Eady problem, baroclinic motion and turublence is generated by the interaction between turbulent motions and a stationary, geostrophically-balanced basic state that is unstable to baroclinic instability. In this example, the baroclinic generation of turbulence due to extraction of energy from the geostrophic basic state is balanced by a bottom boundary condition that extracts momentum from turbulent motions and serves as a crude model for the drag associated with an unresolved and small-scale turbulent bottom boundary layer.

The geostrophic basic state

The geostrophic basic state in the Eady problem is represented by the streamfunction,

\[ψ(y, z) = - α y (z + L_z) \, ,\]

where $α$ is the geostrophic shear and $L_z$ is the depth of the domain. The background buoyancy includes both the geostrophic flow component, $f ∂_z ψ$, where $f$ is the Coriolis parameter, and a background stable stratification component, $N^2 z$, where $N$ is the buoyancy frequency:

\[B(y, z) = f ∂_z ψ + N^2 z = - α f y + N^2 z \, .\]

The background velocity field is related to the geostrophic streamfunction via $U = - ∂_y ψ$ such that

\[U(z) = α (z + L_z) \, .\]

Boundary conditions

All fields are periodic in the horizontal directions. We use "insulating", or zero-flux boundary conditions on the buoyancy perturbation at the top and bottom. We thus implicitly assume that the background vertical density gradient, $N^2 z$, is maintained by a process external to our simulation. We use free-slip, or zero-flux boundary conditions on $u$ and $v$ at the surface where $z=0$. At the bottom, we impose a momentum flux that extracts momentum and energy from the flow.

Bottom boundary condition: quadratic bottom drag

We model the effects of a turbulent bottom boundary layer on the eddy momentum budget with quadratic bottom drag. A quadratic cottom drag is introduced by imposing a vertical flux of horizontal momentum that removes momentum from the layer immediately above: in other words, the flux is negative (downwards) when the velocity at the bottom boundary is positive, and positive (upwards) with the velocity at the bottom boundary is negative. This drag term is "quadratic" because the rate at which momentum is removed is proportional to $\boldsymbol{u} |\boldsymbol{u}|$, where $\boldsymbol{u} = u \boldsymbol{\hat{x}} + v \boldsymbol{\hat{y}}$ is the horizontal velocity.

The $x$-component of the quadratic bottom drag is thus

\[\tau_{xz}(z=L_z) = - c^D u \sqrt{u^2 + v^2} \, ,\]

while the $y$-component is

\[\tau_{yz}(z=L_z) = - c^D v \sqrt{u^2 + v^2} \, ,\]

where $c^D$ is a dimensionless drag coefficient and $\tau_{xz}(z=L_z)$ and $\tau_{yz}(z=L_z)$ denote the flux of $u$ and $v$ momentum at $z = L_z$, the bottom of the domain.

Vertical and horizontal viscosity and diffusivity

Vertical and horizontal viscosities and diffusivities are required to stabilize the Eady problem and can be idealized as modeling the effect of turbulent mixing below the grid scale. For both tracers and velocities we use a Laplacian vertical diffusivity $κ_z ∂_z^2 c$ and a horizontal hyperdiffusivity $ϰ_h (∂_x^4 + ∂_y^4) c$.

Eady problem summary and parameters

To summarize, the Eady problem parameters along with the values we use in this example are

Parameter nameDescriptionValueUnits
$f$Coriolis parameter$10^{-4}$$\mathrm{s^{-1}}$
$N$Buoyancy frequency (square root of $\partial_z B$)$10^{-3}$$\mathrm{s^{-1}}$
$\alpha$Background vertical shear $\partial_z U$$10^{-3}$$\mathrm{s^{-1}}$
$c^D$Bottom quadratic drag coefficient$10^{-4}$none
$κ_z$Laplacian vertical diffusivity$10^{-2}$$\mathrm{m^2 s^{-1}}$
$ϰ_h$Biharmonic horizontal diffusivity$10^{-2} \times \Delta x^4 / \mathrm{day}$$\mathrm{m^4 s^{-1}}$

We start off by importing Oceananigans, Printf, and some convenient constants for specifying dimensional units:

using Printf
using Oceananigans
using Oceananigans.Units: hours, day, days

The grid

We use a three-dimensional grid with a depth of 4000 m and a horizontal extent of 1000 km, appropriate for mesoscale ocean dynamics with characteristic scales of 50-200 km.

grid = RectilinearGrid(size=(48, 48, 16), x=(0, 1e6), y=(0, 1e6), z=(-4e3, 0))
48×48×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── Periodic x ∈ [0.0, 1.0e6)   regularly spaced with Δx=20833.3
├── Periodic y ∈ [0.0, 1.0e6)   regularly spaced with Δy=20833.3
└── Bounded  z ∈ [-4000.0, 0.0] regularly spaced with Δz=250.0

Rotation

The classical Eady problem is posed on an $f$-plane. We use a Coriolis parameter typical to mid-latitudes on Earth,

coriolis = FPlane(f=1e-4) # [s⁻¹]
FPlane{Float64}(f=0.0001)

The background flow

We build a NamedTuple of parameters that describe the background flow,

basic_state_parameters = ( α = 10 * coriolis.f, # s⁻¹, geostrophic shear
                           f = coriolis.f,      # s⁻¹, Coriolis parameter
                           N = 1e-3,            # s⁻¹, buoyancy frequency
                          Lz = grid.Lz)         # m, ocean depth
(α = 0.001, f = 0.0001, N = 0.001, Lz = 4000.0)

and then construct the background fields $U$ and $B$

# Background fields are defined via functions of x, y, z, t, and optional parameters
U(x, y, z, t, p) = + p.α * (z + p.Lz)
B(x, y, z, t, p) = - p.α * p.f * y + p.N^2 * z

U_field = BackgroundField(U, parameters=basic_state_parameters)
B_field = BackgroundField(B, parameters=basic_state_parameters)
BackgroundField{typeof(Main.B), NamedTuple{(:α, :f, :N, :Lz), NTuple{4, Float64}}}
├── func: B (generic function with 1 method)
└── parameters: (α = 0.001, f = 0.0001, N = 0.001, Lz = 4000.0)

Boundary conditions

The boundary conditions prescribe a quadratic drag at the bottom as a flux condition.

cᴰ = 1e-4 # quadratic drag coefficient

@inline drag_u(x, y, t, u, v, cᴰ) = - cᴰ * u * sqrt(u^2 + v^2)
@inline drag_v(x, y, t, u, v, cᴰ) = - cᴰ * v * sqrt(u^2 + v^2)

drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=cᴰ)
drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=cᴰ)

u_bcs = FieldBoundaryConditions(bottom = drag_bc_u)
v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition
├── east: DefaultBoundaryCondition
├── south: DefaultBoundaryCondition
├── north: DefaultBoundaryCondition
├── bottom: FluxBoundaryCondition: ContinuousBoundaryFunction drag_v at (Nothing, Nothing, Nothing)
├── top: DefaultBoundaryCondition
└── immersed: FluxBoundaryCondition: Nothing

Turbulence closures

We use a horizontal hyperdiffusivity and a Laplacian vertical diffusivity to dissipate energy in the Eady problem. To use both of these closures at the same time, we set the keyword argument closure to a tuple of two closures.

κ₂z = 1e-2 # [m² s⁻¹] Laplacian vertical viscosity and diffusivity
κ₄h = 1e-1 / day * grid.Δxᶜᵃᵃ^4 # [m⁴ s⁻¹] horizontal hyperviscosity and hyperdiffusivity

Laplacian_vertical_diffusivity = VerticalScalarDiffusivity(ν=κ₂z, κ=κ₂z)
biharmonic_horizontal_diffusivity = HorizontalScalarBiharmonicDiffusivity(ν=κ₄h, κ=κ₄h)
ScalarBiharmonicDiffusivity{HorizontalFormulation}(ν=2.18033e11, κ=2.18033e11)

Model instantiation

We instantiate the model with the fifth-order WENO advection scheme, a 3rd order Runge-Kutta time-stepping scheme, and a BuoyancyTracer.

model = NonhydrostaticModel(
                   grid = grid,
              advection = WENO5(),
            timestepper = :RungeKutta3,
               coriolis = coriolis,
                tracers = :b,
               buoyancy = BuoyancyTracer(),
      background_fields = (b=B_field, u=U_field),
                closure = (Laplacian_vertical_diffusivity, biharmonic_horizontal_diffusivity),
    boundary_conditions = (u=u_bcs, v=v_bcs)
)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 48×48×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── tracers: b
├── closure: Tuple with 2 closures:
│   ├── ScalarBiharmonicDiffusivity{HorizontalFormulation}(ν=2.18033e11, κ=(b=2.18033e11,))
│   └── VerticalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.01, κ=(b=0.01,))
├── buoyancy: BuoyancyTracer with -ĝ = ZDirection
└── coriolis: FPlane{Float64}(f=0.0001)

Initial conditions

We seed our initial conditions with random noise stimulate the growth of baroclinic instability.

# A noise function, damped at the top and bottom
Ξ(z) = randn() * z/grid.Lz * (z/grid.Lz + 1)

# Scales for the initial velocity and buoyancy
Ũ = 1e-1 * basic_state_parameters.α * grid.Lz
B̃ = 1e-2 * basic_state_parameters.α * coriolis.f

uᵢ(x, y, z) = Ũ * Ξ(z)
vᵢ(x, y, z) = Ũ * Ξ(z)
bᵢ(x, y, z) = B̃ * Ξ(z)

set!(model, u=uᵢ, v=vᵢ, b=bᵢ)

We subtract off any residual mean velocity to avoid exciting domain-scale inertial oscillations. We use a sum over the entire parent arrays or data to ensure this operation is efficient on the GPU (set architecture = GPU() in NonhydrostaticModel constructor to run this problem on the GPU if one is available).

ū = sum(model.velocities.u.data.parent) / (grid.Nx * grid.Ny * grid.Nz)
v̄ = sum(model.velocities.v.data.parent) / (grid.Nx * grid.Ny * grid.Nz)

model.velocities.u.data.parent .-= ū
model.velocities.v.data.parent .-= v̄

Simulation set-up

We set up a simulation that runs for 10 days with a JLD2OutputWriter that saves the vertical vorticity and divergence every 2 hours. We limit the time-step to the maximum allowable due either to diffusion, internal waves, or advection by the background flow.

# Calculate absolute limit on time-step using diffusivities and
# background velocity.
Ū = basic_state_parameters.α * grid.Lz

max_Δt = min(grid.Δxᶜᵃᵃ / Ū, grid.Δxᶜᵃᵃ^4 / κ₄h, grid.Δxᶜᵃᵃ^2 / κ₂z, 1 / basic_state_parameters.N)

simulation = Simulation(model, Δt = max_Δt, stop_time = 8days)
Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 16.667 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN years
├── Stop time: 8 days
├── Stop iteration : Inf
├── Wall time limit: Inf
├── 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
└── Diagnostics: OrderedDict with no entries

The TimeStepWizard

The TimeStepWizard manages the time-step adaptively, keeping the Courant-Freidrichs-Lewy (CFL) number close to 1.0 while ensuring the time-step does not increase beyond the maximum allowable value for numerical stability given the specified background flow, internal wave time scales, and diffusion time scales.

wizard = TimeStepWizard(cfl=0.85, max_change=1.1, max_Δt=max_Δt)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))
Callback of TimeStepWizard(cfl=0.85, max_Δt=1000.0, min_Δt=0.0) on IterationInterval(10)

A progress messenger

We add a callback that prints out a helpful progress message while the simulation runs.

start_time = time_ns()

progress(sim) = @printf("i: % 6d, sim time: % 10s, wall time: % 10s, Δt: % 10s, CFL: %.2e\n",
                        sim.model.clock.iteration,
                        prettytime(sim.model.clock.time),
                        prettytime(1e-9 * (time_ns() - start_time)),
                        prettytime(sim.Δt),
                        AdvectiveCFL(sim.Δt)(sim.model))

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))
Callback of progress on IterationInterval(10)

Output

To visualize the baroclinic turbulence ensuing in the Eady problem, we use computed Fields to diagnose and output vertical vorticity and divergence. Note that computed Fields take "AbstractOperations" on Fields as input:

u, v, w = model.velocities # unpack velocity `Field`s

# Vertical vorticity [s⁻¹]
ζ = ∂x(v) - ∂y(u)

# Horizontal divergence, or ∂x(u) + ∂y(v) [s⁻¹]
δ = -∂z(w)
UnaryOperation at (Center, Center, Center)
├── grid: 48×48×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree: 
    - at (Center, Center, Center) via identity
    └── ∂zᶜᶜᶜ at (Center, Center, Center) via identity
        └── 48×48×17 Field{Center, Center, Face} on RectilinearGrid on CPU

With the vertical vorticity, ζ, and the horizontal divergence, δ in hand, we create a JLD2OutputWriter that saves ζ and δ and add them to simulation.

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; ζ, δ),
                                                      schedule = TimeInterval(4hours),
                                                        prefix = "eady_turbulence",
                                                         force = true)

All that's left is to press the big red button:

run!(simulation)
[ Info: Initializing simulation...
i:      0, sim time:  0 seconds, wall time: 1.146 minutes, Δt: 16.667 minutes, CFL: 4.82e-02
[ Info:     ... simulation initialization complete (9.700 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (6.210 minutes).
i:     10, sim time: 2.778 hours, wall time: 7.526 minutes, Δt: 16.667 minutes, CFL: 1.84e-02
i:     20, sim time: 5.389 hours, wall time: 7.571 minutes, Δt: 16.667 minutes, CFL: 1.68e-02
i:     30, sim time:    8 hours, wall time: 7.610 minutes, Δt: 16.667 minutes, CFL: 1.77e-02
i:     40, sim time: 10.778 hours, wall time: 7.649 minutes, Δt: 16.667 minutes, CFL: 1.95e-02
i:     50, sim time: 13.389 hours, wall time: 7.688 minutes, Δt: 16.667 minutes, CFL: 1.87e-02
i:     60, sim time:   16 hours, wall time: 7.726 minutes, Δt: 16.667 minutes, CFL: 2.23e-02
i:     70, sim time: 18.778 hours, wall time: 7.765 minutes, Δt: 16.667 minutes, CFL: 2.46e-02
i:     80, sim time: 21.389 hours, wall time: 7.803 minutes, Δt: 16.667 minutes, CFL: 2.65e-02
i:     90, sim time:      1 day, wall time: 7.842 minutes, Δt: 16.667 minutes, CFL: 2.83e-02
i:    100, sim time: 1.116 days, wall time: 7.883 minutes, Δt: 16.667 minutes, CFL: 3.15e-02
i:    110, sim time: 1.225 days, wall time: 7.922 minutes, Δt: 16.667 minutes, CFL: 3.59e-02
i:    120, sim time: 1.333 days, wall time: 7.963 minutes, Δt: 16.667 minutes, CFL: 3.89e-02
i:    130, sim time: 1.449 days, wall time: 8.004 minutes, Δt: 16.667 minutes, CFL: 4.44e-02
i:    140, sim time: 1.558 days, wall time: 8.044 minutes, Δt: 16.667 minutes, CFL: 5.76e-02
i:    150, sim time: 1.667 days, wall time: 8.083 minutes, Δt: 16.667 minutes, CFL: 6.02e-02
i:    160, sim time: 1.782 days, wall time: 8.123 minutes, Δt: 16.667 minutes, CFL: 6.38e-02
i:    170, sim time: 1.891 days, wall time: 8.163 minutes, Δt: 16.667 minutes, CFL: 6.72e-02
i:    180, sim time:     2 days, wall time: 8.223 minutes, Δt: 16.667 minutes, CFL: 6.86e-02
i:    190, sim time: 2.116 days, wall time: 8.269 minutes, Δt: 16.667 minutes, CFL: 7.63e-02
i:    200, sim time: 2.225 days, wall time: 8.331 minutes, Δt: 16.667 minutes, CFL: 8.39e-02
i:    210, sim time: 2.333 days, wall time: 8.387 minutes, Δt: 16.667 minutes, CFL: 1.06e-01
i:    220, sim time: 2.449 days, wall time: 8.431 minutes, Δt: 16.667 minutes, CFL: 1.30e-01
i:    230, sim time: 2.558 days, wall time: 8.481 minutes, Δt: 16.667 minutes, CFL: 1.53e-01
i:    240, sim time: 2.667 days, wall time: 8.522 minutes, Δt: 16.667 minutes, CFL: 1.75e-01
i:    250, sim time: 2.782 days, wall time: 8.562 minutes, Δt: 16.667 minutes, CFL: 1.98e-01
i:    260, sim time: 2.891 days, wall time: 8.601 minutes, Δt: 16.667 minutes, CFL: 2.26e-01
i:    270, sim time:     3 days, wall time: 8.640 minutes, Δt: 16.667 minutes, CFL: 2.56e-01
i:    280, sim time: 3.116 days, wall time: 8.680 minutes, Δt: 16.667 minutes, CFL: 2.83e-01
i:    290, sim time: 3.225 days, wall time: 8.722 minutes, Δt: 16.667 minutes, CFL: 3.18e-01
i:    300, sim time: 3.333 days, wall time: 8.765 minutes, Δt: 16.667 minutes, CFL: 3.62e-01
i:    310, sim time: 3.449 days, wall time: 8.803 minutes, Δt: 16.667 minutes, CFL: 4.04e-01
i:    320, sim time: 3.558 days, wall time: 8.842 minutes, Δt: 16.667 minutes, CFL: 4.52e-01
i:    330, sim time: 3.667 days, wall time: 8.881 minutes, Δt: 16.667 minutes, CFL: 5.11e-01
i:    340, sim time: 3.782 days, wall time: 8.920 minutes, Δt: 16.667 minutes, CFL: 5.51e-01
i:    350, sim time: 3.891 days, wall time: 8.959 minutes, Δt: 16.667 minutes, CFL: 6.13e-01
i:    360, sim time:     4 days, wall time: 8.998 minutes, Δt: 16.667 minutes, CFL: 6.29e-01
i:    370, sim time: 4.116 days, wall time: 9.037 minutes, Δt: 16.667 minutes, CFL: 6.00e-01
i:    380, sim time: 4.225 days, wall time: 9.076 minutes, Δt: 16.667 minutes, CFL: 5.46e-01
i:    390, sim time: 4.333 days, wall time: 9.116 minutes, Δt: 16.667 minutes, CFL: 5.68e-01
i:    400, sim time: 4.449 days, wall time: 9.155 minutes, Δt: 16.667 minutes, CFL: 5.61e-01
i:    410, sim time: 4.558 days, wall time: 9.195 minutes, Δt: 16.667 minutes, CFL: 5.38e-01
i:    420, sim time: 4.667 days, wall time: 9.233 minutes, Δt: 16.667 minutes, CFL: 5.63e-01
i:    430, sim time: 4.782 days, wall time: 9.287 minutes, Δt: 16.667 minutes, CFL: 5.45e-01
i:    440, sim time: 4.891 days, wall time: 9.338 minutes, Δt: 16.667 minutes, CFL: 5.62e-01
i:    450, sim time:     5 days, wall time: 9.389 minutes, Δt: 16.667 minutes, CFL: 6.43e-01
i:    460, sim time: 5.116 days, wall time: 9.440 minutes, Δt: 16.667 minutes, CFL: 7.97e-01
i:    470, sim time: 5.225 days, wall time: 9.491 minutes, Δt: 15.853 minutes, CFL: 8.50e-01
i:    480, sim time: 5.333 days, wall time: 9.540 minutes, Δt: 16.213 minutes, CFL: 8.50e-01
i:    490, sim time: 5.446 days, wall time: 9.591 minutes, Δt: 16.667 minutes, CFL: 7.00e-01
i:    500, sim time: 5.558 days, wall time: 9.638 minutes, Δt: 16.667 minutes, CFL: 5.09e-01
i:    510, sim time: 5.667 days, wall time: 9.676 minutes, Δt: 16.667 minutes, CFL: 4.93e-01
i:    520, sim time: 5.782 days, wall time: 9.717 minutes, Δt: 16.667 minutes, CFL: 5.50e-01
i:    530, sim time: 5.891 days, wall time: 9.756 minutes, Δt: 16.667 minutes, CFL: 5.46e-01
i:    540, sim time:     6 days, wall time: 9.797 minutes, Δt: 16.667 minutes, CFL: 5.25e-01
i:    550, sim time: 6.116 days, wall time: 9.837 minutes, Δt: 16.667 minutes, CFL: 4.88e-01
i:    560, sim time: 6.225 days, wall time: 9.877 minutes, Δt: 16.667 minutes, CFL: 4.39e-01
i:    570, sim time: 6.333 days, wall time: 9.918 minutes, Δt: 16.667 minutes, CFL: 4.35e-01
i:    580, sim time: 6.449 days, wall time: 9.959 minutes, Δt: 16.667 minutes, CFL: 4.76e-01
i:    590, sim time: 6.558 days, wall time: 9.999 minutes, Δt: 16.667 minutes, CFL: 5.10e-01
i:    600, sim time: 6.667 days, wall time: 10.039 minutes, Δt: 16.667 minutes, CFL: 6.24e-01
i:    610, sim time: 6.782 days, wall time: 10.079 minutes, Δt: 16.667 minutes, CFL: 5.75e-01
i:    620, sim time: 6.891 days, wall time: 10.118 minutes, Δt: 16.667 minutes, CFL: 6.00e-01
i:    630, sim time:     7 days, wall time: 10.159 minutes, Δt: 16.667 minutes, CFL: 5.32e-01
i:    640, sim time: 7.116 days, wall time: 10.198 minutes, Δt: 16.667 minutes, CFL: 5.25e-01
i:    650, sim time: 7.225 days, wall time: 10.238 minutes, Δt: 16.667 minutes, CFL: 5.22e-01
i:    660, sim time: 7.333 days, wall time: 10.278 minutes, Δt: 16.667 minutes, CFL: 5.25e-01
i:    670, sim time: 7.449 days, wall time: 10.318 minutes, Δt: 16.667 minutes, CFL: 5.23e-01
i:    680, sim time: 7.558 days, wall time: 10.358 minutes, Δt: 16.667 minutes, CFL: 5.30e-01
i:    690, sim time: 7.667 days, wall time: 10.398 minutes, Δt: 16.667 minutes, CFL: 6.89e-01
i:    700, sim time: 7.782 days, wall time: 10.436 minutes, Δt: 16.667 minutes, CFL: 7.16e-01
i:    710, sim time: 7.891 days, wall time: 10.478 minutes, Δt: 16.667 minutes, CFL: 5.81e-01
[ Info: Simulation is stopping. Model time 8 days has hit or exceeded simulation stop time 8 days.
i:    720, sim time:     8 days, wall time: 10.517 minutes, Δt: 16.667 minutes, CFL: 6.25e-01

Visualizing Eady turbulence

We animate slices of FieldTimeSeries of the saved data. To prepare the animation we create coordinate arrays and defining a function for computing colorbar limits:

using Plots

filepath = simulation.output_writers[:fields].filepath

ζ_timeseries = FieldTimeSeries(filepath, "ζ")
δ_timeseries = FieldTimeSeries(filepath, "δ")

# Times and spatial coordinates
times = ζ_timeseries.times
xζ, yζ, zζ = nodes(ζ_timeseries)
xδ, yδ, zδ = nodes(δ_timeseries)
([10416.666666666666, 31250.0, 52083.333333333336, 72916.66666666667, 93750.0, 114583.33333333333, 135416.66666666666, 156250.0, 177083.33333333334, 197916.66666666666  …  802083.3333333334, 822916.6666666666, 843750.0, 864583.3333333334, 885416.6666666666, 906250.0, 927083.3333333334, 947916.6666666666, 968750.0, 989583.3333333334], [10416.666666666666, 31250.0, 52083.333333333336, 72916.66666666667, 93750.0, 114583.33333333333, 135416.66666666666, 156250.0, 177083.33333333334, 197916.66666666666  …  802083.3333333334, 822916.6666666666, 843750.0, 864583.3333333334, 885416.6666666666, 906250.0, 927083.3333333334, 947916.6666666666, 968750.0, 989583.3333333334], [-3875.0, -3625.0, -3375.0, -3125.0, -2875.0, -2625.0, -2375.0, -2125.0, -1875.0, -1625.0, -1375.0, -1125.0, -875.0, -625.0, -375.0, -125.0])

This utility is handy for calculating nice contour intervals:

function nice_divergent_levels(c, clim, nlevels=31)
    levels = range(-clim, stop=clim, length=nlevels)
    cmax = maximum(abs, c)
    clim < cmax && (levels = vcat([-cmax], levels, [cmax]))
    return levels
end

Now we're ready to animate.

@info "Making an animation from saved data..."

anim = @animate for (i, t) in enumerate(times)

    # Load 3D fields from file
    R_snapshot = ζ_timeseries[i] ./ coriolis.f
    δ_snapshot = δ_timeseries[i]

    surface_R = R_snapshot[:, :, grid.Nz]
    surface_δ = δ_snapshot[:, :, grid.Nz]

    slice_R = R_snapshot[:, 1, :]
    slice_δ = δ_snapshot[:, 1, :]

    Rlim = 0.5 * maximum(abs, R_snapshot) + 1e-9
    δlim = 0.5 * maximum(abs, δ_snapshot) + 1e-9

    Rlevels = nice_divergent_levels(R_snapshot, Rlim)
    δlevels = nice_divergent_levels(δ_snapshot, δlim)

    @info @sprintf("Drawing frame %d from iteration %d: max(ζ̃ / f) = %.3f \n",
                   i, iter, maximum(abs, surface_R))

    xy_kwargs = (xlims = (0, grid.Lx), ylims = (0, grid.Lx),
                 xlabel = "x (m)", ylabel = "y (m)",
                 aspectratio = 1,
                 linewidth = 0, color = :balance, legend = false)

    xz_kwargs = (xlims = (0, grid.Lx), ylims = (-grid.Lz, 0),
                 xlabel = "x (m)", ylabel = "z (m)",
                 aspectratio = grid.Lx / grid.Lz * 0.5,
                 linewidth = 0, color = :balance, legend = false)

    R_xy = contourf(xζ, yζ, surface_R'; clims=(-Rlim, Rlim), levels=Rlevels, xy_kwargs...)
    δ_xy = contourf(xδ, yδ, surface_δ'; clims=(-δlim, δlim), levels=δlevels, xy_kwargs...)
    R_xz = contourf(xζ, zζ, slice_R'; clims=(-Rlim, Rlim), levels=Rlevels, xz_kwargs...)
    δ_xz = contourf(xδ, zδ, slice_δ'; clims=(-δlim, δlim), levels=δlevels, xz_kwargs...)

    plot(R_xy, δ_xy, R_xz, δ_xz,
           size = (1000, 800),
           link = :x,
         layout = Plots.grid(2, 2, heights=[0.5, 0.5, 0.2, 0.2]),
          title = [@sprintf("ζ(t=%s) / f", prettytime(t)) @sprintf("δ(t=%s) (s⁻¹)", prettytime(t)) "" ""])
end

This page was generated using Literate.jl.