[Dycore] Implement tendency term for moist static energy following Pauluis (2008)

[glwagner]

Yet another orphaned change from #29; if we want to evolve moist static energy, we need to include the buoyancy flux term in the energy tendency. See equation 15 in Pauluis 2008.

@kaiyuan-cheng probably worth revisiting the validation tests on this branch.

[glwagner]

The thermal bubble example needs work. But it does run.

[glwagner]

@kaiyuan-cheng this is the necessary term

I think a free convection example is sufficient to reveal the necessity of this.

[codecov]

Codecov Report

❌ Patch coverage is 78.20513% with 17 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/AtmosphereModels/set_atmosphere_model.jl	33.33%	4 Missing ⚠️
.../AtmosphereModels/update_atmosphere_model_state.jl	75.00%	4 Missing ⚠️
src/AtmosphereModels/anelastic_formulation.jl	90.00%	2 Missing ⚠️
src/AtmosphereModels/dynamics_kernel_functions.jl	92.00%	2 Missing ⚠️
src/AtmosphereModels/saturation_adjustment.jl	0.00%	2 Missing ⚠️
src/Thermodynamics/reference_states.jl	50.00%	2 Missing ⚠️
src/AtmosphereModels/atmosphere_model.jl	66.66%	1 Missing ⚠️

📢 Thoughts on this report? Let us know!

[glwagner]

For documentation, I'm running this simple free convection case:

using Breeze
using Printf
using GLMakie

Nx = Nz = 128
Lz = 4 * 1024
grid = RectilinearGrid(size=(Nx, Nz), x=(0, 2Lz), z=(0, Lz), topology=(Periodic, Flat, Bounded))

e_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1000))
model = AtmosphereModel(grid, advection=WENO(), boundary_conditions=(; e=e_bcs))

Tₛ = model.formulation.constants.reference_potential_temperature
g = model.thermodynamics.gravitational_acceleration
N² = 1e-6
dθdz = N² * Tₛ / g  # Background potential temperature gradient
θᵢ(x, z) = Tₛ + dθdz * z # background stratification
Ξᵢ(x, z) = 1e-2 * randn()
set!(model, θ=θᵢ, u=Ξᵢ, v=Ξᵢ)

simulation = Simulation(model, Δt=10, stop_iteration=500)
progress(sim) = @info string(iteration(sim), ": t = ", time(sim))
add_callback!(simulation, progress, IterationInterval(100))
run!(simulation)
heatmap(model.velocities.w)

onmain this produces

this is wrong because we should not get convection/flow developing at the top of the domain. Here we are heating from the bottom.

on this PR we get

Unfortunately while this is a bit improved (less convection at the top) it is still not right, so we have a little work to do still.

[glwagner]

err I realized my example above is not doing what I thought because the boundary conditions were on a non-existing variable (we evolve ρe not e). So we need

ρe_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1000))
model = AtmosphereModel(grid, advection=WENO(), boundary_conditions=(; ρe=ρe_bcs))

iterating...

[glwagner]

okkkk... now with

using Breeze
using Oceananigans.Units
using Printf
using GLMakie

Nx = Nz = 64
Δx, Δz = 100, 50
Lx, Lz = Δx * Nx, Δz * Nz
grid = RectilinearGrid(size=(Nx, Nz), x=(0, Lx), z=(0, Lz), topology=(Periodic, Flat, Bounded))

ρe_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1000))
model = AtmosphereModel(grid, advection=WENO(), boundary_conditions=(; ρe=ρe_bcs))

Tₛ = model.formulation.constants.reference_potential_temperature
g = model.thermodynamics.gravitational_acceleration
N² = 1e-6
dθdz = N² * Tₛ / g  # Background potential temperature gradient
θᵢ(x, z) = Tₛ + dθdz * z # background stratification
Ξᵢ(x, z) = 1e-2 * randn()
set!(model, θ=θᵢ, u=Ξᵢ, v=Ξᵢ)

simulation = Simulation(model, Δt=10, stop_time=20minutes)
conjure_time_step_wizard!(simulation, cfl=0.7)
∫ρe = Field(Integral(model.energy))

function progress(sim)
    compute!(∫ρe)
    max_w = maximum(abs, model.velocities.w)

    @info @sprintf("%s: t = %s, ∫ρe dV = %.9e J/kg, max|w| = %.2e m/s",
                   iteration(sim), prettytime(sim), first(∫ρe), max_w)

    return nothing
end

add_callback!(simulation, progress, IterationInterval(5))

∫ρe_series = Float64[]
w_series = []
t = []

function record_data!(sim)
    compute!(∫ρe)
    push!(∫ρe_series, first(∫ρe))
    push!(w_series, deepcopy(model.velocities.w))
    push!(t, sim.model.clock.time)
    return nothing
end
add_callback!(simulation, record_data!)

run!(simulation)

fig = Figure(size=(800, 550), fontsize=12)

Nt = length(∫ρe_series)
axw = Axis(fig[1, 1], xlabel="x (m)", ylabel="z (m)", aspect=DataAspect())
axe = Axis(fig[2, 1], xlabel="Time (s)", ylabel="∫ρe dV (J/kg)", height=Relative(0.2))
slider = Slider(fig[3, 1], range=1:Nt, startvalue=1)
n = slider.value

∫ρe_n = @lift ∫ρe_series[$n] / ∫ρe_series[1]
t_n = @lift t[$n]
w_n = @lift w_series[$n]
w_lim = maximum(abs, model.velocities.w)

lines!(axe, t, ∫ρe_series / ∫ρe_series[1])
scatter!(axe, t_n, ∫ρe_n, markersize=10)

hm = heatmap!(axw, w_n, colorrange=(-w_lim, w_lim), colormap=:balance)
Colorbar(fig[1, 2], hm, label="w (m/s)")

rowgap!(fig.layout, 1, Relative(-0.2))
rowgap!(fig.layout, 2, Relative(-0.2))

display(fig)

we get

I will push this script for future reference.

I have a theory for causes the behavior in the no flux, freely decyaing case that was first posted. I believe this has to do with how we treat pressure / buoyancy. There is an extremely weak instability that occurs (something we have seen over on Oceananigans). I'd like to look into dealing with this in a future PR.

[glwagner]

I also started some simple dynamics docs on this PR.

[glwagner]

thank you @giordano <3

[giordano]

Was it intended to change this from an example which shows the output to a simple code block? This makes building the docs a lot faster, but not sure this was done intentionally.

[glwagner]

that was a mistake while prototyping the docs, shoudl not have been committed

[giordano]

Yeah, that was why I had 6d0c930 in #77. This makes me think that maybe it'd be nice if Documenter had a "draft" feature which doesn't expand @example and doesn't run doctests and turns all errors into warnings, just to get a quick visual feedback, without doing a full build.

[giordano]

Ah, it actually does already! deploydocs(; draft=true)