Simple diffusion example

This is Oceananigans.jl's simplest example: the diffusion of a one-dimensional Gaussian. This example demonstrates

• How to load Oceananigans.jl.

• How to instantiate an Oceananigans.jl model.

• How to create simple Oceananigans.jl output.

• How to set an initial condition with a function.

• How to time-step a model forward.

• How to look at results.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, CairoMakie"

Using Oceananigans.jl

Write

using Oceananigans

to load Oceananigans functions and objects into our script.

Instantiating and configuring a model

A core Oceananigans type is NonhydrostaticModel. We build a NonhydrostaticModel by passing it a grid, plus information about the equations we would like to solve.

Below, we build a rectilinear grid with 128 regularly-spaced grid points in the z-direction, where z spans from z = -0.5 to z = 0.5,

grid = RectilinearGrid(size=128, z=(-0.5, 0.5), topology=(Flat, Flat, Bounded))

1×1×128 RectilinearGrid{Float32, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── Flat x                   
├── Flat y                   
└── Bounded  z ∈ [-0.5, 0.5] regularly spaced with Δz=0.0078125

The default topology is (Periodic, Periodic, Bounded). In this example, we're trying to solve a one-dimensional problem, so we assign Flat to the x and y topologies. We excise halos and avoid interpolation or differencing in Flat directions, saving computation and memory.

We next specify a ScalarDiffusivity with diffusivity , which models either molecular or turbulent diffusion,

κ = 1
closure = ScalarDiffusivity(; κ)

ScalarDiffusivity{ExplicitTimeDiscretization}(ν=nothing, κ=1.0)

We finally pass these two ingredients to NonhydrostaticModel,

model = NonhydrostaticModel(grid; closure, tracers=:T)

NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×128 RectilinearGrid{Float32, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme:
│   ├── momentum: Centered(order=2)
│   └── T: Centered(order=2)
├── tracers: T
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=nothing, κ=(T=1.0,))
├── buoyancy: Nothing
└── coriolis: Nothing

By default, NonhydrostaticModel has no-flux (insulating and stress-free) boundary conditions on all fields.

Next, we set! an initial condition on the temperature field, model.tracers.T. Our objective is to observe the diffusion of a Gaussian.

width = 0.1
initial_temperature(z) = exp(-z^2 / (2width^2))
set!(model, T=initial_temperature)

Visualizing model data

Calling set! above changes the data contained in model.tracers.T, which was initialized as 0's when the model was created. To see the new data in model.tracers.T, we plot it:

using CairoMakie
set_theme!(Theme(fontsize = 20, linewidth=3))

axis = (xlabel = "Temperature (ᵒC)", ylabel = "z")
label = "t = 0"
lines(model.tracers.T; label, axis)

Running a Simulation

Next we set up a Simulation that time-steps the model forward and manages output.

Δz = minimum_zspacing(grid)
Δt = 0.1 * Δz^2 / κ

simulation = Simulation(model; Δt, stop_iteration = 1000)

Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 6.104 μs
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: Inf days
├── stop_iteration: 1000.0
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entries

simulation will run for 1000 iterations with a time-step that resolves the time-scale for diffusion across a grid cell. All that's left is to

run!(simulation)

[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (98.778 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (416.611 ms).
[ Info: Simulation is stopping after running for 554.400 ms.
[ Info: Model iteration 1000 equals or exceeds stop iteration 1000.

Visualizing the results

Let's look at how model.tracers.T changed during the simulation.

label = "t = $(round(model.clock.time, digits=3))"
lines!(model.tracers.T; label)
axislegend()

Very interesting! Next, we run the simulation a bit longer and make an animation. For this, we use the JLD2Writer to write data to disk as the simulation progresses.

simulation.output_writers[:temperature] =
    JLD2Writer(model, model.tracers,
               filename = "one_dimensional_diffusion.jld2",
               schedule = IterationInterval(100),
               overwrite_files = true)

JLD2Writer scheduled on IterationInterval(100):
├── filepath: one_dimensional_diffusion.jld2
├── 1 outputs: T
├── array_type: Array{Float32}
├── including: [:coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

We run the simulation for 10,000 more iterations,

simulation.stop_iteration += 10000
run!(simulation)

[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (833.459 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (67.591 μs).
[ Info: Simulation is stopping after running for 1.464 seconds.
[ Info: Model iteration 11000 equals or exceeds stop iteration 11000.

To animate the results, we load the saved temperature as a FieldTimeSeries and plot the temperature profile at each saved time.

T_timeseries = FieldTimeSeries("one_dimensional_diffusion.jld2", "T")
times = T_timeseries.times

fig = Figure()
ax = Axis(fig[2, 1]; xlabel = "Temperature (ᵒC)", ylabel = "z")
xlims!(ax, 0, 1)

n = Observable(1)

T = @lift T_timeseries[$n]
lines!(ax, T)

label = @lift "t = $(round(times[$n], digits=3))"
Label(fig[1, 1], label, tellwidth=false)

fig

Finally, we record a movie.

frames = 1:length(times)

@info "Making an animation..."

record(fig, "one_dimensional_diffusion.mp4", frames, framerate=24) do i
    n[] = i
end

[ Info: Making an animation...

---

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.13.2
Commit 45b6d15e34b (2026-10-14 20:33 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 128 × AMD EPYC 9374F 32-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.2+0.x64.linux.gnu/local/share/julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.2+0.x64.linux.gnu/share/julia
  JULIA_VERSION_ENZYME = 1.13.1
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  LD_LIBRARY_PATH = 
  JULIA_MAX_NUM_PRECOMPILE_FILES = 24
  JULIA_VERSION = 1.13.2
  JULIA_CUDA_USE_COMPAT = false
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-35368/docs/
  JULIA_DEBUG = Literate

These were the top-level packages installed in the environment:

import Pkg
Pkg.status()

Status `~/Oceananigans.jl-35368/docs/Project.toml`
  [79e6a3ab] Adapt v4.7.4
⌃ [052768ef] CUDA v6.2.2
  [13f3f980] CairoMakie v0.15.16
  [8e50ac2c] ConservativeRegridding v0.2.18
⌅ [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.5.0
  [4710194d] DocumenterVitepress v0.3.9
  [7da242da] Enzyme v0.13.216
  [033835bb] JLD2 v0.6.7
⌅ [63c18a36] KernelAbstractions v0.9.45
  [98b081ad] Literate v2.21.0
  [da04e1cc] MPI v0.20.29
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.116.1 `..`
  [429524aa] Optim v2.3.2
  [f27b6e38] Polynomials v4.1.4
⌅ [3c362404] Reactant v0.2.292
  [6038ab10] Rotations v1.7.1
  [d496a93d] SeawaterPolynomials v0.3.11
  [09ab397b] StructArrays v0.7.4
  [bdfc003b] TimesDates v0.3.3
  [0a941bbe] Zarr v0.11.2
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.13.0
  [44cfe95a] Pkg v1.13.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

---

This page was generated using Literate.jl.