Thermal bubble

This script sets up, runs, and visualizes a "thermal bubble" (just a circular region of warm air) rising through a stably-stratified background.

using Breeze
using Oceananigans.Units
using Statistics
using Printf
using CairoMakie

A simple model on a RectilinearGrid

grid = RectilinearGrid(CPU(); size = (128, 128), halo = (5, 5),
                       x = (-10e3, 10e3), z = (0, 10e3),
                       topology = (Periodic, Flat, Bounded))

advection = WENO(order=9)
model = AtmosphereModel(grid; advection)

AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── formulation: AnelasticFormulation(p₀=101325.0, θ₀=288.0)
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: WENO{5, Float64, Float32}(order=9)
├── tracers: ()
├── coriolis: Nothing
└── microphysics: Nothing

Moist static energy perturbation

We add a localized potential temperature perturbation that translates into a moist static energy anomaly.

r₀ = 2e3
Δθ = 10 # K
N² = 1e-6
θ₀ = model.formulation.reference_state.potential_temperature
g = model.thermodynamics.gravitational_acceleration
dθdz = N² * θ₀ / g

function θᵢ(x, z; x₀=mean(xnodes(grid, Center())), z₀=0.3*grid.Lz)
    θ̄ = θ₀ + dθdz * z
    r = sqrt((x - x₀)^2 + (z - z₀)^2)
    θ′ = Δθ * max(0, 1 - r / r₀)
    return θ̄ + θ′
end

set!(model, θ = θᵢ)

ρe = model.energy_density
ρE = Field(Average(ρe, dims=1))
ρe′ = Field(ρe - ρE)

128×1×128 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: BinaryOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 138×1×138 OffsetArray(::Array{Float64, 3}, -4:133, 1:1, -4:133) with eltype Float64 with indices -4:133×1:1×-4:133
    └── max=7277.39, min=-849.551, mean=6.75016e-13

Initial energy perturbation visualization

Plot the initial moist static energy perturbation to ensure the bubble looks as expected.

fig = Figure()
ax = Axis(fig[1, 1], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Initial energy perturbation ρe′ (J / kg)")
hm = heatmap!(ax, ρe′)
Colorbar(fig[1, 2], hm, label = "ρe′ (J/kg)")
fig

Simulation rising

simulation = Simulation(model; Δt=2, stop_time=15minutes)
conjure_time_step_wizard!(simulation, cfl=0.7)

function progress(sim)
    ρe = sim.model.energy_density
    u, v, w = sim.model.velocities

    msg = @sprintf("Iter: %d, t: %s, Δt: %s, extrema(ρe): (%.2f, %.2f) J/kg, max|u|: %.2f m/s, max|w|: %.2f m/s",
                   iteration(sim), prettytime(sim), prettytime(sim.Δt),
                   minimum(ρe), maximum(ρe),
                   maximum(abs, u), maximum(abs, w))

    @info msg
    return nothing
end

add_callback!(simulation, progress, TimeInterval(1minute))

u, v, w = model.velocities
T = model.temperature

outputs = merge(model.velocities, model.tracers, (; ρe′, ρe, T))

filename = "thermal_bubble.jld2"
writer = JLD2Writer(model, outputs; filename,
                    schedule = TimeInterval(10seconds),
                    overwrite_existing = true)

simulation.output_writers[:jld2] = writer

run!(simulation)

[ Info: Initializing simulation...
[ Info: Iter: 0, t: 0 seconds, Δt: 2.200 seconds, extrema(ρe): (126698.98, 353638.32) J/kg, max|u|: 0.00 m/s, max|w|: 0.00 m/s
[ Info:     ... simulation initialization complete (10.017 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.202 seconds).
[ Info: Iter: 28, t: 1 minute, Δt: 2.662 seconds, extrema(ρe): (126699.02, 353636.65) J/kg, max|u|: 3.22 m/s, max|w|: 8.17 m/s
[ Info: Iter: 52, t: 2 minutes, Δt: 3.543 seconds, extrema(ρe): (126699.18, 353630.57) J/kg, max|u|: 6.37 m/s, max|w|: 14.79 m/s
[ Info: Iter: 73, t: 3 minutes, Δt: 3.103 seconds, extrema(ρe): (126699.46, 353620.44) J/kg, max|u|: 9.25 m/s, max|w|: 17.61 m/s
[ Info: Iter: 98, t: 4 minutes, Δt: 3.040 seconds, extrema(ρe): (126699.83, 353605.91) J/kg, max|u|: 11.61 m/s, max|w|: 17.68 m/s
[ Info: Iter: 123, t: 5 minutes, Δt: 3.261 seconds, extrema(ρe): (126700.70, 353587.35) J/kg, max|u|: 13.02 m/s, max|w|: 16.43 m/s
[ Info: Iter: 144, t: 6 minutes, Δt: 3.294 seconds, extrema(ρe): (126703.37, 353565.57) J/kg, max|u|: 13.42 m/s, max|w|: 13.75 m/s
[ Info: Iter: 168, t: 7 minutes, Δt: 3.327 seconds, extrema(ρe): (126710.01, 353541.57) J/kg, max|u|: 12.15 m/s, max|w|: 13.45 m/s
[ Info: Iter: 188, t: 8 minutes, Δt: 3.744 seconds, extrema(ρe): (126723.41, 353516.38) J/kg, max|u|: 12.23 m/s, max|w|: 12.40 m/s
[ Info: Iter: 208, t: 9 minutes, Δt: 3.332 seconds, extrema(ρe): (126740.76, 353491.12) J/kg, max|u|: 12.90 m/s, max|w|: 13.72 m/s
[ Info: Iter: 232, t: 10 minutes, Δt: 3.018 seconds, extrema(ρe): (126746.88, 353466.95) J/kg, max|u|: 14.15 m/s, max|w|: 15.67 m/s
[ Info: Iter: 256, t: 11.000 minutes, Δt: 3.059 seconds, extrema(ρe): (126757.51, 353447.68) J/kg, max|u|: 15.44 m/s, max|w|: 18.03 m/s
[ Info: Iter: 282, t: 12 minutes, Δt: 2.308 seconds, extrema(ρe): (126778.09, 353439.44) J/kg, max|u|: 14.86 m/s, max|w|: 19.93 m/s
[ Info: Iter: 312, t: 13.000 minutes, Δt: 2.238 seconds, extrema(ρe): (126803.54, 353430.40) J/kg, max|u|: 11.08 m/s, max|w|: 20.41 m/s
[ Info: Iter: 342, t: 14.000 minutes, Δt: 2.568 seconds, extrema(ρe): (126820.62, 353420.35) J/kg, max|u|: 11.96 m/s, max|w|: 19.44 m/s
[ Info: Simulation is stopping after running for 36.215 seconds.
[ Info: Simulation time 15 minutes equals or exceeds stop time 15 minutes.
[ Info: Iter: 366, t: 15 minutes, Δt: 3.031 seconds, extrema(ρe): (126836.05, 353409.38) J/kg, max|u|: 10.77 m/s, max|w|: 17.01 m/s

Visualization

Visualize the moist static energy perturbation and the vertical velocity through time and create an animation.

@info "Creating visualization..."

ρe′t = FieldTimeSeries(filename, "ρe′")
wt = FieldTimeSeries(filename, "w")

times = ρe′t.times
Nt = length(ρe′t)

fig = Figure(size = (800, 800), fontsize = 12)
axρ = Axis(fig[1, 1], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Energy perturbation ρe′ (J / kg)")
axw = Axis(fig[2, 1], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Vertical velocity w (m / s)")

n = Observable(Nt)

ρe′n = @lift ρe′t[$n]
wn = @lift wt[$n]

title = @lift "Thermal bubble evolution — t = $(prettytime(times[$n]))"
fig[0, :] = Label(fig, title, fontsize = 16, tellwidth = false)

ρe′_range = (minimum(ρe′t), maximum(ρe′t))
w_range = maximum(abs, wt)

hmρ = heatmap!(axρ, ρe′n, colorrange = ρe′_range, colormap = :balance)
hmw = heatmap!(axw, wn, colorrange = (-w_range, w_range), colormap = :balance)

Colorbar(fig[1, 2], hmρ, label = "ρe′ (J/kg)", vertical = true)
Colorbar(fig[2, 2], hmw, label = "w (m/s)", vertical = true)

CairoMakie.record(fig, "thermal_bubble.mp4", 1:Nt, framerate = 12) do nn
    n[] = nn
end

[ Info: Creating visualization...

This page was generated using Literate.jl.