# 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, JLD2, Plots"

## 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 an NonhydrostaticModel by passing it a grid, plus information about the equations we would like to solve.

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

grid = RegularRectilinearGrid(size=128, z=(-0.5, 0.5), topology=(Flat, Flat, Bounded))
RegularRectilinearGrid{Float64, Flat, Flat, Bounded}
domain: x ∈ [0.0, 0.0], y ∈ [0.0, 0.0], z ∈ [-0.5, 0.5]
topology: (Flat, Flat, Bounded)
resolution (Nx, Ny, Nz): (1, 1, 128)
halo size (Hx, Hy, Hz): (0, 0, 1)
grid spacing (Δx, Δy, Δz): (0.0, 0.0, 0.0078125)

The default topology is (Periodic, Periodic, Bounded) but since we only want to solve a one-dimensional problem, we assign the x and y dimensions to Flat. We could specify each of them to be either Periodic or Bounded but that will define a halo in each of those directions, and that is numerically more costly. Note that we only specify the extent and size for the Bounded dimension.

We next specify a model with an IsotropicDiffusivity, which models either molecular or turbulent diffusion,

closure = IsotropicDiffusivity(κ=1.0)
IsotropicDiffusivity: ν=0.0, κ=1.0

We finally pass these two ingredients to NonhydrostaticModel,

model = NonhydrostaticModel(grid=grid, closure=closure)
NonhydrostaticModel{CPU, Float64}(time = 0 seconds, iteration = 0)
├── grid: RegularRectilinearGrid{Float64, Flat, Flat, Bounded}(Nx=1, Ny=1, Nz=128)
├── tracers: (:T, :S)
├── closure: IsotropicDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Float64, NamedTuple{(:T, :S), Tuple{Float64, Float64}}}
├── buoyancy: SeawaterBuoyancy{Float64, LinearEquationOfState{Float64}, Nothing, Nothing}
└── coriolis: Nothing

Our simple grid and model use a number of defaults:

• The default grid topology is periodic in x, y and bounded in z.
• The default Model has no-flux (insulating and stress-free) boundary conditions on non-periodic boundaries for velocities u, v, w and tracers.
• The default Model has two tracers: temperature T, and salinity S.
• The default Model uses a SeawaterBuoyancy model with a LinearEquationOfState. However, buoyancy is not active in the simulation we run below.

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(x, y, 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 Plots

z = znodes(model.tracers.T)

T_plot = plot(interior(model.tracers.T)[1, 1, :], z,
linewidth = 2,
label = "t = 0",
xlabel = "Temperature (ᵒC)",
ylabel = "z")

The function interior above extracts a view of the physical interior points of model.tracers.T. This is useful because model.tracers.T also contains "halo" points that lie outside the physical domain (halo points are used to set boundary conditions during time-stepping).

## Running a Simulation

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

# Time-scale for diffusion across a grid cell
diffusion_time_scale = model.grid.Δz^2 / model.closure.κ.T

simulation = Simulation(model, Δt = 0.1 * diffusion_time_scale, stop_iteration = 1000)
Simulation{typename(NonhydrostaticModel){typename(CPU), Float64}}
├── Model clock: time = 0 seconds, iteration = 0
├── Next time step (Float64): 6.104 μs
├── Iteration interval: 1
├── Stop criteria: Any[Oceananigans.Simulations.iteration_limit_exceeded, Oceananigans.Simulations.stop_time_exceeded, Oceananigans.Simulations.wall_time_limit_exceeded]
├── Run time: 0 seconds, wall time limit: Inf
├── Stop time: Inf years, stop iteration: 1000
├── Diagnostics: typename(OrderedCollections.OrderedDict) with 1 entry:
│   └── nan_checker => typename(NaNChecker)
└── Output writers: typename(OrderedCollections.OrderedDict) with no entries

simulation will run for 1000 iterations with a time-step that resolves the time-scale at which our temperature field diffuses. All that's left is to

run!(simulation)
[ Info: Updating model auxiliary state before the first time step...
[ Info:     ... updated in 1.199 ms.
[ Info: Executing first time step...
[ Info: Simulation is stopping. Model iteration 1000 has hit or exceeded simulation stop iteration 1000.

## Visualizing the results

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

using Printf

plot!(T_plot, interior(model.tracers.T)[1, 1, :], z, linewidth=2,
label=@sprintf("t = %.3f", model.clock.time))

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

using Oceananigans.OutputWriters: JLD2OutputWriter, IterationInterval

simulation.output_writers[:temperature] =
JLD2OutputWriter(model, model.tracers, prefix = "one_dimensional_diffusion",
schedule=IterationInterval(100), force = true)
JLD2OutputWriter scheduled on IterationInterval(100):
├── filepath: ./one_dimensional_diffusion.jld2
├── 2 outputs: (:T, :S)
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB

We run the simulation for 10,000 more iterations,

simulation.stop_iteration += 10000

run!(simulation)
[ Info: Updating model auxiliary state before the first time step...
[ Info:     ... updated in 1.250 ms.
[ Info: Simulation is stopping. Model iteration 11000 has hit or exceeded simulation stop iteration 11000.

Finally, we animate the results by opening the JLD2 file, extract the iterations we ended up saving at, and plot the evolution of the temperature profile in a loop over the iterations.

using JLD2

file = jldopen(simulation.output_writers[:temperature].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

anim = @animate for (i, iter) in enumerate(iterations)

T = file["timeseries/T/$iter"][1, 1, :] t = file["timeseries/t/$iter"]

plot(T, z, linewidth=2, title=@sprintf("t = %.3f", t),
label="", xlabel="Temperature", ylabel="z", xlims=(0, 1))
end