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: 3.260 s
progress: 0.5 %, i: 0200, t: 3.850 min, Δt: 1.210 s, wall time: 552.928 ms
progress: 0.8 %, i: 0300, t: 6.068 min, Δt: 1.331 s, wall time: 545.391 ms
progress: 1.2 %, i: 0400, t: 8.508 min, Δt: 1.464 s, wall time: 562.755 ms
progress: 1.6 %, i: 0500, t: 11.193 min, Δt: 1.611 s, wall time: 568.661 ms
progress: 2.0 %, i: 0600, t: 14.145 min, Δt: 1.772 s, wall time: 569.234 ms
progress: 2.4 %, i: 0700, t: 17.393 min, Δt: 1.949 s, wall time: 547.880 ms
progress: 2.9 %, i: 0800, t: 20.966 min, Δt: 2.144 s, wall time: 592.567 ms
progress: 3.5 %, i: 0900, t: 24.896 min, Δt: 2.358 s, wall time: 634.789 ms
progress: 4.1 %, i: 1000, t: 29.219 min, Δt: 2.594 s, wall time: 547.103 ms
progress: 4.7 %, i: 1100, t: 33.974 min, Δt: 2.853 s, wall time: 553.472 ms
progress: 5.4 %, i: 1200, t: 39.205 min, Δt: 3.138 s, wall time: 555.010 ms
progress: 6.2 %, i: 1300, t: 44.958 min, Δt: 3.452 s, wall time: 544.561 ms
progress: 7.1 %, i: 1400, t: 51.287 min, Δt: 3.797 s, wall time: 551.862 ms
progress: 8.1 %, i: 1500, t: 58.250 min, Δt: 4.177 s, wall time: 608.767 ms
progress: 9.2 %, i: 1600, t: 1.098 hr, Δt: 4.595 s, wall time: 635.164 ms
progress: 10.3 %, i: 1700, t: 1.239 hr, Δt: 5.054 s, wall time: 647.083 ms
progress: 11.6 %, i: 1800, t: 1.393 hr, Δt: 5.560 s, wall time: 602.473 ms
progress: 13.0 %, i: 1900, t: 1.563 hr, Δt: 6.116 s, wall time: 551.847 ms
progress: 14.6 %, i: 2000, t: 1.750 hr, Δt: 6.727 s, wall time: 592.436 ms
progress: 16.2 %, i: 2100, t: 1.948 hr, Δt: 7.142 s, wall time: 619.121 ms
progress: 17.2 %, i: 2200, t: 2.061 hr, Δt: 4.053 s, wall time: 612.942 ms
progress: 18.2 %, i: 2300, t: 2.185 hr, Δt: 4.458 s, wall time: 601.268 ms
progress: 19.3 %, i: 2400, t: 2.321 hr, Δt: 4.904 s, wall time: 549.283 ms
progress: 20.6 %, i: 2500, t: 2.471 hr, Δt: 5.394 s, wall time: 542.492 ms
progress: 22.0 %, i: 2600, t: 2.636 hr, Δt: 5.934 s, wall time: 545.569 ms
progress: 23.5 %, i: 2700, t: 2.817 hr, Δt: 6.527 s, wall time: 544.103 ms
progress: 25.0 %, i: 2800, t: 2.996 hr, Δt: 6.455 s, wall time: 542.619 ms
progress: 26.5 %, i: 2900, t: 3.179 hr, Δt: 6.562 s, wall time: 533.468 ms
progress: 27.6 %, i: 3000, t: 3.316 hr, Δt: 4.931 s, wall time: 545.617 ms
progress: 28.9 %, i: 3100, t: 3.466 hr, Δt: 5.424 s, wall time: 536.661 ms
progress: 30.2 %, i: 3200, t: 3.630 hr, Δt: 5.886 s, wall time: 538.802 ms
progress: 31.7 %, i: 3300, t: 3.799 hr, Δt: 6.093 s, wall time: 541.059 ms
progress: 32.9 %, i: 3400, t: 3.945 hr, Δt: 5.255 s, wall time: 539.929 ms
progress: 34.2 %, i: 3500, t: 4.104 hr, Δt: 5.724 s, wall time: 548.488 ms
progress: 35.2 %, i: 3600, t: 4.221 hr, Δt: 4.216 s, wall time: 555.840 ms
progress: 36.2 %, i: 3700, t: 4.350 hr, Δt: 4.638 s, wall time: 550.023 ms
progress: 37.3 %, i: 3800, t: 4.477 hr, Δt: 4.587 s, wall time: 552.256 ms
progress: 38.5 %, i: 3900, t: 4.617 hr, Δt: 5.045 s, wall time: 556.991 ms
progress: 39.8 %, i: 4000, t: 4.772 hr, Δt: 5.550 s, wall time: 577.960 ms
progress: 40.9 %, i: 4100, t: 4.904 hr, Δt: 4.773 s, wall time: 566.014 ms
progress: 42.0 %, i: 4200, t: 5.046 hr, Δt: 5.092 s, wall time: 577.381 ms
progress: 43.2 %, i: 4300, t: 5.185 hr, Δt: 5.008 s, wall time: 621.970 ms
progress: 44.2 %, i: 4400, t: 5.307 hr, Δt: 4.397 s, wall time: 699.774 ms
progress: 45.3 %, i: 4500, t: 5.437 hr, Δt: 4.688 s, wall time: 692.543 ms
progress: 46.4 %, i: 4600, t: 5.565 hr, Δt: 4.588 s, wall time: 718.147 ms
progress: 47.4 %, i: 4700, t: 5.692 hr, Δt: 4.583 s, wall time: 703.314 ms
progress: 48.4 %, i: 4800, t: 5.813 hr, Δt: 4.368 s, wall time: 601.112 ms
progress: 49.6 %, i: 4900, t: 5.947 hr, Δt: 4.805 s, wall time: 586.110 ms
progress: 50.7 %, i: 5000, t: 6.084 hr, Δt: 4.949 s, wall time: 574.507 ms
progress: 51.7 %, i: 5100, t: 6.209 hr, Δt: 4.500 s, wall time: 567.387 ms
progress: 52.8 %, i: 5200, t: 6.331 hr, Δt: 4.376 s, wall time: 563.935 ms
progress: 53.6 %, i: 5300, t: 6.436 hr, Δt: 3.795 s, wall time: 558.283 ms
progress: 54.5 %, i: 5400, t: 6.543 hr, Δt: 3.830 s, wall time: 569.132 ms
progress: 55.5 %, i: 5500, t: 6.660 hr, Δt: 4.213 s, wall time: 565.763 ms
progress: 56.6 %, i: 5600, t: 6.787 hr, Δt: 4.581 s, wall time: 559.894 ms
progress: 57.6 %, i: 5700, t: 6.916 hr, Δt: 4.665 s, wall time: 558.962 ms
progress: 58.6 %, i: 5800, t: 7.027 hr, Δt: 3.995 s, wall time: 566.860 ms
progress: 59.3 %, i: 5900, t: 7.118 hr, Δt: 3.258 s, wall time: 563.388 ms
progress: 60.1 %, i: 6000, t: 7.217 hr, Δt: 3.584 s, wall time: 559.407 ms
progress: 61.0 %, i: 6100, t: 7.323 hr, Δt: 3.798 s, wall time: 578.016 ms
progress: 61.8 %, i: 6200, t: 7.419 hr, Δt: 3.461 s, wall time: 563.982 ms
progress: 62.7 %, i: 6300, t: 7.525 hr, Δt: 3.807 s, wall time: 563.491 ms
progress: 63.5 %, i: 6400, t: 7.616 hr, Δt: 3.266 s, wall time: 561.257 ms
progress: 64.3 %, i: 6500, t: 7.715 hr, Δt: 3.592 s, wall time: 543.480 ms
progress: 65.2 %, i: 6600, t: 7.825 hr, Δt: 3.952 s, wall time: 552.605 ms
progress: 66.2 %, i: 6700, t: 7.946 hr, Δt: 4.347 s, wall time: 557.295 ms
progress: 67.1 %, i: 6800, t: 8.049 hr, Δt: 3.727 s, wall time: 562.500 ms
progress: 68.0 %, i: 6900, t: 8.163 hr, Δt: 4.099 s, wall time: 565.996 ms
progress: 68.9 %, i: 7000, t: 8.268 hr, Δt: 3.770 s, wall time: 551.094 ms
progress: 69.7 %, i: 7100, t: 8.367 hr, Δt: 3.570 s, wall time: 542.261 ms
progress: 70.6 %, i: 7200, t: 8.469 hr, Δt: 3.682 s, wall time: 586.601 ms
progress: 71.5 %, i: 7300, t: 8.575 hr, Δt: 3.819 s, wall time: 685.178 ms
progress: 72.4 %, i: 7400, t: 8.686 hr, Δt: 3.966 s, wall time: 636.003 ms
progress: 73.3 %, i: 7500, t: 8.800 hr, Δt: 4.102 s, wall time: 618.849 ms
progress: 74.1 %, i: 7600, t: 8.898 hr, Δt: 3.540 s, wall time: 612.424 ms
progress: 74.9 %, i: 7700, t: 8.989 hr, Δt: 3.283 s, wall time: 560.381 ms
progress: 75.7 %, i: 7800, t: 9.081 hr, Δt: 3.289 s, wall time: 550.115 ms
progress: 76.5 %, i: 7900, t: 9.181 hr, Δt: 3.618 s, wall time: 597.361 ms
progress: 77.4 %, i: 8000, t: 9.292 hr, Δt: 3.980 s, wall time: 617.865 ms
progress: 78.2 %, i: 8100, t: 9.386 hr, Δt: 3.390 s, wall time: 608.874 ms
progress: 79.1 %, i: 8200, t: 9.489 hr, Δt: 3.729 s, wall time: 609.668 ms
progress: 79.9 %, i: 8300, t: 9.593 hr, Δt: 3.728 s, wall time: 632.496 ms
progress: 80.8 %, i: 8400, t: 9.692 hr, Δt: 3.561 s, wall time: 654.600 ms
progress: 81.4 %, i: 8500, t: 9.774 hr, Δt: 2.949 s, wall time: 600.997 ms
progress: 82.1 %, i: 8600, t: 9.853 hr, Δt: 2.851 s, wall time: 559.754 ms
progress: 82.8 %, i: 8700, t: 9.940 hr, Δt: 3.136 s, wall time: 546.435 ms
progress: 83.6 %, i: 8800, t: 10.036 hr, Δt: 3.450 s, wall time: 554.814 ms
progress: 84.5 %, i: 8900, t: 10.141 hr, Δt: 3.795 s, wall time: 552.486 ms
progress: 85.3 %, i: 9000, t: 10.231 hr, Δt: 3.250 s, wall time: 555.312 ms
progress: 86.0 %, i: 9100, t: 10.316 hr, Δt: 3.031 s, wall time: 549.306 ms
progress: 86.6 %, i: 9200, t: 10.394 hr, Δt: 2.832 s, wall time: 553.112 ms
progress: 87.3 %, i: 9300, t: 10.481 hr, Δt: 3.116 s, wall time: 549.075 ms
progress: 88.1 %, i: 9400, t: 10.571 hr, Δt: 3.226 s, wall time: 551.217 ms
progress: 88.9 %, i: 9500, t: 10.666 hr, Δt: 3.453 s, wall time: 544.440 ms
progress: 89.8 %, i: 9600, t: 10.772 hr, Δt: 3.798 s, wall time: 549.004 ms
progress: 90.7 %, i: 9700, t: 10.886 hr, Δt: 4.102 s, wall time: 556.782 ms
progress: 91.5 %, i: 9800, t: 10.981 hr, Δt: 3.429 s, wall time: 622.234 ms
progress: 92.4 %, i: 9900, t: 11.086 hr, Δt: 3.772 s, wall time: 636.248 ms
progress: 93.3 %, i: 10000, t: 11.193 hr, Δt: 3.853 s, wall time: 565.544 ms
progress: 94.2 %, i: 10100, t: 11.308 hr, Δt: 4.149 s, wall time: 562.664 ms
progress: 95.1 %, i: 10200, t: 11.410 hr, Δt: 3.656 s, wall time: 549.880 ms
progress: 96.0 %, i: 10300, t: 11.521 hr, Δt: 3.999 s, wall time: 551.828 ms
progress: 97.0 %, i: 10400, t: 11.635 hr, Δt: 4.118 s, wall time: 555.138 ms
progress: 97.8 %, i: 10500, t: 11.734 hr, Δt: 3.538 s, wall time: 555.662 ms
progress: 98.7 %, i: 10600, t: 11.842 hr, Δt: 3.892 s, wall time: 555.079 ms
progress: 99.6 %, i: 10700, t: 11.950 hr, Δt: 3.885 s, wall time: 550.222 ms
progress: 100.4 %, i: 10800, t: 12.051 hr, Δt: 3.666 s, wall time: 579.058 ms

Plot the result at the end

makeplot!(axs, model)
gcf()

This page was generated using Literate.jl.