### Minimal Simulation Setup

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/simulations/simulations_overview.md

A basic example demonstrating the creation of a `Simulation` object, setting a time-step and stop iteration, and running the simulation.

```julia
using Oceananigans

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid)
simulation = Simulation(model; Δt=7, stop_iteration=6)
run!(simulation)

simulation
```

--------------------------------

### Set up a simple simulation

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/simulations/schedules.md

Initializes a basic simulation with a grid and model, and defines a dummy callback function for demonstration purposes. This setup is reused across examples.

```julia
using Oceananigans
grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
model = NonhydrostaticModel(grid)
simulation = Simulation(model, Δt=0.1, stop_time=2.5, verbose=false)
dummy(sim) = @info string("Iter: ", iteration(sim), " -- I was called at t = ", time(sim), 
                          " and wall time = ", prettytime(sim.run_wall_time))

```

--------------------------------

### Install Documentation Environment

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/developer_docs/contributing.md

Installs the documentation environment and develops the package locally. This is a prerequisite for building the documentation.

```julia
julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
```

--------------------------------

### Install Pre-commit Hooks (Original)

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/developer_docs/contributing.md

Installs the original pre-commit hooks for the repository. This ensures code formatting consistency.

```bash
pre-commit install
```

--------------------------------

### Instantiate Development Environment

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/developer_docs/contributing.md

Set up the project's development environment by installing all necessary dependencies using Julia's Pkg manager.

```julia
julia --project
] instantiate
```

--------------------------------

### Install Pre-commit Hooks (Prek)

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/developer_docs/contributing.md

Installs pre-commit hooks using the prek manager. This is an alternative to the original pre-commit tool for ensuring code style.

```bash
prek install
```

--------------------------------

### Example of using `launch!` for GPU kernels

Source: https://github.com/clima/oceananigans.jl/blob/main/AGENTS.md

Demonstrates the correct way to launch a kernel for grid point computations on a GPU. Avoid looping over grid points outside of kernels.

```julia
using Oceananigans

model = NonhydrostaticModel(grid)

# Compute tendencies for velocity and tracer fields
t = model.clock.time

compute_tendencies!(model, model.tracers, model.velocities)

```

--------------------------------

### Setup CPU Simulation

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/quick_start.md

Sets up a 2D RectilinearGrid and a NonhydrostaticModel for CPU execution. Initializes velocity fields with random noise and configures a simulation to run for 100 iterations.

```julia
using CairoMakie
CairoMakie.activate!(type = "png")
```

```julia
using Oceananigans

grid = RectilinearGrid(size = (128, 128),
                       x = (0, 2π),
                       y = (0, 2π),
                       topology = (Periodic, Periodic, Flat))

model = NonhydrostaticModel(grid; advection=WENO())

ϵ(x, y) = 2rand() - 1
set!(model, u=ϵ, v=ϵ)

simulation = Simulation(model; Δt=0.01, stop_iteration=100)
run!(simulation)
```

--------------------------------

### Create GPU-enabled Grid and Model

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/simulation_tips.md

Demonstrates how to create a RectilinearGrid on a GPU and initialize a NonhydrostaticModel. This setup is necessary for GPU computations.

```julia
using Oceananigans, CUDA

grid = RectilinearGrid(GPU(); size=(1, 1, 1), extent=(1, 1, 1), halo=(1, 1, 1))
model = NonhydrostaticModel(grid)
typeof(model.velocities.u.data)
```

--------------------------------

### Create a Rectilinear Grid on a GPU

Source: https://github.com/clima/oceananigans.jl/wiki/Accessing-GPUs-and-using-Oceananigans-on-GPUs

This snippet demonstrates how to initialize a RectilinearGrid on a GPU. Ensure CUDA is installed and accessible.

```julia
using Oceananigans
grid = RectilinearGrid(GPU(), size=(128, 128, 64), x=(0, 1), y=(0, 1), z=(0, 1))
```

--------------------------------

### Instantiate Convergence Test Environment

Source: https://github.com/clima/oceananigans.jl/blob/main/validation/convergence_tests/README.md

Activates the project environment, installs dependencies, and develops the local Oceananigans package. Run this once before executing any convergence tests.

```julia
julia -e 'using Pkg; Pkg.activate(pwd()); Pkg.instantiate(); Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..", "..")))'
```

--------------------------------

### Setting up an Immersed Boundary Grid

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/boundary_conditions.md

This example demonstrates how to create an ImmersedBoundaryGrid with a fitted bottom boundary. It also shows how to define no-slip boundary conditions for velocity fields on the immersed boundary.

```jldoctest
hill(x, y) = 0.1 + 0.1 * exp(-x^2 - y^2)
underlying_grid = RectilinearGrid(size=(32, 32, 16), x=(-3, 3), y=(-3, 3), z=(0, 1), topology=(Periodic, Periodic, Bounded))
grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(hill))

# Create a no-slip boundary condition for velocity fields.
# Note that the no-slip boundary condition is _only_ applied on immersed boundaries.
velocity_bcs = FieldBoundaryConditions(immersed=ValueBoundaryCondition(0))
model = NonhydrostaticModel(grid; boundary_conditions=(u=velocity_bcs, v=velocity_bcs, w=velocity_bcs))
```

--------------------------------

### RectilinearGrid Topology Examples

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Illustrates different topology configurations for a RectilinearGrid, specifying periodicity and boundedness for each dimension (x, y, z).

```julia
topology = (Periodic, Periodic, Periodic) # triply periodic
```

```julia
topology = (Periodic, Periodic, Bounded)  # periodic in x, y, bounded in z
```

```julia
topology = (Periodic, Bounded, Bounded)   # periodic in x, but bounded in y, z (a "channel")
```

```julia
topology = (Bounded, Bounded, Bounded)    # bounded in x, y, z (a closed box)
```

```julia
topology = (Periodic, Periodic, Flat)     # two-dimensional, doubly-periodic in x, y (a torus)
```

```julia
topology = (Periodic, Flat, Flat)         # one-dimensional, periodic in x (a line)
```

```julia
topology = (Flat, Flat, Bounded)          # one-dimensional and bounded in z (a single column)
```

--------------------------------

### Setting Up Immersed Boundary Conditions (No-Slip)

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/boundary_conditions.md

This example demonstrates setting up an `ImmersedBoundaryGrid` with a no-slip boundary condition for velocity fields. The no-slip condition is applied only on the immersed boundary.

```jldoctest
# Generate a simple ImmersedBoundaryGrid
hill(x, y) = 0.1 + 0.1 * exp(-x^2 - y^2)
underlying_grid = RectilinearGrid(size=(32, 32, 16), x=(-3, 3), y=(-3, 3), z=(0, 1), topology=(Periodic, Periodic, Bounded))
grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(hill))

# Create a no-slip boundary condition for velocity fields.
# Note that the no-slip boundary condition is _only_ applied on immersed boundaries.
velocity_bcs = FieldBoundaryConditions(immersed=ValueBoundaryCondition(0))
model = NonhydrostaticModel(grid; boundary_conditions=(u=velocity_bcs, v=velocity_bcs, w=velocity_bcs))
```

--------------------------------

### Simulation Setup for LorenzModel

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/developer_docs/model_interface.md

Sets up a Simulation with the LorenzModel, records its state trajectory, and runs the simulation. This demonstrates how to integrate a custom model into the Oceananigans simulation framework.

```julia
lorenz = LorenzModel()
set!(lorenz, x=1)
simulation = Simulation(lorenz; Δt=0.01, stop_time=100, verbose=false)
trajectory = NTuple{3, Float64}[]

function record_trajectory!(sim)
    push!(trajectory, values(sim.model.state))
end

add_callback!(simulation, record_trajectory!)
run!(simulation)
nothing # hide
```

--------------------------------

### Field Location Example

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/operations.md

Demonstrates the default location of a field and its derivatives on a staggered grid. Note the shift in x-location for x-derivatives.

```julia
location(c) = (Center, Center, Center)
location(dx_c) = (Face, Center, Center)
```

--------------------------------

### Running a Distributed Grid Example with 2 Ranks

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Executes the previously defined distributed grid creation script using 2 MPI ranks. This demonstrates how the global grid size is divided among the ranks.

```bash
using MPI
run(`$(mpiexec()) -n 2 $(Base.julia_cmd()) --project distributed_grid_example.jl`)
nothing # hide
```

--------------------------------

### Setup GPU Simulation with Enhanced Features

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/quick_start.md

Sets up a high-resolution 2D simulation on the GPU, including a passive tracer and adaptive time-stepping using `TimeStepWizard`. Initializes velocity and tracer fields and runs the simulation until a specified stop time.

```julia
using CairoMakie
CairoMakie.activate!(type = "png")
```

```julia
using Oceananigans
using CairoMakie
using CUDA

grid = RectilinearGrid(GPU(),
                       size = (1024, 1024),
                       x = (-π, π),
                       y = (-π, π),
                       topology = (Periodic, Periodic, Flat))

model = NonhydrostaticModel(grid; advection=WENO(), tracers=:c)

δ = 0.5
cᵢ(x, y) = exp(-(x^2 + y^2) / 2δ^2)
ϵ(x, y) = 2rand() - 1
set!(model, u=ϵ, v=ϵ, c=cᵢ)

simulation = Simulation(model; Δt=1e-3, stop_time=10)
conjure_time_step_wizard!(simulation, cfl = 0.2, IterationInterval(10))
run!(simulation)
```

--------------------------------

### Setting Field Values with an Array

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/fields.md

Shows how to initialize a Field with values from a multi-dimensional array using `set!`. This is useful for setting complex initial conditions from pre-computed data.

```julia
using Random
Random.seed!(123)
random_stuff = rand(size(c)...)
set!(c, random_stuff)

heatmap(view(c, :, :, 1))
```

--------------------------------

### Running a Distributed Architecture Example

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

This snippet demonstrates how to execute a Julia script using MPI for distributed computing from the terminal. It shows how to launch a script with multiple MPI ranks and manage its execution.

```bash
using MPI
run(`$(mpiexec()) -n 2 $(Base.julia_cmd()) --project distributed_arch_example.jl`)
rm("distributed_arch_example.jl")
nothing # hide
```

--------------------------------

### Running a Custom Partition Example with 2 Ranks

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Executes the script that defines a custom partition for distributed grids, using 2 MPI ranks. This demonstrates the application of the custom partition strategy.

```bash
using MPI
run(`$(mpiexec()) -n 2 $(Base.julia_cmd()) --project partition_example.jl`)
nothing # hide
```

--------------------------------

### Define and Run Distributed Grid Example

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Sets up a Julia program to demonstrate distributed grid creation using `Oceananigans.DistributedComputations` and `MPI`. The program initializes MPI, creates a `Distributed()` architecture, and prints it on ranks 0 and 1. It's written to a file and can be executed via `mpiexec`.

```julia
using Oceananigans
using Oceananigans.DistributedComputations
using MPI; MPI.Init()
architecture = Distributed()
@onrank 0 @show architecture
@onrank 1 @show architecture

```

--------------------------------

### Inspecting Applied Boundary Conditions

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/boundary_conditions.md

After creating the model, you can inspect the applied boundary conditions on specific fields like velocity ('u') and tracers ('c') to verify the setup.

```jldoctest
model.velocities.u
```

```jldoctest
model.tracers.c
```

--------------------------------

### Create New Project Environment and Add Oceananigans

Source: https://github.com/clima/oceananigans.jl/wiki/Productive-Oceananigans-workflows-and-Julia-environments

Creates a new directory for a project, initializes a Project.toml file, and adds the latest tagged version of Oceananigans to the environment. This ensures all necessary dependencies are installed.

```julia
mkdir WowCoolSciences
cd WowCoolSciences
touch Project.toml
julia --project -e 'using Pkg; Pkg.add("Oceananigans")'
```

--------------------------------

### Running a Distributed Grid Example with 3 Ranks

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Executes the distributed grid creation script using 3 MPI ranks. This shows how the global grid size is further divided, resulting in smaller local grid sizes per rank.

```bash
run(`$(mpiexec()) -n 3 $(Base.julia_cmd()) --project distributed_grid_example.jl`)
nothing # hide
```

--------------------------------

### Initialize Model with Discrete Forcing

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/forcing_functions.md

Initializes a NonhydrostaticModel with defined discrete forcing functions for velocity and tracers. This setup allows for complex physics to be incorporated into the simulation.

```julia
grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; tracers=:b, buoyancy=BuoyancyTracer(), forcing=(u=u_forcing, b=b_forcing))
```

--------------------------------

### Run a 2D Turbulence Simulation

Source: https://github.com/clima/oceananigans.jl/blob/main/README.md

Sets up and runs a basic two-dimensional, horizontally-periodic simulation of turbulence using Oceananigans. This example uses a RectilinearGrid on the CPU.

```julia
using Oceananigans
grid = RectilinearGrid(CPU(), size=(128, 128), x=(0, 2π), y=(0, 2π), topology=(Periodic, Periodic, Flat))
model = NonhydrostaticModel(grid; advection=WENO())
ϵ(x, y) = 2rand() - 1
set!(model, u=ϵ, v=ϵ)
simulation = Simulation(model; Δt=0.01, stop_time=4)
run!(simulation)
```

--------------------------------

### Example of FieldBoundaryConditions with custom top boundary

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/boundary_conditions.md

Demonstrates how to set custom boundary conditions for a field, specifically overriding the default for the top boundary while keeping defaults for others. Useful for scenarios like free-slip surfaces.

```julia
julia> no_slip_bc = ValueBoundaryCondition(0)
ValueBoundaryCondition: 0

julia> free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing))
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
├── east: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
├── south: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
├── north: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
├── bottom: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
├── top: FluxBoundaryCondition: Nothing
└── immersed: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
```

--------------------------------

### Evenly Partition Grid Across X and Y

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Allocate half the ranks to the y-dimension and the rest to the x-dimension using `Partition(x=Rx, y=Equal())`. This example demonstrates how to create and display grids on each rank when partitioning is applied.

```julia
using Oceananigans
using Oceananigans.DistributedComputations: Equal, barrier
using MPI
MPI.Init()

# Total number of ranks
Nr = MPI.Comm_size(MPI.COMM_WORLD)

# Allocate half the ranks to y, and the rest to x
Rx = Nr ÷ 2
partition = Partition(x=Rx, y=Equal())
arch = Distributed(CPU(); partition)

grid = RectilinearGrid(arch,
                       size = (48, 48, 16),
                       x = (0, 64),
                       y = (0, 64),
                       z = (0, 16),
                       topology = (Periodic, Periodic, Bounded))

# Let's see all the grids this time.
for r in 0:Nr-1
    if r == arch.local_rank
        msg = string("On rank ", r, ":", '\n', '\n',
                     arch, '\n',
                     grid)
        @info msg
    end

    barrier(arch)
end
```

--------------------------------

### Testing Closure Construction and Basic Model Behavior

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/developer_docs/turbulence_closures.md

Verify that a closure can be constructed with default parameters and that a model using this closure behaves as expected during a time step. This example checks model type and clock time after a single time step.

```julia
using Test

closure = PacanowskiPhilanderVerticalDiffusivity()
grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; closure=closure, buoyancy=BuoyancyTracer(), tracers=:b)

@test model isa NonhydrostaticModel

time_step!(model, 1)
@test model.clock.time == 1
```

--------------------------------

### Configure NonhydrostaticModel with WENO Advection, Buoyancy, Coriolis, and Flux Boundary Conditions

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/models_overview.md

Configures a NonhydrostaticModel with advanced settings including a ninth-order WENO advection scheme, temperature/salinity-based buoyancy, FPlane Coriolis force, and flux boundary conditions for momentum and a passive tracer. This setup is suitable for more complex simulations.

```julia
using Oceananigans

grid = RectilinearGrid(size=(128, 128), halo=(5, 5), x=(0, 256), z=(-128, 0),
                       topology = (Periodic, Flat, Bounded))

# Numerical method and physics choices
advection = WENO(order=9) # ninth‑order upwind for momentum and tracers
buoyancy = BuoyancyTracer()
coriolis = FPlane(f=1e-4)

τx = - 8e-5 # m² s⁻² (τₓ < 0 ⟹ eastward wind stress)
u_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(τx))

@inline Jc(x, t, Lx) = cos(2π / Lx * x)
c_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(Jc, parameters=grid.Lx))

model = NonhydrostaticModel(grid; advection, buoyancy, coriolis,
                            tracers = (:b, :c),
                            boundary_conditions = (; u=u_bcs, c=c_bcs))
```

--------------------------------

### Install MPICH_jll before Oceananigans on M1

Source: https://github.com/clima/oceananigans.jl/wiki/Installation-and-getting-started-with-Oceananigans

Installs a specific version of MPICH_jll before installing Oceananigans. This is a workaround for issues on Apple M1 chips with Julia v1.7.

```julia
julia> using Pkg; Pkg.add(name="MPICH_jll", version="4.0.1"); Pkg.build()

julia> Pkg.add("Oceananigans")
```

--------------------------------

### Initialize Distributed Grid

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Demonstrates how to initialize a distributed grid on a CPU architecture. This is useful for utilizing multiple CPUs for computation.

```julia
using Oceananigans

child_architecture = CPU()
architecture = Distributed(child_architecture)
```

--------------------------------

### AndSchedule Example for Combined Conditions

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/simulations/schedules.md

Triggers a callback only if all child schedules actuate in the same iteration. This example combines an iteration interval with a time-based condition.

```julia
Oceananigans.Simulations.reset!(simulation)
simulation.stop_time = 2.5

after_one_point_seven(model) = model.clock.time > 1.7
schedule = AndSchedule(IterationInterval(2), after_one_point_seven)
add_callback!(simulation, dummy, schedule, name=:dummy)
run!(simulation)
```

--------------------------------

### Simulations

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/appendix/library.md

Module for managing simulation setups in Oceananigans.jl.

```APIDOC
## Simulations

```@autodocs
Modules = [Oceananigans.Simulations]
Private = false
```
```

--------------------------------

### Setting Up and Running a Simulation

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/simulations/simulations_overview.md

This snippet demonstrates how to initialize a NonhydrostaticModel, set up a Simulation with adaptive time-stepping, add a progress callback, configure a JLD2 output writer, and run the simulation.

```julia
using Oceananigans.Units: hours, minutes

model = NonhydrostaticModel(grid; tracers=:T)
simulation = Simulation(model, Δt=20, stop_time=2hours)

progress(sim) = @info "t = $(prettytime(sim)), Δt = $(prettytime(sim.Δt))"
add_callback!(simulation, progress, IterationInterval(10))
conjure_time_step_wizard!(simulation, cfl=0.8)

simulation.output_writers[:snapshots] = JLD2Writer(simulation.model, simulation.model.velocities;
                                                   filename = "snapshots.jld2",
                                                   schedule = TimeInterval(30minutes))

run!(simulation)
```

--------------------------------

### Simulation with Multiple Callbacks

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/simulations/simulations_overview.md

Illustrates setting up a `Simulation` with multiple callbacks, including one for printing progress at iteration intervals and another for reporting wall time.

```julia
using Oceananigans

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid)
simulation = Simulation(model; Δt=7, stop_time=14)

print_progress(sim) = @info string("Iteration: ", iteration(sim), ", time: ", time(sim))
add_callback!(simulation, print_progress, IterationInterval(2))

declare_time(sim) = @info string("The simulation has been running for ", prettytime(sim.run_wall_time), "!")
add_callback!(simulation, declare_time, TimeInterval(10))

run!(simulation)
```

--------------------------------

### Check Oceananigans Package Status

Source: https://github.com/clima/oceananigans.jl/blob/main/README.md

Verifies the installed version of the Oceananigans package.

```julia
using Pkg

pkg> Pkg.status("Oceananigans")
```

--------------------------------

### FluxBoundaryCondition with BulkDragFunction

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/boundary_conditions.md

Example of a FluxBoundaryCondition using BulkDragFunction for quadratic drag formulation. The direction is inferred during model construction.

```jldoctest
using Oceananigans

FluxBoundaryCondition(BulkDragFunction(QuadraticFormulation(), XDirection(), Cᴰ=0.002))
```

--------------------------------

### Constructing and Setting a CenterField

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/operations.md

Creates a RectilinearGrid and a CenterField, then sets its values using a function. This is a common setup for simulations.

```jldoctest
using Oceananigans

grid = RectilinearGrid(topology = (Periodic, Flat, Bounded),
                       size = (4, 4),
                       x = (0, 2π),
                       z = (-4, 0))

c = CenterField(grid)

periodic_but_decaying(x, z) = sin(x) * exp(z)
set!(c, periodic_but_decaying)

# output
4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 10×1×10 OffsetArray(::Array{Float64, 3}, -2:7, 1:1, -2:7) with eltype Float64 with indices -2:7×1:1×-2:7
    └── max=0.428882, min=-0.428882, mean=1.04083e-17

```

--------------------------------

### Add Oceananigans Package

Source: https://github.com/clima/oceananigans.jl/wiki/Installation-and-getting-started-with-Oceananigans

Installs the Oceananigans package using Julia's package manager. This is the standard method for most users.

```julia
julia> using Pkg

julia> Pkg.add("Oceananigans")
```

--------------------------------

### Build and Preview Documentation

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/developer_docs/contributing.md

Builds the documentation locally and provides debug information. The built documentation can be previewed in a web browser.

```julia
JULIA_DEBUG=Documenter julia --project=docs/ docs/make.jl
```

--------------------------------

### Incorporating ImmersedBoundaryCondition into FieldBoundaryConditions

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/boundary_conditions.md

This example demonstrates how to use a pre-defined ImmersedBoundaryCondition, like the bottom drag condition, within the FieldBoundaryConditions for a model.

```jldoctest
velocity_bcs = FieldBoundaryConditions(immersed=bottom_drag_bc)
```

--------------------------------

### Tupled Closures for Summed Contributions

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/turbulence_closures.md

Combine multiple closures in a tuple to sum their contributions. This example shows horizontal and vertical viscosity.

```jldoctest
julia> using Oceananigans

julia> closure = (HorizontalScalarDiffusivity(ν=1e-3), VerticalScalarDiffusivity(ν=1e-4))
(HorizontalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.001, κ=0.0), VerticalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0001, κ=0.0))
```

--------------------------------

### Setting up a Rectilinear Grid

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/boundary_conditions.md

Demonstrates the creation of a 3D rectilinear grid with mixed topology (Periodic, Bounded, Bounded). This is a prerequisite for defining boundary conditions in specific directions.

```julia
using Oceananigans

using Random
Random.seed!(1234)

grid = RectilinearGrid(size=(16, 16, 16), x=(0, 2π), y=(0, 1), z=(0, 1), topology=(Periodic, Bounded, Bounded))
```

--------------------------------

### Write Lagrangian Particles to JLD2

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/lagrangian_particles.md

Saves Lagrangian particle data to a JLD2 file. Pass `model.particles` as part of the output tuple to the `JLD2Writer`.

```julia
JLD2Writer(model, (; model.particles), filename="particles", schedule=TimeInterval(15))
```

--------------------------------

### Visualize ImmersedBoundaryGrid Bottom Height

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Visualizes the bottom height of an immersed boundary grid using a heatmap. This helps in understanding the shape of the topography defined.

```julia
using CairoMakie

h = mountain_grid.immersed_boundary.bottom_height

fig = Figure()
ax = Axis(fig[2, 1], xlabel="x (m)", ylabel="y (m)", aspect=1)
hm = heatmap!(ax, h)
Colorbar(fig[1, 1], hm, vertical=false, label="Bottom height (m)")

fig
```

--------------------------------

### Get z-spacings at Face location

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/fields.md

Retrieves the spacings of grid cells along the z-dimension when located at the faces (interfaces) of the cells. This is relevant for staggered grids.

```jldoctest
zspacings(grid, Face())[:, :, 1:5]

# output
5-element Vector{Float64}:
 0.1
 0.15000000000000002
 0.24999999999999994
 0.3500000000000001
 0.3999999999999999

```

--------------------------------

### Get z-spacings at Center location

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/fields.md

Retrieves the spacings of grid cells along the z-dimension when located at the center of the cells. This is useful for understanding grid resolution.

```jldoctest
zspacings(grid, Center())[:, :, 1:4]

# output
4-element Vector{Float64}:
 0.1
 0.19999999999999998
 0.3
 0.4

```

--------------------------------

### Build Vertical Grids with ReferenceToStretchedDiscretization

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Demonstrates building three single-column vertical grids with varying power-law stretching factors and a maximum stretching extent. Useful for controlling spacing at the ocean surface.

```julia
using Oceananigans
using CairoMakie
set_theme!(Theme(fontsize=16))
```

```julia
bias = :right
bias_edge = 0
extent = 800
constant_spacing = 25
constant_spacing_extent = 160

z = ReferenceToStretchedDiscretization(; extent, bias, bias_edge,
                                       constant_spacing,  constant_spacing_extent,
                                       stretching =  PowerLawStretching(1.06))
grid = RectilinearGrid(; size=length(z), z, topology=(Flat, Flat, Bounded))
zf = znodes(grid, Face())
zc = znodes(grid, Center())
Δz = zspacings(grid, Center())
Δz = view(Δz, 1, 1, :)  # for plotting

fig = Figure(size=(800, 550), colgap = 5)

axΔz1 = Axis(fig[1, 1];
             xlabel = "z-spacing (m)",
             ylabel = "z (m)",
             title = "PowerLawStretching(1.06)\n $(length(zf)) cells\n bottom @ z = $(zf[1]) m\n ")

axz1 = Axis(fig[1, 2])

ldepth = hlines!(axΔz1, bias_edge - extent, color = :salmon, linestyle=:dash)
lzbottom = hlines!(axΔz1, zf[1], color = :grey)
scatter!(axΔz1, Δz, zc)
hidespines!(axΔz1, :t, :r)

lines!(axz1, [0, 0], [zf[1], 0], color=:gray)
scatter!(axz1, 0 * zf, zf, marker=:hline, color=:gray, markersize=20)
scatter!(axz1, 0 * zc, zc)
hidedecorations!(axz1)
hidespines!(axz1)


z = ReferenceToStretchedDiscretization(; extent, bias, bias_edge,
                                       constant_spacing,  constant_spacing_extent,
                                       stretching =  PowerLawStretching(1.03))
grid = RectilinearGrid(; size=length(z), z, topology=(Flat, Flat, Bounded))
zf = znodes(grid, Face())
zc = znodes(grid, Center())
Δz = zspacings(grid, Center())
Δz = view(Δz, 1, 1, :)  # for plotting

axΔz2 = Axis(fig[1, 3];
             xlabel = "z-spacing (m)",
             ylabel = "z (m)",
             title = "PowerLawStretching(1.03)\n $(length(zf)) cells\n bottom @ z = $(zf[1]) m\n ")
axz2 = Axis(fig[1, 4])

ldepth = hlines!(axΔz2, bias_edge - extent, color = :salmon, linestyle=:dash)
lzbottom = hlines!(axΔz2, zf[1], color = :grey)
scatter!(axΔz2, Δz, zc)
hidespines!(axΔz2, :t, :r)

lines!(axz2, [0, 0], [zf[1], 0], color=:gray)
scatter!(axz2, 0 * zf, zf, marker=:hline, color=:gray, markersize=20)
scatter!(axz2, 0 * zc, zc)
hidedecorations!(axz2)
hidespines!(axz2)

z = ReferenceToStretchedDiscretization(; extent, bias, bias_edge,
                                       constant_spacing,  constant_spacing_extent,
                                       stretching =  PowerLawStretching(1.03),
                                       maximum_stretching_extent =  500)

grid = RectilinearGrid(; size=length(z), z, topology=(Flat, Flat, Bounded))
zf = znodes(grid, Face())
zc = znodes(grid, Center())
Δz = zspacings(grid, Center())
Δz = view(Δz, 1, 1, :)  # for plotting

axΔz3 = Axis(fig[1, 5];
             xlabel = "z-spacing (m)",
             ylabel = "z (m)",
             title = "PowerLawStretching(1.03)\n $(length(zf)) cells\n bottom @ z = $(zf[1]) m\n maximum_stretching_extent = 500")
axz3 = Axis(fig[1, 6])

ldepth = hlines!(axΔz3, bias_edge - extent, color = :salmon, linestyle=:dash)
lzbottom = hlines!(axΔz3, zf[1], color = :grey)
scatter!(axΔz3, Δz, zc)

hidespines!(axΔz3, :t, :r)

lines!(axz3, [0, 0], [zf[1], 0], color=:gray)
scatter!(axz3, 0 * zf, zf, marker=:hline, color=:gray, markersize=20)
scatter!(axz3, 0 * zc, zc)
hidedecorations!(axz3)
hidespines!(axz3)


linkaxes!(axΔz1, axz1, axΔz2, axz2, axΔz3, axz3)

Legend(fig[2, :], [ldepth, lzbottom], ["prescribed extent", "bottom z interface"], orientation = :horizontal)

colsize!(fig.layout, 2, Relative(0.1))

```

--------------------------------

### Import Libraries and Check Versions

Source: https://github.com/clima/oceananigans.jl/blob/main/validation/regridding/xesmf_python_regridding.ipynb

Imports necessary libraries like numpy, xesmf, xarray, and matplotlib. It also checks the installed versions of xesmf and esmpy.

```python
import numpy as np
import xesmf as xe
import xarray as xr
import matplotlib.pyplot as plt
xe.__version__
```

```python
import esmpy
esmpy.__version__
```

--------------------------------

### Custom Callback at Tendency Callsite

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/simulations/simulations_overview.md

Implement a custom callback that modifies model tendencies before a time-step. This example uses `TendencyCallsite` to hook into the tendency calculation phase.

```julia
using Oceananigans
using Oceananigans: TendencyCallsite

grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
model = NonhydrostaticModel(grid)
simulation = Simulation(model, Δt=1, stop_iteration=10)

function modify_tendency!(model, params)
    model.timestepper.Gⁿ[params.c] .+= params.δ
    return nothing
end

simulation.callbacks[:modify_u] = Callback(modify_tendency!, IterationInterval(1),
                                           callsite = TendencyCallsite(),
                                           parameters = (c = :u, δ = 1))

run!(simulation)
```

--------------------------------

### Create RectilinearGrid on CPU

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Demonstrates creating a simple rectilinear grid on the CPU. This is the default architecture if not specified.

```jldoctest
grid     = RectilinearGrid(size=3, z=(0, 1), topology=(Flat, Flat, Bounded))
cpu_grid = RectilinearGrid(CPU(), size=3, z=(0, 1), topology=(Flat, Flat, Bounded))

grid == cpu_grid
```

--------------------------------

### Basic Forcing Function Example

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/forcing_functions.md

Defines a simple forcing function for the 'u' velocity component and initializes a model with it. The forcing depends on spatial coordinates and time.

```jldoctest
u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t)

grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; forcing=(u=u_forcing,))

model.forcing.u

# output
ContinuousForcing{Nothing} at (Face, Center, Center)
├── func: u_forcing (generic function with 1 method)
├── parameters: nothing
└── field dependencies: ()

```

--------------------------------

### Model constructor with positional arguments

Source: https://github.com/clima/oceananigans.jl/blob/main/AGENTS.md

Shows the correct syntax for instantiating a `ShallowWaterModel`, where `grid` and `gravitational_acceleration` are positional arguments.

```julia
using Oceananigans

grid = RectilinearGrid(size=(10, 10, 10), extent=(100, 100, 100))
g = 9.81
model = ShallowWaterModel(grid, g)
```

--------------------------------

### Initialize and Run Kuramoto-Sivashinsky Simulation

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/developer_docs/model_interface.md

Sets up a simulation for the Kuramoto-Sivashinsky equation. It initializes the grid, the model, perturbs the initial solution, and configures a JLD2 writer to save the solution every time unit.

```julia
using Oceananigans
using Oceananigans.OutputWriters: JLD2Writer
using Oceananigans.Simulations: Simulation

grid = RectilinearGrid(size=256, x=(0, 32π), topology=(Periodic, Flat, Flat), halo=4)
ks_model = KuramotoSivashinskyModel(grid)

# Initialize with a combination of sinusoidal modes
set!(ks_model.solution, x -> cos(x/16) * (1 + sin(x/16)))
simulation = Simulation(ks_model; Δt=0.002, stop_time=60)

simulation.output_writers[:solution] = JLD2Writer(ks_model, (; u=ks_model.solution),
                                                  filename = "ks_solution.jld2",
                                                  schedule = TimeInterval(1),
                                                  overwrite_existing = true)

run!(simulation)
nothing # hide
```

--------------------------------

### ConsecutiveIterations Schedule Example

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/simulations/schedules.md

Triggers a callback for a specified number of iterations after a parent schedule actuates. Useful for averaging callbacks that need data at scheduled times and immediately after.

```julia
Oceananigans.Simulations.reset!(simulation)
simulation.stop_time = 2.5

times = SpecifiedTimes(0.55, 1.5, 2.12)
schedule = ConsecutiveIterations(times)
add_callback!(simulation, dummy, schedule, name=:dummy)
run!(simulation)
```

--------------------------------

### Set Up Simulation Parameters

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/developer_docs/turbulence_closures.md

Defines constants and parameters for the boundary layer simulation, such as domain depth, resolution, stratification, surface stress, and Coriolis parameter.

```julia
using Oceananigans
using Oceananigans.Units

Lz = 256       # Domain depth (m)
Nz = 64        # Vertical resolution
N² = 1e-5      # Background stratification (s⁻²)
τˣ = -1e-4     # Surface kinematic stress (m² s⁻²), i.e. wind stress / density
Jᵇ = 0         # Surface buoyancy flux (m² s⁻³)
f  = 1e-4      # Coriolis parameter (s⁻¹)

stop_time = 2days
nothing # hide
```

--------------------------------

### Initialize Advanced HydrostaticFreeSurfaceModel

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/models_overview.md

Initializes a complex `HydrostaticFreeSurfaceModel` with a spherical grid, custom physics, and boundary conditions. Requires specific packages like `SeawaterPolynomials`.

```julia
using Oceananigans
using SeawaterPolynomials: TEOS10EquationOfState

grid = LatitudeLongitudeGrid(size = (180, 80, 10),
                             longitude = (0, 360),
                             latitude = (-80, 80),
                             z = (-1000, 0),
                             halo = (6, 6, 3))

momentum_advection = WENOVectorInvariant()
coriolis = HydrostaticSphericalCoriolis()
equation_of_state = TEOS10EquationOfState()
buoyancy = SeawaterBuoyancy(; equation_of_state)
closure = CATKEVerticalDiffusivity()

# Generate a zonal wind stress that mimics Earth's mean winds
# with westerlies in mid-latitudes and easterlies near equator and poles
function zonal_wind_stress(λ, φ, t)
    # Parameters
    τ₀ = 1e-4  # Maximum wind stress magnitude (N/m²)
    φ₀ = 30   # Latitude of maximum westerlies (degrees)
    dφ = 10

    # Approximate wind stress pattern
    return - τ₀ * (+ exp(-(φ - φ₀)^2 / 2dφ^2)
                   - exp(-(φ + φ₀)^2 / 2dφ^2)
                   - 0.3 * exp(-φ^2 / dφ^2))
end

u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(zonal_wind_stress))

model = HydrostaticFreeSurfaceModel(grid; momentum_advection, coriolis, closure, buoyancy,
                                    boundary_conditions = (; u=u_bcs), tracers=(:T, :S))
```

--------------------------------

### Forced Free-Slip Flow Streamfunction and Velocity

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/appendix/convergence_tests.md

Defines the streamfunction and velocity components for a forced flow with free-slip conditions at y=0 and y=pi. This setup is used for convergence tests.

```math
\psi(x, y, t) = - \cos [x - \xi(t)] \sin (y) \quad \text{with} \quad g(y) = \sin y
```

```math
u = \cos (x - \xi) \cos y \quad \text{and} \quad v = \sin (x - \xi) \sin y
```

--------------------------------

### Initialize SeawaterBuoyancy with an idealized nonlinear equation of state

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/buoyancy_and_equation_of_state.md

Use an idealized nonlinear equation of state from `SeawaterPolynomials.jl` by first creating the desired polynomial (e.g., `RoquetEquationOfState`) and then passing it to `SeawaterBuoyancy`.

```jldoctest
julia> using SeawaterPolynomials.SecondOrderSeawaterPolynomials

julia> eos = RoquetEquationOfState(:Freezing)
BoussinesqEquationOfState{Float64}:
├── seawater_polynomial: SecondOrderSeawaterPolynomial{Float64}
└── reference_density: 1024.6

julia> eos.seawater_polynomial # the density anomaly
ρ' = 0.7718 Sᴬ - 0.0491 Θ - 0.005027 Θ² - 2.5681e-5 Θ Z + 0.0 Sᴬ² + 0.0 Sᴬ Z + 0.0 Sᴬ Θ

julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: BoussinesqEquationOfState{Float64}
```

--------------------------------

### Get z-nodes at Center location

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/fields.md

Retrieves the z-coordinates of the nodes (cell centers) of the primary mesh, excluding halo cells. This shows the centers of the active grid cells.

```jldoctest
znodes(grid, Center())

# output
4-element view(::Vector{Float64}, 2:5) with eltype Float64:
 0.05
 0.2
 0.44999999999999996
 0.8

```

--------------------------------

### Initialize NonhydrostaticModel with SeawaterBuoyancy

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/models/buoyancy_and_equation_of_state.md

Use `SeawaterBuoyancy()` for a linear equation of state with standard gravity. Requires `:T` and `:S` as tracers.

```jldoctest
julia> model = NonhydrostaticModel(grid; buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection()
└── coriolis: Nothing
```

--------------------------------

### Instantiate HydrostaticFreeSurfaceModel with QuasiAdamsBashforth2

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/numerical_implementation/time_stepping.md

Demonstrates how to manually construct and specify the QuasiAdamsBashforth2 time stepper for a HydrostaticFreeSurfaceModel. The \(\\chi\) parameter can be customized, with 0.125 being an optimal value suggested by Ascher et al. (1995).

```jldoctest
using Oceananigans
using Oceananigans.TimeSteppers: QuasiAdamsBashforth2TimeStepper

# Use QuasiAdamsBashforth2 time stepper (default)
grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
timestepper = QuasiAdamsBashforth2TimeStepper(χ = 0.125)
model = HydrostaticFreeSurfaceModel(grid; timestepper)
model.timestepper


```

--------------------------------

### Creating an Exponential Discretization

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/grids.md

Constructs a coordinate with 10 cells spanning the range [-700, 300] using the default right-biased ExponentialDiscretization. This is a common setup for oceanographic vertical coordinates.

```julia
using Oceananigans

N = 10
l = -700
r = 300

x = ExponentialDiscretization(N, l, r)
```

--------------------------------

### Modify Velocity Tendency with Callback

Source: https://github.com/clima/oceananigans.jl/blob/main/docs/src/simulations/callbacks.md

This example shows how to use a callback with `TendencyCallsite` to directly modify a model's tendency field, in this case, adding a value to the 'u' velocity tendency.

```julia
using Oceananigans
using Oceananigans: TendencyCallsite

grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
model = NonhydrostaticModel(grid)
simulation = Simulation(model, Δt=1, stop_iteration=10)

function modify_tendency!(model, params)
    Gⁿc = model.timestepper.Gⁿ[params.c]
    parent(Gⁿc) .+= params.δ
    return nothing
end

add_callback!(simulation,
              modify_tendency!,
              callsite = TendencyCallsite(),
              parameters = (c = :u, δ = 1))


run!(simulation)
```