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.120 s
progress: 0.5 %, i: 0200, t: 3.850 min, Δt: 1.210 s, wall time: 545.484 ms
progress: 0.8 %, i: 0300, t: 6.068 min, Δt: 1.331 s, wall time: 542.975 ms
progress: 1.2 %, i: 0400, t: 8.508 min, Δt: 1.464 s, wall time: 534.907 ms
progress: 1.6 %, i: 0500, t: 11.193 min, Δt: 1.611 s, wall time: 536.352 ms
progress: 2.0 %, i: 0600, t: 14.145 min, Δt: 1.772 s, wall time: 526.220 ms
progress: 2.4 %, i: 0700, t: 17.393 min, Δt: 1.949 s, wall time: 526.593 ms
progress: 2.9 %, i: 0800, t: 20.966 min, Δt: 2.144 s, wall time: 539.509 ms
progress: 3.5 %, i: 0900, t: 24.896 min, Δt: 2.358 s, wall time: 573.279 ms
progress: 4.1 %, i: 1000, t: 29.219 min, Δt: 2.594 s, wall time: 595.527 ms
progress: 4.7 %, i: 1100, t: 33.974 min, Δt: 2.853 s, wall time: 533.274 ms
progress: 5.4 %, i: 1200, t: 39.205 min, Δt: 3.138 s, wall time: 527.293 ms
progress: 6.2 %, i: 1300, t: 44.958 min, Δt: 3.452 s, wall time: 523.358 ms
progress: 7.1 %, i: 1400, t: 51.287 min, Δt: 3.797 s, wall time: 557.421 ms
progress: 8.1 %, i: 1500, t: 58.250 min, Δt: 4.177 s, wall time: 596.004 ms
progress: 9.2 %, i: 1600, t: 1.098 hr, Δt: 4.595 s, wall time: 592.992 ms
progress: 10.3 %, i: 1700, t: 1.239 hr, Δt: 5.054 s, wall time: 583.177 ms
progress: 11.6 %, i: 1800, t: 1.393 hr, Δt: 5.560 s, wall time: 583.240 ms
progress: 13.0 %, i: 1900, t: 1.563 hr, Δt: 6.116 s, wall time: 541.288 ms
progress: 14.6 %, i: 2000, t: 1.750 hr, Δt: 6.727 s, wall time: 544.386 ms
progress: 15.7 %, i: 2100, t: 1.889 hr, Δt: 5.015 s, wall time: 543.625 ms
progress: 16.6 %, i: 2200, t: 1.997 hr, Δt: 3.890 s, wall time: 530.329 ms
progress: 17.6 %, i: 2300, t: 2.116 hr, Δt: 4.279 s, wall time: 525.149 ms
progress: 18.7 %, i: 2400, t: 2.247 hr, Δt: 4.707 s, wall time: 529.681 ms
progress: 19.9 %, i: 2500, t: 2.391 hr, Δt: 5.178 s, wall time: 530.464 ms
progress: 21.2 %, i: 2600, t: 2.549 hr, Δt: 5.696 s, wall time: 514.038 ms
progress: 22.7 %, i: 2700, t: 2.723 hr, Δt: 6.265 s, wall time: 511.900 ms
progress: 24.2 %, i: 2800, t: 2.906 hr, Δt: 6.599 s, wall time: 511.407 ms
progress: 25.4 %, i: 2900, t: 3.044 hr, Δt: 4.953 s, wall time: 505.686 ms
progress: 26.6 %, i: 3000, t: 3.195 hr, Δt: 5.449 s, wall time: 512.650 ms
progress: 27.7 %, i: 3100, t: 3.326 hr, Δt: 4.684 s, wall time: 512.931 ms
progress: 28.9 %, i: 3200, t: 3.469 hr, Δt: 5.152 s, wall time: 515.950 ms
progress: 30.0 %, i: 3300, t: 3.603 hr, Δt: 4.849 s, wall time: 516.938 ms
progress: 31.3 %, i: 3400, t: 3.752 hr, Δt: 5.334 s, wall time: 517.621 ms
progress: 32.6 %, i: 3500, t: 3.915 hr, Δt: 5.868 s, wall time: 513.061 ms
progress: 34.0 %, i: 3600, t: 4.086 hr, Δt: 6.156 s, wall time: 521.218 ms
progress: 35.4 %, i: 3700, t: 4.250 hr, Δt: 5.923 s, wall time: 552.280 ms
progress: 36.7 %, i: 3800, t: 4.403 hr, Δt: 5.518 s, wall time: 566.293 ms
progress: 37.9 %, i: 3900, t: 4.543 hr, Δt: 5.045 s, wall time: 569.580 ms
progress: 39.1 %, i: 4000, t: 4.688 hr, Δt: 5.194 s, wall time: 568.831 ms
progress: 40.1 %, i: 4100, t: 4.815 hr, Δt: 4.590 s, wall time: 525.387 ms
progress: 41.1 %, i: 4200, t: 4.933 hr, Δt: 4.231 s, wall time: 515.364 ms
progress: 42.2 %, i: 4300, t: 5.062 hr, Δt: 4.654 s, wall time: 533.932 ms
progress: 43.3 %, i: 4400, t: 5.194 hr, Δt: 4.736 s, wall time: 516.040 ms
progress: 44.5 %, i: 4500, t: 5.338 hr, Δt: 5.210 s, wall time: 512.500 ms
progress: 45.6 %, i: 4600, t: 5.471 hr, Δt: 4.786 s, wall time: 514.787 ms
progress: 46.6 %, i: 4700, t: 5.596 hr, Δt: 4.494 s, wall time: 514.279 ms
progress: 47.6 %, i: 4800, t: 5.718 hr, Δt: 4.372 s, wall time: 527.701 ms
progress: 48.5 %, i: 4900, t: 5.816 hr, Δt: 3.541 s, wall time: 524.619 ms
progress: 49.4 %, i: 5000, t: 5.924 hr, Δt: 3.896 s, wall time: 516.408 ms
progress: 50.4 %, i: 5100, t: 6.043 hr, Δt: 4.285 s, wall time: 524.516 ms
progress: 51.3 %, i: 5200, t: 6.153 hr, Δt: 3.951 s, wall time: 521.735 ms
progress: 52.3 %, i: 5300, t: 6.274 hr, Δt: 4.346 s, wall time: 519.098 ms
progress: 53.2 %, i: 5400, t: 6.381 hr, Δt: 3.852 s, wall time: 521.408 ms
progress: 54.1 %, i: 5500, t: 6.497 hr, Δt: 4.198 s, wall time: 534.817 ms
progress: 55.0 %, i: 5600, t: 6.602 hr, Δt: 3.754 s, wall time: 567.335 ms
progress: 56.0 %, i: 5700, t: 6.716 hr, Δt: 4.129 s, wall time: 514.364 ms
progress: 57.0 %, i: 5800, t: 6.842 hr, Δt: 4.542 s, wall time: 517.410 ms
progress: 58.2 %, i: 5900, t: 6.981 hr, Δt: 4.996 s, wall time: 515.026 ms
progress: 59.2 %, i: 6000, t: 7.106 hr, Δt: 4.509 s, wall time: 515.519 ms
progress: 60.1 %, i: 6100, t: 7.210 hr, Δt: 3.742 s, wall time: 516.501 ms
progress: 61.0 %, i: 6200, t: 7.317 hr, Δt: 3.850 s, wall time: 505.407 ms
progress: 61.8 %, i: 6300, t: 7.415 hr, Δt: 3.507 s, wall time: 509.152 ms
progress: 62.6 %, i: 6400, t: 7.511 hr, Δt: 3.467 s, wall time: 521.234 ms
progress: 63.5 %, i: 6500, t: 7.617 hr, Δt: 3.814 s, wall time: 513.745 ms
progress: 64.4 %, i: 6600, t: 7.724 hr, Δt: 3.861 s, wall time: 515.975 ms
progress: 65.4 %, i: 6700, t: 7.842 hr, Δt: 4.247 s, wall time: 510.901 ms
progress: 66.4 %, i: 6800, t: 7.968 hr, Δt: 4.537 s, wall time: 507.539 ms
progress: 67.3 %, i: 6900, t: 8.076 hr, Δt: 3.887 s, wall time: 519.103 ms
progress: 68.2 %, i: 7000, t: 8.183 hr, Δt: 3.836 s, wall time: 558.657 ms
progress: 69.1 %, i: 7100, t: 8.293 hr, Δt: 3.959 s, wall time: 553.355 ms
progress: 70.0 %, i: 7200, t: 8.406 hr, Δt: 4.074 s, wall time: 516.462 ms
progress: 70.9 %, i: 7300, t: 8.506 hr, Δt: 3.616 s, wall time: 513.044 ms
progress: 71.8 %, i: 7400, t: 8.617 hr, Δt: 3.978 s, wall time: 516.386 ms
progress: 72.7 %, i: 7500, t: 8.723 hr, Δt: 3.833 s, wall time: 514.710 ms
progress: 73.5 %, i: 7600, t: 8.824 hr, Δt: 3.616 s, wall time: 513.269 ms
progress: 74.5 %, i: 7700, t: 8.934 hr, Δt: 3.978 s, wall time: 522.792 ms
progress: 75.5 %, i: 7800, t: 9.056 hr, Δt: 4.376 s, wall time: 528.163 ms
progress: 76.3 %, i: 7900, t: 9.154 hr, Δt: 3.544 s, wall time: 517.885 ms
progress: 77.2 %, i: 8000, t: 9.263 hr, Δt: 3.898 s, wall time: 511.082 ms
progress: 78.2 %, i: 8100, t: 9.381 hr, Δt: 4.250 s, wall time: 509.580 ms
progress: 79.0 %, i: 8200, t: 9.479 hr, Δt: 3.550 s, wall time: 518.818 ms
progress: 79.9 %, i: 8300, t: 9.585 hr, Δt: 3.808 s, wall time: 513.863 ms
progress: 80.6 %, i: 8400, t: 9.674 hr, Δt: 3.198 s, wall time: 508.373 ms
progress: 81.4 %, i: 8500, t: 9.772 hr, Δt: 3.517 s, wall time: 506.400 ms
progress: 82.3 %, i: 8600, t: 9.877 hr, Δt: 3.781 s, wall time: 504.785 ms
progress: 83.1 %, i: 8700, t: 9.973 hr, Δt: 3.460 s, wall time: 519.656 ms
progress: 83.8 %, i: 8800, t: 10.055 hr, Δt: 2.951 s, wall time: 517.048 ms
progress: 84.5 %, i: 8900, t: 10.145 hr, Δt: 3.246 s, wall time: 516.796 ms
progress: 85.3 %, i: 9000, t: 10.233 hr, Δt: 3.160 s, wall time: 506.565 ms
progress: 86.0 %, i: 9100, t: 10.317 hr, Δt: 3.041 s, wall time: 507.250 ms
progress: 86.7 %, i: 9200, t: 10.406 hr, Δt: 3.185 s, wall time: 509.145 ms
progress: 87.4 %, i: 9300, t: 10.491 hr, Δt: 3.062 s, wall time: 514.238 ms
progress: 88.2 %, i: 9400, t: 10.584 hr, Δt: 3.368 s, wall time: 506.276 ms
progress: 89.0 %, i: 9500, t: 10.675 hr, Δt: 3.284 s, wall time: 511.262 ms
progress: 89.8 %, i: 9600, t: 10.776 hr, Δt: 3.613 s, wall time: 505.916 ms
progress: 90.7 %, i: 9700, t: 10.886 hr, Δt: 3.974 s, wall time: 508.234 ms
progress: 91.5 %, i: 9800, t: 10.977 hr, Δt: 3.257 s, wall time: 514.095 ms
progress: 92.2 %, i: 9900, t: 11.068 hr, Δt: 3.284 s, wall time: 518.321 ms
progress: 93.1 %, i: 10000, t: 11.168 hr, Δt: 3.612 s, wall time: 510.630 ms
progress: 93.8 %, i: 10100, t: 11.260 hr, Δt: 3.304 s, wall time: 510.933 ms
progress: 94.6 %, i: 10200, t: 11.350 hr, Δt: 3.240 s, wall time: 521.292 ms
progress: 95.4 %, i: 10300, t: 11.449 hr, Δt: 3.564 s, wall time: 526.562 ms
progress: 96.2 %, i: 10400, t: 11.540 hr, Δt: 3.276 s, wall time: 519.948 ms
progress: 97.0 %, i: 10500, t: 11.640 hr, Δt: 3.598 s, wall time: 529.684 ms
progress: 97.7 %, i: 10600, t: 11.721 hr, Δt: 2.907 s, wall time: 521.156 ms
progress: 98.4 %, i: 10700, t: 11.807 hr, Δt: 3.102 s, wall time: 519.423 ms
progress: 99.2 %, i: 10800, t: 11.902 hr, Δt: 3.412 s, wall time: 507.918 ms
progress: 100.0 %, i: 10900, t: 11.995 hr, Δt: 3.372 s, wall time: 504.856 ms
progress: 100.7 %, i: 11000, t: 12.085 hr, Δt: 3.246 s, wall time: 514.133 ms

Plot the result at the end

makeplot!(axs, model)
gcf()

This page was generated using Literate.jl.