diff --git a/docs/src/microphysics/saturation_adjustment.md b/docs/src/microphysics/saturation_adjustment.md index 166834aa..e88fb7a1 100644 --- a/docs/src/microphysics/saturation_adjustment.md +++ b/docs/src/microphysics/saturation_adjustment.md @@ -34,7 +34,7 @@ The saturation specific humidity is then using Breeze using Breeze.MoistAirBuoyancies: saturation_specific_humidity, HeightReferenceThermodynamicState -thermo = AtmosphereThermodynamics() +thermo = ThermodynamicConstants() ref = ReferenceStateConstants(base_pressure=101325, potential_temperature=288) z = 0.0 # [m] height diff --git a/docs/src/thermodynamics.md b/docs/src/thermodynamics.md index 731639e3..e9317c28 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 -- @@ -148,7 +148,7 @@ using Breeze using Breeze.Thermodynamics: reference_pressure, reference_density using CairoMakie -thermo = AtmosphereThermodynamics() +thermo = ThermodynamicConstants() constants = ReferenceStateConstants(base_pressure=101325, potential_temperature=288) grid = RectilinearGrid(size=160, z=(0, 12_000), topology=(Flat, Flat, Bounded)) @@ -212,7 +212,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: @@ -298,7 +298,7 @@ The saturation vapor pressure is using Breeze using Breeze.Thermodynamics: saturation_vapor_pressure -thermo = AtmosphereThermodynamics() +thermo = ThermodynamicConstants() T = collect(200:0.1:320) pᵛˡ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, thermo.liquid) for Tⁱ in T] @@ -327,7 +327,7 @@ and this is what it looks like: using Breeze using Breeze.MoistAirBuoyancies: saturation_specific_humidity -thermo = AtmosphereThermodynamics() +thermo = ThermodynamicConstants() ref = ReferenceStateConstants(base_pressure=101325, potential_temperature=288) z = 0 diff --git a/examples/JRA55_saturation_specific_humidity.jl b/examples/JRA55_saturation_specific_humidity.jl index d22355d5..a584f449 100644 --- a/examples/JRA55_saturation_specific_humidity.jl +++ b/examples/JRA55_saturation_specific_humidity.jl @@ -1,11 +1,11 @@ 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)) diff --git a/examples/large_yeager_saturation_specific_humidity.jl b/examples/large_yeager_saturation_specific_humidity.jl index e0c50d65..d8a6b521 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 = Breeze.ThermodynamicConstants() saturation_specific_humidity_large_yeager(T, ρ) = 640380 * exp(-5107.4 / T) / ρ diff --git a/src/AtmosphereModels/anelastic_formulation.jl b/src/AtmosphereModels/anelastic_formulation.jl index 4a1ee40c..7ff6e114 100644 --- a/src/AtmosphereModels/anelastic_formulation.jl +++ b/src/AtmosphereModels/anelastic_formulation.jl @@ -1,5 +1,5 @@ using ..Thermodynamics: - AtmosphereThermodynamics, + ThermodynamicConstants, ReferenceStateConstants, reference_pressure, reference_density, diff --git a/src/AtmosphereModels/atmosphere_model.jl b/src/AtmosphereModels/atmosphere_model.jl index b03238e5..93dc9ba6 100644 --- a/src/AtmosphereModels/atmosphere_model.jl +++ b/src/AtmosphereModels/atmosphere_model.jl @@ -1,5 +1,5 @@ using ..Thermodynamics: - AtmosphereThermodynamics, + ThermodynamicConstants, ReferenceStateConstants, reference_pressure, reference_density, @@ -69,7 +69,7 @@ end """ AtmosphereModel(grid; clock = Clock(grid), - thermodynamics = AtmosphereThermodynamics(eltype(grid)), + thermodynamics = ThermodynamicConstants(eltype(grid)), formulation = default_formulation(grid, thermodynamics), absolute_humidity = DefaultValue(), tracers = tuple(), @@ -109,7 +109,7 @@ Pauluis, O. (2008). Thermodynamic consistency of the anelastic approximation for """ function AtmosphereModel(grid; clock = Clock(grid), - thermodynamics = AtmosphereThermodynamics(eltype(grid)), + thermodynamics = ThermodynamicConstants(eltype(grid)), formulation = default_formulation(grid, thermodynamics), absolute_humidity = DefaultValue(), tracers = tuple(), diff --git a/src/Breeze.jl b/src/Breeze.jl index 548a9180..6c517c43 100644 --- a/src/Breeze.jl +++ b/src/Breeze.jl @@ -7,7 +7,7 @@ module Breeze export MoistAirBuoyancy, - AtmosphereThermodynamics, + ThermodynamicConstants, ReferenceStateConstants, AnelasticFormulation, AtmosphereModel, diff --git a/src/MoistAirBuoyancies.jl b/src/MoistAirBuoyancies.jl index 0b579373..cc7a9ab4 100644 --- a/src/MoistAirBuoyancies.jl +++ b/src/MoistAirBuoyancies.jl @@ -19,7 +19,7 @@ import Oceananigans.BuoyancyFormulations: AbstractBuoyancyFormulation, required_tracers using ..Thermodynamics: - AtmosphereThermodynamics, + ThermodynamicConstants, ReferenceStateConstants, mixture_heat_capacity, mixture_gas_constant, @@ -38,7 +38,7 @@ end """ MoistAirBuoyancy(FT=Oceananigans.defaults.FloatType; - thermodynamics = AtmosphereThermodynamics(FT), + thermodynamics = ThermodynamicConstants(FT), reference_constants = ReferenceStateConstants{FT}(101325, 290)) Return a MoistAirBuoyancy formulation that can be provided as input to an `AtmosphereModel` @@ -56,7 +56,7 @@ julia> using Breeze, Oceananigans julia> buoyancy = MoistAirBuoyancy() MoistAirBuoyancy ├── reference_constants: Breeze.Thermodynamics.ReferenceStateConstants{Float64} -└── thermodynamics: AtmosphereThermodynamics +└── thermodynamics: ThermodynamicConstants julia> model = NonhydrostaticModel(; grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 2, 3)), buoyancy, tracers = (:θ, :q)) @@ -71,7 +71,7 @@ NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ``` """ function MoistAirBuoyancy(FT=Oceananigans.defaults.FloatType; - thermodynamics = AtmosphereThermodynamics(FT), + thermodynamics = ThermodynamicConstants(FT), reference_constants = ReferenceStateConstants{FT}(101325, 290)) AT = typeof(thermodynamics) diff --git a/src/Thermodynamics/Thermodynamics.jl b/src/Thermodynamics/Thermodynamics.jl index 16d6fc5e..c0ad354c 100644 --- a/src/Thermodynamics/Thermodynamics.jl +++ b/src/Thermodynamics/Thermodynamics.jl @@ -1,6 +1,6 @@ module Thermodynamics -export AtmosphereThermodynamics, ReferenceStateConstants, IdealGas, +export ThermodynamicConstants, ReferenceStateConstants, IdealGas, PhaseTransitionConstants, CondensedPhase, mixture_gas_constant, mixture_heat_capacity diff --git a/src/Thermodynamics/atmosphere_thermodynamics.jl b/src/Thermodynamics/atmosphere_thermodynamics.jl index 4b706b03..38d37919 100644 --- a/src/Thermodynamics/atmosphere_thermodynamics.jl +++ b/src/Thermodynamics/atmosphere_thermodynamics.jl @@ -91,7 +91,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 @@ -103,9 +103,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", @@ -118,9 +118,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) @@ -133,38 +133,37 @@ 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, - gravitational_acceleration, - energy_reference_temperature, - triple_point_temperature, - triple_point_pressure, - dry_air, - vapor, - liquid, - solid) + return ThermodynamicConstants{FT, C, S}(molar_gas_constant, + gravitational_acceleration, + energy_reference_temperature, + triple_point_temperature, + triple_point_pressure, + dry_air, + vapor, + liquid, + solid) end """ - AtmosphereThermodynamics(FT = Oceananigans.defaults.FloatType; - gravitational_acceleration = 9.81, - molar_gas_constant = 8.314462618, - energy_reference_temperature = 273.16, - triple_point_temperature = 273.16, - triple_point_pressure = 611.657, - dry_air_molar_mass = 0.02897, - dry_air_heat_capacity = 1005, - vapor_molar_mass = 0.018015, - vapor_heat_capacity = 1850, - liquid = liquid_water(FT), - solid = water_ice(FT), - condensed_phases = nothing) - -Create `AtmosphereThermodynamics` with parameters that represent gaseous mixture of dry "air" + ThermodynamicConstants(FT = Oceananigans.defaults.FloatType; + molar_gas_constant = 8.314462618, + gravitational_acceleration = 9.81, + energy_reference_temperature = 273.16, + triple_point_temperature = 273.16, + triple_point_pressure = 611.657, + dry_air_molar_mass = 0.02897, + dry_air_heat_capacity = 1005, + vapor_molar_mass = 0.018015, + vapor_heat_capacity = 1850, + liquid = liquid_water(FT), + solid = water_ice(FT)) + +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 at the boundary between vapor and a homogeneous sample of the condensed phase. -The `gravitational_acceleration` parameter is included to compute reference_state +The `gravitational_acceleration` parameter is included to compute `reference_state` quantities associated with hydrostatic balance. The Clausius-Clapeyron relation describes the pressure-temperature relationship during phase @@ -198,18 +197,18 @@ Note: any reference values for pressure and temperature can be used in principle The advantage of using reference values at the triple point is that the same values can then be used for both condensation (vapor → liquid) and deposition (vapor → ice). """ -function AtmosphereThermodynamics(FT = Oceananigans.defaults.FloatType; - molar_gas_constant = 8.314462618, - gravitational_acceleration = 9.81, - energy_reference_temperature = 273.16, - triple_point_temperature = 273.16, - triple_point_pressure = 611.657, - dry_air_molar_mass = 0.02897, - dry_air_heat_capacity = 1005, - vapor_molar_mass = 0.018015, - vapor_heat_capacity = 1850, - liquid = liquid_water(FT), - solid = water_ice(FT)) +function ThermodynamicConstants(FT = Oceananigans.defaults.FloatType; + molar_gas_constant = 8.314462618, + gravitational_acceleration = 9.81, + energy_reference_temperature = 273.16, + triple_point_temperature = 273.16, + triple_point_pressure = 611.657, + dry_air_molar_mass = 0.02897, + dry_air_heat_capacity = 1005, + vapor_molar_mass = 0.018015, + vapor_heat_capacity = 1850, + liquid = liquid_water(FT), + solid = water_ice(FT)) dry_air = IdealGas(FT; molar_mass = dry_air_molar_mass, heat_capacity = dry_air_heat_capacity) @@ -217,24 +216,24 @@ 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 TC = 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 +@inline vapor_gas_constant(thermo::TC) = thermo.molar_gas_constant / thermo.vapor.molar_mass +@inline dry_air_gas_constant(thermo::TC) = thermo.molar_gas_constant / thermo.dry_air.molar_mass -const NonCondensingAtmosphereThermodynamics{FT} = AtmosphereThermodynamics{FT, Nothing, Nothing} +const NonCondensingThermodynamicConstants{FT} = ThermodynamicConstants{FT, Nothing, Nothing} """ mixture_gas_constant(q, thermo) @@ -256,12 +255,12 @@ where: # Arguments - `q`: Specific humidity (dimensionless) -- `thermo`: `AtmosphereThermodynamics` instance containing gas constants +- `thermo`: `ThermodynamicConstants` instance containing gas constants # 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, thermo::TC) Rᵈ = dry_air_gas_constant(thermo) Rᵛ = vapor_gas_constant(thermo) return Rᵈ * (1 - q) + Rᵛ * q @@ -274,7 +273,7 @@ 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, thermo::TC) cᵖᵈ = thermo.dry_air.heat_capacity cᵖᵛ = thermo.vapor.heat_capacity return cᵖᵈ * (1 - q) + cᵖᵛ * q diff --git a/test/runtests.jl b/test/runtests.jl index 441210c1..74b53db3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,7 @@ using Oceananigans @testset "Breeze.jl" begin @testset "Thermodynamics" begin - thermo = AtmosphereThermodynamics() + thermo = ThermodynamicConstants() # Test Saturation specific humidity calculation T = 293.15 # 20°C @@ -16,7 +16,7 @@ using Oceananigans @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) + thermo = ThermodynamicConstants(FT) for p₀ in (101325, 100000) for θ₀ in (288, 300) 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..c267846b --- /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_constants = ReferenceStateConstants(FT; base_pressure=101325, potential_temperature=288) + + formulation = AnelasticFormulation(grid, reference_constants, 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_constants = ReferenceStateConstants(FT; base_pressure=101325.0, potential_temperature=288.0) + formulation = AnelasticFormulation(grid, reference_constants, 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, 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..1950e7cd --- /dev/null +++ b/test/test_atmosphere_model_unit.jl @@ -0,0 +1,26 @@ +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_constants = ReferenceStateConstants(FT, base_pressure=p₀, potential_temperature=θ₀) + formulation = AnelasticFormulation(grid, reference_constants, 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 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..3ffe6ac2 --- /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_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) + + success = try + time_step!(model, 1e-2) + true + catch + false + end + + @test success +end