Edit on GitHub

Ocean convection example

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

  • to set boundary conditions;
  • to defined and insert a user-defined forcing function into a simulation.
  • to use the TimeStepWizard to manage and adapt the simulation time-step.

To begin, we load Oceananigans, a plotting package, and a few miscellaneous useful packages.

using Oceananigans, PyPlot, Random, Printf

Parameters

We choose a modest two-dimensional resolution of 128² in a 64² m² domain , implying a resolution of 0.5 m. Our fluid is initially stratified with a squared buoyancy frequency

$ N² = 10⁻⁵ \rm{s⁻²} $

and a surface buoyancy flux

$ Q_b = 10⁻⁸ \rm{m³ s⁻²} $

Because we use the physics-based convection whereby buoyancy flux by a positive vertical velocity implies positive flux, a positive buoyancy flux at the top of the domain carries buoyancy out of the fluid and causes convection. Finally, we end the simulation after 1 day.

Nz = 128
Lz = 64.0
N² = 1e-5
Qb = 1e-8
end_time = day / 2
43200.0

Creating boundary conditions

Create boundary conditions. Note that temperature is buoyancy in our problem.

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

Define a forcing function

Our forcing function roughly corresponds to the growth of phytoplankton in light (with a penetration depth of 16 meters here), and death due to natural mortality at a rate of 1 phytoplankton unit per second.

growth_and_decay = SimpleForcing((x, y, z, t) -> exp(z/16) - 1)

# Instantiate the model
model = Model(
                   grid = RegularCartesianGrid(size = (Nz, 1, Nz), length = (Lz, Lz, Lz)),
                closure = ConstantIsotropicDiffusivity(ν=1e-4, κ=1e-4),
               coriolis = FPlane(f=1e-4),
                tracers = (:b, :plankton),
               buoyancy = BuoyancyTracer(),
                forcing = ModelForcing(plankton=growth_and_decay),
    boundary_conditions = BoundaryConditions(b=buoyancy_bcs)
)
Oceananigans.Model on a CPU architecture (time = 0.000 ns, iteration = 0) 
├── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
├── tracers: (:b, :plankton)
├── closure: ConstantIsotropicDiffusivity{Float64,NamedTuple{(:b, :plankton),Tuple{Float64,Float64}}}
├── buoyancy: BuoyancyTracer
├── coriolis: FPlane{Float64}
├── output writers: OrderedCollections.OrderedDict{Symbol,Oceananigans.AbstractOutputWriter} with no entries
└── diagnostics: OrderedCollections.OrderedDict{Symbol,Oceananigans.AbstractDiagnostic} with no entries

Set makeplot = true to live-update a plot of vertical velocity, temperature, and salinity as the simulation runs.

makeplot = false

# Set initial condition. Initial velocity and salinity fluctuations needed for AMD.
Ξ(z) = randn() * z / Lz * (1 + z / Lz) # noise
b₀(x, y, z) = N² * z + N² * Lz * 1e-6 * Ξ(z)
set!(model, b=b₀)

# A wizard for managing the simulation time-step.
wizard = TimeStepWizard(cfl=0.1, Δt=1.0, max_change=1.1, max_Δt=90.0)

# Create a plot
fig, axs = subplots(ncols=2, figsize=(10, 6))

"""
    makeplot!(axs, model)

Make side-by-side x-z slices of vertical velocity and plankton associated with `model` in `axs`.
"""
function makeplot!(axs, model)
    xC = repeat(model.grid.xC, 1, model.grid.Nz)
    zF = repeat(reshape(model.grid.zF[1:end-1], 1, model.grid.Nz), model.grid.Nx, 1)
    zC = repeat(reshape(model.grid.zC, 1, model.grid.Nz), model.grid.Nx, 1)

    sca(axs[1]); cla()
    # Calling the Array() constructor here allows plots to be made for GPU models:
    pcolormesh(xC, zF, Array(interior(model.velocities.w))[:, 1, :])
    title("Vertical velocity")
    xlabel("\$ x \$ (m)")
    ylabel("\$ z \$ (m)")

    sca(axs[2]); cla()
    # Calling the Array() constructor here allows plots to be made for GPU models:
    pcolormesh(xC, zC, Array(interior(model.tracers.plankton))[:, 1, :])
    title("Phytoplankton concentration")
    xlabel("\$ x \$ (m)")
    axs[2].tick_params(left=false, labelleft=false)

    suptitle(@sprintf("\$ t = %.2f\$ hours", model.clock.time / hour))
    [ax.set_aspect(1) for ax in axs]
    pause(0.01)
    gcf()
    return nothing
end
Main.ex-ocean_convection_with_plankton.makeplot!

Run the model:

while model.clock.time < end_time
    update_Δt!(wizard, model)
    walltime = @elapsed time_step!(model, 100, wizard.Δt)

    # Print a progress message
    @printf("progress: %.1f %%, i: %04d, t: %s, Δt: %s, wall time: %s\n",
            model.clock.time / end_time * 100, model.clock.iteration,
            prettytime(model.clock.time), prettytime(wizard.Δt), prettytime(walltime))

    makeplot && makeplot!(axs, model)
end
progress: 0.3 %, i: 0100, t: 1.833 min, Δt: 1.100 s, wall time: 2.916 s
progress: 0.5 %, i: 0200, t: 3.850 min, Δt: 1.210 s, wall time: 574.060 ms
progress: 0.8 %, i: 0300, t: 6.068 min, Δt: 1.331 s, wall time: 566.691 ms
progress: 1.2 %, i: 0400, t: 8.508 min, Δt: 1.464 s, wall time: 578.322 ms
progress: 1.6 %, i: 0500, t: 11.193 min, Δt: 1.611 s, wall time: 571.777 ms
progress: 2.0 %, i: 0600, t: 14.145 min, Δt: 1.772 s, wall time: 581.540 ms
progress: 2.4 %, i: 0700, t: 17.393 min, Δt: 1.949 s, wall time: 568.769 ms
progress: 2.9 %, i: 0800, t: 20.966 min, Δt: 2.144 s, wall time: 572.755 ms
progress: 3.5 %, i: 0900, t: 24.896 min, Δt: 2.358 s, wall time: 570.563 ms
progress: 4.1 %, i: 1000, t: 29.219 min, Δt: 2.594 s, wall time: 572.504 ms
progress: 4.7 %, i: 1100, t: 33.974 min, Δt: 2.853 s, wall time: 572.720 ms
progress: 5.4 %, i: 1200, t: 39.205 min, Δt: 3.138 s, wall time: 576.781 ms
progress: 6.2 %, i: 1300, t: 44.958 min, Δt: 3.452 s, wall time: 572.665 ms
progress: 7.1 %, i: 1400, t: 51.287 min, Δt: 3.797 s, wall time: 581.402 ms
progress: 8.1 %, i: 1500, t: 58.250 min, Δt: 4.177 s, wall time: 575.389 ms
progress: 9.2 %, i: 1600, t: 1.098 hr, Δt: 4.595 s, wall time: 573.795 ms
progress: 10.3 %, i: 1700, t: 1.239 hr, Δt: 5.054 s, wall time: 576.901 ms
progress: 11.6 %, i: 1800, t: 1.393 hr, Δt: 5.560 s, wall time: 575.366 ms
progress: 13.0 %, i: 1900, t: 1.563 hr, Δt: 6.116 s, wall time: 572.539 ms
progress: 14.6 %, i: 2000, t: 1.750 hr, Δt: 6.727 s, wall time: 577.592 ms
progress: 16.1 %, i: 2100, t: 1.933 hr, Δt: 6.584 s, wall time: 572.369 ms
progress: 17.1 %, i: 2200, t: 2.046 hr, Δt: 4.083 s, wall time: 574.423 ms
progress: 18.1 %, i: 2300, t: 2.171 hr, Δt: 4.491 s, wall time: 582.559 ms
progress: 19.2 %, i: 2400, t: 2.308 hr, Δt: 4.940 s, wall time: 580.212 ms
progress: 20.5 %, i: 2500, t: 2.459 hr, Δt: 5.434 s, wall time: 572.975 ms
progress: 21.9 %, i: 2600, t: 2.625 hr, Δt: 5.978 s, wall time: 569.404 ms
progress: 23.3 %, i: 2700, t: 2.795 hr, Δt: 6.102 s, wall time: 581.090 ms
progress: 24.8 %, i: 2800, t: 2.981 hr, Δt: 6.712 s, wall time: 569.299 ms
progress: 26.5 %, i: 2900, t: 3.174 hr, Δt: 6.939 s, wall time: 573.664 ms
progress: 27.9 %, i: 3000, t: 3.346 hr, Δt: 6.199 s, wall time: 574.701 ms
progress: 29.3 %, i: 3100, t: 3.520 hr, Δt: 6.242 s, wall time: 579.462 ms
progress: 30.5 %, i: 3200, t: 3.660 hr, Δt: 5.062 s, wall time: 587.369 ms
progress: 31.7 %, i: 3300, t: 3.809 hr, Δt: 5.342 s, wall time: 575.019 ms
progress: 32.9 %, i: 3400, t: 3.947 hr, Δt: 4.993 s, wall time: 570.140 ms
progress: 33.8 %, i: 3500, t: 4.054 hr, Δt: 3.835 s, wall time: 568.006 ms
progress: 34.8 %, i: 3600, t: 4.171 hr, Δt: 4.219 s, wall time: 574.759 ms
progress: 35.8 %, i: 3700, t: 4.300 hr, Δt: 4.641 s, wall time: 578.457 ms
progress: 37.0 %, i: 3800, t: 4.442 hr, Δt: 5.105 s, wall time: 575.476 ms
progress: 38.1 %, i: 3900, t: 4.573 hr, Δt: 4.722 s, wall time: 576.498 ms
progress: 39.3 %, i: 4000, t: 4.717 hr, Δt: 5.194 s, wall time: 579.088 ms
progress: 40.3 %, i: 4100, t: 4.839 hr, Δt: 4.396 s, wall time: 573.907 ms
progress: 41.4 %, i: 4200, t: 4.974 hr, Δt: 4.836 s, wall time: 569.840 ms
progress: 42.5 %, i: 4300, t: 5.101 hr, Δt: 4.571 s, wall time: 571.921 ms
progress: 43.5 %, i: 4400, t: 5.219 hr, Δt: 4.253 s, wall time: 571.170 ms
progress: 44.5 %, i: 4500, t: 5.343 hr, Δt: 4.477 s, wall time: 581.387 ms
progress: 45.7 %, i: 4600, t: 5.480 hr, Δt: 4.925 s, wall time: 589.789 ms
progress: 46.9 %, i: 4700, t: 5.627 hr, Δt: 5.297 s, wall time: 574.060 ms
progress: 47.9 %, i: 4800, t: 5.752 hr, Δt: 4.510 s, wall time: 576.693 ms
progress: 49.0 %, i: 4900, t: 5.879 hr, Δt: 4.552 s, wall time: 575.618 ms
progress: 50.0 %, i: 5000, t: 5.998 hr, Δt: 4.275 s, wall time: 573.157 ms
progress: 51.1 %, i: 5100, t: 6.128 hr, Δt: 4.702 s, wall time: 569.660 ms
progress: 52.1 %, i: 5200, t: 6.258 hr, Δt: 4.672 s, wall time: 581.782 ms
progress: 53.3 %, i: 5300, t: 6.391 hr, Δt: 4.779 s, wall time: 571.247 ms
progress: 54.5 %, i: 5400, t: 6.537 hr, Δt: 5.256 s, wall time: 575.198 ms
progress: 55.5 %, i: 5500, t: 6.657 hr, Δt: 4.318 s, wall time: 573.094 ms
progress: 56.6 %, i: 5600, t: 6.789 hr, Δt: 4.750 s, wall time: 573.495 ms
progress: 57.6 %, i: 5700, t: 6.909 hr, Δt: 4.337 s, wall time: 577.922 ms
progress: 58.5 %, i: 5800, t: 7.026 hr, Δt: 4.194 s, wall time: 578.461 ms
progress: 59.6 %, i: 5900, t: 7.152 hr, Δt: 4.548 s, wall time: 569.852 ms
progress: 60.7 %, i: 6000, t: 7.282 hr, Δt: 4.673 s, wall time: 560.567 ms
progress: 61.6 %, i: 6100, t: 7.396 hr, Δt: 4.104 s, wall time: 566.741 ms
progress: 62.4 %, i: 6200, t: 7.493 hr, Δt: 3.513 s, wall time: 562.104 ms
progress: 63.3 %, i: 6300, t: 7.601 hr, Δt: 3.865 s, wall time: 568.637 ms
progress: 64.1 %, i: 6400, t: 7.693 hr, Δt: 3.336 s, wall time: 571.361 ms
progress: 64.9 %, i: 6500, t: 7.783 hr, Δt: 3.236 s, wall time: 564.290 ms
progress: 65.7 %, i: 6600, t: 7.882 hr, Δt: 3.560 s, wall time: 574.634 ms
progress: 66.6 %, i: 6700, t: 7.986 hr, Δt: 3.743 s, wall time: 575.088 ms
progress: 67.4 %, i: 6800, t: 8.092 hr, Δt: 3.822 s, wall time: 575.247 ms
progress: 68.3 %, i: 6900, t: 8.201 hr, Δt: 3.933 s, wall time: 581.216 ms
progress: 69.2 %, i: 7000, t: 8.306 hr, Δt: 3.762 s, wall time: 570.387 ms
progress: 70.0 %, i: 7100, t: 8.403 hr, Δt: 3.496 s, wall time: 573.519 ms
progress: 70.9 %, i: 7200, t: 8.502 hr, Δt: 3.576 s, wall time: 573.012 ms
progress: 71.6 %, i: 7300, t: 8.595 hr, Δt: 3.350 s, wall time: 576.716 ms
progress: 72.5 %, i: 7400, t: 8.698 hr, Δt: 3.685 s, wall time: 582.030 ms
progress: 73.4 %, i: 7500, t: 8.805 hr, Δt: 3.869 s, wall time: 572.133 ms
progress: 74.2 %, i: 7600, t: 8.900 hr, Δt: 3.412 s, wall time: 579.174 ms
progress: 75.0 %, i: 7700, t: 9.004 hr, Δt: 3.753 s, wall time: 577.830 ms
progress: 75.9 %, i: 7800, t: 9.111 hr, Δt: 3.842 s, wall time: 575.939 ms
progress: 76.8 %, i: 7900, t: 9.213 hr, Δt: 3.685 s, wall time: 574.961 ms
progress: 77.5 %, i: 8000, t: 9.296 hr, Δt: 2.961 s, wall time: 579.632 ms
progress: 78.2 %, i: 8100, t: 9.383 hr, Δt: 3.134 s, wall time: 583.811 ms
progress: 79.0 %, i: 8200, t: 9.477 hr, Δt: 3.411 s, wall time: 580.744 ms
progress: 79.8 %, i: 8300, t: 9.582 hr, Δt: 3.752 s, wall time: 578.203 ms
progress: 80.7 %, i: 8400, t: 9.683 hr, Δt: 3.636 s, wall time: 579.478 ms
progress: 81.5 %, i: 8500, t: 9.779 hr, Δt: 3.472 s, wall time: 569.854 ms
progress: 82.4 %, i: 8600, t: 9.885 hr, Δt: 3.806 s, wall time: 566.402 ms
progress: 83.3 %, i: 8700, t: 9.997 hr, Δt: 4.033 s, wall time: 570.211 ms
progress: 84.1 %, i: 8800, t: 10.095 hr, Δt: 3.524 s, wall time: 572.459 ms
progress: 84.9 %, i: 8900, t: 10.191 hr, Δt: 3.479 s, wall time: 568.711 ms
progress: 85.7 %, i: 9000, t: 10.282 hr, Δt: 3.259 s, wall time: 573.348 ms
progress: 86.5 %, i: 9100, t: 10.382 hr, Δt: 3.584 s, wall time: 571.701 ms
progress: 87.4 %, i: 9200, t: 10.487 hr, Δt: 3.803 s, wall time: 567.029 ms
progress: 88.3 %, i: 9300, t: 10.592 hr, Δt: 3.778 s, wall time: 569.000 ms
progress: 89.2 %, i: 9400, t: 10.708 hr, Δt: 4.156 s, wall time: 570.381 ms
progress: 90.1 %, i: 9500, t: 10.814 hr, Δt: 3.825 s, wall time: 576.333 ms
progress: 90.9 %, i: 9600, t: 10.905 hr, Δt: 3.287 s, wall time: 565.002 ms
progress: 91.7 %, i: 9700, t: 11.000 hr, Δt: 3.423 s, wall time: 569.213 ms
progress: 92.4 %, i: 9800, t: 11.085 hr, Δt: 3.049 s, wall time: 587.702 ms
progress: 93.0 %, i: 9900, t: 11.166 hr, Δt: 2.907 s, wall time: 578.383 ms
progress: 93.8 %, i: 10000, t: 11.254 hr, Δt: 3.197 s, wall time: 570.464 ms
progress: 94.6 %, i: 10100, t: 11.352 hr, Δt: 3.517 s, wall time: 567.707 ms
progress: 95.5 %, i: 10200, t: 11.460 hr, Δt: 3.869 s, wall time: 568.092 ms
progress: 96.4 %, i: 10300, t: 11.569 hr, Δt: 3.925 s, wall time: 565.681 ms
progress: 97.2 %, i: 10400, t: 11.667 hr, Δt: 3.536 s, wall time: 574.957 ms
progress: 98.0 %, i: 10500, t: 11.756 hr, Δt: 3.219 s, wall time: 565.267 ms
progress: 98.7 %, i: 10600, t: 11.849 hr, Δt: 3.353 s, wall time: 570.275 ms
progress: 99.6 %, i: 10700, t: 11.952 hr, Δt: 3.688 s, wall time: 569.564 ms
progress: 100.3 %, i: 10800, t: 12.038 hr, Δt: 3.113 s, wall time: 576.939 ms

Plot the result at the end

makeplot!(axs, model)
gcf()

This page was generated using Literate.jl.