Geostrophic adjustment in the hydrostatic Boussinesq equations

This example simulates a one-dimensional geostrophic adjustement problem using the ClimateMachine.Ocean subcomponent to solve the hydrostatic Boussinesq equations.

First we ClimateMachine.init().

using ClimateMachine

ClimateMachine.init()

Domain setup

We formulate our problem in a Cartesian domain 100 km in $x, y$ and 400 m deep, and discretized on a grid with 100 fourth-order elements in $x$, and 1 fourth-order element in $y, z$,

using ClimateMachine.CartesianDomains

domain = RectangularDomain(
    Ne = (25, 1, 1),
    Np = 4,
    x = (0, 1e6),
    y = (0, 1e6),
    z = (-400, 0),
    periodicity = (false, true, false),
)

RectangularDomain{Float64}
    Np = 4, Ne = (x = 25, y = 1, z = 1)
    L = (x = 1.0e6, y = 1.0e6, z = 400.0)

Physical parameters

We use a Coriolis parameter appropriate for mid-latitudes,

f = 1e-4 # s⁻¹, Coriolis parameter
nothing # hide

and Earth's gravitational acceleration,

using CLIMAParameters: AbstractEarthParameterSet, Planet
struct EarthParameters <: AbstractEarthParameterSet end

g = Planet.grav(EarthParameters()) # m s⁻²

9.81

An unbalanced initial state

We use a Gaussian, partially-balanced initial condition with parameters

U = 0.1              # geostrophic velocity (m s⁻¹)
L = domain.L.x / 40  # Gaussian width (m)
a = f * U * L / g    # amplitude of the geostrophic surface displacement (m)
x₀ = domain.L.x / 4  # Gaussian origin (m, recall that x ∈ [0, Lx])

250000.0

and functional form

Gaussian(x, L) = exp(-x^2 / (2 * L^2))

# Geostrophic ``y``-velocity: f V = g ∂_x η
vᵍ(x, y, z) = -U * (x - x₀) / L * Gaussian(x - x₀, L)

# Geostrophic surface displacement
ηᵍ(x, y, z) = a * Gaussian(x - x₀, L)

ηᵍ (generic function with 1 method)

We double the initial surface displacement so that the surface is half-balanced, half unbalanced,

ηⁱ(x, y, z) = 2 * ηᵍ(x, y, z)

ηⁱ (generic function with 1 method)

In summary,

using ClimateMachine.Ocean.OceanProblems: InitialConditions

initial_conditions = InitialConditions(v = vᵍ, η = ηⁱ)

@info """ Parameters for the Geostrophic adjustment problem are...

    Coriolis parameter:                            $f s⁻¹
    Gravitational acceleration:                    $g m s⁻²
    Geostrophic velocity:                          $U m s⁻¹
    Width of the initial geostrophic perturbation: $L m
    Amplitude of the initial surface perturbation: $a m
    Rossby number (U / f L):                       $(U / (f * L))

"""

┌ Info:  Parameters for the Geostrophic adjustment problem are...
│ 
│     Coriolis parameter:                            0.0001 s⁻¹
│     Gravitational acceleration:                    9.81 m s⁻²
│     Geostrophic velocity:                          0.1 m s⁻¹
│     Width of the initial geostrophic perturbation: 25000.0 m
│     Amplitude of the initial surface perturbation: 0.0254841997961264 m
│     Rossby number (U / f L):                       0.04
│ 
└ @ Main.##382 string:5

Boundary conditions, Driver configuration, and Solver configuration

Next, we configure the HydrostaticBoussinesqModel and build the DriverConfiguration. We configure our model in a domain which is bounded in the $x$ direction. Both the boundary conditions in $x$ and in $z$ require boundary conditions, which we define:

using ClimateMachine.Ocean:
    Impenetrable, Penetrable, FreeSlip, Insulating, OceanBC

solid_surface_boundary_conditions = OceanBC(
    Impenetrable(FreeSlip()), # Velocity boundary conditions
    Insulating(),             # Temperature boundary conditions
)

free_surface_boundary_conditions = OceanBC(
    Penetrable(FreeSlip()),   # Velocity boundary conditions
    Insulating(),             # Temperature boundary conditions
)

boundary_conditions =
    (solid_surface_boundary_conditions, free_surface_boundary_conditions)

(ClimateMachine.Ocean.OceanBC{ClimateMachine.Ocean.Impenetrable{ClimateMachine.Ocean.FreeSlip},ClimateMachine.Ocean.Insulating}(ClimateMachine.Ocean.Impenetrable{ClimateMachine.Ocean.FreeSlip}(ClimateMachine.Ocean.FreeSlip()), ClimateMachine.Ocean.Insulating()), ClimateMachine.Ocean.OceanBC{ClimateMachine.Ocean.Penetrable{ClimateMachine.Ocean.FreeSlip},ClimateMachine.Ocean.Insulating}(ClimateMachine.Ocean.Penetrable{ClimateMachine.Ocean.FreeSlip}(ClimateMachine.Ocean.FreeSlip()), ClimateMachine.Ocean.Insulating()))

We refer to these boundary conditions by their indices in the boundary_tags tuple when specifying the boundary conditions for the state; in other words, "1" corresponds to solid_surface_boundary_conditions, while 2 corresponds to free_surface_boundary_conditions,

boundary_tags = (
    (1, 1), # (west, east) boundary conditions
    (0, 0), # (south, north) boundary conditions
    (1, 2), # (bottom, top) boundary conditions
)

((1, 1), (0, 0), (1, 2))

We're now ready to build the model.

using ClimateMachine.Ocean

model = Ocean.HydrostaticBoussinesqSuperModel(
    domain = domain,
    time_step = 2.0,
    initial_conditions = initial_conditions,
    parameters = EarthParameters(),
    turbulence_closure = (νʰ = 0, κʰ = 0, νᶻ = 0, κᶻ = 0),
    coriolis = (f₀ = f, β = 0),
    boundary_tags = boundary_tags,
    boundary_conditions = boundary_conditions,
);
nothing

┌ Info: Initializing 
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/solver_configs.jl:185

Animating the solution

To animate the ClimateMachine.Ocean solution, we'll create a callback that draws a plot and stores it in an array. When the simulation is finished, we'll string together the plotted frames into an animation.

using Printf
using Plots
using ClimateMachine.GenericCallbacks: EveryXSimulationSteps
using ClimateMachine.CartesianFields: assemble
using ClimateMachine.Ocean: current_step, current_time

u, v, η, θ = model.fields

# Container to hold the plotted frames
movie_plots = []

plot_every = 200 # iterations

plot_maker = EveryXSimulationSteps(plot_every) do

    @info "Steps: $(current_step(model)), time: $(current_time(model))"

    assembled_u = assemble(u.elements)
    assembled_v = assemble(v.elements)
    assembled_η = assemble(η.elements)

    umax = 0.5 * max(maximum(abs, u), maximum(abs, v))
    ulim = (-umax, umax)

    u_plot = plot(
        assembled_u.x[:, 1, 1],
        [assembled_u.data[:, 1, 1] assembled_v.data[:, 1, 1]],
        xlim = domain.x,
        ylim = (-0.7U, 0.7U),
        label = ["u" "v"],
        linewidth = 2,
        xlabel = "x (m)",
        ylabel = "Velocities (m s⁻¹)",
    )

    η_plot = plot(
        assembled_η.x[:, 1, 1],
        assembled_η.data[:, 1, 1],
        xlim = domain.x,
        ylim = (-0.01a, 1.2a),
        linewidth = 2,
        label = nothing,
        xlabel = "x (m)",
        ylabel = "η (m)",
    )

    push!(movie_plots, (u = u_plot, η = η_plot, time = current_time(model)))

    return nothing
end

ClimateMachine.GenericCallbacks.EveryXSimulationSteps(Main.##382.var"#1#2"(), 200, 0)

Running the simulation and animating the results

Finally, we run the simulation,

hours = 3600.0

model.solver_configuration.timeend = 2hours

result = ClimateMachine.invoke!(
    model.solver_configuration;
    user_callbacks = [plot_maker],
)

0.9092270436985049

and animate the results,

animation = @animate for p in movie_plots
    title = @sprintf("Geostrophic adjustment at t = %.2f hours", p.time / hours)
    frame = plot(
        p.u,
        p.η,
        layout = (2, 1),
        size = (800, 600),
        title = [title ""],
    )
end

gif(animation, "geostrophic_adjustment.mp4", fps = 5) # hide

Plots.AnimatedGif("/central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/docs/src/generated/Ocean/geostrophic_adjustment.mp4")

This page was generated using Literate.jl.