Refactor tests; implement anelastic pressure solver tests

Project.toml

@@ -18,7 +18,7 @@ Adapt = "4.3.0"
 AtmosphericProfilesLibrary = "0.1.7"
 CloudMicrophysics = "0.22.13"
 JLD2 = "0.5.13"
-Oceananigans = "0.99"
+Oceananigans = "0.100"
 Printf = "1"
 RootSolvers = "0.4.4"
 julia = "1.10.10"

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+Breeze = "660aa2fb-d4c8-4359-a52c-9c057bc511da"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"

docs/src/index.md

@@ -41,8 +41,8 @@ Nx = Nz = 64
 Lz = 4 * 1024
 grid = RectilinearGrid(CPU(), size=(Nx, Nz), x=(0, 2Lz), z=(0, Lz), topology=(Periodic, Flat, Bounded))
 
-reference_constants = Breeze.Thermodynamics.ReferenceStateConstants(base_pressure=1e5, potential_temperature=288)
-buoyancy = Breeze.MoistAirBuoyancy(; reference_constants)
+reference_state = Breeze.Thermodynamics.ReferenceState(base_pressure=1e5, potential_temperature=288)
+buoyancy = Breeze.MoistAirBuoyancy(; reference_state)
 
 Q₀ = 1000 # heat flux in W / m²
 ρ₀ = Breeze.MoistAirBuoyancies.base_density(buoyancy) # air density at z=0
@@ -57,7 +57,7 @@ model = NonhydrostaticModel(; grid, advection, buoyancy,
                             boundary_conditions = (θ=θ_bcs, q=q_bcs))
 
 Δθ = 2 # ᵒK
-Tₛ = reference_constants.reference_potential_temperature # K
+Tₛ = reference_state.potential_temperature # K
 θᵢ(x, z) = Tₛ + Δθ * z / grid.Lz + 1e-2 * Δθ * randn()
 qᵢ(x, z) = 0 # 1e-2 + 1e-5 * rand()
 set!(model, θ=θᵢ, q=qᵢ)

docs/src/thermodynamics.md

@@ -2,7 +2,7 @@
 
 ```@setup thermo
 using Breeze
-thermo = AtmosphereThermodynamics()
+thermo = ThermodynamicConstants()
 ```
 
 Breeze implements thermodynamic relations for moist atmospheres ---
@@ -106,7 +106,7 @@ and
 ```
 
 ```@example thermo
-thermo = AtmosphereThermodynamics()
+thermo = ThermodynamicConstants()
 ```
 
 We can visualise the hydrostatic reference column implied by `thermo` by
@@ -119,20 +119,20 @@ using Breeze
 using Breeze.Thermodynamics: reference_pressure, reference_density
 using CairoMakie
 
-thermo = AtmosphereThermodynamics()
-constants = ReferenceStateConstants(base_pressure=101325, potential_temperature=288)
+thermo = ThermodynamicConstants()
+thermo = ReferenceState(base_pressure=101325, potential_temperature=288)
 grid = RectilinearGrid(size=160, z=(1, 12_000), topology=(Flat, Flat, Bounded))
 
 pᵣ = CenterField(grid)
 ρᵣ = CenterField(grid)
 
-set!(pᵣ, z -> reference_pressure(z, constants, thermo))
-set!(ρᵣ, z -> reference_density(z, constants, thermo))
+set!(pᵣ, z -> reference_pressure(z, thermo, thermo))
+set!(ρᵣ, z -> reference_density(z, thermo, thermo))
 
 Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo)
 cᵖᵈ = thermo.dry_air.heat_capacity
-p₀ = constants.base_pressure
-θ₀ = constants.reference_potential_temperature
+p₀ = thermo.base_pressure
+θ₀ = thermo.potential_temperature
 g = thermo.gravitational_acceleration
 
 # Verify that Tᵣ = θ₀ (1 - g z / (cᵖᵈ θ₀))
@@ -179,7 +179,7 @@ Central to Breeze's implementation of moist thermodynamics is a struct that
 holds parameters like the molar gas constant and molar masses,
 
 ```@example thermo
-thermo = AtmosphereThermodynamics()
+thermo = ThermodynamicConstants()
 ```
 
 The default parameter evince basic facts about water vapor air typical to Earth's atmosphere:
@@ -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/JRA55_saturation_specific_humidity.jl

@@ -1,14 +1,14 @@
 using JLD2
 using Breeze
-using Breeze: saturation_specific_humidity, AtmosphereThermodynamics
+using Breeze: saturation_specific_humidity, ThermodynamicConstants
 using GLMakie
 
 @load "JRA55_atmospheric_state_Jan_1_1991.jld2" q T p
 
-thermo = AtmosphereThermodynamics()
+thermo = ThermodynamicConstants()
 
 ρ = 1.2
-qᵛ★ = saturation_specific_humidity.(T, 1.2, Ref(thermo))
+qᵛ⁺ = saturation_specific_humidity.(T, 1.2, Ref(thermo))
 qˡ = @. max(0, q - q★)
 
 fig = Figure(size=(1200, 600))
@@ -17,23 +17,23 @@ ax2 = Axis(fig[2, 1], title="Saturation specific humidity")
 ax3 = Axis(fig[1, 1], title="Total specific humidity")
 ax4 = Axis(fig[2, 2], title="Liquid specific humidity")
 # heatmap!(ax1, T, colormap=:magma)
-# heatmap!(ax2, qᵛ★)
+# heatmap!(ax2, qᵛ⁺)
 # heatmap!(ax3, q, colormap=:grays)
 # heatmap!(ax4, qˡ, colormap=:grays, colorrange=(0, 0.001))
 
 # Compute cloudiness for instantaneous drop
 θ⁻ = T .- 10
 Ψ = Breeze.ThermodynamicState{Float64}.(θ⁻, q, 0)
-ℛ = Breeze.ReferenceStateConstants{Float64}(101325, 20)
+ℛ = Breeze.ReferenceState{Float64}(101325, 20)
 T⁻ = Breeze.temperature.(Ψ, Ref(ℛ), Ref(thermo))
-qᵛ★ = saturation_specific_humidity.(T⁻, 1.2, Ref(thermo))
-qˡ = @. max(0, q - qᵛ★)
+qᵛ⁺ = saturation_specific_humidity.(T⁻, 1.2, Ref(thermo))
+qˡ = @. max(0, q - qᵛ⁺)
 Π = Breeze.exner_function.(Ψ, Ref(ℛ), Ref(thermo))
 Tu = @. Π * θ⁻
 ΔT = T⁻ - Tu
 
 heatmap!(ax1, ΔT, colormap=:magma)
-heatmap!(ax2, qᵛ★)
+heatmap!(ax2, qᵛ⁺)
 heatmap!(ax3, q, colormap=:grays)
 heatmap!(ax4, qˡ, colormap=:grays, colorrange=(0, 0.01))
 display(fig)

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.thermo.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=5, 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/bomex.jl

@@ -42,8 +42,8 @@ u_bomex = AtmosphericProfilesLibrary.Bomex_u(FT)
 
 p₀ = 101500 # Pa
 θ₀ = 299.1 # K
-reference_constants = Breeze.Thermodynamics.ReferenceStateConstants(base_pressure=p₀, potential_temperature=θ₀)
-buoyancy = Breeze.MoistAirBuoyancy(; reference_constants) #, microphysics)
+reference_state = Breeze.Thermodynamics.ReferenceState(base_pressure=p₀, potential_temperature=θ₀)
+buoyancy = Breeze.MoistAirBuoyancy(; reference_state) #, microphysics)
 
 # Simple precipitation scheme from CloudMicrophysics    
 FT = eltype(grid)
@@ -197,8 +197,8 @@ add_callback!(simulation, compute_averages!)
 
 T = Breeze.TemperatureField(model)
 qˡ = Breeze.CondensateField(model, T)
-qᵛ★ = Breeze.SaturationField(model, T)
-rh = Field(model.tracers.q / qᵛ★) # relative humidity
+qᵛ⁺ = Breeze.SaturationField(model, T)
+rh = Field(model.tracers.q / qᵛ⁺) # relative humidity
 
 function progress(sim)
     compute!(T)
@@ -235,8 +235,8 @@ add_callback!(simulation, progress, IterationInterval(10))
 # The commented out lines below diagnose the forcing applied to model.tracers.q
 # using Oceananigans.Models: ForcingOperation
 # Sʳ = ForcingOperation(:q, model)
-# outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★, Sʳ))
-outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★))
+# outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ⁺, Sʳ))
+outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ⁺))
 averaged_outputs = NamedTuple(name => Average(outputs[name], dims=(1, 2)) for name in keys(outputs))
 
 filename = string("bomex_", Nx, "_", Ny, "_", Nz, ".jld2")

examples/free_convection.jl

@@ -11,17 +11,8 @@ grid = RectilinearGrid(arch, size=(Nx, Nz), x=(0, 2Lz), z=(0, Lz), topology=(Per
 
 p₀ = 101325 # Pa
 θ₀ = 288 # K
-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)
+reference_state = Breeze.Thermodynamics.ReferenceState(base_pressure=p₀, potential_temperature=θ₀)
+buoyancy = Breeze.MoistAirBuoyancy(; reference_state)
 
 ρ₀ = Breeze.MoistAirBuoyancies.base_density(buoyancy) # air density at z=0
 cₚ = buoyancy.thermodynamics.dry_air.heat_capacity
@@ -41,7 +32,7 @@ model = NonhydrostaticModel(; grid, advection, buoyancy,
 
 Lz = grid.Lz
 Δθ = 5 # K
-Tₛ = reference_constants.reference_potential_temperature # K
+Tₛ = reference_state.potential_temperature # K
 θᵢ(x, z) = Tₛ + Δθ * z / Lz + 1e-2 * Δθ * randn()
 qᵢ(x, z) = 0 # 1e-2 + 1e-5 * rand()
 set!(model, θ=θᵢ, q=qᵢ)
@@ -51,8 +42,8 @@ conjure_time_step_wizard!(simulation, cfl=0.7)
 
 T = Breeze.TemperatureField(model)
 qˡ = Breeze.CondensateField(model, T)
-qᵛ★ = Breeze.SaturationField(model, T)
-δ = Field(model.tracers.q - qᵛ★)
+qᵛ⁺ = Breeze.SaturationField(model, T)
+δ = Field(model.tracers.q - qᵛ⁺)
 
 function progress(sim)
     compute!(T)
@@ -89,7 +80,7 @@ add_callback!(simulation, progress, IterationInterval(10))
 
 using Oceananigans.Models: ForcingOperation
 Sʳ = ForcingOperation(:q, model)
-outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★, Sʳ))
+outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ⁺, Sʳ))
 
 ow = JLD2Writer(model, outputs,
                 filename = "free_convection.jld2",

examples/gate.jl

@@ -155,8 +155,8 @@ conjure_time_step_wizard!(simulation, cfl=0.7)
 
 T = Breeze.TemperatureField(model)
 qˡ = Breeze.CondensateField(model, T)
-qᵛ★ = Breeze.SaturationField(model, T)
-δ = Field(model.tracers.q - qᵛ★)
+qᵛ⁺ = Breeze.SaturationField(model, T)
+δ = Field(model.tracers.q - qᵛ⁺)
 
 function progress(sim)
     compute!(T)
@@ -193,8 +193,8 @@ add_callback!(simulation, progress, IterationInterval(10))
 
 # using Oceananigans.Models: ForcingOperation
 # Sʳ = ForcingOperation(:q, model)
-# outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★, Sʳ))
-outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★))
+# outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ⁺, Sʳ))
+outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ⁺))
 
 ow = JLD2Writer(model, outputs,
                 filename = "bomex.jld2",

examples/large_yeager_saturation_specific_humidity.jl

@@ -1,7 +1,7 @@
 using GLMakie
 using Breeze
 
-thermo = Breeze.AtmosphereThermodynamics()
+thermo = ThermodynamicConstants()
 
 saturation_specific_humidity_large_yeager(T, ρ) = 640380 * exp(-5107.4 / T) / ρ
 

examples/thermal_bubble.jl

@@ -28,8 +28,8 @@ grid = RectilinearGrid(arch,
 # Thermodynamic setup
 p₀ = 101325  # Pa - standard atmospheric pressure
 θ₀ = 300.0   # K - reference potential temperature
-reference_constants = ReferenceStateConstants(base_pressure=p₀, potential_temperature=θ₀)
-buoyancy = MoistAirBuoyancy(; reference_constants)
+reference_state = ReferenceState(base_pressure=p₀, potential_temperature=θ₀)
+buoyancy = MoistAirBuoyancy(; reference_state)
 
 # Advection scheme - WENO for high-order accuracy
 advection = WENO(order=5)