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.950 s
progress: 0.5 %, i: 0200, t: 3.850 min, Δt: 1.210 s, wall time: 657.197 ms
progress: 0.8 %, i: 0300, t: 6.068 min, Δt: 1.331 s, wall time: 646.143 ms
progress: 1.2 %, i: 0400, t: 8.508 min, Δt: 1.464 s, wall time: 637.609 ms
progress: 1.6 %, i: 0500, t: 11.193 min, Δt: 1.611 s, wall time: 637.539 ms
progress: 2.0 %, i: 0600, t: 14.145 min, Δt: 1.772 s, wall time: 659.771 ms
progress: 2.4 %, i: 0700, t: 17.393 min, Δt: 1.949 s, wall time: 644.049 ms
progress: 2.9 %, i: 0800, t: 20.966 min, Δt: 2.144 s, wall time: 664.023 ms
progress: 3.5 %, i: 0900, t: 24.896 min, Δt: 2.358 s, wall time: 648.265 ms
progress: 4.1 %, i: 1000, t: 29.219 min, Δt: 2.594 s, wall time: 658.792 ms
progress: 4.7 %, i: 1100, t: 33.974 min, Δt: 2.853 s, wall time: 654.423 ms
progress: 5.4 %, i: 1200, t: 39.205 min, Δt: 3.138 s, wall time: 660.702 ms
progress: 6.2 %, i: 1300, t: 44.958 min, Δt: 3.452 s, wall time: 645.074 ms
progress: 7.1 %, i: 1400, t: 51.287 min, Δt: 3.797 s, wall time: 655.633 ms
progress: 8.1 %, i: 1500, t: 58.250 min, Δt: 4.177 s, wall time: 664.161 ms
progress: 9.2 %, i: 1600, t: 1.098 hr, Δt: 4.595 s, wall time: 651.306 ms
progress: 10.3 %, i: 1700, t: 1.239 hr, Δt: 5.054 s, wall time: 652.581 ms
progress: 11.6 %, i: 1800, t: 1.393 hr, Δt: 5.560 s, wall time: 657.264 ms
progress: 13.0 %, i: 1900, t: 1.563 hr, Δt: 6.116 s, wall time: 651.743 ms
progress: 14.6 %, i: 2000, t: 1.750 hr, Δt: 6.727 s, wall time: 651.324 ms
progress: 16.2 %, i: 2100, t: 1.944 hr, Δt: 6.976 s, wall time: 664.054 ms
progress: 17.1 %, i: 2200, t: 2.047 hr, Δt: 3.705 s, wall time: 794.303 ms
progress: 18.0 %, i: 2300, t: 2.160 hr, Δt: 4.076 s, wall time: 717.391 ms
progress: 19.0 %, i: 2400, t: 2.285 hr, Δt: 4.484 s, wall time: 705.959 ms
progress: 20.2 %, i: 2500, t: 2.422 hr, Δt: 4.932 s, wall time: 694.088 ms
progress: 21.4 %, i: 2600, t: 2.572 hr, Δt: 5.425 s, wall time: 676.130 ms
progress: 22.8 %, i: 2700, t: 2.738 hr, Δt: 5.968 s, wall time: 662.178 ms
progress: 24.1 %, i: 2800, t: 2.894 hr, Δt: 5.633 s, wall time: 656.812 ms
progress: 25.6 %, i: 2900, t: 3.067 hr, Δt: 6.197 s, wall time: 659.816 ms
progress: 27.1 %, i: 3000, t: 3.250 hr, Δt: 6.599 s, wall time: 655.171 ms
progress: 28.6 %, i: 3100, t: 3.435 hr, Δt: 6.645 s, wall time: 658.980 ms
progress: 29.9 %, i: 3200, t: 3.583 hr, Δt: 5.337 s, wall time: 665.661 ms
progress: 31.1 %, i: 3300, t: 3.730 hr, Δt: 5.312 s, wall time: 682.003 ms
progress: 32.4 %, i: 3400, t: 3.883 hr, Δt: 5.479 s, wall time: 641.533 ms
progress: 33.4 %, i: 3500, t: 4.007 hr, Δt: 4.472 s, wall time: 629.381 ms
progress: 34.5 %, i: 3600, t: 4.143 hr, Δt: 4.920 s, wall time: 643.778 ms
progress: 35.7 %, i: 3700, t: 4.279 hr, Δt: 4.876 s, wall time: 636.688 ms
progress: 36.8 %, i: 3800, t: 4.413 hr, Δt: 4.843 s, wall time: 648.637 ms
progress: 38.0 %, i: 3900, t: 4.558 hr, Δt: 5.220 s, wall time: 649.205 ms
progress: 39.2 %, i: 4000, t: 4.710 hr, Δt: 5.445 s, wall time: 640.877 ms
progress: 40.5 %, i: 4100, t: 4.860 hr, Δt: 5.430 s, wall time: 647.273 ms
progress: 41.6 %, i: 4200, t: 4.997 hr, Δt: 4.918 s, wall time: 639.861 ms
progress: 42.8 %, i: 4300, t: 5.142 hr, Δt: 5.203 s, wall time: 635.089 ms
progress: 44.1 %, i: 4400, t: 5.296 hr, Δt: 5.564 s, wall time: 658.496 ms
progress: 45.1 %, i: 4500, t: 5.414 hr, Δt: 4.242 s, wall time: 645.006 ms
progress: 46.2 %, i: 4600, t: 5.544 hr, Δt: 4.666 s, wall time: 652.166 ms
progress: 47.4 %, i: 4700, t: 5.683 hr, Δt: 5.026 s, wall time: 652.683 ms
progress: 48.5 %, i: 4800, t: 5.826 hr, Δt: 5.127 s, wall time: 657.132 ms
progress: 49.7 %, i: 4900, t: 5.964 hr, Δt: 4.976 s, wall time: 655.491 ms
progress: 50.8 %, i: 5000, t: 6.094 hr, Δt: 4.681 s, wall time: 626.366 ms
progress: 51.7 %, i: 5100, t: 6.202 hr, Δt: 3.891 s, wall time: 654.503 ms
progress: 52.6 %, i: 5200, t: 6.317 hr, Δt: 4.160 s, wall time: 647.718 ms
progress: 53.7 %, i: 5300, t: 6.445 hr, Δt: 4.576 s, wall time: 645.440 ms
progress: 54.8 %, i: 5400, t: 6.574 hr, Δt: 4.662 s, wall time: 648.169 ms
progress: 55.9 %, i: 5500, t: 6.703 hr, Δt: 4.627 s, wall time: 638.133 ms
progress: 56.9 %, i: 5600, t: 6.832 hr, Δt: 4.661 s, wall time: 635.331 ms
progress: 57.9 %, i: 5700, t: 6.943 hr, Δt: 3.991 s, wall time: 638.023 ms
progress: 58.7 %, i: 5800, t: 7.049 hr, Δt: 3.817 s, wall time: 655.679 ms
progress: 59.7 %, i: 5900, t: 7.166 hr, Δt: 4.198 s, wall time: 638.579 ms
progress: 60.7 %, i: 6000, t: 7.279 hr, Δt: 4.083 s, wall time: 654.386 ms
progress: 61.5 %, i: 6100, t: 7.385 hr, Δt: 3.827 s, wall time: 648.350 ms
progress: 62.5 %, i: 6200, t: 7.496 hr, Δt: 3.967 s, wall time: 647.580 ms
progress: 63.3 %, i: 6300, t: 7.600 hr, Δt: 3.769 s, wall time: 632.759 ms
progress: 64.2 %, i: 6400, t: 7.707 hr, Δt: 3.827 s, wall time: 665.735 ms
progress: 65.1 %, i: 6500, t: 7.809 hr, Δt: 3.694 s, wall time: 650.852 ms
progress: 66.0 %, i: 6600, t: 7.922 hr, Δt: 4.056 s, wall time: 639.361 ms
progress: 67.0 %, i: 6700, t: 8.039 hr, Δt: 4.235 s, wall time: 639.133 ms
progress: 67.9 %, i: 6800, t: 8.143 hr, Δt: 3.721 s, wall time: 653.002 ms
progress: 68.8 %, i: 6900, t: 8.256 hr, Δt: 4.091 s, wall time: 639.474 ms
progress: 69.8 %, i: 7000, t: 8.373 hr, Δt: 4.185 s, wall time: 644.346 ms
progress: 70.8 %, i: 7100, t: 8.494 hr, Δt: 4.348 s, wall time: 644.694 ms
progress: 71.9 %, i: 7200, t: 8.623 hr, Δt: 4.644 s, wall time: 642.954 ms
progress: 72.7 %, i: 7300, t: 8.725 hr, Δt: 3.680 s, wall time: 632.170 ms
progress: 73.6 %, i: 7400, t: 8.837 hr, Δt: 4.048 s, wall time: 621.087 ms
progress: 74.5 %, i: 7500, t: 8.936 hr, Δt: 3.552 s, wall time: 642.349 ms
progress: 75.4 %, i: 7600, t: 9.044 hr, Δt: 3.907 s, wall time: 654.065 ms
progress: 76.3 %, i: 7700, t: 9.156 hr, Δt: 4.024 s, wall time: 653.217 ms
progress: 77.1 %, i: 7800, t: 9.257 hr, Δt: 3.630 s, wall time: 639.862 ms
progress: 77.9 %, i: 7900, t: 9.354 hr, Δt: 3.477 s, wall time: 631.578 ms
progress: 78.8 %, i: 8000, t: 9.451 hr, Δt: 3.495 s, wall time: 626.412 ms
progress: 79.5 %, i: 8100, t: 9.545 hr, Δt: 3.409 s, wall time: 635.222 ms
progress: 80.3 %, i: 8200, t: 9.641 hr, Δt: 3.438 s, wall time: 636.815 ms
progress: 81.2 %, i: 8300, t: 9.745 hr, Δt: 3.764 s, wall time: 646.578 ms
progress: 82.1 %, i: 8400, t: 9.854 hr, Δt: 3.900 s, wall time: 656.866 ms
progress: 83.0 %, i: 8500, t: 9.959 hr, Δt: 3.787 s, wall time: 650.627 ms
progress: 83.9 %, i: 8600, t: 10.068 hr, Δt: 3.939 s, wall time: 650.272 ms
progress: 84.7 %, i: 8700, t: 10.163 hr, Δt: 3.394 s, wall time: 643.055 ms
progress: 85.4 %, i: 8800, t: 10.253 hr, Δt: 3.269 s, wall time: 645.025 ms
progress: 86.3 %, i: 8900, t: 10.353 hr, Δt: 3.596 s, wall time: 619.382 ms
progress: 87.2 %, i: 9000, t: 10.460 hr, Δt: 3.858 s, wall time: 637.257 ms
progress: 88.1 %, i: 9100, t: 10.570 hr, Δt: 3.949 s, wall time: 647.092 ms
progress: 89.0 %, i: 9200, t: 10.675 hr, Δt: 3.779 s, wall time: 654.140 ms
progress: 89.8 %, i: 9300, t: 10.781 hr, Δt: 3.802 s, wall time: 636.021 ms
progress: 90.7 %, i: 9400, t: 10.881 hr, Δt: 3.591 s, wall time: 639.581 ms
progress: 91.5 %, i: 9500, t: 10.981 hr, Δt: 3.626 s, wall time: 647.625 ms
progress: 92.4 %, i: 9600, t: 11.090 hr, Δt: 3.917 s, wall time: 640.461 ms
progress: 93.4 %, i: 9700, t: 11.205 hr, Δt: 4.147 s, wall time: 643.316 ms
progress: 94.3 %, i: 9800, t: 11.320 hr, Δt: 4.115 s, wall time: 645.427 ms
progress: 95.2 %, i: 9900, t: 11.427 hr, Δt: 3.873 s, wall time: 637.281 ms
progress: 96.0 %, i: 10000, t: 11.523 hr, Δt: 3.441 s, wall time: 650.339 ms
progress: 96.9 %, i: 10100, t: 11.623 hr, Δt: 3.613 s, wall time: 654.138 ms
progress: 97.7 %, i: 10200, t: 11.726 hr, Δt: 3.692 s, wall time: 630.517 ms
progress: 98.6 %, i: 10300, t: 11.830 hr, Δt: 3.746 s, wall time: 652.241 ms
progress: 99.5 %, i: 10400, t: 11.939 hr, Δt: 3.933 s, wall time: 661.528 ms
progress: 100.3 %, i: 10500, t: 12.041 hr, Δt: 3.688 s, wall time: 650.323 ms

Plot the result at the end

makeplot!(axs, model)
gcf()

This page was generated using Literate.jl.