Mode-1 internal wave reflection

This example simulates the propagation of a mode-1 internal wave using the ClimateMachine.Ocean subcomponent to solve the hydrostatic Boussinesq equations.

First we ClimateMachine.init().

using ClimateMachine

ClimateMachine.init()

Domain setup

We formulate a non-dimension problem in a Cartesian domain with oceanic anisotropy,

using ClimateMachine.CartesianDomains

domain = RectangularDomain(
    Ne = (32, 1, 4),
    Np = 4,
    x = (-128, 128),
    y = (-128, 128),
    z = (-1, 0),
    periodicity = (false, false, false),
)

RectangularDomain{Float64}
    Np = 4, Ne = (x = 32, y = 1, z = 4)
    L = (x = 256.0, y = 256.0, z = 1.0)

Parameters

We choose parameters appropriate for a hydrostatic internal wave,

Non-dimensional internal wave parameters

f = 1  # Coriolis
N = 10 # Buoyancy frequency

10

Note that the validity of the hydrostatic approximation requires small aspect ratio motions with $k / m \\ll 1$. The hydrostatic dispersion relation for inertia-gravity waves then implies that

λ = 8      # horizontal wave-length
k = 2π / λ # horizontal wavenumber
m = π

ω² = f^2 + N^2 * k^2 / m^2 # and

ω = √(ω²)

2.692582403567252

Internal wave initial condition

We impose modest gravitational acceleration to render time-stepping feasible,

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

Planet.grav(::NonDimensionalParameters) = 256.0

we'd like to use θ as a buoyancy variable, which requires setting the thermal expansion coefficient $αᵀ$ to

g = Planet.grav(NonDimensionalParameters())

αᵀ = 1 / g

0.00390625

We then use the "polarization relations" for vertically-standing, horizontally- propagating hydrostatic internal waves to initialze two wave packets. The hydrostatic polarization relations require

\[\begin{gather} (∂_t^2 + f^2) u = - ∂_x ∂_t p ∂_t v = - f u b = ∂_z p \end{gather}\]

Thus given $p = \cos (k x - ω t) \cos (m z)$, we find

δ = domain.L.x / 15
a(x) = 1e-6 * exp(-x^2 / 2 * δ^2)

ũ(x, z, t) = +a(x) * ω * sin(k * x - ω * t) * cos(m * z)
ṽ(x, z, t) = -a(x) * f * cos(k * x - ω * t) * cos(m * z)
θ̃(x, z, t) = -a(x) * m / k * (ω^2 - f^2) * sin(k * x - ω * t) * sin(m * z)

uᵢ(x, y, z) = ũ(x, z, 0)
vᵢ(x, y, z) = ṽ(x, z, 0)
θᵢ(x, y, z) = θ̃(x, z, 0) + N^2 * z

using ClimateMachine.Ocean.OceanProblems: InitialConditions

initial_conditions = InitialConditions(u = uᵢ, v = vᵢ, θ = θᵢ)

ClimateMachine.Ocean.OceanProblems.InitialConditions{typeof(Main.##382.uᵢ),typeof(Main.##382.vᵢ),typeof(Main.##382.θᵢ),typeof(ClimateMachine.Ocean.OceanProblems.resting)}(Main.##382.uᵢ, Main.##382.vᵢ, Main.##382.θᵢ, ClimateMachine.Ocean.OceanProblems.resting)

Model configuration

We choose a time-step that resolves the gravity wave phase speed,

time_step = 0.005 # close to Δx / c = 0.5 * 1/16, where Δx is nominal resolution

0.005

and build a model with a smidgeon of viscosity and diffusion,

using ClimateMachine.Ocean: HydrostaticBoussinesqSuperModel

model = HydrostaticBoussinesqSuperModel(
    domain = domain,
    time_step = time_step,
    initial_conditions = initial_conditions,
    parameters = NonDimensionalParameters(),
    turbulence_closure = (νʰ = 1e-6, νᶻ = 1e-6, κʰ = 1e-6, κᶻ = 1e-6),
    coriolis = (f₀ = f, β = 0),
    buoyancy = (αᵀ = αᵀ,),
    boundary_tags = ((1, 1), (1, 1), (1, 2)),
);
nothing

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

Fetching data for an animation

To animate the ClimateMachine.Ocean solution, we assemble and cache the horizontal velocity $u$ at periodic intervals:

using ClimateMachine.Ocean: current_time
using ClimateMachine.CartesianFields: assemble
using ClimateMachine.GenericCallbacks: EveryXSimulationTime

fetched_states = []
fetch_every = 0.2 * 2π / ω # time

data_fetcher = EveryXSimulationTime(fetch_every) do
    push!(
        fetched_states,
        (
            u = assemble(model.fields.u.elements),
            θ = assemble(model.fields.θ.elements),
            η = assemble(model.fields.η.elements),
            time = current_time(model),
        ),
    )
    return nothing
end

ClimateMachine.GenericCallbacks.EveryXSimulationTime(Main.##382.var"#1#2"(), 0.46670328817831874, 0)

We also build a callback to log the progress of our simulation,

using Printf
using ClimateMachine.GenericCallbacks: EveryXSimulationSteps
using ClimateMachine.Ocean: current_time, current_step, Δt

print_every = 100 # iterations
wall_clock = [time_ns()]

tiny_progress_printer = EveryXSimulationSteps(print_every) do

    @info(@sprintf(
        "Steps: %d, time: %.2f, Δt: %.2f, max(|u|): %.2e, elapsed time: %.2f secs",
        current_step(model),
        current_time(model),
        Δt(model),
        maximum(abs, model.fields.u),
        1e-9 * (time_ns() - wall_clock[1])
    ))

    wall_clock[1] = time_ns()
end

ClimateMachine.GenericCallbacks.EveryXSimulationSteps(Main.##382.var"#3#4"(), 100, 0)

Running the simulation and animating the results

We're ready to launch.

model.solver_configuration.timeend = 6 * 2π / ω
# model.solver.dt = 0.05 # make this work

@info """ Simulating a hydrostatic Gaussian wave packet with parameters

    f (Coriolis parameter):       $f
    N (buoyancy frequency):       $N
    Internal wave frequency:      $(abs(ω))
    Surface wave frequency:       $(k * sqrt(g * domain.L.z))
    Surface wave group velocity:  $(sqrt(g * domain.L.z))
    Internal wave group velocity: $(N^2 * k / (ω * m))
    Domain width:                 $(domain.L.x)
    Domain height:                $(domain.L.z)

"""

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

0.9999575475529562

Animating the result

We first analye the results to generate plotting limits and contour levels

ηmax = maximum([maximum(abs, state.η.data) for state in fetched_states])
umax = maximum([maximum(abs, state.u.data) for state in fetched_states])

ηlim = (-ηmax, ηmax)
ulim = (-umax, umax)
ulevels = range(ulim[1], ulim[2], length = 31)

-2.2110987101884194e-7:1.4740658067922796e-8:2.2110987101884194e-7

and then animate both fields in a loop,

using Plots

animation = @animate for (i, state) in enumerate(fetched_states)
    @info "Plotting frame $i of $(length(fetched_states))..."

    η_plot = plot(
        state.u.x[:, 1, 1],
        state.η.data[:, 1, 1],
        ylim = ηlim,
        label = nothing,
        title = @sprintf("η at t = %.2f", state.time),
    )

    u_plot = contourf(
        state.u.x[:, 1, 1],
        state.u.z[1, 1, :],
        clamp.(state.u.data[:, 1, :], ulim[1], ulim[2])';
        aspectratio = 64,
        linewidth = 0,
        xlim = domain.x,
        ylim = domain.z,
        xlabel = "x",
        ylabel = "z",
        color = :balance,
        colorbar = false,
        clim = ulim,
        levels = ulevels,
        title = @sprintf("u at t = %.2f", state.time),
    )

    plot(
        η_plot,
        u_plot,
        layout = Plots.grid(2, 1, heights = (0.3, 0.7)),
        link = :x,
        size = (600, 300),
    )
end

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

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

This page was generated using Literate.jl.