diff --git a/Project.toml b/Project.toml index 49c6bd0e..2ba6f94c 100644 --- a/Project.toml +++ b/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" diff --git a/docs/Project.toml b/docs/Project.toml index 0ffdc57e..b25e47a0 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/src/index.md b/docs/src/index.md index 37ab5280..518155f3 100644 --- a/docs/src/index.md +++ b/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ᵢ) diff --git a/docs/src/thermodynamics.md b/docs/src/thermodynamics.md index 8d8bad7c..b272e284 100644 --- a/docs/src/thermodynamics.md +++ b/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 diff --git a/examples/JRA55_saturation_specific_humidity.jl b/examples/JRA55_saturation_specific_humidity.jl index d22355d5..787df337 100644 --- a/examples/JRA55_saturation_specific_humidity.jl +++ b/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) \ No newline at end of file diff --git a/examples/anelastic_free_convection.jl b/examples/anelastic_free_convection.jl new file mode 100644 index 00000000..f7114210 --- /dev/null +++ b/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 diff --git a/examples/bomex.jl b/examples/bomex.jl index 7c5f72d3..028e69f2 100644 --- a/examples/bomex.jl +++ b/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") diff --git a/examples/free_convection.jl b/examples/free_convection.jl index b64e4835..411b2e4d 100644 --- a/examples/free_convection.jl +++ b/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", diff --git a/examples/gate.jl b/examples/gate.jl index e261ef7c..3369a81e 100644 --- a/examples/gate.jl +++ b/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", diff --git a/examples/large_yeager_saturation_specific_humidity.jl b/examples/large_yeager_saturation_specific_humidity.jl index 7ccf849a..24eb3b7d 100644 --- a/examples/large_yeager_saturation_specific_humidity.jl +++ b/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) / ρ diff --git a/examples/test/runtests.jl b/examples/test/runtests.jl deleted file mode 100644 index 59918ade..00000000 --- a/examples/test/runtests.jl +++ /dev/null @@ -1,3 +0,0 @@ -using Test - -@test 1==1 \ No newline at end of file diff --git a/examples/thermal_bubble.jl b/examples/thermal_bubble.jl index de3fa431..0cdda6fa 100644 --- a/examples/thermal_bubble.jl +++ b/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) diff --git a/examples/time_step_atmos_model.jl b/examples/time_step_atmos_model.jl index 58b97adb..5a9edeae 100644 --- a/examples/time_step_atmos_model.jl +++ b/examples/time_step_atmos_model.jl @@ -15,7 +15,7 @@ model = AtmosphereModel(grid, advection=WENO(), boundary_conditions=(; e=e_bcs)) Lz = grid.Lz Δθ = 5 # K -Tₛ = model.formulation.constants.reference_potential_temperature +Tₛ = model.formulation.thermo.potential_temperature # θᵢ(x, y, z) = Tₛ + Δθ * z / Lz qᵢ(x, y, z) = 0 Ξᵢ(x, y, z) = 1e-2 * randn() @@ -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" -=# \ No newline at end of file +=# diff --git a/src/AtmosphereModels/AtmosphereModels.jl b/src/AtmosphereModels/AtmosphereModels.jl index 68752a8f..bbc0a381 100644 --- a/src/AtmosphereModels/AtmosphereModels.jl +++ b/src/AtmosphereModels/AtmosphereModels.jl @@ -1,11 +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") include("set_atmosphere_model.jl") diff --git a/src/AtmosphereModels/anelastic_pressure_solver.jl b/src/AtmosphereModels/anelastic_pressure_solver.jl deleted file mode 100644 index f1e97ede..00000000 --- a/src/AtmosphereModels/anelastic_pressure_solver.jl +++ /dev/null @@ -1,195 +0,0 @@ -using Oceananigans.Operators: Δxᶜᵃᵃ, Δxᶠᵃᵃ, Δyᵃᶜᵃ, Δyᵃᶠᵃ, Δzᵃᵃᶜ, Δzᵃᵃᶠ -using Oceananigans.Grids: XYRegularRG, XZRegularRG, YZRegularRG, XYZRegularRG, stretched_dimensions - -import Oceananigans.Architectures: architecture - -struct FourierTridiagonalPoissonSolver{G, B, R, S, β, T} - grid :: G - batched_tridiagonal_solver :: B - source_term :: R - storage :: S - buffer :: β - transforms :: T -end - -architecture(solver::FourierTridiagonalPoissonSolver) = architecture(solver.grid) - -@kernel function compute_main_diagonal!(D, grid, λy, λz, ::XDirection) - j, k = @index(Global, NTuple) - Nx = size(grid, 1) - - # Using a homogeneous Neumann (zero Gradient) boundary condition: - @inbounds D[1, j, k] = -1 / Δxᶠᵃᵃ(2, j, k, grid) - Δxᶜᵃᵃ(1, j, k, grid) * (λy[j] + λz[k]) - for i in 2:Nx-1 - @inbounds D[i, j, k] = - (1 / Δxᶠᵃᵃ(i+1, j, k, grid) + 1 / Δxᶠᵃᵃ(i, j, k, grid)) - Δxᶜᵃᵃ(i, j, k, grid) * (λy[j] + λz[k]) - end - @inbounds D[Nx, j, k] = -1 / Δxᶠᵃᵃ(Nx, j, k, grid) - Δxᶜᵃᵃ(Nx, j, k, grid) * (λy[j] + λz[k]) -end - -@kernel function compute_main_diagonal!(D, grid, λx, λz, ::YDirection) - i, k = @index(Global, NTuple) - Ny = size(grid, 2) - - # Using a homogeneous Neumann (zero Gradient) boundary condition: - @inbounds D[i, 1, k] = -1 / Δyᵃᶠᵃ(i, 2, k, grid) - Δyᵃᶜᵃ(i, 1, k, grid) * (λx[i] + λz[k]) - for j in 2:Ny-1 - @inbounds D[i, j, k] = - (1 / Δyᵃᶠᵃ(i, j+1, k, grid) + 1 / Δyᵃᶠᵃ(i, j, k, grid)) - Δyᵃᶜᵃ(i, j, k, grid) * (λx[i] + λz[k]) - end - @inbounds D[i, Ny, k] = -1 / Δyᵃᶠᵃ(i, Ny, k, grid) - Δyᵃᶜᵃ(i, Ny, k, grid) * (λx[i] + λz[k]) -end - -@kernel function compute_main_diagonal!(D, grid, λx, λy, ::ZDirection) - i, j = @index(Global, NTuple) - Nz = size(grid, 3) - - # Using a homogeneous Neumann (zero Gradient) boundary condition: - @inbounds D[i, j, 1] = -1 / Δzᵃᵃᶠ(i, j, 2, grid) - Δzᵃᵃᶜ(i, j, 1, grid) * (λx[i] + λy[j]) - for k in 2:Nz-1 - @inbounds D[i, j, k] = - (1 / Δzᵃᵃᶠ(i, j, k+1, grid) + 1 / Δzᵃᵃᶠ(i, j, k, grid)) - Δzᵃᵃᶜ(i, j, k, grid) * (λx[i] + λy[j]) - end - @inbounds D[i, j, Nz] = -1 / Δzᵃᵃᶠ(i, j, Nz, grid) - Δzᵃᵃᶜ(i, j, Nz, grid) * (λx[i] + λy[j]) -end - -stretched_direction(::YZRegularRG) = XDirection() -stretched_direction(::XZRegularRG) = YDirection() -stretched_direction(::XYRegularRG) = ZDirection() - -dimension(::XDirection) = 1 -dimension(::YDirection) = 2 -dimension(::ZDirection) = 3 - -infer_launch_configuration(::XDirection) = :yz -infer_launch_configuration(::YDirection) = :xz -infer_launch_configuration(::ZDirection) = :xy - -Δξᶠ(i, grid, ::XDirection) = Δxᶠᵃᵃ(i, 1, 1, grid) -Δξᶠ(j, grid, ::YDirection) = Δyᵃᶠᵃ(1, j, 1, grid) -Δξᶠ(k, grid, ::ZDirection) = Δzᵃᵃᶠ(1, 1, k, grid) - -extent(grid) = (grid.Lx, grid.Ly, grid.Lz) - -""" - FourierTridiagonalPoissonSolver(grid, planner_flag = FFTW.PATIENT; - tridiagonal_direction = stretched_direction(grid)) - -Construct a `FourierTridiagonalPoissonSolver` on `grid` with `tridiagonal_direction` either -`XDirection()`, `YDirection()`, or `ZDirection()`. By default, the `tridiagonal_direction` will -be selected as `stretched_direction(grid)`, or `ZDirection()` if no directions are stretched. -variably spaced, the tridiagonal direction is -selected to be the direction of stretched grid spacing. -The Poisson equation is solved with an FFT-based eigenfunction expansion in the non-tridiagonal-directions -augmented by a tridiagonal solve in `tridiagonal_direction`. -The non-tridiagonal-directions must be uniformly spaced. -""" -function FourierTridiagonalPoissonSolver(grid, planner_flag=FFTW.PATIENT; tridiagonal_direction=nothing) - - # Try to guess what direction should be tridiagonal - if isnothing(tridiagonal_direction) - tridiagonal_direction = grid isa XYZRegularRG ? ZDirection() : stretched_direction(grid) - end - - tridiagonal_dim = dimension(tridiagonal_direction) - if topology(grid, tridiagonal_dim) != Bounded - msg = "`FourierTridiagonalPoissonSolver` can only be used \ - when the stretched direction's topology is `Bounded`." - throw(ArgumentError(msg)) - end - - # Compute discrete Poisson eigenvalues - N1, N2 = Tuple(el for (i, el) in enumerate(size(grid)) if i ≠ tridiagonal_dim) - T1, T2 = Tuple(el for (i, el) in enumerate(topology(grid)) if i ≠ tridiagonal_dim) - L1, L2 = Tuple(el for (i, el) in enumerate(extent(grid)) if i ≠ tridiagonal_dim) - - λ1 = poisson_eigenvalues(N1, L1, 1, T1()) - λ2 = poisson_eigenvalues(N2, L2, 2, T2()) - - arch = architecture(grid) - λ1 = on_architecture(arch, λ1) - λ2 = on_architecture(arch, λ2) - - # Plan required transforms for x and y - sol_storage = on_architecture(arch, zeros(complex(eltype(grid)), size(grid)...)) - transforms = plan_transforms(grid, sol_storage, planner_flag) - - # Lower and upper diagonals are the same - lower_diagonal = CUDA.@allowscalar [1 / Δξᶠ(q, grid, tridiagonal_direction) for q in 2:size(grid, tridiagonal_dim)] - lower_diagonal = on_architecture(arch, lower_diagonal) - upper_diagonal = lower_diagonal - - # Compute diagonal coefficients for each grid point - diagonal = on_architecture(arch, zeros(size(grid)...)) - launch_config = infer_launch_configuration(tridiagonal_direction) - launch!(arch, grid, launch_config, compute_main_diagonal!, diagonal, grid, λ1, λ2, tridiagonal_direction) - - # Set up batched tridiagonal solver - btsolver = BatchedTridiagonalSolver(grid; lower_diagonal, diagonal, upper_diagonal, tridiagonal_direction) - - # Need buffer for index permutations and transposes. - buffer_needed = arch isa GPU && Bounded in (T1, T2) - buffer = buffer_needed ? similar(sol_storage) : nothing - - # Storage space for right hand side of Poisson equation - rhs = on_architecture(arch, zeros(complex(eltype(grid)), size(grid)...)) - - return FourierTridiagonalPoissonSolver(grid, btsolver, rhs, sol_storage, buffer, transforms) -end - -function solve!(x, solver::FourierTridiagonalPoissonSolver, b=nothing) - !isnothing(b) && set_source_term!(solver, b) # otherwise, assume source term is set correctly - - arch = architecture(solver) - ϕ = solver.storage - - # Apply forward transforms in order - for transform! in solver.transforms.forward - transform!(solver.source_term, solver.buffer) - end - - # Solve tridiagonal system of linear equations at every column. - solve!(ϕ, solver.batched_tridiagonal_solver, solver.source_term) - - # Apply backward transforms in order - for transform! in solver.transforms.backward - transform!(ϕ, solver.buffer) - end - - # Set the volume mean of the solution to be zero. - # Solutions to Poisson's equation are only unique up to a constant (the global mean - # of the solution), so we need to pick a constant. We choose the constant to be zero - # so that the solution has zero-mean. - ϕ .= ϕ .- mean(ϕ) - - launch!(arch, solver.grid, :xyz, copy_real_component!, x, ϕ, indices(x)) - - return nothing -end - -""" - set_source_term!(solver, source_term) - -Sets the source term in the discrete Poisson equation `solver` to `source_term` by -multiplying it by the vertical grid spacing at cell centers in the stretched direction. -""" -function set_source_term!(solver::FourierTridiagonalPoissonSolver, source_term) - grid = solver.grid - arch = architecture(solver) - solver.source_term .= source_term - launch!(arch, grid, :xyz, multiply_by_stretched_spacing!, solver.source_term, grid) - return nothing -end - -@kernel function multiply_by_stretched_spacing!(a, grid::YZRegularRG) - i, j, k = @index(Global, NTuple) - @inbounds a[i, j, k] *= Δxᶜᵃᵃ(i, j, k, grid) -end - -@kernel function multiply_by_stretched_spacing!(a, grid::XZRegularRG) - i, j, k = @index(Global, NTuple) - @inbounds a[i, j, k] *= Δyᵃᶜᵃ(i, j, k, grid) -end - -@kernel function multiply_by_stretched_spacing!(a, grid::XYRegularRG) - i, j, k = @index(Global, NTuple) - @inbounds a[i, j, k] *= Δzᵃᵃᶜ(i, j, k, grid) -end - diff --git a/src/AtmosphereModels/atmosphere_model.jl b/src/AtmosphereModels/atmosphere_model.jl index 95912506..7d9e7ef8 100644 --- a/src/AtmosphereModels/atmosphere_model.jl +++ b/src/AtmosphereModels/atmosphere_model.jl @@ -1,11 +1,15 @@ using ..Thermodynamics: - AtmosphereThermodynamics, - ReferenceStateConstants, + ThermodynamicConstants, + ReferenceState, reference_pressure, reference_density, mixture_gas_constant, mixture_heat_capacity, - dry_air_gas_constant + dry_air_gas_constant, + field_names, + materialize_momentum_and_velocities, + collect_prognostic_fields, + formulation_pressure_solver using Oceananigans using Oceananigans.Advection: Centered, adapt_advection_order @@ -28,8 +32,6 @@ struct DefaultValue end tupleit(t::Tuple) = t tupleit(t) = tuple(t) -formulation_pressure_solver(formulation, grid) = nothing - mutable struct AtmosphereModel{Frm, Arc, Tst, Grd, Clk, Thm, Den, Mom, Eng, Wat, Hum, Tmp, Prs, Ppa, Sol, Vel, Trc, Adv, Cor, Frc, Mic, Cnd, Cls, Dif} <: AbstractModel{Tst, Arc} architecture :: Arc @@ -58,18 +60,19 @@ mutable struct AtmosphereModel{Frm, Arc, Tst, Grd, Clk, Thm, Den, Mom, Eng, Wat, diffusivity_fields :: Dif end -function default_formulation(grid, thermo) +function default_formulation(grid, thermodynamics) FT = eltype(grid) base_pressure = convert(FT, 101325) potential_temperature = convert(FT, 288) - constants = ReferenceStateConstants(base_pressure, potential_temperature) - return AnelasticFormulation(grid, constants, thermo) + reference_state = ReferenceState(base_pressure=base_pressure, + potential_temperature=potential_temperature) + return AnelasticFormulation(grid, reference_state, thermodynamics) end """ AtmosphereModel(grid; clock = Clock(grid), - thermodynamics = AtmosphereThermodynamics(eltype(grid)), + thermodynamics = ThermodynamicConstants(eltype(grid)), formulation = default_formulation(grid, thermodynamics), absolute_humidity = DefaultValue(), tracers = tuple(), @@ -77,7 +80,7 @@ end boundary_conditions = NamedTuple(), forcing = NamedTuple(), advection = WENO(order=5), - microphysics = WarmPhaseSaturationAdjustment(), + microphysics = nothing, timestepper = :RungeKutta3) Return an AtmosphereModel that uses the anelastic approximation following @@ -102,7 +105,7 @@ AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) """ function AtmosphereModel(grid; clock = Clock(grid), - thermodynamics = AtmosphereThermodynamics(eltype(grid)), + thermodynamics = ThermodynamicConstants(eltype(grid)), formulation = default_formulation(grid, thermodynamics), absolute_humidity = DefaultValue(), tracers = tuple(), @@ -110,7 +113,7 @@ function AtmosphereModel(grid; boundary_conditions = NamedTuple(), forcing = NamedTuple(), advection = WENO(order=5), - microphysics = WarmPhaseSaturationAdjustment(), + microphysics = nothing, timestepper = :RungeKutta3) arch = grid.architecture @@ -129,7 +132,7 @@ function AtmosphereModel(grid; density = materialize_density(formulation, grid) velocities, momentum = materialize_momentum_and_velocities(formulation, grid, boundary_conditions) - tracers = NamedTuple(n => CenterField(grid, boundary_conditions=boundary_conditions[n]) for name in tracers) + tracers = NamedTuple(name => CenterField(grid, boundary_conditions=boundary_conditions[name]) for name in tracers) condensates = materialize_condenstates(microphysics, grid) advection = adapt_advection_order(advection, grid) @@ -181,7 +184,10 @@ function AtmosphereModel(grid; closure, diffusivity_fields) - update_state!(model) + + # Provide a sensible default initial state (assumes anelastic formulation) + Tₛ = formulation.reference_state_constants.potential_temperature # K + set!(model, θ=Tₛ, enforce_mass_conservation=false) # consistent resting state return model end diff --git a/src/AtmosphereModels/saturation_adjustment.jl b/src/AtmosphereModels/saturation_adjustment.jl deleted file mode 100644 index dc0a6282..00000000 --- a/src/AtmosphereModels/saturation_adjustment.jl +++ /dev/null @@ -1,62 +0,0 @@ -##### -##### Saturation adjustment -##### - -function condensate_specific_humidity(T, state::AnelasticThermodynamicState, thermo) - qᵛ★ = saturation_specific_humidity(T, state.reference_density, thermo, thermo.liquid) - q = state.specific_humidity - return max(0, q - qᵛ★) -end - -@inline function compute_temperature(state::AnelasticThermodynamicState{FT}, thermo) where FT - θ = state.potential_temperature - θ == 0 && return zero(FT) - - # Generate guess for unsaturated conditions - Π = state.exner_function - T₁ = Π * θ - qˡ₁ = condensate_specific_humidity(T₁, state, thermo) - qˡ₁ <= 0 && return T₁ - - # If we made it this far, we have condensation - r₁ = saturation_adjustment_residual(T₁, qˡ₁, state, thermo) - - q = state.specific_humidity - ℒ = thermo.liquid.latent_heat - cᵖᵐ = mixture_heat_capacity(q, thermo) - T₂ = (T₁ + sqrt(T₁^2 + 4 * ℒ * qˡ₁ / cᵖᵐ)) / 2 - qˡ₂ = condensate_specific_humidity(T₂, state, thermo) - r₂ = saturation_adjustment_residual(T₂, qˡ₂, state, thermo) - - # Saturation adjustment - R = sqrt(max(T₂, T₁)) - ϵ = convert(FT, 1e-4) - δ = ϵ * R - iter = 0 - - while abs(r₂ - r₁) > δ - # Compute slope - ΔTΔr = (T₂ - T₁) / (r₂ - r₁) - - # Store previous values - r₁ = r₂ - T₁ = T₂ - - # Update - T₂ -= r₂ * ΔTΔr - qˡ₂ = condensate_specific_humidity(T₂, state, thermo) - r₂ = saturation_adjustment_residual(T₂, qˡ₂, state, thermo) - iter += 1 - end - - return T₂ -end - -@inline function saturation_adjustment_residual(T, qˡ, state, thermo) - ℒᵛ = thermo.liquid.latent_heat - q = state.specific_humidity - θ = state.potential_temperature - Π = state.exner_function - cᵖᵐ = mixture_heat_capacity(q, thermo) - return T^2 - ℒᵛ * qˡ / cᵖᵐ - Π * θ * T -end diff --git a/src/AtmosphereModels/set_atmosphere_model.jl b/src/AtmosphereModels/set_atmosphere_model.jl index 91f14e8b..834d6d9f 100644 --- a/src/AtmosphereModels/set_atmosphere_model.jl +++ b/src/AtmosphereModels/set_atmosphere_model.jl @@ -1,21 +1,44 @@ -import Oceananigans.Fields: set! +using Oceananigans.Grids: znode, Center using Oceananigans.TimeSteppers: update_state! using Oceananigans.BoundaryConditions: fill_halo_regions! using Oceananigans.Models.NonhydrostaticModels: compute_pressure_correction!, make_pressure_correction! +using ..Thermodynamics: exner_function, SpecificHumidities, mixture_heat_capacity, PotentialTemperatureState, temperature + +import Oceananigans.Fields: set! + +const c = Center() + +move_to_front(names, name) = tuple(name, filter(n -> n != name, names)...) + +function prioritize_names(names) + for n in (:w, :ρw, :v, :ρv, :u, :ρu, :qᵗ, :ρqᵗ) + if n ∈ names + names = move_to_front(names, n) + end + end + + return names +end + function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) - for (name, value) in kw + names = collect(keys(kw)) + prioritized = prioritize_names(names) + + for name in prioritized + value = kw[name] # Prognostic variables if name ∈ propertynames(model.momentum) - ϕ = getproperty(model.momentum, name) - set!(ϕ, value) + ρu = getproperty(model.momentum, name) + set!(ρu, value) elseif name ∈ propertynames(model.tracers) - ϕ = getproperty(model.tracers, name) + c = getproperty(model.tracers, name) + set!(c, value) elseif name == :ρe - ϕ = model.energy - elseif name == :ρq - ϕ = model.absolute_humidity + set!(model.energy, value) + elseif name == :ρqᵗ + set!(model.absolute_humidity, value) end # Setting diagnostic variables @@ -23,17 +46,22 @@ function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) θ = model.temperature # use scratch set!(θ, value) - ρʳ = model.formulation.reference_density - cᵖᵈ = model.thermodynamics.dry_air.heat_capacity - ϕ = model.energy - value = ρʳ * cᵖᵈ * θ - elseif name == :q - q = model.specific_humidity - set!(q, value) + grid = model.grid + arch = grid.architecture + thermo = model.thermodynamics + formulation = model.formulation + energy = model.energy + specific_humidity = model.specific_humidity + launch!(arch, grid, :xyz, _energy_from_potential_temperature!, energy, grid, + θ, specific_humidity, formulation, thermo) + elseif name == :qᵗ + qᵗ = model.specific_humidity + set!(qᵗ, value) ρʳ = model.formulation.reference_density - ϕ = model.absolute_humidity - value = ρʳ * q + ρqᵗ = model.absolute_humidity + set!(ρqᵗ, ρʳ * qᵗ) + elseif name ∈ (:u, :v, :w) u = model.velocities[name] set!(u, value) @@ -41,14 +69,12 @@ function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) ρʳ = model.formulation.reference_density ϕ = model.momentum[Symbol(:ρ, name)] value = ρʳ * u + set!(ϕ, value) end - - set!(ϕ, value) - fill_halo_regions!(ϕ, model.clock, fields(model)) end # Apply a mask - # foreach(mask_immersed_field!, prognostic_fields(model)) + foreach(mask_immersed_field!, prognostic_fields(model)) update_state!(model, compute_tendencies=false) if enforce_mass_conservation @@ -59,5 +85,33 @@ function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) update_state!(model, compute_tendencies=false) end + fill_halo_regions!(model.energy) + return nothing end + +@kernel function _energy_from_potential_temperature!(moist_static_energy, grid, + potential_temperature, + specific_humidity, + formulation, + thermo) + i, j, k = @index(Global, NTuple) + + @inbounds begin + ρʳ = formulation.reference_density[i, j, k] + qᵗ = specific_humidity[i, j, k] + θ = potential_temperature[i, j, k] + z = znode(i, j, k, grid, c, c, c) + end + + # Assume non-condensed state + # TODO: relax this assumption + q = SpecificHumidities(qᵗ, zero(qᵗ), zero(qᵗ)) + ref = formulation.reference_state_constants + 𝒰 = PotentialTemperatureState(θ, q, z, ref) + T = temperature(𝒰, thermo) + ℒ₀ = thermo.liquid.latent_heat + g = thermo.gravitational_acceleration + cᵖᵐ = mixture_heat_capacity(q, thermo) + @inbounds moist_static_energy[i, j, k] = ρʳ * (cᵖᵐ * T + g * z + qᵗ * ℒ₀) +end diff --git a/src/AtmosphereModels/update_atmosphere_model_state.jl b/src/AtmosphereModels/update_atmosphere_model_state.jl index a80ee5df..b288878c 100644 --- a/src/AtmosphereModels/update_atmosphere_model_state.jl +++ b/src/AtmosphereModels/update_atmosphere_model_state.jl @@ -1,17 +1,42 @@ using ..Thermodynamics: saturation_specific_humidity, mixture_heat_capacity, - mixture_gas_constant + mixture_gas_constant, + AnelasticFormulation, + thermodynamic_state, + solve_for_anelastic_pressure!, + formulation_pressure_solver, + total_specific_humidity, + _pressure_correct_momentum! using Oceananigans.BoundaryConditions: fill_halo_regions!, compute_x_bcs!, compute_y_bcs!, compute_z_bcs! using Oceananigans.ImmersedBoundaries: mask_immersed_field! using Oceananigans.Architectures: architecture +using Oceananigans.Operators: ℑzᵃᵃᶠ, ℑzᵃᵃᶜ -import Oceananigans.TimeSteppers: update_state!, compute_flux_bc_tendencies! +import Oceananigans.TimeSteppers: update_state!, compute_flux_bc_tendencies!, compute_pressure_correction!, make_pressure_correction! import Oceananigans: fields, prognostic_fields const AnelasticModel = AtmosphereModel{<:AnelasticFormulation} +function compute_pressure_correction!(model::AnelasticModel, Δt) + foreach(mask_immersed_field!, model.momentum) + fill_halo_regions!(model.momentum, model.clock, fields(model)) + + ρŨ = model.momentum + solver = model.pressure_solver + if solver === nothing + solver = formulation_pressure_solver(model.formulation, model.grid) + model.pressure_solver = solver + end + pₙ = model.nonhydrostatic_pressure + solve_for_anelastic_pressure!(pₙ, solver, ρŨ, Δt) + + fill_halo_regions!(pₙ) + + return nothing +end + function prognostic_fields(model::AnelasticModel) thermodynamic_fields = (ρe=model.energy, ρq=model.absolute_humidity) return merge(model.momentum, thermodynamic_fields, model.condensates, model.tracers) @@ -22,11 +47,22 @@ fields(model::AnelasticModel) = prognostic_fields(model) function update_state!(model::AnelasticModel, callbacks=[]; compute_tendencies=true) fill_halo_regions!(prognostic_fields(model), model.clock, fields(model), async=true) compute_auxiliary_variables!(model) - update_hydrostatic_pressure!(model) compute_tendencies && compute_tendencies!(model) return nothing end +function make_pressure_correction!(model::AnelasticModel, Δt) + launch!(model.architecture, model.grid, :xyz, + _pressure_correct_momentum!, + model.momentum, + model.grid, + Δt, + model.nonhydrostatic_pressure, + model.formulation.reference_density) + + return nothing +end + function compute_auxiliary_variables!(model) grid = model.grid arch = grid.architecture @@ -43,6 +79,7 @@ function compute_auxiliary_variables!(model) model.temperature, model.specific_humidity, grid, + model.microphysics, model.thermodynamics, formulation, model.energy, @@ -70,31 +107,22 @@ end end end -@kernel function _compute_auxiliary_thermodynamic_variables!(temperature, +@kernel function _compute_auxiliary_thermodynamic_variables!(temperature_field, specific_humidity, grid, + microphysics, thermo, formulation, energy, absolute_humidity) i, j, k = @index(Global, NTuple) - 𝒰 = thermodynamic_state(i, j, k, grid, formulation, thermo, energy, absolute_humidity) - @inbounds specific_humidity[i, j, k] = 𝒰.specific_humidity - - # Saturation adjustment - T = compute_temperature(𝒰, thermo) - @inbounds temperature[i, j, k] = T -end + 𝒰 = thermodynamic_state(i, j, k, grid, formulation, energy, absolute_humidity) + @inbounds specific_humidity[i, j, k] = total_specific_humidity(𝒰) -#= -@inline function specific_volume(state, ref, thermo) - T = temperature(state, ref, thermo) - Rᵐ = mixture_gas_constant(state.q, thermo) - pᵣ = reference_pressure(state.z, ref, thermo) - return Rᵐ * T / pᵣ + T = temperature(𝒰, thermo) + @inbounds temperature_field[i, j, k] = T end -=# using Oceananigans.Advection: div_𝐯u, div_𝐯v, div_𝐯w, div_Uc using Oceananigans.Coriolis: x_f_cross_U, y_f_cross_U, z_f_cross_U @@ -131,8 +159,10 @@ function compute_tendencies!(model::AnelasticModel) Gρe = model.timestepper.Gⁿ.ρe ρe = model.energy Fρe = model.forcing.ρe - ρe_args = tuple(ρe, Fρe, scalar_args...) - launch!(arch, grid, :xyz, compute_scalar_tendency!, Gρe, grid, ρe_args) + ρe_args = tuple(ρe, Fρe, scalar_args..., ρᵣ, + model.formulation, model.temperature, + model.specific_humidity, model.thermodynamics, model.condensates, model.microphysics) + launch!(arch, grid, :xyz, compute_energy_tendency!, Gρe, grid, ρe_args) ρq = model.absolute_humidity Gρq = model.timestepper.Gⁿ.ρq @@ -146,11 +176,17 @@ end hydrostatic_pressure_gradient_x(i, j, k, grid, pₕ′) = ∂xᶠᶜᶜ(i, j, k, grid, pₕ′) hydrostatic_pressure_gradient_y(i, j, k, grid, pₕ′) = ∂yᶜᶠᶜ(i, j, k, grid, pₕ′) -@kernel function compute_scalar_tendency!(Gc, grid, args) +@kernel function compute_scalar_tendency!(Gρc, grid, args) + i, j, k = @index(Global, NTuple) + @inbounds Gρc[i, j, k] = scalar_tendency(i, j, k, grid, args...) +end + +@kernel function compute_energy_tendency!(Gρe, grid, args) i, j, k = @index(Global, NTuple) - @inbounds Gc[i, j, k] = scalar_tendency(i, j, k, grid, args...) + @inbounds Gρe[i, j, k] = energy_tendency(i, j, k, grid, args...) end + @kernel function compute_x_momentum_tendency!(Gρu, grid, args) i, j, k = @index(Global, NTuple) @inbounds Gρu[i, j, k] = x_momentum_tendency(i, j, k, grid, args...) @@ -206,6 +242,19 @@ end + forcing(i, j, k, grid, clock, model_fields)) end +@inline function ρᵣbᶜᶜᶠ(i, j, k, grid, ρᵣ, T, q, formulation, thermo) + + ρᵣᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, ρᵣ) + bᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, buoyancy, formulation, T, q, thermo) + + return ρᵣᶜᶜᶠ * bᶜᶜᶠ +end + +@inline function ρᵣwbᶜᶜᶠ(i, j, k, grid, w, ρᵣ, T, q, formulation, thermo) + ρᵣb = ρᵣbᶜᶜᶠ(i, j, k, grid, ρᵣ, T, q, formulation, thermo) + return @inbounds ρᵣb * w[i, j, k] +end + @inline function z_momentum_tendency(i, j, k, grid, advection, velocities, @@ -220,14 +269,12 @@ end specific_humidity, thermo) - ρᵣᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, reference_density) - bᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, buoyancy, - formulation, temperature, specific_humidity, thermo) return ( - div_𝐯w(i, j, k, grid, advection, velocities, momentum.ρw) - z_f_cross_U(i, j, k, grid, coriolis, momentum) - + ρᵣᶜᶜᶠ * bᶜᶜᶠ - + forcing(i, j, k, grid, clock, model_fields)) + + ρᵣbᶜᶜᶠ(i, j, k, grid, reference_density, temperature, specific_humidity, formulation, thermo) + + forcing(i, j, k, grid, clock, model_fields) + ) end @inline function scalar_tendency(i, j, k, grid, @@ -242,23 +289,30 @@ end + forcing(i, j, k, grid, clock, model_fields)) end -#= @inline function energy_tendency(i, j, k, grid, - formulation, energy, forcing, advection, velocities, - condensates, - microphysics clock, - model_fields) + model_fields, + reference_density, + formulation, + temperature, + specific_humidity, + thermo, + condensates, + microphysics) + + + ρᵣwbᶜᶜᶜ = ℑzᵃᵃᶜ(i, j, k, grid, ρᵣwbᶜᶜᶠ, velocities.w, reference_density, + temperature, specific_humidity, formulation, thermo) return ( - div_Uc(i, j, k, grid, advection, velocities, energy) - + microphysical_energy_tendency(i, j, k, grid, formulation, microphysics, condensates) + + ρᵣwbᶜᶜᶜ + # + microphysical_energy_tendency(i, j, k, grid, formulation, microphysics, condensates) + forcing(i, j, k, grid, clock, model_fields)) end -=# """ Apply boundary conditions by adding flux divergences to the right-hand-side. """ function compute_flux_bc_tendencies!(model::AtmosphereModel) diff --git a/src/AtmosphereModels/update_hydrostatic_pressure.jl b/src/AtmosphereModels/update_hydrostatic_pressure.jl index 84257f70..687034d3 100644 --- a/src/AtmosphereModels/update_hydrostatic_pressure.jl +++ b/src/AtmosphereModels/update_hydrostatic_pressure.jl @@ -7,6 +7,9 @@ using Oceananigans.ImmersedBoundaries: PartialCellBottom, ImmersedBoundaryGrid using Oceananigans.Grids: topology using Oceananigans.Grids: XFlatGrid, YFlatGrid using Oceananigans.Utils: KernelParameters +using Oceananigans.BoundaryConditions: fill_halo_regions! + +using ..Thermodynamics: specific_volume, reference_specific_volume const c = Center() const f = Face() diff --git a/src/Breeze.jl b/src/Breeze.jl index 548a9180..7cd87a24 100644 --- a/src/Breeze.jl +++ b/src/Breeze.jl @@ -7,8 +7,8 @@ module Breeze export MoistAirBuoyancy, - AtmosphereThermodynamics, - ReferenceStateConstants, + ThermodynamicConstants, + ReferenceState, AnelasticFormulation, AtmosphereModel, TemperatureField, @@ -16,6 +16,7 @@ export PhaseTransitionConstants, CondensedPhase, mixture_gas_constant, + prognostic_fields, mixture_heat_capacity using Oceananigans @@ -51,6 +52,10 @@ export include("Thermodynamics/Thermodynamics.jl") using .Thermodynamics +import .Thermodynamics: AnelasticFormulation + +include("Microphysics/Microphysics.jl") +using .Microphysics include("MoistAirBuoyancies.jl") using .MoistAirBuoyancies diff --git a/src/Microphysics/Microphysics.jl b/src/Microphysics/Microphysics.jl new file mode 100644 index 00000000..78b1f8e6 --- /dev/null +++ b/src/Microphysics/Microphysics.jl @@ -0,0 +1,16 @@ +module Microphysics + +export temperature, specific_volume + +using ..Thermodynamics: + mixture_heat_capacity, + mixture_gas_constant, + exner_function, + reference_pressure + +import ..Thermodynamics: condensate_specific_humidity, temperature + +include("nothing_microphysics.jl") +include("saturation_adjustment.jl") + +end # module diff --git a/src/Microphysics/nothing_microphysics.jl b/src/Microphysics/nothing_microphysics.jl new file mode 100644 index 00000000..90240a66 --- /dev/null +++ b/src/Microphysics/nothing_microphysics.jl @@ -0,0 +1,26 @@ +using ..Thermodynamics: + ThermodynamicConstants, + ReferenceState, + SpecificHumidities, + exner_function, + reference_pressure, + mixture_heat_capacity, + mixture_gas_constant + + +#= +# fully compressible case +@inline function temperature(e, U, q, z, thermo::ThermodynamicConstants) + Π = exner_function(q, z, ref, thermo) + return Π * θ +end +=# + +@inline function specific_volume(T, q::SpecificHumidities, z, + ref::ReferenceState, + thermo::ThermodynamicConstants) + + Rᵐ = mixture_gas_constant(q, thermo) + pᵣ = reference_pressure(z, ref, thermo) + return Rᵐ * T / pᵣ +end diff --git a/src/Microphysics/saturation_adjustment.jl b/src/Microphysics/saturation_adjustment.jl new file mode 100644 index 00000000..71e25edd --- /dev/null +++ b/src/Microphysics/saturation_adjustment.jl @@ -0,0 +1,148 @@ +using ..Thermodynamics: total_specific_humidity, saturation_specific_humidity + +struct SaturationAdjustmentMicrophysics end + +##### +##### Microphysics saturation adjustment utilities +##### + +# Solve θ = T/Π (1 - ℒ qˡ / (cᵖᵐ T)) for temperature T with qˡ = max(0, q - qᵛ⁺) +# by iterating on the root of f(T) = T² - Π θ T - ℒ qˡ / cᵖᵐ. +@inline function adjust_temperature_and_humidities(𝒰, thermo) + θ = 𝒰.potential_temperature + θ == 0 && return zero(θ), 𝒰.humidities + + # qᵈ = dry_air_mass_fraction(𝒰.humidities)) + q = 𝒰.humidities + qᵗ = total_specific_humidity(q) + z = 𝒰.height + ref = 𝒰.reference_state + Π = exner_function(𝒰.humidities, 𝒰.height, 𝒰.reference_state, thermo) + T₁ = Π * θ + qˡ₁ = adjusted_condensate_specific_humidity(T₁, qᵗ, z, ref, thermo) + + if qˡ₁ <= 0 + qᵛ = total_specific_humidity(𝒰.humidities) + qˡ = zero(qᵛ) + qˢ = zero(qᵛ) + q = SpecificHumidities(qᵛ, qˡ, qˢ) + return T₁, q + end + + qᵛ₁ = qᵗ - qˡ₁ + q₁ = SpecificHumidities(qᵛ₁, qˡ₁, zero(qᵗ)) + r₁ = saturation_adjustment_residual(T₁, Π, q₁, θ, thermo) + + ℒ = thermo.liquid.latent_heat + cᵖᵐ = mixture_heat_capacity(q, thermo) + T₂ = (T₁ + sqrt(T₁^2 + 4 * ℒ * qˡ₁ / cᵖᵐ)) / 2 + qˡ₂ = adjusted_condensate_specific_humidity(T₂, qᵗ, z, ref, thermo) + qᵛ₂ = qᵗ - qˡ₂ + q₂ = SpecificHumidities(qᵛ₂, qˡ₂, zero(qᵗ)) + r₂ = saturation_adjustment_residual(T₂, Π, q₂, θ, thermo) + + R = sqrt(max(T₂, T₁)) + ϵ = convert(typeof(T₂), 1e-4) + δ = ϵ * R + + while abs(r₂ - r₁) > δ + ΔTΔr = (T₂ - T₁) / (r₂ - r₁) + + r₁ = r₂ + T₁ = T₂ + + T₂ -= r₂ * ΔTΔr + qˡ₂ = adjusted_condensate_specific_humidity(T₂, qᵗ, z, ref, thermo) + q₂ = SpecificHumidities(qᵛ₂, qˡ₂, zero(qᵗ)) + r₂ = saturation_adjustment_residual(T₂, Π, q₂, θ, thermo) + end + + qᵗ = total_specific_humidity(𝒰.humidities) + qᵛ = qᵗ - qˡ₂ + qˢ = zero(qˡ₂) + adjusted_q = SpecificHumidities(qᵛ, qˡ₂, qˢ) + + return T₂, adjusted_q +end + +function adjusted_condensate_specific_humidity(T, qᵗ, z, ref::ReferenceState, thermo) + qᵛ⁺ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid) + return max(0, qᵗ - qᵛ⁺) +end + +function adjusted_ice_specific_humidity(T, qᵗ, z, ref::ReferenceState, thermo) + qˢ⁺ = saturation_specific_humidity(T, z, ref, thermo, thermo.solid) + return max(0, qᵗ - qˢ⁺) +end + +@inline function saturation_adjustment_residual(T, Π, q, θ, thermo) + ℒᵛ = thermo.liquid.latent_heat + cᵖᵐ = mixture_heat_capacity(q, thermo) + qˡ = q.liquid + return T^2 - ℒᵛ * qˡ / cᵖᵐ - Π * θ * T +end + + +##### +##### Microphysics schemes +##### + +#= +# Solve +# θ = T/Π ( 1 - ℒ qˡ / (cᵖᵐ T)) +# for temperature T with qˡ = max(0, q - qᵛ⁺). +# root of: f(T) = T² - Π θ T - ℒ qˡ / cᵖᵐ +@inline function temperature(state::ThermodynamicState{FT}, ref, thermo) where FT + state.θ == 0 && return zero(FT) + + qᵛ = state.q + qᵈ = one(qᵛ) - qᵛ + + # Generate guess for unsaturated conditions + Π = exner_function(state, ref, thermo) + T₁ = Π * state.θ + qˡ₁ = condensate_specific_humidity(T₁, state, ref, thermo) + qˡ₁ <= 0 && return T₁ + + # If we made it this far, we have condensation + r₁ = saturation_adjustment_residual(T₁, Π, qˡ₁, state, thermo) + + ℒ = thermo.liquid.latent_heat + cᵖᵐ = mixture_heat_capacity(qᵈ, qᵛ, thermo) + T₂ = (T₁ + sqrt(T₁^2 + 4 * ℒ * qˡ₁ / cᵖᵐ)) / 2 + qˡ₂ = condensate_specific_humidity(T₂, state, ref, thermo) + r₂ = saturation_adjustment_residual(T₂, Π, qˡ₂, state, thermo) + + # Saturation adjustment + R = sqrt(max(T₂, T₁)) + ϵ = convert(FT, 1e-4) + δ = ϵ * R + iter = 0 + + while abs(r₂ - r₁) > δ + # Compute slope + ΔTΔr = (T₂ - T₁) / (r₂ - r₁) + + # Store previous values + r₁ = r₂ + T₁ = T₂ + + # Update + T₂ -= r₂ * ΔTΔr + qˡ₂ = condensate_specific_humidity(T₂, state, ref, thermo) + r₂ = saturation_adjustment_residual(T₂, Π, qˡ₂, state, thermo) + iter += 1 + end + + return T₂ +end + +@inline function specific_volume(state, ref, thermo) + T = temperature(state, ref, thermo) + qᵛ = state.q + qᵈ = one(qᵛ) - qᵛ + Rᵐ = mixture_gas_constant(qᵈ, qᵛ, thermo) + pᵣ = reference_pressure(state.z, ref, thermo) + return Rᵐ * T / pᵣ +end +=# diff --git a/src/MoistAirBuoyancies.jl b/src/MoistAirBuoyancies.jl index c6843f12..52c4dd52 100644 --- a/src/MoistAirBuoyancies.jl +++ b/src/MoistAirBuoyancies.jl @@ -19,27 +19,36 @@ import Oceananigans.BuoyancyFormulations: AbstractBuoyancyFormulation, required_tracers using ..Thermodynamics: - AtmosphereThermodynamics, - ReferenceStateConstants, - mixture_heat_capacity, - mixture_gas_constant, + ThermodynamicConstants, + ReferenceState, + SpecificHumidities, reference_specific_volume, - reference_pressure + PotentialTemperatureState + +using ..Microphysics: + SaturationAdjustmentMicrophysics, + adjust_temperature_and_humidities import ..Thermodynamics: base_density, - saturation_specific_humidity, - condensate_specific_humidity + saturation_specific_humidity + +import ..Microphysics: + temperature, + specific_volume, + adjusted_condensate_specific_humidity -struct MoistAirBuoyancy{FT, AT} <: AbstractBuoyancyFormulation{Nothing} - reference_constants :: ReferenceStateConstants{FT} +struct MoistAirBuoyancy{FT, AT, M} <: AbstractBuoyancyFormulation{Nothing} + reference_state :: ReferenceState{FT} thermodynamics :: AT + microphysics :: M end """ MoistAirBuoyancy(FT=Oceananigans.defaults.FloatType; - thermodynamics = AtmosphereThermodynamics(FT), - reference_constants = ReferenceStateConstants{FT}(101325, 290)) + thermodynamics = ThermodynamicConstants(FT), + reference_state = ReferenceState{FT}(101325, 290)) + microphysics = nothing) Return a MoistAirBuoyancy formulation that can be provided as input to an [`AtmosphereModel`](@ref Breeze.AtmosphereModels.AtmosphereModel) or an @@ -56,8 +65,8 @@ julia> using Breeze, Oceananigans julia> buoyancy = MoistAirBuoyancy() MoistAirBuoyancy -├── reference_constants: Breeze.Thermodynamics.ReferenceStateConstants{Float64} -└── thermodynamics: AtmosphereThermodynamics +├── reference_state: Breeze.Thermodynamics.ReferenceState{Float64} +└── thermodynamics: ThermodynamicConstants julia> model = NonhydrostaticModel(; grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 2, 3)), buoyancy, tracers = (:θ, :q)) @@ -72,56 +81,67 @@ NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ``` """ function MoistAirBuoyancy(FT=Oceananigans.defaults.FloatType; - thermodynamics = AtmosphereThermodynamics(FT), - reference_constants = ReferenceStateConstants{FT}(101325, 290)) + thermodynamics = ThermodynamicConstants(FT), + reference_state = ReferenceState{FT}(101325, 290), + microphysics = SaturationAdjustmentMicrophysics()) AT = typeof(thermodynamics) - return MoistAirBuoyancy{FT, AT}(reference_constants, thermodynamics) + MT = typeof(microphysics) + return MoistAirBuoyancy{FT, AT, MT}(reference_state, thermodynamics, microphysics) end Base.summary(b::MoistAirBuoyancy) = "MoistAirBuoyancy" function Base.show(io::IO, b::MoistAirBuoyancy) print(io, summary(b), "\n", - "├── reference_constants: ", summary(b.reference_constants), "\n", - "└── thermodynamics: ", summary(b.thermodynamics)) + "├── reference_state: ", summary(b.reference_state), "\n", + "├── thermodynamics: ", summary(b.thermodynamics), "\n", + "└── microphysics: ", summary(b.microphysics)) end required_tracers(::MoistAirBuoyancy) = (:θ, :q) -reference_density(z, mb::MoistAirBuoyancy) = reference_density(z, mb.reference_constants, mb.thermodynamics) -base_density(mb::MoistAirBuoyancy) = base_density(mb.reference_constants, mb.thermodynamics) +reference_density(z, mb::MoistAirBuoyancy) = reference_density(z, mb.reference_state, mb.thermodynamics) +base_density(mb::MoistAirBuoyancy) = base_density(mb.reference_state, mb.thermodynamics) ##### -##### +##### buoyancy ##### const c = Center() +# Nothing microphysics: no condensates +function compute_temperature_and_humidities(::Nothing, θ, qᵗ, z, mb) + q = SpecificHumidities(qᵗ, zero(qᵗ), zero(qᵗ)) + 𝒰 = PotentialTemperatureState(θ, q, z, mb.reference_state) + T = temperature(𝒰, mb.thermodynamics) + return T, q +end + +function compute_temperature_and_humidities(::SaturationAdjustmentMicrophysics, θ, qᵗ, z, mb) + q = SpecificHumidities(qᵗ, zero(qᵗ), zero(qᵗ)) + 𝒰 = PotentialTemperatureState(θ, q, z, mb.reference_state) + T, q = adjust_temperature_and_humidities(𝒰, mb.thermodynamics) + return T, q +end + @inline function buoyancy_perturbationᶜᶜᶜ(i, j, k, grid, mb::MoistAirBuoyancy, tracers) z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) θ = @inbounds tracers.θ[i, j, k] - q = @inbounds tracers.q[i, j, k] - 𝒰 = HeightReferenceThermodynamicState(θ, q, z) + qᵗ = @inbounds tracers.q[i, j, k] - # Perform Saturation adjustment - α = specific_volume(𝒰, mb.reference_constants, mb.thermodynamics) + T, q = compute_temperature_and_humidities(mb.microphysics, θ, qᵗ, z, mb) + α = specific_volume(T, q, z, mb.reference_state, mb.thermodynamics) - # Compute reference specific volume - αʳ = reference_specific_volume(z, mb.reference_constants, mb.thermodynamics) + # Compute buoyancy + αʳ = reference_specific_volume(z, mb.reference_state, mb.thermodynamics) g = mb.thermodynamics.gravitational_acceleration - # Formulation in terms of base density: - # ρ₀ = base_density(mb.reference_constants, mb.thermodynamics) - # return ρ₀ * g * (α - αʳ) - return g * (α - αʳ) / αʳ end @inline ∂z_b(i, j, k, grid, mb::MoistAirBuoyancy, tracers) = ∂zᶜᶜᶠ(i, j, k, grid, buoyancy_perturbationᶜᶜᶜ, mb, tracers) -const c = Center() - ##### ##### Temperature ##### @@ -130,8 +150,9 @@ function temperature(i, j, k, grid::AbstractGrid, mb::MoistAirBuoyancy, θ, q) z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) θi = @inbounds θ[i, j, k] qi = @inbounds q[i, j, k] - 𝒰 = HeightReferenceThermodynamicState(θi, qi, z) - return temperature(𝒰, mb.reference_constants, mb.thermodynamics) + q = SpecificHumidities(qi, zero(qi), zero(qi)) + 𝒰 = PotentialTemperatureState(θi, q, z, mb.reference_state) + return temperature(𝒰, mb.thermodynamics) end struct TemperatureKernelFunction end @@ -156,7 +177,7 @@ end @inline function saturation_specific_humidity(i, j, k, grid, mb::MoistAirBuoyancy, T, condensed_phase) z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) Ti = @inbounds T[i, j, k] - return saturation_specific_humidity(Ti, z, mb.reference_constants, mb.thermodynamics, condensed_phase) + return saturation_specific_humidity(Ti, z, mb.reference_state, mb.thermodynamics, condensed_phase) end struct PhaseTransitionConstantsKernel{T, P} @@ -176,6 +197,7 @@ end function SaturationField(model, T = TemperatureField(model); condensed_phase = model.buoyancy.formulation.thermodynamics.liquid) + func = PhaseTransitionConstantsKernel(condensed_phase, T) grid = model.grid buoyancy = model.buoyancy.formulation @@ -196,8 +218,8 @@ Adapt.adapt_structure(to, ck::CondensateKernel) = CondensateKernel(adapt(to, ck. @inline function condensate_specific_humidity(i, j, k, grid, mb::MoistAirBuoyancy, T, q) z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) Ti = @inbounds T[i, j, k] - qi = @inbounds q[i, j, k] - qˡ = condensate_specific_humidity(Ti, qi, z, mb.reference_constants, mb.thermodynamics) + qᵗ = @inbounds q[i, j, k] + qˡ = adjusted_condensate_specific_humidity(Ti, qᵗ, z, mb.reference_state, mb.thermodynamics) return qˡ end @@ -215,120 +237,4 @@ function CondensateField(model, T=TemperatureField(model)) return Field(op) end -##### -##### Saturation adjustment -##### - -# Organizing information about the state is a WIP -struct HeightReferenceThermodynamicState{FT} - θ :: FT - q :: FT - z :: FT -end - -condensate_specific_humidity(T, state::HeightReferenceThermodynamicState, ref, thermo) = - condensate_specific_humidity(T, state.q, state.z, ref, thermo) - -# Solve -# θ = T/Π ( 1 - ℒ qˡ / (cᵖᵐ T)) -# for temperature T with qˡ = max(0, q - qᵛ★). -# root of: f(T) = T² - Π θ T - ℒ qˡ / cᵖᵐ -@inline function temperature(state::HeightReferenceThermodynamicState{FT}, ref, thermo) where FT - state.θ == 0 && return zero(FT) - - # Generate guess for unsaturated conditions - Π = exner_function(state, ref, thermo) - T₁ = Π * state.θ - qˡ₁ = condensate_specific_humidity(T₁, state, ref, thermo) - qˡ₁ <= 0 && return T₁ - - # If we made it this far, we have condensation - r₁ = saturation_adjustment_residual(T₁, Π, qˡ₁, state, thermo) - - ℒ = thermo.liquid.latent_heat - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - T₂ = (T₁ + sqrt(T₁^2 + 4 * ℒ * qˡ₁ / cᵖᵐ)) / 2 - qˡ₂ = condensate_specific_humidity(T₂, state, ref, thermo) - r₂ = saturation_adjustment_residual(T₂, Π, qˡ₂, state, thermo) - - # Saturation adjustment - R = sqrt(max(T₂, T₁)) - ϵ = convert(FT, 1e-4) - δ = ϵ * R - iter = 0 - - while abs(r₂ - r₁) > δ - # Compute slope - ΔTΔr = (T₂ - T₁) / (r₂ - r₁) - - # Store previous values - r₁ = r₂ - T₁ = T₂ - - # Update - T₂ -= r₂ * ΔTΔr - qˡ₂ = condensate_specific_humidity(T₂, state, ref, thermo) - r₂ = saturation_adjustment_residual(T₂, Π, qˡ₂, state, thermo) - iter += 1 - end - - return T₂ -end - -@inline function saturation_adjustment_residual(T, Π, qˡ, state::HeightReferenceThermodynamicState, thermo) - ℒᵛ = thermo.liquid.latent_heat - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - return T^2 - ℒᵛ * qˡ / cᵖᵐ - Π * state.θ * T -end - -@inline function specific_volume(state::HeightReferenceThermodynamicState, ref, thermo) - T = temperature(state, ref, thermo) - Rᵐ = mixture_gas_constant(state.q, thermo) - pᵣ = reference_pressure(state.z, ref, thermo) - return Rᵐ * T / pᵣ -end - -@inline function exner_function(state::HeightReferenceThermodynamicState, ref, thermo) - Rᵐ = mixture_gas_constant(state.q, thermo) - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - inv_ϰᵐ = Rᵐ / cᵖᵐ - pᵣ = reference_pressure(state.z, ref, thermo) - p₀ = ref.base_pressure - return (pᵣ / p₀)^inv_ϰᵐ -end - -##### -##### Reference implementation of an "unsaturated" moist air buoyancy model, -##### which assumes unsaturated air -##### - -struct UnsaturatedMoistAirBuoyancy{FT} <: AbstractBuoyancyFormulation{Nothing} - expansion_coefficient :: FT - reference_potential_temperature :: FT - gas_constant_ratio :: FT -end - -function UnsaturatedMoistAirBuoyancy(FT=Oceananigans.defaults.FloatType; - expansion_coefficient = 3.27e-2, - reference_potential_temperature = 0, - gas_constant_ratio = 1.61) - - return UnsaturatedMoistAirBuoyancy{FT}(expansion_coefficient, - reference_potential_temperature, - gas_constant_ratio) -end - -required_tracers(::UnsaturatedMoistAirBuoyancy) = (:θ, :q) - -@inline function buoyancy_perturbationᶜᶜᶜ(i, j, k, grid, mb::UnsaturatedMoistAirBuoyancy, tracers) - β = mb.expansion_coefficient - θ₀ = mb.reference_potential_temperature - ϵᵥ = mb.gas_constant_ratio - δ = ϵᵥ - 1 - θ = @inbounds tracers.θ[i, j, k] - q = @inbounds tracers.q[i, j, k] - θᵥ = θ * (1 + δ * q) - return β * (θᵥ - θ₀) -end - end # module diff --git a/src/Thermodynamics/Thermodynamics.jl b/src/Thermodynamics/Thermodynamics.jl index 1ab6ef6c..72d8165c 100644 --- a/src/Thermodynamics/Thermodynamics.jl +++ b/src/Thermodynamics/Thermodynamics.jl @@ -1,9 +1,12 @@ module Thermodynamics -export AtmosphereThermodynamics, ReferenceStateConstants, IdealGas, PhaseTransitionConstants, CondensedPhase, mixture_gas_constant, mixture_heat_capacity +export ThermodynamicConstants, ReferenceState, IdealGas, PhaseTransitionConstants, CondensedPhase, mixture_gas_constant, mixture_heat_capacity -include("atmosphere_thermodynamics.jl") +include("thermodynamics_constants.jl") include("vapor_saturation.jl") include("reference_states.jl") +include("dynamic_states.jl") +include("anelastic_formulation.jl") +include("boussinesq_formulation.jl") -end # module \ No newline at end of file +end # module diff --git a/src/AtmosphereModels/anelastic_formulation.jl b/src/Thermodynamics/anelastic_formulation.jl similarity index 69% rename from src/AtmosphereModels/anelastic_formulation.jl rename to src/Thermodynamics/anelastic_formulation.jl index 5d983b5b..22165157 100644 --- a/src/AtmosphereModels/anelastic_formulation.jl +++ b/src/Thermodynamics/anelastic_formulation.jl @@ -1,12 +1,14 @@ using Oceananigans.Architectures: architecture -using Oceananigans.Grids: inactive_cell -using Oceananigans.Utils: prettysummary -using Oceananigans.Operators: Δzᵃᵃᶜ, Δzᵃᵃᶠ, divᶜᶜᶜ -using Oceananigans.Solvers: solve! +using Oceananigans.Grids: inactive_cell, Center, znode, ZDirection +using Oceananigans.Utils: prettysummary, launch! +using Oceananigans.Operators: Δzᵃᵃᶜ, Δzᵃᵃᶠ, Δzᶜᶜᶜ, ℑzᵃᵃᶠ, ∂xᶠᶜᶜ, ∂yᶜᶠᶜ, ∂zᶜᶜᶠ, divᶜᶜᶜ +using Oceananigans.Solvers: solve!, FourierTridiagonalPoissonSolver +using Oceananigans.Fields: Field, XFaceField, YFaceField, ZFaceField, set! +using Oceananigans.BoundaryConditions: fill_halo_regions!, FieldBoundaryConditions, regularize_field_boundary_conditions using ..Thermodynamics: - AtmosphereThermodynamics, - ReferenceStateConstants, + ThermodynamicConstants, + ReferenceState, reference_pressure, reference_density, mixture_gas_constant, @@ -16,72 +18,63 @@ using ..Thermodynamics: using KernelAbstractions: @kernel, @index import Oceananigans.Solvers: tridiagonal_direction, compute_main_diagonal!, compute_lower_diagonal! -import Oceananigans.TimeSteppers: compute_pressure_correction!, make_pressure_correction! ##### ##### Formulation definition ##### struct AnelasticFormulation{FT, F} - constants :: ReferenceStateConstants{FT} + reference_state_constants :: ReferenceState{FT} reference_pressure :: F reference_density :: F end -const AnelasticModel = AtmosphereModel{<:AnelasticFormulation} - -function Base.summary(formulation::AnelasticFormulation) - p₀ = formulation.constants.base_pressure - θᵣ = formulation.constants.reference_potential_temperature - return string("AnelasticFormulation(p₀=", prettysummary(p₀), - ", θᵣ=", prettysummary(θᵣ), ")") -end +Base.summary(formulation::AnelasticFormulation) = + string("AnelasticFormulation with ", summary(formulation.reference_state_constants)) Base.show(io::IO, formulation::AnelasticFormulation) = print(io, "AnelasticFormulation") field_names(::AnelasticFormulation, tracer_names) = (:ρu, :ρv, :ρw, :ρe, :ρq, tracer_names...) -struct AnelasticThermodynamicState{FT} - potential_temperature :: FT - specific_humidity :: FT - reference_density :: FT - reference_pressure :: FT - exner_function :: FT -end - -function AnelasticFormulation(grid, state_constants, thermo) +function AnelasticFormulation(grid, ref, thermo) pᵣ = Field{Nothing, Nothing, Center}(grid) ρᵣ = Field{Nothing, Nothing, Center}(grid) - set!(pᵣ, z -> reference_pressure(z, state_constants, thermo)) - set!(ρᵣ, z -> reference_density(z, state_constants, thermo)) + set!(pᵣ, z -> reference_pressure(z, ref, thermo)) + set!(ρᵣ, z -> reference_density(z, ref, thermo)) fill_halo_regions!(pᵣ) fill_halo_regions!(ρᵣ) - return AnelasticFormulation(state_constants, pᵣ, ρᵣ) + return AnelasticFormulation(ref, pᵣ, ρᵣ) end ##### ##### Thermodynamic state ##### -function thermodynamic_state(i, j, k, grid, formulation::AnelasticFormulation, thermo, energy, absolute_humidity) +struct AnelasticThermodynamicState{FT} + specific_moist_static_energy :: FT + specific_humidity :: FT + height :: FT +end + +const c = Center() + +function thermodynamic_state(i, j, k, grid, + formulation::AnelasticFormulation, + energy, + absolute_humidity) @inbounds begin - e = energy[i, j, k] - pᵣ = formulation.reference_pressure[i, j, k] + ρe = energy[i, j, k] ρᵣ = formulation.reference_density[i, j, k] - ρq = absolute_humidity[i, j, k] + ρqᵗ = absolute_humidity[i, j, k] end - cᵖᵈ = thermo.dry_air.heat_capacity - θ = e / (cᵖᵈ * ρᵣ) - - q = ρq / ρᵣ - Rᵐ = mixture_gas_constant(q, thermo) - cᵖᵐ = mixture_heat_capacity(q, thermo) - - p₀ = formulation.constants.base_pressure - Π = (pᵣ / p₀)^(Rᵐ / cᵖᵐ) + e = ρe / ρᵣ + z = znode(i, j, k, grid, c, c, c) - return AnelasticThermodynamicState(θ, q, ρᵣ, pᵣ, Π) + qᵗ = ρqᵗ / ρᵣ + q = SpecificHumidities(qᵗ, zero(qᵗ), zero(qᵗ)) + + return MoistStaticEnergyState(e, q, z) end @inline function specific_volume(i, j, k, grid, formulation, temperature, specific_humidity, thermo) @@ -91,7 +84,9 @@ end T = temperature[i, j, k] end - Rᵐ = mixture_gas_constant(q, thermo) + qᵛ = q + qᵈ = one(qᵛ) - qᵛ + Rᵐ = mixture_gas_constant(qᵈ, qᵛ, thermo) return Rᵐ * T / pᵣ end @@ -99,7 +94,7 @@ end @inline function reference_specific_volume(i, j, k, grid, formulation, thermo) Rᵈ = dry_air_gas_constant(thermo) pᵣ = @inbounds formulation.reference_pressure[i, j, k] - θᵣ = formulation.constants.reference_potential_temperature + θᵣ = formulation.thermo.potential_temperature return Rᵈ * θᵣ / pᵣ end @@ -146,6 +141,7 @@ function formulation_pressure_solver(anelastic_formulation::AnelasticFormulation reference_density = anelastic_formulation.reference_density tridiagonal_formulation = AnelasticTridiagonalSolverFormulation(reference_density) solver = FourierTridiagonalPoissonSolver(grid; tridiagonal_formulation) + # solver = FourierTridiagonalPoissonSolver(grid) #; tridiagonal_formulation) return solver end @@ -199,22 +195,6 @@ end end end -function compute_pressure_correction!(model::AnelasticModel, Δt) - # Mask immersed velocities - foreach(mask_immersed_field!, model.momentum) - fill_halo_regions!(model.momentum, model.clock, fields(model)) - - ρʳ = model.formulation.reference_density - ρŨ = model.momentum - solver = model.pressure_solver - pₙ = model.nonhydrostatic_pressure - solve_for_anelastic_pressure!(pₙ, solver, ρŨ, Δt) - - fill_halo_regions!(pₙ) - - return nothing -end - function solve_for_anelastic_pressure!(pₙ, solver, ρŨ, Δt) compute_anelastic_source_term!(solver, ρŨ, Δt) solve!(pₙ, solver) @@ -237,17 +217,6 @@ end @inbounds rhs[i, j, k] = active * Δzᶜᶜᶜ(i, j, k, grid) * δ / Δt end -#= -function compute_source_term!(solver::DistributedFourierTridiagonalPoissonSolver, Ũ) - rhs = solver.storage.zfield - arch = architecture(solver) - grid = solver.local_grid - tdir = solver.batched_tridiagonal_solver.tridiagonal_direction - launch!(arch, grid, :xyz, _fourier_tridiagonal_source_term!, rhs, tdir, grid, Ũ) - return nothing -end -=# - ##### ##### Fractional and time stepping ##### @@ -267,16 +236,3 @@ Update the predictor momentum (ρu, ρv, ρw) with the non-hydrostatic pressure @inbounds M.ρv[i, j, k] -= ρᶜ * Δt * ∂yᶜᶠᶜ(i, j, k, grid, αʳ_pₙ) @inbounds M.ρw[i, j, k] -= ρᶠ * Δt * ∂zᶜᶜᶠ(i, j, k, grid, αʳ_pₙ) end - -function make_pressure_correction!(model::AnelasticModel, Δt) - - launch!(model.architecture, model.grid, :xyz, - _pressure_correct_momentum!, - model.momentum, - model.grid, - Δt, - model.nonhydrostatic_pressure, - model.formulation.reference_density) - - return nothing -end \ No newline at end of file diff --git a/src/Thermodynamics/boussinesq_formulation.jl b/src/Thermodynamics/boussinesq_formulation.jl new file mode 100644 index 00000000..fb671448 --- /dev/null +++ b/src/Thermodynamics/boussinesq_formulation.jl @@ -0,0 +1,12 @@ +struct BoussinesqThermodynamicState{FT} + potential_temperature :: FT + specific_humidity :: FT + height :: FT +end + +@inline function specific_volume(state, microphysics, ref, thermo) + T = temperature(state, microphysics, ref, thermo) + Rᵐ = mixture_gas_constant(q, thermo) + pᵣ = reference_pressure(z, ref, thermo) + return Rᵐ * T / pᵣ +end diff --git a/src/Thermodynamics/dynamic_states.jl b/src/Thermodynamics/dynamic_states.jl new file mode 100644 index 00000000..2d581089 --- /dev/null +++ b/src/Thermodynamics/dynamic_states.jl @@ -0,0 +1,47 @@ +struct PotentialTemperatureState{FT, H, R} + potential_temperature :: FT + humidities :: H + height :: FT + reference_state :: R +end + +@inline function exner_function(𝒰::PotentialTemperatureState, thermo::ThermodynamicConstants) + q = 𝒰.humidities + z = 𝒰.height + ref = 𝒰.reference_state + return exner_function(q, z, ref, thermo) +end + +@inline total_specific_humidity(state::PotentialTemperatureState) = + total_specific_humidity(state.humidities) + +@inline function temperature(state::PotentialTemperatureState, thermo::ThermodynamicConstants) + θ = state.potential_temperature + q = state.humidities + z = state.height + ref = state.reference_state + Π = exner_function(q, z, ref, thermo) + return Π * θ +end + +struct MoistStaticEnergyState{FT, H} + moist_static_energy :: FT + humidities :: H + height :: FT +end + +@inline total_specific_humidity(state::MoistStaticEnergyState) = + total_specific_humidity(state.humidities) + + # No microphysics: no liquid, only vapor +@inline function temperature(state::MoistStaticEnergyState, thermo::ThermodynamicConstants) + e = state.moist_static_energy + q = state.humidities + z = state.height + g = thermo.gravitational_acceleration + cᵖᵐ = mixture_heat_capacity(q, thermo) + qᵗ = total_specific_humidity(q) + ℒ₀ = thermo.liquid.latent_heat + cᵖᵐ_T = e - g * z - qᵗ * ℒ₀ + return cᵖᵐ_T / cᵖᵐ +end diff --git a/src/Thermodynamics/reference_states.jl b/src/Thermodynamics/reference_states.jl index 457d45e2..02a74467 100644 --- a/src/Thermodynamics/reference_states.jl +++ b/src/Thermodynamics/reference_states.jl @@ -4,17 +4,25 @@ using Oceananigans ##### Reference state computations for Boussinesq and Anelastic models ##### -struct ReferenceStateConstants{FT} +struct ReferenceState{FT} base_pressure :: FT # base pressure: reference pressure at z=0 - reference_potential_temperature :: FT # constant reference potential temperature + potential_temperature :: FT # constant reference potential temperature end -function ReferenceStateConstants(FT = Oceananigans.defaults.FloatType; - base_pressure = 101325, - potential_temperature = 288) +function Base.summary(ref::ReferenceState) + return string("ReferenceState(", + "p₀=", prettysummary(ref.base_pressure), ", ", + "θᵣ=", prettysummary(ref.potential_temperature), ")") +end + +Base.show(io::IO, ref::ReferenceState) = print(io, summary(ref)) + +function ReferenceState(FT = Oceananigans.defaults.FloatType; + base_pressure = 101325, + potential_temperature = 288) - return ReferenceStateConstants{FT}(convert(FT, base_pressure), - convert(FT, potential_temperature)) + return ReferenceState{FT}(convert(FT, base_pressure), + convert(FT, potential_temperature)) end """ @@ -23,7 +31,7 @@ end Compute the reference density associated with the reference pressure and potential temperature. The reference density is defined as the density of dry air at the reference pressure and temperature. """ -@inline function reference_density(z, ref::ReferenceStateConstants, thermo) +@inline function reference_density(z, ref::ReferenceState, thermo) Rᵈ = dry_air_gas_constant(thermo) cᵖᵈ = thermo.dry_air.heat_capacity pᵣ = reference_pressure(z, ref, thermo) @@ -32,40 +40,39 @@ The reference density is defined as the density of dry air at the reference pres return ρ₀ * (pᵣ / p₀)^(1 - Rᵈ / cᵖᵈ) end -@inline function base_density(ref::ReferenceStateConstants, thermo) +@inline function base_density(ref::ReferenceState, thermo) Rᵈ = dry_air_gas_constant(thermo) p₀ = ref.base_pressure - θᵣ = ref.reference_potential_temperature + θᵣ = ref.potential_temperature return p₀ / (Rᵈ * θᵣ) end -@inline function reference_specific_volume(z, ref::ReferenceStateConstants, thermo) +@inline function reference_specific_volume(z, ref::ReferenceState, thermo) Rᵈ = dry_air_gas_constant(thermo) pᵣ = reference_pressure(z, ref, thermo) - θᵣ = ref.reference_potential_temperature + θᵣ = ref.potential_temperature return Rᵈ * θᵣ / pᵣ end -@inline function reference_pressure(z, ref::ReferenceStateConstants, thermo) +@inline function reference_pressure(z, ref::ReferenceState, thermo) cᵖᵈ = thermo.dry_air.heat_capacity Rᵈ = dry_air_gas_constant(thermo) g = thermo.gravitational_acceleration - θᵣ = ref.reference_potential_temperature + θᵣ = ref.potential_temperature p₀ = ref.base_pressure return p₀ * (1 - g * z / (cᵖᵈ * θᵣ))^(cᵖᵈ / Rᵈ) end -@inline function saturation_specific_humidity(T, z, ref::ReferenceStateConstants, thermo, condensed_phase) - ρ = reference_density(z, ref, thermo) - return saturation_specific_humidity(T, ρ, thermo, condensed_phase) -end - -function condensate_specific_humidity(T, q, z, ref::ReferenceStateConstants, thermo) - qᵛ★ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid) - return max(0, q - qᵛ★) +@inline function exner_function(q, z, ref, thermo) + Rᵐ = mixture_gas_constant(q, thermo) + cᵖᵐ = mixture_heat_capacity(q, thermo) + inv_ϰᵐ = Rᵐ / cᵖᵐ + pᵣ = reference_pressure(z, ref, thermo) + p₀ = ref.base_pressure + return (pᵣ / p₀)^inv_ϰᵐ end -function ice_specific_humidity(T, q, z, ref::ReferenceStateConstants, thermo) - qi★ = saturation_specific_humidity(T, z, ref, thermo, thermo.solid) - return max(0, q - qi★) +@inline function saturation_specific_humidity(T, z, ref::ReferenceState, thermo, condensed_phase) + ρ = reference_density(z, ref, thermo) + return saturation_specific_humidity(T, ρ, thermo, condensed_phase) end diff --git a/src/Thermodynamics/atmosphere_thermodynamics.jl b/src/Thermodynamics/thermodynamics_constants.jl similarity index 66% rename from src/Thermodynamics/atmosphere_thermodynamics.jl rename to src/Thermodynamics/thermodynamics_constants.jl index 672ea284..0ad2e6a4 100644 --- a/src/Thermodynamics/atmosphere_thermodynamics.jl +++ b/src/Thermodynamics/thermodynamics_constants.jl @@ -90,7 +90,7 @@ end liquid_water(FT) = CondensedPhase(FT; latent_heat=2500800, heat_capacity=4181) water_ice(FT) = CondensedPhase(FT; latent_heat=2834000, heat_capacity=2108) -struct AtmosphereThermodynamics{FT, C, S} +struct ThermodynamicConstants{FT, C, S} molar_gas_constant :: FT gravitational_acceleration :: FT energy_reference_temperature :: FT @@ -102,9 +102,9 @@ struct AtmosphereThermodynamics{FT, C, S} solid :: S end -Base.summary(at::AtmosphereThermodynamics{FT}) where FT = "AtmosphereThermodynamics{$FT}" +Base.summary(at::ThermodynamicConstants{FT}) where FT = "ThermodynamicConstants{$FT}" -function Base.show(io::IO, at::AtmosphereThermodynamics) +function Base.show(io::IO, at::ThermodynamicConstants) print(io, summary(at), ":", '\n', "├── molar_gas_constant: ", at.molar_gas_constant, "\n", "├── gravitational_acceleration: ", at.gravitational_acceleration, "\n", @@ -117,9 +117,9 @@ function Base.show(io::IO, at::AtmosphereThermodynamics) "└── solid: ", at.solid) end -Base.eltype(::AtmosphereThermodynamics{FT}) where FT = FT +Base.eltype(::ThermodynamicConstants{FT}) where FT = FT -function Adapt.adapt_structure(to, thermo::AtmosphereThermodynamics) +function Adapt.adapt_structure(to, thermo::ThermodynamicConstants) molar_gas_constant = adapt(to, thermo.molar_gas_constant) gravitational_acceleration = adapt(to, thermo.gravitational_acceleration) dry_air = adapt(to, thermo.dry_air) @@ -132,7 +132,7 @@ function Adapt.adapt_structure(to, thermo::AtmosphereThermodynamics) FT = typeof(molar_gas_constant) C = typeof(liquid) S = typeof(solid) - return AtmosphereThermodynamics{FT, C, S}(molar_gas_constant, + return ThermodynamicConstants{FT, C, S}(molar_gas_constant, gravitational_acceleration, energy_reference_temperature, triple_point_temperature, @@ -144,7 +144,7 @@ function Adapt.adapt_structure(to, thermo::AtmosphereThermodynamics) end """ - AtmosphereThermodynamics(FT = Oceananigans.defaults.FloatType; + ThermodynamicConstants(FT = Oceananigans.defaults.FloatType; gravitational_acceleration = 9.81, molar_gas_constant = 8.314462618, energy_reference_temperature = 273.16, @@ -158,7 +158,7 @@ end solid = water_ice(FT), condensed_phases = nothing) -Create `AtmosphereThermodynamics` with parameters that represent gaseous mixture of dry "air" +Create `ThermodynamicConstants` with parameters that represent gaseous mixture of dry "air" and vapor, as well as condensed liquid and solid phases. The `triple_point_temperature` and `triple_point_pressure` may be combined with internal energy parameters for condensed phases to compute the vapor pressure @@ -193,7 +193,7 @@ The advantage of using reference values at the triple point is that the same val can then be used for both condensation (vapor → liquid) and deposition (vapor → ice). """ -function AtmosphereThermodynamics(FT = Oceananigans.defaults.FloatType; +function ThermodynamicConstants(FT = Oceananigans.defaults.FloatType; molar_gas_constant = 8.314462618, gravitational_acceleration = 9.81, energy_reference_temperature = 273.16, @@ -212,146 +212,78 @@ function AtmosphereThermodynamics(FT = Oceananigans.defaults.FloatType; vapor = IdealGas(FT; molar_mass = vapor_molar_mass, heat_capacity = vapor_heat_capacity) - return AtmosphereThermodynamics(convert(FT, molar_gas_constant), - convert(FT, gravitational_acceleration), - convert(FT, energy_reference_temperature), - convert(FT, triple_point_temperature), - convert(FT, triple_point_pressure), - dry_air, - vapor, - liquid, - solid) + return ThermodynamicConstants(convert(FT, molar_gas_constant), + convert(FT, gravitational_acceleration), + convert(FT, energy_reference_temperature), + convert(FT, triple_point_temperature), + convert(FT, triple_point_pressure), + dry_air, + vapor, + liquid, + solid) end -const AT = AtmosphereThermodynamics +const AT = ThermodynamicConstants const IG = IdealGas @inline vapor_gas_constant(thermo::AT) = thermo.molar_gas_constant / thermo.vapor.molar_mass @inline dry_air_gas_constant(thermo::AT) = thermo.molar_gas_constant / thermo.dry_air.molar_mass -const NonCondensingAtmosphereThermodynamics{FT} = AtmosphereThermodynamics{FT, Nothing, Nothing} +struct SpecificHumidities{FT} + vapor :: FT + liquid :: FT + ice :: FT +end + +@inline dry_air_mass_fraction(q::SpecificHumidities) = 1 - q.vapor - q.liquid - q.ice +@inline total_specific_humidity(q::SpecificHumidities) = q.vapor + q.liquid + q.ice """ - mixture_gas_constant(q, thermo) + mixture_gas_constant(qᵈ, qᵛ, thermo) Compute the gas constant of moist air given the specific humidity `q` and thermodynamic parameters `thermo`. The mixture gas constant is calculated as a weighted average of the dry air -and water vapor gas constants: +and water vapor gas thermo: ```math -R_m = R_d (1 - q) + R_v q +Rᵐ = qᵈ * Rᵈ + qᵛ * Rᵛ ``` where: -- `R_d` is the dry air gas constant -- `R_v` is the water vapor gas constant -- `q` is the specific humidity (mass fraction of water vapor) +- `Rᵈ` is the dry air gas constant +- `Rᵛ` is the water vapor gas constant +- `qᵈ` is the mass fraction of dry air +- `qᵛ` is the mass fraction of water vapor # Arguments -- `q`: Specific humidity (dimensionless) -- `thermo`: `AtmosphereThermodynamics` instance containing gas constants +- `qᵈ`: Mass fraction of dry air (dimensionless) +- `qᵛ`: Mass fraction of water vapor (dimensionless) +- `thermo`: `ThermodynamicConstants` instance containing gas thermo # Returns - Gas constant of the moist air mixture in J/(kg·K) """ -@inline function mixture_gas_constant(q, thermo::AT) +@inline function mixture_gas_constant(q::SpecificHumidities, thermo::AT) + qᵈ = dry_air_mass_fraction(q) + qᵛ = q.vapor Rᵈ = dry_air_gas_constant(thermo) Rᵛ = vapor_gas_constant(thermo) - return Rᵈ * (1 - q) + Rᵛ * q + return qᵈ * Rᵈ + qᵛ * Rᵛ end """ - mixture_heat_capacity(q, thermo) + mixture_heat_capacity(qᵈ, qᵛ, thermo) Compute the heat capacity of state air given the total specific humidity q and assuming that condensate mass ratio qᶜ ≪ q, where qℓ is the mass ratio of liquid condensate. """ -@inline function mixture_heat_capacity(q, thermo::AT) +@inline function mixture_heat_capacity(q::SpecificHumidities, thermo::AT) + qᵈ = dry_air_mass_fraction(q) + qᵛ = q.vapor cᵖᵈ = thermo.dry_air.heat_capacity cᵖᵛ = thermo.vapor.heat_capacity - return cᵖᵈ * (1 - q) + cᵖᵛ * q -end - -##### -##### state thermodynamics for a Boussinesq model -##### - -# Organizing information about the state is a WIP -struct ThermodynamicState{FT} - θ :: FT - q :: FT - z :: FT -end - -struct ReferenceState{FT} - p₀ :: FT # base pressure: reference pressure at z=0 - θ :: FT # constant reference potential temperature -end - -Adapt.adapt_structure(to, ref::ReferenceState) = - ReferenceState(adapt(to, ref.p₀), - adapt(to, ref.θ)) - -function ReferenceState(FT = Oceananigans.defaults.FloatType; - base_pressure = 101325, - potential_temperature = 288) - - return ReferenceState{FT}(convert(FT, base_pressure), - convert(FT, potential_temperature)) -end - -""" - reference_density(z, ref, thermo) - -Compute the reference density associated with the reference pressure and potential temperature. -The reference density is defined as the density of dry air at the reference pressure and temperature. -""" -@inline function reference_density(z, ref, thermo) - Rᵈ = dry_air_gas_constant(thermo) - p = reference_pressure(z, ref, thermo) - return p / (Rᵈ * ref.θ) -end - -@inline function base_density(ref, thermo) - Rᵈ = dry_air_gas_constant(thermo) - return ref.p₀ / (Rᵈ * ref.θ) -end - -@inline function reference_specific_volume(z, ref, thermo) - Rᵈ = dry_air_gas_constant(thermo) - p = reference_pressure(z, ref, thermo) - return Rᵈ * ref.θ / p -end - -@inline function reference_pressure(z, ref, thermo) - cᵖᵈ = thermo.dry_air.heat_capacity - Rᵈ = dry_air_gas_constant(thermo) - inv_ϰᵈ = Rᵈ / cᵖᵈ - g = thermo.gravitational_acceleration - return ref.p₀ * (1 - g * z / (cᵖᵈ * ref.θ))^inv_ϰᵈ -end - -@inline function saturation_specific_humidity(T, z, ref::ReferenceState, thermo, condensed_phase) - ρ = reference_density(z, ref, thermo) - return saturation_specific_humidity(T, ρ, thermo, condensed_phase) -end - -@inline function exner_function(state, ref, thermo) - Rᵐ = mixture_gas_constant(state.q, thermo) - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - inv_ϰᵐ = Rᵐ / cᵖᵐ - pᵣ = reference_pressure(state.z, ref, thermo) - p₀ = ref.base_pressure - return (pᵣ / p₀)^inv_ϰᵐ -end - -condensate_specific_humidity(T, state, ref, thermo) = - condensate_specific_humidity(T, state.q, state.z, ref, thermo) - -function condensate_specific_humidity(T, q, z, ref, thermo) - qᵛ★ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid) - return max(0, q - qᵛ★) + return qᵈ * cᵖᵈ + qᵛ * cᵖᵛ end diff --git a/src/microphysics.jl b/src/microphysics.jl deleted file mode 100644 index 0c807df7..00000000 --- a/src/microphysics.jl +++ /dev/null @@ -1,62 +0,0 @@ -##### -##### Microphysics schemes -##### - -# Solve -# θ = T/Π ( 1 - ℒ qˡ / (cᵖᵐ T)) -# for temperature T with qˡ = max(0, q - qᵛ★). -# root of: f(T) = T² - Π θ T - ℒ qˡ / cᵖᵐ -@inline function temperature(state::ThermodynamicState{FT}, ref, thermo) where FT - state.θ == 0 && return zero(FT) - - # Generate guess for unsaturated conditions - Π = exner_function(state, ref, thermo) - T₁ = Π * state.θ - qˡ₁ = condensate_specific_humidity(T₁, state, ref, thermo) - qˡ₁ <= 0 && return T₁ - - # If we made it this far, we have condensation - r₁ = saturation_adjustment_residual(T₁, Π, qˡ₁, state, thermo) - - ℒ = thermo.liquid.latent_heat - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - T₂ = (T₁ + sqrt(T₁^2 + 4 * ℒ * qˡ₁ / cᵖᵐ)) / 2 - qˡ₂ = condensate_specific_humidity(T₂, state, ref, thermo) - r₂ = saturation_adjustment_residual(T₂, Π, qˡ₂, state, thermo) - - # Saturation adjustment - R = sqrt(max(T₂, T₁)) - ϵ = convert(FT, 1e-4) - δ = ϵ * R - iter = 0 - - while abs(r₂ - r₁) > δ - # Compute slope - ΔTΔr = (T₂ - T₁) / (r₂ - r₁) - - # Store previous values - r₁ = r₂ - T₁ = T₂ - - # Update - T₂ -= r₂ * ΔTΔr - qˡ₂ = condensate_specific_humidity(T₂, state, ref, thermo) - r₂ = saturation_adjustment_residual(T₂, Π, qˡ₂, state, thermo) - iter += 1 - end - - return T₂ -end - -@inline function saturation_adjustment_residual(T, Π, qˡ, state, thermo) - ℒᵛ = thermo.liquid.latent_heat - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - return T^2 - ℒᵛ * qˡ / cᵖᵐ - Π * state.θ * T -end - -@inline function specific_volume(state, ref, thermo) - T = temperature(state, ref, thermo) - Rᵐ = mixture_gas_constant(state.q, thermo) - pᵣ = reference_pressure(state.z, ref, thermo) - return Rᵐ * T / pᵣ -end diff --git a/test/runtests.jl b/test/runtests.jl index 441210c1..c29430e4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,66 +1,23 @@ +include("runtests_setup.jl") + using Test using Breeze using Oceananigans @testset "Breeze.jl" begin - @testset "Thermodynamics" begin - thermo = AtmosphereThermodynamics() - - # Test Saturation specific humidity calculation - T = 293.15 # 20°C - ρ = 1.2 # kg/m³ - q★ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, thermo, thermo.liquid) - @test q★ > 0 + @testset "Atmosphere Thermodynamics" begin + include("test_atmosphere_thermodynamics.jl") end - @testset "AtmosphereModel" begin - for FT in (Float32, Float64) - grid = RectilinearGrid(FT, size=(8, 8, 8), x=(0, 1_000), y=(0, 1_000), z=(0, 1_000)) - thermo = AtmosphereThermodynamics(FT) - - for p₀ in (101325, 100000) - for θ₀ in (288, 300) - constants = Breeze.Thermodynamics.ReferenceStateConstants(p₀, θ₀) - formulation = AnelasticFormulation(grid, constants, thermo) - model = AtmosphereModel(grid; formulation) - - # test set! - ρᵣ = model.formulation.reference_density - cᵖᵈ = model.thermodynamics.dry_air.heat_capacity - ρeᵢ = ρᵣ * cᵖᵈ * θ₀ - - set!(model; θ = θ₀) - ρe₁ = deepcopy(model.energy) - - set!(model; ρe = ρeᵢ) - @test model.energy == ρe₁ - end - end - end + @testset "Atmosphere Model Unit" begin + include("test_atmosphere_model_unit.jl") end - @testset "NonhydrostaticModel with MoistAirBuoyancy" begin - reference_constants = ReferenceStateConstants(potential_temperature=300) - buoyancy = MoistAirBuoyancy(; reference_constants) - - grid = RectilinearGrid(size=(8, 8, 8), x=(0, 400), y=(0, 400), z=(0, 400)) - model = NonhydrostaticModel(; grid, buoyancy, tracers = (:θ, :q)) - - θ₀ = reference_constants.reference_potential_temperature - Δθ = 2 - Lz = grid.Lz - - θᵢ(x, y, z) = θ₀ + Δθ * z / Lz - set!(model; θ = θᵢ, q = 0) - - # Can time-step - success = try - time_step!(model, 1e-2) - true - catch - false - end + @testset "Moist Air Buoyancy" begin + include("test_moist_air_buoyancy.jl") + end - @test success + @testset "Anelastic Pressure Solver" begin + include("test_anelastic_pressure_solver.jl") end end diff --git a/test/runtests_setup.jl b/test/runtests_setup.jl new file mode 100644 index 00000000..152812f6 --- /dev/null +++ b/test/runtests_setup.jl @@ -0,0 +1,3 @@ +using Test +using Breeze +using Oceananigans diff --git a/test/test_anelastic_pressure_solver.jl b/test/test_anelastic_pressure_solver.jl new file mode 100644 index 00000000..8ecdecd6 --- /dev/null +++ b/test/test_anelastic_pressure_solver.jl @@ -0,0 +1,87 @@ +include(joinpath(@__DIR__, "runtests_setup.jl")) + +using Test +using Breeze +using Oceananigans +using Oceananigans.Fields: fill_halo_regions! + +for FT in (Float32, Float64) + @testset "Pressure solver matches NonhydrostaticModel with ρᵣ == 1 [$FT]" begin + @info "Testing the anelastic pressure solver against NonhydrostaticModel [$FT]..." + Nx = Ny = Nz = 32 + z = 0:(1/Nz):1 + grid = RectilinearGrid(FT; size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z) + thermodynamics = ThermodynamicConstants(FT) + reference_state = ReferenceState(FT; base_pressure=101325, potential_temperature=288) + + formulation = AnelasticFormulation(grid, reference_state, thermodynamics) + parent(formulation.reference_density) .= 1 + + anelastic = AtmosphereModel(grid; thermodynamics=thermodynamics, formulation) + boussinesq = NonhydrostaticModel(; grid) + + uᵢ = rand(size(grid)...) + vᵢ = rand(size(grid)...) + wᵢ = rand(size(grid)...) + + set!(anelastic, ρu=uᵢ, ρv=vᵢ, ρw=wᵢ) + set!(boussinesq, u=uᵢ, v=vᵢ, w=wᵢ) + + ρu = anelastic.momentum.ρu + ρv = anelastic.momentum.ρv + ρw = anelastic.momentum.ρw + δᵃ = Field(∂x(ρu) + ∂y(ρv) + ∂z(ρw)) + + u = boussinesq.velocities.u + v = boussinesq.velocities.v + w = boussinesq.velocities.w + δᵇ = Field(∂x(u) + ∂y(v) + ∂z(w)) + + boussinesq_solver = boussinesq.pressure_solver + anelastic_solver = anelastic.pressure_solver + @test anelastic_solver.batched_tridiagonal_solver.a == boussinesq_solver.batched_tridiagonal_solver.a + @test anelastic_solver.batched_tridiagonal_solver.b == boussinesq_solver.batched_tridiagonal_solver.b + @test anelastic_solver.batched_tridiagonal_solver.c == boussinesq_solver.batched_tridiagonal_solver.c + @test anelastic_solver.source_term == boussinesq_solver.source_term + + @test maximum(abs, δᵃ) < prod(size(grid)) * eps(FT) + @test maximum(abs, δᵇ) < prod(size(grid)) * eps(FT) + @test anelastic.nonhydrostatic_pressure == boussinesq.pressures.pNHS + end + + @testset "Anelastic pressure solver recovers analytic solution [$FT]" begin + @info "Test that anelastic pressure solver recovers analytic solution [$FT]..." + grid = RectilinearGrid(FT; size=48, z=(0, 1), topology=(Flat, Flat, Bounded)) + thermodynamics = ThermodynamicConstants(FT) + reference_state = ReferenceState(FT; base_pressure=101325.0, potential_temperature=288.0) + formulation = AnelasticFormulation(grid, reference_state, thermodynamics) + + #= + ρᵣ = 2 + cos(π z / 2) + ∂z ρᵣ ∂z ϕ = ? + + ϕ = cos(π z) + ⟹ ∂z ϕ = -π sin(π z) + ⟹ (2 + cos(π z)) ∂z ϕ = -π (2 sin(π z) + cos(π z) sin(π z)) + ⟹ ∂z (1 + cos(π z / 2)) ∂z ϕ = -π² (2 cos(π z) + 2 cos²(π z) - 1) + + ϕ = z² / 2 - z³ / 3 = z² (1/2 - z/3) + ∂z ϕ = z (1 - z) = z - z² + ∂z² ϕ = 1 - 2z + ⟹ z ∂z ϕ = z² - z³ + ⟹ ∂z (z ∂z ϕ) = 2 z - 3 z² + + R = ∂z ρw = 2 z - 3 z² + ⟹ ρw = z² - z³ + =# + + set!(formulation.reference_density, z -> z) + fill_halo_regions!(formulation.reference_density) + model = AtmosphereModel(grid; thermodynamics=thermodynamics, formulation) + set!(model, ρw = z -> z^2 - z^3) + + ϕ_exact = CenterField(grid) + set!(ϕ_exact, z -> z^2 / 2 - z^3 / 3 - 1 / 12) + @test isapprox(ϕ_exact, model.nonhydrostatic_pressure, rtol=1e-3) + end +end diff --git a/test/test_atmosphere_model_unit.jl b/test/test_atmosphere_model_unit.jl new file mode 100644 index 00000000..05612773 --- /dev/null +++ b/test/test_atmosphere_model_unit.jl @@ -0,0 +1,27 @@ +include(joinpath(@__DIR__, "runtests_setup.jl")) + +@testset "AtmosphereModel" begin + for FT in (Float32, Float64) + grid = RectilinearGrid(FT, size=(8, 8, 8), x=(0, 1_000), y=(0, 1_000), z=(0, 1_000)) + thermodynamics = ThermodynamicConstants(FT) + + for p₀ in (101325, 100000) + for θ₀ in (288, 300) + reference_state = Breeze.Thermodynamics.ReferenceState(base_pressure=p₀, + potential_temperature=θ₀) + formulation = AnelasticFormulation(grid, reference_state, thermodynamics) + model = AtmosphereModel(grid; formulation) + + ρᵣ = model.formulation.reference_density + cᵖᵈ = model.thermodynamics.dry_air.heat_capacity + ρeᵢ = ρᵣ * cᵖᵈ * θ₀ + + set!(model; θ = θ₀) + ρe₁ = deepcopy(model.energy) + + set!(model; ρe = ρeᵢ) + @test_broken model.energy == ρe₁ + end + end + end +end diff --git a/test/test_atmosphere_thermodynamics.jl b/test/test_atmosphere_thermodynamics.jl new file mode 100644 index 00000000..78d03a9b --- /dev/null +++ b/test/test_atmosphere_thermodynamics.jl @@ -0,0 +1,10 @@ +include(joinpath(@__DIR__, "runtests_setup.jl")) + +@testset "Thermodynamics" begin + thermo = ThermodynamicConstants() + + T = 293.15 # 20°C + ρ = 1.2 # kg/m³ + q★ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, thermo, thermo.liquid) + @test q★ > 0 +end diff --git a/test/test_moist_air_buoyancy.jl b/test/test_moist_air_buoyancy.jl new file mode 100644 index 00000000..e1f83bf3 --- /dev/null +++ b/test/test_moist_air_buoyancy.jl @@ -0,0 +1,25 @@ +include(joinpath(@__DIR__, "runtests_setup.jl")) + +@testset "NonhydrostaticModel with MoistAirBuoyancy" begin + reference_state = ReferenceState(potential_temperature=300) + buoyancy = MoistAirBuoyancy(; reference_state) + + grid = RectilinearGrid(size=(8, 8, 8), x=(0, 400), y=(0, 400), z=(0, 400)) + model = NonhydrostaticModel(; grid, buoyancy, tracers = (:θ, :q)) + + θ₀ = reference_state.potential_temperature + Δθ = 2 + Lz = grid.Lz + + θᵢ(x, y, z) = θ₀ + Δθ * z / Lz + set!(model; θ = θᵢ, q = 0) + + success = try + time_step!(model, 1e-2) + true + catch + false + end + + @test success +end