Refactor tests; implement anelastic pressure solver tests

docs/src/thermodynamics.md

@@ -207,15 +207,17 @@ ratio of condensed species. In most situations on Earth, ``qᶜ ≪ qᵛ``.
 ```@example thermo
 # Compute mixture properties for air with 0.01 specific humidity
 qᵛ = 0.01 # 1% water vapor by mass
-Rᵐ = mixture_gas_constant(qᵛ, thermo)
+qᵈ = 1 - qᵛ
+Rᵐ = mixture_gas_constant(qᵈ, qᵛ, thermo)
 ```
 
 We likewise define a mixture heat capacity via ``cᵖᵐ = qᵈ cᵖᵈ + qᵛ cᵖᵛ``,
 
 
 ```@example thermo
-q = 0.01 # 1% water vapor by mass
-cᵖᵐ = mixture_heat_capacity(qᵛ, thermo)
+qᵛ = 0.01 # 1% water vapor by mass
+qᵈ = 1 - qᵛ
+cᵖᵐ = mixture_heat_capacity(qᵈ, qᵛ, thermo)
 ```
 
 ## The Clasuius-Claperyon relation and saturation specific humidity

examples/anelastic_free_convection.jl

@@ -0,0 +1,98 @@
+using Breeze
+using Oceananigans.Units
+using Printf
+
+arch = CPU()
+Nx = Nz = 128
+Lz = 10e3
+grid = RectilinearGrid(arch, size=(Nx, Nz), x=(0, 2Lz), z=(0, Lz), topology=(Periodic, Flat, Bounded))
+
+ρe_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1000))
+advection = WENO()
+model = AtmosphereModel(grid; advection, boundary_conditions=(; ρe=ρe_bcs))
+
+Lz = grid.Lz
+Δθ = 2 # K
+Tₛ = model.formulation.constants.reference_potential_temperature # K
+θᵢ(x, z) = Tₛ + Δθ * z / Lz + 1e-2 * Δθ * randn()
+Ξ(x, z) = 1e-2 * randn()
+set!(model, θ=θᵢ, ρu=Ξ, ρv=Ξ, ρw=Ξ)
+
+simulation = Simulation(model, Δt=10, stop_iteration=1000) #0stop_time=4hours)
+
+# TODO make this work
+# conjure_time_step_wizard!(simulation, cfl=0.7)
+
+ρu, ρv, ρw = model.momentum
+δ = Field(∂x(ρu) + ∂y(ρv) + ∂z(ρw))
+
+function progress(sim)
+    compute!(δ)
+    u, v, w = sim.model.velocities
+    T = sim.model.temperature
+    umax = maximum(abs, u)
+    vmax = maximum(abs, v)
+    wmax = maximum(abs, w)
+    δmax = maximum(abs, δ)
+    Tmin = minimum(T)
+    Tmax = maximum(T)
+
+    msg = @sprintf("Iter: %d, t: %s, max|u|: (%.2e, %.2e, %.2e), max|δ|: %.2e, extrema(T): (%.2e, %.2e)",
+                    iteration(sim), prettytime(sim), umax, vmax, wmax, δmax, Tmin, Tmax)
+
+    @info msg
+
+    return nothing
+end
+
+add_callback!(simulation, progress, IterationInterval(10))
+
+ow = JLD2Writer(model, prognostic_fields(model),
+                filename = "free_convection.jld2",
+                schedule = TimeInterval(1minutes),
+                overwrite_existing = true)
+
+simulation.output_writers[:jld2] = ow
+
+run!(simulation)
+
+ρwt = FieldTimeSeries("free_convection.jld2", "ρw")
+ρet = FieldTimeSeries("free_convection.jld2", "ρe")
+times = ρet.times
+Nt = length(ρet)
+
+using GLMakie, Printf
+
+fig = Figure(size=(1200, 600), fontsize=12)
+axw = Axis(fig[1, 1], xlabel="x (m)", ylabel="z (m)")
+axe = Axis(fig[1, 2], xlabel="x (m)", ylabel="z (m)")
+slider = Slider(fig[2, 1:2], range=1:Nt, startvalue=1)
+
+n = slider.value
+ρwn = @lift interior(ρwt[$n], :, 1, :)
+ρen = @lift interior(ρet[$n], :, 1, :)
+title = @lift "t = $(prettytime(times[$n]))"
+
+fig[0, :] = Label(fig, title, fontsize=22, tellwidth=false)
+
+Tmin = minimum(ρet)
+Tmax = maximum(ρet)
+wlim = maximum(abs, ρwt) / 2
+
+hmw = heatmap!(axw, ρwn, colorrange=(-wlim, wlim), colormap=:balance)
+hme = heatmap!(axe, ρen, colorrange=(Tmin, Tmax), colormap=:balance)
+
+# Label(fig[0, 1], "θ", tellwidth=false)
+# Label(fig[0, 2], "q", tellwidth=false)
+# Label(fig[0, 1], "θ", tellwidth=false)
+# Label(fig[0, 2], "q", tellwidth=false)
+
+Colorbar(fig[1, 0], hmw, label = "w", vertical=true, flipaxis=true)
+Colorbar(fig[1, 3], hme, label = "ρe", vertical=true)
+
+fig
+
+record(fig, "free_convection.mp4", 1:Nt, framerate=12) do nn
+    @info "Drawing frame $nn of $Nt..."
+    n[] = nn
+end

examples/free_convection.jl

@@ -14,15 +14,6 @@ p₀ = 101325 # Pa
 reference_constants = Breeze.Thermodynamics.ReferenceStateConstants(base_pressure=p₀, potential_temperature=θ₀)
 buoyancy = Breeze.MoistAirBuoyancy(; reference_constants)
 
-# Simple precipitation scheme from CloudMicrophysics
-using CloudMicrophysics
-using CloudMicrophysics.Microphysics0M: remove_precipitation
-
-FT = eltype(grid)
-microphysics = CloudMicrophysics.Parameters.Parameters0M{FT}(τ_precip=600, S_0=0, qc_0=0.02)
-@inline precipitation(x, z, t, q, params) = remove_precipitation(params, q, 0)
-q_forcing = Forcing(precipitation, field_dependencies=:q, parameters=microphysics)
-
 ρ₀ = Breeze.MoistAirBuoyancies.base_density(buoyancy) # air density at z=0
 cₚ = buoyancy.thermodynamics.dry_air.heat_capacity
 Q₀ = 1000 # heat flux in W / m²

examples/time_step_atmos_model.jl

@@ -42,7 +42,7 @@ set!(model, θ=θᵢ, q=qᵢ, u=Ξᵢ, v=Ξᵢ)
 compute!(δ)
 
 stop_time = 5minutes
-simulation = Simulation(model, Δt=0.1; stop_iteration=1000)
+simulation = Simulation(model, Δt=0.1, stop_iteration=10)
 # conjure_time_step_wizard!(simulation, cfl=0.7)
 
 using Printf
@@ -94,10 +94,7 @@ writer = JLD2Writer(model, outputs; filename,
 
 simulation.output_writers[:jld2] = writer
 
-try
-    run!(simulation)
-catch
-end
+run!(simulation)
 
 using GLMakie
 
@@ -164,4 +161,4 @@ GLMakie.record(fig, "thermal_bubble.mp4", 1:Nt, framerate=10) do nn
     n[] = nn
 end
 @info "Saved animation to thermal_bubble.mp4"
-=#
+=#

src/AtmosphereModels/AtmosphereModels.jl

@@ -1,9 +1,8 @@
 module AtmosphereModels
 
-export AtmosphereModel, AnelasticFormulation
+export AtmosphereModel, prognostic_fields   
 
 include("atmosphere_model.jl")
-include("anelastic_formulation.jl")
 include("saturation_adjustment.jl")
 include("update_hydrostatic_pressure.jl")
 include("update_atmosphere_model_state.jl")