Simple diffusion example
This script provides our simplest example of Oceananigans.jl functionality: 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, and finally
- how to look at results.
Using Oceananigans.jl
To use Oceananigans.jl
after it has been installed, we bring Oceananigans.jl
functions and names into our 'namespace' by writing
using Oceananigans
In addition, we import the submodule Grids
.
using Oceananigans.Grids
Instantiating and configuring a model
To begin using Oceananigans, we instantiate an incompressible model by calling the IncompressibleModel
constructor:
model = IncompressibleModel(
grid = RegularCartesianGrid(size = (1, 1, 128), x = (0, 1), y = (0, 1), z = (-0.5, 0.5)),
closure = IsotropicDiffusivity(κ = 1.0)
)
The keyword arguments grid
and closure
indicate that our model grid is Cartesian with uniform grid spacing, that our diffusive stress and tracer fluxes are determined by diffusion with a constant diffusivity κ
(note that we do not use viscosity in this example).
Note that by default, a Model
has no-flux boundary condition on all variables. Next, we set an initial condition on our "passive tracer", temperature. Our objective is to observe the diffusion of a Gaussian.
# Build a Gaussian initial condition function with width `δ`:
δ = 0.1
Tᵢ(x, y, z) = exp(-z^2 / (2δ^2))
# Set `model.tracers.T` to the function `Tᵢ`:
set!(model, T=Tᵢ)
Running your first Model
Finally, we time-step the model forward using the function time_step!
, with a time-step size that ensures numerical stability.
# Time-scale for diffusion across a grid cell
cell_diffusion_time_scale = model.grid.Δz^2 / model.closure.κ.T
# We create a `Simulation` which will handle time stepping the model. It will
# execute `Nt` time steps with step size `Δt` using a second-order Adams-Bashforth method.
simulation = Simulation(model, Δt = 0.1 * cell_diffusion_time_scale, stop_iteration = 1000)
run!(simulation)
[ Info: Simulation is stopping. Model iteration 1000 has hit or exceeded simulation stop iteration 1000.
Visualizing the results
We use Plots.jl
to look at the results. Fields are stored as 3D arrays in Oceananigans so we plot interior(T)[1, 1, :]
which will return a 1D array.
using Plots, Printf
# A convenient function for generating a label with the current model time
tracer_label(time) = @sprintf("t = %.3f", time)
# Plot initial condition
T = model.tracers.T
z = znodes(T)[:]
p = plot(Tᵢ.(0, 0, z), z, linewidth=2, label="t = 0",
xlabel="Temperature", ylabel="z")
# Plot current solution
plot!(p, interior(T)[1, 1, :], z, linewidth=2, label=tracer_label(model.clock.time))
Interesting! Next, we add an output writer that saves the temperature field in JLD2 files, and run the simulation for longer so that we can animate the results.
using Oceananigans.OutputWriters: JLD2OutputWriter
simulation.output_writers[:temperature] =
JLD2OutputWriter(model, model.tracers, prefix = "one_dimensional_diffusion",
iteration_interval = 100, force = true)
# Run simulation for 10,000 more iterations
simulation.stop_iteration += 10000
run!(simulation)
nothing
[ 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=tracer_label(t),
label="", xlabel="Temperature", ylabel="z", xlims=(0, 1))
end
This page was generated using Literate.jl.