Edit on GitHub

Convecting plankton

In this example, two-dimensional convection into a stratified fluid mixes a phytoplankton-like tracer. This example demonstrates:

  • How to set boundary conditions.
  • How to insert a user-defined forcing function into a model.
  • How to use the TimeStepWizard to adapt the simulation time-step.
  • How to use AveragedField to diagnose spatial averages of model fields.

The grid

We use a two-dimensional grid with 128² points and 2 m grid spacing:

using Oceananigans

grid = RegularCartesianGrid(size=(128, 1, 128), extent=(64, 1, 64))
RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}
                   domain: x ∈ [0.0, 64.0], y ∈ [0.0, 1.0], z ∈ [-64.0, 0.0]
                 topology: (Periodic, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (128, 1, 128)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (0.5, 1.0, 0.5)

Boundary conditions

We impose buoyancy loss at the surface with the buoyancy flux

Qb = 1e-8 # m³ s⁻²

Note that a positive flux at the top boundary means that buoyancy is carried upwards, out of the fluid. This reduces the fluid's buoyancy near the surface, causing convection.

The initial condition consists of the constant buoyancy gradient

N² = 1e-6 # s⁻²

which we also impose as a boundary condition at the bottom.

buoyancy_bcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, Qb),
                                              bottom = BoundaryCondition(Gradient, N²))
Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions
├── x: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
├── y: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
└── z: CoordinateBoundaryConditions{BoundaryCondition{Gradient,Float64},BoundaryCondition{Flux,Float64}}

Forcing: the growing and grazing of phytoplankton

We add a forcing term to the plankton equation that crudely models:

  • the growth of phytoplankton at $γ = 10$ "phytoplankton units" per day, in sunlight with a penetration depth of $λ = 16$ meters;

  • death due to viruses and grazing by zooplankton at a rate of $μ = 1$ phytoplankton unit per day.

using Oceananigans.Utils

growing_and_grazing(x, y, z, t, p) = p.γ * exp(z / p.λ) - p.μ

plankton_forcing = Forcing(growing_and_grazing, parameters=(λ=16, γ=10/day, μ=1/day))
ContinuousForcing{NamedTuple{(:λ, :γ, :μ),Tuple{Int64,Float64,Float64}}}
├── func: growing_and_grazing
├── parameters: (λ = 16, γ = 0.00011574074074074075, μ = 1.1574074074074073e-5)
└── field dependencies: ()

Finally, we're ready to build an IncompressibleModel. We use a third-order advection scheme, third-order Runge-Kutta time-stepping, and add Coriolis forces appropriate for planktonic convection at mid-latitudes on Earth.

using Oceananigans.Advection: UpwindBiasedThirdOrder

model = IncompressibleModel(
                   grid = grid,
              advection = UpwindBiasedThirdOrder(),
            timestepper = :RungeKutta3,
                closure = IsotropicDiffusivity(ν=1e-4, κ=1e-4),
               coriolis = FPlane(f=1e-4),
                tracers = (:b, :plankton),
               buoyancy = BuoyancyTracer(),
                forcing = (plankton=plankton_forcing,),
    boundary_conditions = (b=buoyancy_bcs,)
)
IncompressibleModel{CPU, Float64}(time = 0 seconds, iteration = 0) 
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=1, Nz=128)
├── tracers: (:b, :plankton)
├── closure: IsotropicDiffusivity{Float64,NamedTuple{(:b, :plankton),Tuple{Float64,Float64}}}
├── buoyancy: BuoyancyTracer
└── coriolis: FPlane{Float64}

Initial condition

Our initial condition consists of a linear buoyancy gradient superposed with random noise.

Ξ(z) = randn() * z / grid.Lz * (1 + z / grid.Lz) # noise

initial_buoyancy(x, y, z) = N² * z + N² * grid.Lz * 1e-6 * Ξ(z)

set!(model, b=initial_buoyancy)

Simulation setup

We use a TimeStepWizard that limits the time-step to 1 minute, and adapts the time-step such that CFL (Courant-Freidrichs-Lewy) number remains below 1.0,

wizard = TimeStepWizard(cfl=1.0, Δt=2minutes, max_change=1.1, max_Δt=2minutes)
TimeStepWizard{Float64}(1.0, Inf, 1.1, 0.5, 120.0, 120.0)

We also write a function that prints the progress of the simulation

using Printf

progress(sim) = @printf("Iteration: %d, time: %s, Δt: %s\n",
                        sim.model.clock.iteration,
                        prettytime(sim.model.clock.time),
                        prettytime(sim.Δt.Δt))

simulation = Simulation(model, Δt=wizard, stop_time=4hour,
                        iteration_interval=10, progress=progress)
Simulation{IncompressibleModel{CPU, Float64}}
├── Model clock: time = 0 seconds, iteration = 0 
├── Next time step (TimeStepWizard{Float64}): 2 minutes 
├── Iteration interval: 10
├── Stop criteria: Any[Oceananigans.Simulations.iteration_limit_exceeded, Oceananigans.Simulations.stop_time_exceeded, Oceananigans.Simulations.wall_time_limit_exceeded]
├── Run time: 0 seconds, wall time limit: Inf
├── Stop time: 4 hours, stop iteration: Inf
├── Diagnostics: OrderedCollections.OrderedDict with no entries
└── Output writers: OrderedCollections.OrderedDict with no entries

We add a basic JLD2OutputWriter that writes velocities, tracers, and the horizontally-averaged plankton:

using Oceananigans.OutputWriters, Oceananigans.Fields

averaged_plankton = AveragedField(model.tracers.plankton, dims=(1, 2))

outputs = (w = model.velocities.w,
           plankton = model.tracers.plankton,
           averaged_plankton = averaged_plankton)

simulation.output_writers[:fields] =
    JLD2OutputWriter(model, outputs,
                     schedule = TimeInterval(10minute),
                     prefix = "convecting_plankton",
                     force = true)
JLD2OutputWriter scheduled on TimeInterval(10 minutes):
├── filepath: ./convecting_plankton.jld2
├── 3 outputs: (:w, :plankton, :averaged_plankton)
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB

Note that it often makes sense to define different output writers for two- or three-dimensional fields and AveragedFields (since averages take up so much less disk space, it's usually possible to output them a lot more frequently than full fields without blowing up your hard drive).

The simulation is setup. Let there be plankton:

run!(simulation)
Iteration: 10, time: 20 minutes, Δt: 2 minutes
Iteration: 20, time: 40 minutes, Δt: 2 minutes
Iteration: 30, time: 1 hour, Δt: 2 minutes
Iteration: 40, time: 1.333 hours, Δt: 2 minutes
Iteration: 50, time: 1.667 hours, Δt: 2 minutes
Iteration: 60, time: 1.839 hours, Δt: 1.033 minutes
Iteration: 70, time: 1.960 hours, Δt: 43.618 seconds
Iteration: 80, time: 2.093 hours, Δt: 47.980 seconds
Iteration: 90, time: 2.240 hours, Δt: 52.778 seconds
Iteration: 100, time: 2.367 hours, Δt: 45.612 seconds
Iteration: 110, time: 2.498 hours, Δt: 47.404 seconds
Iteration: 120, time: 2.621 hours, Δt: 44.281 seconds
Iteration: 130, time: 2.732 hours, Δt: 39.784 seconds
Iteration: 140, time: 2.838 hours, Δt: 38.093 seconds
Iteration: 150, time: 2.935 hours, Δt: 35.040 seconds
Iteration: 160, time: 3.025 hours, Δt: 32.549 seconds
Iteration: 170, time: 3.113 hours, Δt: 31.538 seconds
Iteration: 180, time: 3.203 hours, Δt: 32.310 seconds
Iteration: 190, time: 3.289 hours, Δt: 31.009 seconds
Iteration: 200, time: 3.363 hours, Δt: 26.657 seconds
Iteration: 210, time: 3.433 hours, Δt: 25.080 seconds
Iteration: 220, time: 3.505 hours, Δt: 26.047 seconds
Iteration: 230, time: 3.584 hours, Δt: 28.652 seconds
Iteration: 240, time: 3.667 hours, Δt: 29.775 seconds
Iteration: 250, time: 3.749 hours, Δt: 29.358 seconds
Iteration: 260, time: 3.830 hours, Δt: 29.092 seconds
Iteration: 270, time: 3.917 hours, Δt: 31.450 seconds
Iteration: 280, time: 4.008 hours, Δt: 32.880 seconds
[ Info: Simulation is stopping. Model time 4.008 hours has hit or exceeded simulation stop time 4 hours.

Visualizing the solution

We'd like to a make a plankton movie. First we load the output file,

using JLD2

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))
25-element Array{Int64,1}:
   0
   5
  10
  15
  20
  25
  30
  35
  40
  45
   ⋮
 125
 140
 158
 176
 197
 220
 240
 261
 280

Next we construct the $x, z$ grid for plotting purposes,

using Oceananigans.Grids

xw, yw, zw = nodes(model.velocities.w)
xp, yp, zp = nodes(model.tracers.plankton)

Finally, we animate the convective plumes and plankton swirls,

using Plots

@info "Making a movie about plankton..."

anim = @animate for (i, iteration) in enumerate(iterations)

    @info "Plotting frame $i from iteration $iteration..."

    t = file["timeseries/t/$iteration"]
    w = file["timeseries/w/$iteration"][:, 1, :]
    p = file["timeseries/plankton/$iteration"][:, 1, :]
    P = file["timeseries/averaged_plankton/$iteration"][1, 1, :]

    w_max = maximum(abs, w) + 1e-9
    w_lim = 0.8 * w_max

    p_min = minimum(p) - 1e-9
    p_max = maximum(p) + 1e-9
    p_lim = 1

    w_levels = vcat([-w_max], range(-w_lim, stop=w_lim, length=21), [w_max])
    p_levels = collect(range(p_min, stop=p_lim, length=20))
    p_max > p_lim && push!(p_levels, p_max)

    kwargs = (xlabel="x", ylabel="y", aspectratio=1, linewidth=0, colorbar=true,
              xlims=(0, model.grid.Lx), ylims=(-model.grid.Lz, 0))

    w_plot = contourf(xw, zw, w';
                       color = :balance,
                      levels = w_levels,
                       clims = (-w_lim, w_lim),
                      kwargs...)

    p_plot = contourf(xp, zp, p';
                       color = :matter,
                      levels = p_levels,
                       clims = (p_min, p_lim),
                      kwargs...)

    P_plot = plot(P, zp[:],
                  linewidth = 2,
                  label = nothing,
                  xlims = (-0.2, 1),
                  ylabel = "Averaged plankton",
                  xlabel = "Plankton concentration")

    plot(w_plot, p_plot, P_plot,
         title=["Vertical velocity" "Plankton" "Averaged plankton"],
         link = :y,
         layout=(1, 3), size=(1700, 400))
end

This page was generated using Literate.jl.