Geostrophic adjustment using Oceananigans.HydrostaticFreeSurfaceModel
This example demonstrates how to simulate the one-dimensional geostrophic adjustment of a free surface using Oceananigans.HydrostaticFreeSurfaceModel. Here, we solve the hydrostatic Boussinesq equations beneath a free surface with a small-amplitude about rest $z = 0$, with boundary conditions expanded around $z = 0$, and free surface dynamics linearized under the assumption $η / H \ll 1$, where $η$ is the free surface displacement, and $H$ is the total depth of the fluid.
Install dependencies
First let's make sure we have all required packages installed.
using Pkg
pkg"add Oceananigans, JLD2, Plots"The HydrostaticFreeSurfaceModel is still "experimental". This means some of its features are not exported, such as the ImplicitFreeSurface, and must be brought into scope manually:
using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ImplicitFreeSurfaceA one-dimensional domain
We use a one-dimensional domain of geophysical proportions,
grid = RectilinearGrid(size = (128, 1),
x = (0, 1000kilometers), z = (-400meters, 0),
topology = (Bounded, Flat, Bounded))128×1×1 RectilinearGrid{Float64, Bounded, Flat, Bounded} on CPU with 1×0×1 halo
├── Bounded x ∈ [0.0, 1.0e6] regularly spaced with Δx=7812.5
├── Flat y
└── Bounded z ∈ [-400.0, 0.0] regularly spaced with Δz=400.0We always have to include the z-direction for HydrostaticFreeSurfaceModel, even if the model is "barotropic" with just one grid point in z.
We deploy a Coriolis parameter appropriate for the mid-latitudes on Earth,
coriolis = FPlane(f=1e-4)FPlane{Float64}(f=0.0001)Building a HydrostaticFreeSurfaceModel
We use grid and coriolis to build a simple HydrostaticFreeSurfaceModel,
model = HydrostaticFreeSurfaceModel(; grid, coriolis, free_surface = ImplicitFreeSurface())HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×1 RectilinearGrid{Float64, Bounded, Flat, Bounded} on CPU with 1×0×1 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with -ĝ = ZDirection
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│ └── solver: FFTImplicitFreeSurfaceSolver
└── coriolis: FPlane{Float64}A geostrophic adjustment initial value problem
We pose a geostrophic adjustment problem that consists of a partially-geostrophic Gaussian height field complemented by a geostrophic $y$-velocity,
gaussian(x, L) = exp(-x^2 / 2L^2)
U = 0.1 # geostrophic velocity
L = grid.Lx / 40 # Gaussian width
x₀ = grid.Lx / 4 # Gaussian center
vᵍ(x, y, z) = - U * (x - x₀) / L * gaussian(x - x₀, L)
g = model.free_surface.gravitational_acceleration
η₀ = coriolis.f * U * L / g # geostrohpic free surface amplitude
ηᵍ(x) = η₀ * gaussian(x - x₀, L)ηᵍ (generic function with 1 method)We use an initial height field that's twice the geostrophic solution, thus superimposing a geostrophic and ageostrophic component in the free surface displacement field:
ηⁱ(x, y) = 2 * ηᵍ(x)ηⁱ (generic function with 1 method)We set the initial condition to $vᵍ$ and $ηⁱ$,
set!(model, v=vᵍ, η=ηⁱ)Running a Simulation
We pick a time-step that resolves the surface dynamics,
gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed
wave_propagation_time_scale = model.grid.Δxᶜᵃᵃ / gravity_wave_speed
Δt = 0.1wave_propagation_time_scale
simulation = Simulation(model; Δt, stop_iteration = 1000)Simulation of HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 12.474 seconds
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN years
├── Stop time: Inf years
├── Stop iteration : 1000.0
├── Wall time limit: Inf
├── 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
└── Diagnostics: OrderedDict with no entriesOutput
We output the velocity field and free surface displacement,
output_fields = merge(model.velocities, (η=model.free_surface.η,))
simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
schedule = IterationInterval(10),
prefix = "geostrophic_adjustment",
force = true)
run!(simulation)[ Info: Initializing simulation...
[ Info: ... simulation initialization complete (2.212 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (1.210 minutes).
[ Info: Simulation is stopping. Model iteration 1000 has hit or exceeded simulation stop iteration 1000.Visualizing the results
using Oceananigans.OutputReaders: FieldTimeSeries
using Plots, Printf
u_timeseries = FieldTimeSeries("geostrophic_adjustment.jld2", "u")
v_timeseries = FieldTimeSeries("geostrophic_adjustment.jld2", "v")
η_timeseries = FieldTimeSeries("geostrophic_adjustment.jld2", "η")
xη = xw = xv = xnodes(v_timeseries)
xu = xnodes(u_timeseries)
t = u_timeseries.times
anim = @animate for i = 1:length(t)
u = interior(u_timeseries[i], :, 1, 1)
v = interior(v_timeseries[i], :, 1, 1)
η = interior(η_timeseries[i], :, 1, 1)
titlestr = @sprintf("Geostrophic adjustment at t = %.1f hours", t[i] / hours)
label = ""
linewidth = 2
u_plot = plot(xu / kilometers, u; linewidth, label,
xlabel = "x (km)", ylabel = "u (m s⁻¹)", ylims = (-2e-3, 2e-3))
v_plot = plot(xv / kilometers, v; linewidth, label, title = titlestr,
xlabel = "x (km)", ylabel = "v (m s⁻¹)", ylims = (-U, U))
η_plot = plot(xη / kilometers, η; linewidth, label,
xlabel = "x (km)", ylabel = "η (m)", ylims = (-η₀/10, 2η₀))
plot(v_plot, u_plot, η_plot, layout = (3, 1), size = (800, 600))
endThis page was generated using Literate.jl.