diff --git a/docs/make.jl b/docs/make.jl index cb1152ef..b34ba2ec 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -18,6 +18,7 @@ mkpath(literated_dir) example_scripts = [ "thermal_bubble.jl", + "moist_thermal_bubble.jl", "cloudy_kelvin_helmholtz.jl", # "prescribed_sst.jl", # this is a WIP ] @@ -31,6 +32,7 @@ end example_pages = Any[ "Thermal bubble" => "literated/thermal_bubble.md", + "Moist and dry thermal bubbles" => "literated/moist_thermal_bubble.md", "Cloudy Kelvin-Helmholtz instability" => "literated/cloudy_kelvin_helmholtz.md", # "Prescribed SST" => "literated/prescribed_sst.md", ] diff --git a/docs/src/appendix/notation.md b/docs/src/appendix/notation.md index 5f7e8085..be243424 100644 --- a/docs/src/appendix/notation.md +++ b/docs/src/appendix/notation.md @@ -17,7 +17,7 @@ A few notes about the following table: a "reference state", which is an adiabatic, hydrostatic solution to the equations of motion. But there is also an "energy reference temperature" and "reference latent heat", which are thermodynamic constants required to define the internal energy of moist atmospheric constituents. -* Mapping to AM fields: `ρe` corresponds to `model.energy_density`, `ρqᵗ` to `model.moisture_density`, and `qᵗ` to `model.specific_moisture`. +* Mapping to AM fields: `ρe` corresponds to `energy_density(model)`, `ρqᵗ` to `model.moisture_density`, and `qᵗ` to `model.specific_moisture`. The following table also uses a few conventions that suffuse the source code and which are internalized by wise developers: diff --git a/docs/src/microphysics/mixed_phase_saturation_adjustment.md b/docs/src/microphysics/mixed_phase_saturation_adjustment.md index ed0af874..044e0f19 100644 --- a/docs/src/microphysics/mixed_phase_saturation_adjustment.md +++ b/docs/src/microphysics/mixed_phase_saturation_adjustment.md @@ -47,7 +47,7 @@ We can compute the saturation specific humidity across the entire range from hom ice nucleation up to freezing using the equilibrium model above: ```@example mixed_phase -using Breeze.Microphysics: adjustment_saturation_specific_humidity +using Breeze.Microphysics: equilibrium_saturation_specific_humidity p = 101325.0 qᵗ = 0.012 @@ -56,7 +56,7 @@ Tᶠ = eq.freezing_temperature Tʰ = eq.homogeneous_ice_nucleation_temperature T = range(Tʰ - 10, Tᶠ + 10; length=81) # slightly beyond the mixed-phase range -qᵛ⁺ = [adjustment_saturation_specific_humidity(Tʲ, p, qᵗ, thermo, eq) for Tʲ in T] +qᵛ⁺ = [equilibrium_saturation_specific_humidity(Tʲ, p, qᵗ, thermo, eq) for Tʲ in T] ``` Optionally, we can visualize how `qᵛ⁺` varies with temperature: diff --git a/docs/src/microphysics/warm_phase_saturation_adjustment.md b/docs/src/microphysics/warm_phase_saturation_adjustment.md index 7bc8f6b5..a434b6d4 100644 --- a/docs/src/microphysics/warm_phase_saturation_adjustment.md +++ b/docs/src/microphysics/warm_phase_saturation_adjustment.md @@ -101,10 +101,10 @@ Next, we compute the saturation specific humidity for moist air with a carefully chosen moist air mass fraction, ```@example microphysics -using Breeze.Microphysics: adjustment_saturation_specific_humidity, WarmPhaseEquilibrium +using Breeze.Microphysics: equilibrium_saturation_specific_humidity, WarmPhaseEquilibrium qᵗ = 0.012 # [kg kg⁻¹] total specific humidity -qᵛ⁺ = Breeze.Microphysics.adjustment_saturation_specific_humidity(T, p, qᵗ, thermo, WarmPhaseEquilibrium()) +qᵛ⁺ = Breeze.Microphysics.equilibrium_saturation_specific_humidity(T, p, qᵗ, thermo, WarmPhaseEquilibrium()) ``` There are two facts of note. First is that we have identified a situation in which ``qᵗ > qᵛ⁺``, @@ -164,7 +164,7 @@ To generate a second guess for the secant solver, we start by estimating the liquid mass fraction using the guess ``T = T₁``, ```@example microphysics -qᵛ⁺₂ = adjustment_saturation_specific_humidity(T₁, p, qᵗ, thermo, WarmPhaseEquilibrium()) +qᵛ⁺₂ = equilibrium_saturation_specific_humidity(T₁, p, qᵗ, thermo, WarmPhaseEquilibrium()) qˡ₁ = qᵗ - qᵛ⁺₂ ``` @@ -189,7 +189,7 @@ using CairoMakie equilibrium = WarmPhaseEquilibrium() T = 230:0.5:320 r = [saturation_adjustment_residual(Tʲ, 𝒰, thermo, equilibrium) for Tʲ in T] -qᵛ⁺ = [adjustment_saturation_specific_humidity(Tʲ, p, qᵗ, thermo, equilibrium) for Tʲ in T] +qᵛ⁺ = [equilibrium_saturation_specific_humidity(Tʲ, p, qᵗ, thermo, equilibrium) for Tʲ in T] fig = Figure() axr = Axis(fig[1, 1], xlabel="Temperature (K)", ylabel="Saturation adjustment residual (K)") @@ -228,7 +228,7 @@ for (i, qᵗⁱ) in enumerate(qᵗ) q = MoistureMassFractions(qᵗⁱ) 𝒰 = StaticEnergyState(e₀, q, z, p) T[i] = compute_temperature(𝒰, microphysics, thermo) - qᵛ⁺ = Breeze.Microphysics.adjustment_saturation_specific_humidity(T[i], p, qᵗⁱ, thermo, WarmPhaseEquilibrium()) + qᵛ⁺ = Breeze.Microphysics.equilibrium_saturation_specific_humidity(T[i], p, qᵗⁱ, thermo, WarmPhaseEquilibrium()) qˡ[i] = max(0, qᵗⁱ - qᵛ⁺) end @@ -295,7 +295,7 @@ for k = 1:grid.Nz T[k] = compute_temperature(𝒰, microphysics, thermo) # Saturation specific humidity via adjustment formula using T[k], pᵣ, and qᵗ - qᵛ⁺[k] = Breeze.Microphysics.adjustment_saturation_specific_humidity(T[k], pᵣ, qᵗ, thermo, WarmPhaseEquilibrium()) + qᵛ⁺[k] = Breeze.Microphysics.equilibrium_saturation_specific_humidity(T[k], pᵣ, qᵗ, thermo, WarmPhaseEquilibrium()) qˡ[k] = max(0, qᵗ - qᵛ⁺[k]) rh[k] = 100 * min(qᵗ, qᵛ⁺[k]) / qᵛ⁺[k] end diff --git a/examples/anelastic_bomex.jl b/examples/anelastic_bomex.jl index 7524d62e..11b457e8 100644 --- a/examples/anelastic_bomex.jl +++ b/examples/anelastic_bomex.jl @@ -32,7 +32,8 @@ u_bomex = AtmosphericProfilesLibrary.Bomex_u(FT) p₀, θ₀ = 101500, 299.1 constants = ThermodynamicConstants() reference_state = ReferenceState(grid, constants, base_pressure=p₀, potential_temperature=θ₀) -formulation = AnelasticFormulation(reference_state) +formulation = AnelasticFormulation(reference_state, thermodynamics=:LiquidIcePotentialTemperature) +# formulation = AnelasticFormulation(reference_state, thermodynamics=:StaticEnergy) q₀ = Breeze.Thermodynamics.MoistureMassFractions{eltype(grid)} |> zero ρ₀ = Breeze.Thermodynamics.density(p₀, θ₀, q₀, constants) @@ -41,7 +42,9 @@ Lˡ = constants.liquid.reference_latent_heat w′T′, w′q′ = 8e-3, 5.2e-5 Q = ρ₀ * cᵖᵈ * w′T′ F = ρ₀ * w′q′ + ρe_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(Q)) +ρθ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(ρ₀ * w′T′)) ρqᵗ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(F)) u★ = 0.28 # m/s @@ -69,6 +72,11 @@ end return @inbounds - p.ρᵣ[i, j, k] * w_dz_E end +@inline function Fρθ_subsidence(i, j, k, grid, clock, fields, p) + w_dz_T = ℑzᵃᵃᶜ(i, j, k, grid, w_dz_ϕ, p.wˢ, p.θ_avg) + return @inbounds - p.ρᵣ[i, j, k] * w_dz_T +end + @inline function Fρqᵗ_subsidence(i, j, k, grid, clock, fields, p) w_dz_Qᵗ = ℑzᵃᵃᶜ(i, j, k, grid, w_dz_ϕ, p.wˢ, p.qᵗ_avg) return @inbounds - p.ρᵣ[i, j, k] * w_dz_Qᵗ @@ -78,6 +86,7 @@ end u_avg_f = Field{Nothing, Nothing, Center}(grid) v_avg_f = Field{Nothing, Nothing, Center}(grid) e_avg_f = Field{Nothing, Nothing, Center}(grid) +θ_avg_f = Field{Nothing, Nothing, Center}(grid) qᵗ_avg_f = Field{Nothing, Nothing, Center}(grid) wˢ = Field{Nothing, Nothing, Face}(grid) @@ -88,17 +97,24 @@ set!(wˢ, z -> w_bomex(z)) ρu_subsidence_forcing = Forcing(Fρu_subsidence, discrete_form=true, parameters=(; u_avg=u_avg_f, wˢ, ρᵣ)) ρv_subsidence_forcing = Forcing(Fρv_subsidence, discrete_form=true, parameters=(; v_avg=v_avg_f, wˢ, ρᵣ)) ρe_subsidence_forcing = Forcing(Fρe_subsidence, discrete_form=true, parameters=(; e_avg=e_avg_f, wˢ, ρᵣ)) +ρθ_subsidence_forcing = Forcing(Fρθ_subsidence, discrete_form=true, parameters=(; θ_avg=θ_avg_f, wˢ, ρᵣ)) ρqᵗ_subsidence_forcing = Forcing(Fρqᵗ_subsidence, discrete_form=true, parameters=(; qᵗ_avg=qᵗ_avg_f, wˢ, ρᵣ)) coriolis = FPlane(f=3.76e-5) -ρuᵍ = Field{Nothing, Nothing, Center}(grid) -ρvᵍ = Field{Nothing, Nothing, Center}(grid) +# ρuᵍ = Field{Nothing, Nothing, Center}(grid) +# ρvᵍ = Field{Nothing, Nothing, Center}(grid) +uᵍ = Field{Nothing, Nothing, Center}(grid) +vᵍ = Field{Nothing, Nothing, Center}(grid) uᵍ_bomex = AtmosphericProfilesLibrary.Bomex_geostrophic_u(FT) vᵍ_bomex = AtmosphericProfilesLibrary.Bomex_geostrophic_v(FT) -set!(ρuᵍ, z -> uᵍ_bomex(z)) -set!(ρvᵍ, z -> vᵍ_bomex(z)) -set!(ρuᵍ, ρᵣ * ρuᵍ) -set!(ρvᵍ, ρᵣ * ρvᵍ) +set!(uᵍ, z -> uᵍ_bomex(z)) +set!(vᵍ, z -> vᵍ_bomex(z)) + +ρuᵍ = Field(ρᵣ * uᵍ) +ρvᵍ = Field(ρᵣ * vᵍ) + +# set!(ρuᵍ, ρᵣ * ρuᵍ) +# set!(ρvᵍ, ρᵣ * ρvᵍ) @inline Fρu_geostrophic(i, j, k, grid, clock, fields, p) = @inbounds - p.f * p.ρvᵍ[i, j, k] @inline Fρv_geostrophic(i, j, k, grid, clock, fields, p) = @inbounds p.f * p.ρuᵍ[i, j, k] @@ -117,13 +133,19 @@ set!(drying, ρᵣ * drying) ρqᵗ_forcing = (ρqᵗ_drying_forcing, ρqᵗ_subsidence_forcing) Fρe_field = Field{Nothing, Nothing, Center}(grid) +Fρθ_field = Field{Nothing, Nothing, Center}(grid) dTdt_bomex = AtmosphericProfilesLibrary.Bomex_dTdt(FT) set!(Fρe_field, z -> dTdt_bomex(1, z)) set!(Fρe_field, ρᵣ * cᵖᵈ * Fρe_field) +set!(Fρθ_field, z -> dTdt_bomex(1, z)) +set!(Fρθ_field, ρᵣ * Fρe_field) ρe_radiation_forcing = Forcing(Fρe_field) ρe_forcing = (ρe_radiation_forcing, ρe_subsidence_forcing) +ρθ_radiation_forcing = Forcing(Fρθ_field) +ρθ_forcing = (ρθ_radiation_forcing, ρθ_subsidence_forcing) + fig = Figure() axe = Axis(fig[1, 1], xlabel="z (m)", ylabel="Fρe (K/s)") axq = Axis(fig[1, 2], xlabel="z (m)", ylabel="Fρqᵗ (1/s)") @@ -144,8 +166,10 @@ closure = AnisotropicMinimumDissipation() model = AtmosphereModel(grid; formulation, coriolis, microphysics, closure, scalar_advection = Centered(order=2), momentum_advection = WENO(order=9), - forcing = (ρqᵗ=ρqᵗ_forcing, ρu=ρu_forcing, ρv=ρv_forcing, ρe=ρe_forcing), - boundary_conditions = (ρe=ρe_bcs, ρqᵗ=ρqᵗ_bcs, ρu=ρu_bcs, ρv=ρv_bcs)) + forcing = (ρqᵗ=ρqᵗ_forcing, ρu=ρu_forcing, ρv=ρv_forcing, ρθ=ρθ_forcing), + boundary_conditions = (ρθ=ρθ_bcs, ρqᵗ=ρqᵗ_bcs, ρu=ρu_bcs, ρv=ρv_bcs)) + # forcing = (ρqᵗ=ρqᵗ_forcing, ρu=ρu_forcing, ρv=ρv_forcing, ρe=ρe_forcing), + # boundary_conditions = (ρe=ρe_bcs, ρqᵗ=ρqᵗ_bcs, ρu=ρu_bcs, ρv=ρv_bcs)) # Values for the initial perturbations can be found in Appendix B # of Siebesma et al 2003, 3rd paragraph @@ -159,32 +183,34 @@ simulation = Simulation(model; Δt=10, stop_time) conjure_time_step_wizard!(simulation, cfl=0.7) # Write a callback to compute *_avg_f -ρe = energy_density(model) +e = static_energy(model) +θ = liquid_ice_potential_temperature(model) u_avg = Field(Average(model.velocities.u, dims=(1, 2))) v_avg = Field(Average(model.velocities.v, dims=(1, 2))) -e_avg = Field(Average(ρe / ρᵣ, dims=(1, 2))) +e_avg = Field(Average(e, dims=(1, 2))) +θ_avg = Field(Average(θ, dims=(1, 2))) qᵗ_avg = Field(Average(model.specific_moisture, dims=(1, 2))) function compute_averages!(sim) compute!(u_avg) compute!(v_avg) compute!(e_avg) + compute!(θ_avg) compute!(qᵗ_avg) parent(u_avg_f) .= parent(u_avg) parent(v_avg_f) .= parent(v_avg) parent(e_avg_f) .= parent(e_avg) + parent(θ_avg_f) .= parent(θ_avg) parent(qᵗ_avg_f) .= parent(qᵗ_avg) return nothing end add_callback!(simulation, compute_averages!) -θ = Breeze.AtmosphereModels.PotentialTemperatureField(model) qˡ = model.microphysical_fields.qˡ qᵛ = model.microphysical_fields.qᵛ qᵗ = model.specific_moisture -qᵛ⁺ = Breeze.AtmosphereModels.SaturationSpecificHumidityField(model) -θ_avg = Average(θ, dims=(1, 2)) |> Field +qᵛ⁺ = SaturationSpecificHumidityField(model) qˡ_avg = Average(qˡ, dims=(1, 2)) |> Field fig = Figure() @@ -212,7 +238,8 @@ ax_ρe = Axis(fig_lower[1, 1], xlabel="Energy density", ylabel="z (m)") ax_e = Axis(fig_lower[1, 2], xlabel="Specific energy", ylabel="z (m)") ax_θ = Axis(fig_lower[1, 3], xlabel="Potential temperature (K)", ylabel="z (m)") -e = specific_energy(model) +e = static_energy(model) +ρe = ρᵣ * e ρe_avg = Average(ρe, dims=(1, 2)) |> Field e_avg = Average(e, dims=(1, 2)) |> Field @@ -247,7 +274,7 @@ function progress(sim) qᵗ = sim.model.specific_moisture qᵗmax = maximum(qᵗ) - ρe = energy_density(sim.model) + ρe = static_energy_density(sim.model) ρemin = minimum(ρe) ρemax = maximum(ρe) diff --git a/examples/atmos_convection.jl b/examples/atmos_convection.jl index 2da9a9d8..85ed250d 100644 --- a/examples/atmos_convection.jl +++ b/examples/atmos_convection.jl @@ -21,7 +21,7 @@ set!(model, θ=θᵢ, u=Ξᵢ, v=Ξᵢ) simulation = Simulation(model, Δt=10, stop_time=20minutes) conjure_time_step_wizard!(simulation, cfl=0.7) -ρe = energy_density(model) +ρe = static_energy_density(model) ∫ρe = Field(Integral(ρe)) function progress(sim) diff --git a/examples/cloudy_kelvin_helmholtz.jl b/examples/cloudy_kelvin_helmholtz.jl index 2c5bdcaf..07b2ae09 100644 --- a/examples/cloudy_kelvin_helmholtz.jl +++ b/examples/cloudy_kelvin_helmholtz.jl @@ -28,11 +28,8 @@ using Printf # Grid resolution is modest but enough to clearly resolve the Kelvin-Helmholtz billows and # rolled-up moisture filament. -Nx = 384 # horizontal resolution -Nz = 128 # vertical resolution - -Lx = 10e3 # domain length -Lz = 3e3 # domain height +Nx, Nz = 384, 128 # resolution +Lx, Lz = 10e3, 3e3 # domain extent grid = RectilinearGrid(; size = (Nx, Nz), x = (0, Lx), z = (0, Lz), topology = (Periodic, Flat, Bounded)) @@ -53,14 +50,14 @@ model = AtmosphereModel(grid; advection=WENO(order=5), microphysics) # N² = \frac{g}{θ₀} \frac{∂θ}{∂z} , # ``` # -# We initialize the potential temperature that gives constant Brunt–Väisälä frequency, +# We initialize with a potential temperature that gives constant Brunt–Väisälä frequency, # representative of mid-tropospheric stability. The (dry) Brunt–Väisälä frequency is # # ```math # N² = \frac{g}{θ} \frac{∂θ}{∂z} # ``` # -# and thus, for constant ``N²`` the above implies ``θ = θ₀ \exp{(N² z / g)}``. +# and thus, for constant ``N²``, the above implies that ``θ = θ₀ \exp{(N² z / g)}``. thermo = ThermodynamicConstants() g = thermo.gravitational_acceleration @@ -81,17 +78,29 @@ N = 0.01 # target dry Brunt–Väisälä frequency (s⁻¹) # First, we set up the shear layer using a ``\tanh`` profile: -z₀ = 1e3 # center of shear & moist layer (m) -Δzᶸ = 150 # shear layer half-thickness (m) -U_top = 25 # upper-layer wind (m/s) -U_bot = 5 # lower-layer wind (m/s) -uᵇ(z) = U_bot + (U_top - U_bot) * (1 + tanh((z - z₀) / Δzᶸ)) / 2 +z₀ = 1e3 # center of shear & moist layer (m) +Δzᵘ = 150 # shear layer half-thickness (m) +U₀ = 5 # base wind speed (m/s) +ΔU = 20 # upper-layer wind (m/s) +uᵇ(z) = U₀ + ΔU * (1 + tanh((z - z₀) / Δzᵘ)) / 2 # For the moisture layer, we use a Gaussian in ``z`` centered at ``z₀``: -q_max = 0.012 # peak specific humidity (kg/kg) -Δz_q = 200 # moist layer half-width (m) -qᵇ(z) = q_max * exp(-(z - z₀)^2 / 2Δz_q^2) +qᵗ₀ = 0.012 # peak specific humidity (kg/kg) +Δzᵗ = 200 # moist layer half-width (m) +qᵇ(z) = qᵗ₀ * exp(-(z - z₀)^2 / 2Δzᵗ^2) + +# We initialize the model via Oceananigans `set!`, adding also a bit of random noise. + +δθ = 0.01 +δu = 1e-3 +δq = 0.05 * qᵗ₀ + +θᵢ(x, z) = θᵇ(z) + δθ * rand() +qᵗᵢ(x, z) = qᵇ(z) + δq * rand() +uᵢ(x, z) = uᵇ(z) + δu * rand() + +set!(model; u=uᵢ, qᵗ=qᵗᵢ, θ=θᵢ) # ## The Kelvin-Helmholtz instability # @@ -107,12 +116,11 @@ qᵇ(z) = q_max * exp(-(z - z₀)^2 / 2Δz_q^2) # # Let's plot the initial state as well as the Richardson number. -z = znodes(grid, Center()) - -dudz = @. (U_top - U_bot) * sech((z - z₀) / Δzᶸ)^2 / 2Δzᶸ -Ri = N^2 ./ dudz.^2 +U = Field(Average(model.velocities.u, dims=(1, 2))) +Ri = N^2 / ∂z(U)^2 -using CairoMakie +Qᵗ = Field(Average(model.specific_moisture, dims=1)) +θ = Field(Average(liquid_ice_potential_temperature(model), dims=1)) fig = Figure(size=(800, 500)) @@ -121,11 +129,12 @@ axq = Axis(fig[1, 2], xlabel = "qᵇ (kg/kg)", title="Total liquid") axθ = Axis(fig[1, 3], xlabel = "θᵇ (K)", title="Potential temperature") axR = Axis(fig[1, 4], xlabel = "Ri", ylabel="z (m)", title="Richardson number") -lines!(axu, uᵇ.(z), z) -lines!(axq, qᵇ.(z), z) -lines!(axθ, θᵇ.(z), z) -lines!(axR, Ri, z) +lines!(axu, U) +lines!(axq, Qᵗ) +lines!(axθ, θ) +lines!(axR, Ri) lines!(axR, [1/4, 1/4], [0, Lz], linestyle = :dash, color = :black) + xlims!(axR, 0, 0.8) axR.xticks = 0:0.25:1 @@ -137,20 +146,6 @@ end fig -# ## Define initial conditions -# -# We initialize the model via Oceananigans `set!`, adding also a bit of random noise. - -δθ = 0.01 -δu = 1e-3 -δq = 0.05 * q_max - -θᵢ(x, z) = θᵇ(z) + δθ * rand() -qᵗᵢ(x, z) = qᵇ(z) + δq * rand() -uᵢ(x, z) = uᵇ(z) + δu * rand() - -set!(model; u=uᵢ, qᵗ=qᵗᵢ, θ=θᵢ) - # ## Set up and run the simulation # # We construct a simulation and use the time-step wizard to keep the CFL number under control. @@ -176,7 +171,7 @@ add_callback!(simulation, progress, TimeInterval(1minute)) # the potential temperatures and the specific humidities (vapour, liquid, ice). u, v, w = model.velocities ξ = ∂z(u) - ∂x(w) -θ = PotentialTemperatureField(model) +θ = liquid_ice_potential_temperature(model) outputs = merge(model.velocities, model.microphysical_fields, (; ξ, θ)) filename = "wave_clouds.jld2" diff --git a/examples/moist_thermal_bubble.jl b/examples/moist_thermal_bubble.jl new file mode 100644 index 00000000..5f493044 --- /dev/null +++ b/examples/moist_thermal_bubble.jl @@ -0,0 +1,282 @@ +# # Dry and cloudy thermal bubbles +# +# This example sets up, runs, and visualizes simulations of "thermal bubbles" +# (just circular regions of warm air) rising through a neutral background. +# We run both a dry simulation and a "cloudy" simulation. In the cloudy case, +# we simulate a pocket of warm air rising in a saturated, condensate-laden environment. + +using Breeze +using Oceananigans.Units +using Statistics +using Printf +using CairoMakie + +# ## Dry thermal bubble +# +# We first set up a dry thermal bubble simulation without moisture processes. +# This serves as a baseline for comparison with the moist case. + +grid = RectilinearGrid(CPU(); + size = (128, 128), halo = (5, 5), + x = (-10e3, 10e3), + z = (0, 10e3), + topology = (Bounded, Flat, Bounded)) + +thermodynamic_constants = ThermodynamicConstants() +reference_state = ReferenceState(grid, thermodynamic_constants, base_pressure=1e5, potential_temperature=300) +#formulation = AnelasticFormulation(reference_state, thermodynamics=:LiquidIcePotentialTemperature) +formulation = AnelasticFormulation(reference_state, thermodynamics=:StaticEnergy) +advection = Centered(order=2) #WENO(order=9) +model = AtmosphereModel(grid; formulation, thermodynamic_constants, advection) + +# ## Potential temperature perturbation +# +# We add a localized potential temperature perturbation for the dry bubble. +# In the dry case, this perturbation directly affects buoyancy without any +# moisture-related effects. + +r₀ = 2e3 +z₀ = 2e3 +Δθ = 2 # K +θ₀ = model.formulation.reference_state.potential_temperature +g = model.thermodynamic_constants.gravitational_acceleration + +function θᵢ(x, z) + r = sqrt((x / r₀)^2 + ((z - z₀) / r₀)^2) + return θ₀ + Δθ * cos(π * min(1, r) / 2)^2 +end + +set!(model, θ=θᵢ) + +# ## Initial dry bubble visualization +# +# Plot the initial potential temperature to visualize the dry thermal bubble. + +θ = liquid_ice_potential_temperature(model) +E = total_energy(model) +∫E = Integral(E) |> Field + +fig = Figure() +ax = Axis(fig[1, 1], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Initial potential temperature θ (K)") +hm = heatmap!(ax, θ) +Colorbar(fig[1, 2], hm, label = "ρe′ (J/kg)") +fig + +# ## Simulation rising + +simulation = Simulation(model; Δt=2, stop_time=1000) +conjure_time_step_wizard!(simulation, cfl=0.7) +θ = liquid_ice_potential_temperature(model) + +function progress(sim) + u, v, w = sim.model.velocities + + msg = @sprintf("Iter: % 4d, t: % 14s, Δt: % 14s, ∫E: %.8e J, extrema(θ): (%.2f, %.2f) K, max|w|: %.2f m/s", + iteration(sim), prettytime(sim), prettytime(sim.Δt), ∫E[], extrema(θ)..., maximum(abs, w)) + + @info msg + return nothing +end + +add_callback!(simulation, progress, TimeInterval(100)) + +u, v, w = model.velocities +outputs = (; θ, w) + +filename = "dry_thermal_bubble.jld2" +writer = JLD2Writer(model, outputs; filename, + schedule = TimeInterval(10seconds), + overwrite_existing = true) + +simulation.output_writers[:jld2] = writer + +run!(simulation) + +fig = Figure() +axθ = Axis(fig[1, 1], aspect=2, xlabel="x (m)", ylabel="z (m)") +axw = Axis(fig[2, 1], aspect=2, xlabel="x (m)", ylabel="z (m)") + +hmθ = heatmap!(axθ, θ) +hmw = heatmap!(axw, w) + +Colorbar(fig[1, 2], hmθ, label = "θ (K) at t = $(prettytime(simulation.model.clock.time))") +Colorbar(fig[2, 2], hmw, label = "w (m/s) at t = $(prettytime(simulation.model.clock.time))") + +fig + +# Just running to t=1000 is pretty boring, Let's run the simulation for a longer time, just for fun! + +# simulation.stop_time = 30minutes +# run!(simulation) + +# ## Visualization +# +# Visualize the potential temperature and the vertical velocity through +# time and create an animation. + +θt = FieldTimeSeries(filename, "θ") +wt = FieldTimeSeries(filename, "w") + +times = θt.times +fig = Figure(size = (800, 800), fontsize = 12) +axθ = Axis(fig[1, 1], aspect=2, xlabel="x (m)", ylabel="z (m)") +axw = Axis(fig[2, 1], aspect=2, xlabel="x (m)", ylabel="z (m)") + +n = Observable(length(θt)) +θn = @lift θt[$n] +wn = @lift wt[$n] + +title = @lift "Dry thermal bubble evolution — t = $(prettytime(times[$n]))" +fig[0, :] = Label(fig, title, fontsize = 16, tellwidth = false) + +θ_range = (minimum(θt), maximum(θt)) +w_range = maximum(abs, wt) + +hmθ = heatmap!(axθ, θn, colorrange = θ_range, colormap = :thermal) +hmw = heatmap!(axw, wn, colorrange = (-w_range, w_range), colormap = :balance) + +Colorbar(fig[1, 2], hmθ, label = "θ (K)", vertical = true) +Colorbar(fig[2, 2], hmw, label = "w (m/s)", vertical = true) + +CairoMakie.record(fig, "dry_thermal_bubble.mp4", 1:length(θt), framerate = 12) do nn + n[] = nn +end + +nothing #hide + +# ![](dry_thermal_bubble.mp4) + +# ## Moist thermal bubble with warm-phase saturation adjustment +# +# Now we set up a moist thermal bubble simulation with warm-phase saturation adjustment, +# following the methodology described by Bryan and Fritsch (2002). This simulation +# includes moisture processes, where excess water vapor condenses to liquid water, +# releasing latent heat that enhances the buoyancy of the rising bubble. +# +# For pedagogical purposes, we build a new model with warm-phase saturation adjustment microphysics. +# (We coudl have also used this model for the dry simulation): + +microphysics = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium()) +moist_model = AtmosphereModel(grid; formulation, thermodynamic_constants, advection, microphysics) + +# ## Moist thermal bubble initial conditions +# +# For the moist bubble, we initialize both temperature and moisture perturbations. +# The bubble is warm and moist, leading to condensation and latent heat release +# as it rises and cools. First, we set the potential temperature to match the dry case, +# then we use the diagnostic saturation specific humidity field to set the moisture. + +# Set potential temperature to match the dry bubble initially +set!(moist_model, θ=θᵢ, qᵗ=0.025) + +# Compute saturation specific humidity using the diagnostic field +using Breeze.Thermodynamics: dry_air_gas_constant, vapor_gas_constant + +qᵛ⁺ = SaturationSpecificHumidityField(moist_model, :equilibrium) +θᵈ = liquid_ice_potential_temperature(moist_model) # note, current state is dry +Rᵈ = dry_air_gas_constant(thermodynamic_constants) +Rᵛ = vapor_gas_constant(thermodynamic_constants) +Rᵐ = Rᵈ * (1 - qᵛ⁺) + Rᵛ * qᵛ⁺ +θᵐ = θᵈ * Rᵈ / Rᵐ + +set!(moist_model, θ=θᵐ) + +# ## Simulation + +moist_simulation = Simulation(moist_model; Δt=2, stop_time=30minutes) +conjure_time_step_wizard!(moist_simulation, cfl=0.7) + +E = total_energy(moist_model) +∫E = Integral(E) |> Field +θ = liquid_ice_potential_temperature(moist_model) + +function progress_moist(sim) + compute!(∫E) + ρqᵗ = sim.model.moisture_density + u, v, w = sim.model.velocities + + msg = @sprintf("Iter: % 4d, t: % 14s, Δt: % 14s, ∫E: %.8e J, extrema(θ): (%.2f, %.2f) K \n", + iteration(sim), prettytime(sim), prettytime(sim.Δt), ∫E[], extrema(θ)...) + + msg *= @sprintf(" extrema(qᵗ): (%.2e, %.2e), max(qˡ): %.2e, max|w|: %.2f m/s, mean(qᵗ): %.2e", + extrema(ρqᵗ)..., maximum(qˡ), maximum(abs, w), mean(ρqᵗ)) + + @info msg + return nothing +end + +add_callback!(moist_simulation, progress_moist, TimeInterval(100)) + +θ = liquid_ice_potential_temperature(moist_model) +u, v, w = moist_model.velocities +qᵗ = moist_model.specific_moisture +qˡ = moist_model.microphysical_fields.qˡ +qˡ′ = qˡ - Field(Average(qˡ, dims=1)) +moist_outputs = (; θ, w, qˡ′) + +moist_filename = "moist_thermal_bubble.jld2" +moist_writer = JLD2Writer(moist_model, moist_outputs; filename=moist_filename, + schedule = TimeInterval(10seconds), + overwrite_existing = true) + +moist_simulation.output_writers[:jld2] = moist_writer + +run!(moist_simulation) + +fig = Figure(size=(1200, 600)) + +axθ = Axis(fig[1, 2], aspect=2, xlabel="x (m)", ylabel="z (m)") +axw = Axis(fig[1, 3], aspect=2, xlabel="x (m)", ylabel="z (m)") +axl = Axis(fig[2, 2:3], aspect=2, xlabel="x (m)", ylabel="z (m)") + +qˡ = moist_model.microphysical_fields.qˡ + +hmθ = heatmap!(axθ, θ) +hmw = heatmap!(axw, w) +hmqˡ = heatmap!(axl, qˡ′) + +t_str = @sprintf("t = %s", prettytime(moist_simulation.model.clock.time)) +Colorbar(fig[1, 1], hmθ, label = "θ (K) at $t_str") +Colorbar(fig[1, 4], hmw, label = "w (m/s) at $t_str") +Colorbar(fig[2, 4], hmqˡ, label = "qˡ (kg/kg) at $t_str") + +fig + +# simulation_moist.stop_time = 30minutes +# run!(simulation_moist) + +# ## Visualization of moist thermal bubble + +θt = FieldTimeSeries(moist_filename, "θ") +wt = FieldTimeSeries(moist_filename, "w") +qˡ′t = FieldTimeSeries(moist_filename, "qˡ′") + +times = θt.times +fig = Figure(size = (1200, 800), fontsize = 12) +axθ = Axis(fig[1, 2], aspect=2, xlabel="x (m)", ylabel="z (m)") +axw = Axis(fig[1, 3], aspect=2, xlabel="x (m)", ylabel="z (m)") +axl = Axis(fig[2, 2:3], aspect=2, xlabel="x (m)", ylabel="z (m)") + +θ_range = (minimum(θt), maximum(θt)) +w_range = maximum(abs, wt) +qˡ′_range = (minimum(qˡ′t), maximum(qˡ′t)) + +n = Observable(length(θt)) +θn = @lift θt[$n] +wn = @lift wt[$n] +qˡ′n = @lift qˡ′t[$n] + +hmθ = heatmap!(axθ, θn, colorrange = θ_range, colormap = :thermal) +hmw = heatmap!(axw, wn, colorrange = (-w_range, w_range), colormap = :balance) +hml = heatmap!(axl, qˡ′n, colorrange = qˡ′_range, colormap = :balance) + +Colorbar(fig[1, 1], hmθ, label = "θ (K)", vertical = true) +Colorbar(fig[1, 4], hmw, label = "w (m/s)", vertical = true) +Colorbar(fig[2, 4], hml, label = "qˡ (kg/kg)", vertical = true) + +CairoMakie.record(fig, "moist_thermal_bubble.mp4", 1:length(θt), framerate = 12) do nn + n[] = nn +end +nothing #hide + +# ![](moist_thermal_bubble.mp4) \ No newline at end of file diff --git a/examples/prescribed_sst.jl b/examples/prescribed_sst.jl index 72c8bbb6..3b2fd983 100644 --- a/examples/prescribed_sst.jl +++ b/examples/prescribed_sst.jl @@ -202,7 +202,7 @@ conjure_time_step_wizard!(simulation, cfl=0.7) # liquid condensate mass fraction qˡ, and supersaturation δ = qᵗ - qᵛ⁺. T = model.temperature -θ = Breeze.AtmosphereModels.PotentialTemperatureField(model) +θ = liquid_ice_potential_temperature(model) qᵛ⁺ = Breeze.AtmosphereModels.SaturationSpecificHumidityField(model) # Liquid condensate mass fraction from microphysical fields diff --git a/examples/thermal_bubble.jl b/examples/thermal_bubble.jl index ea9bb9d8..9a06bcfa 100644 --- a/examples/thermal_bubble.jl +++ b/examples/thermal_bubble.jl @@ -43,7 +43,7 @@ end set!(model, θ = θᵢ) -ρe = energy_density(model) +ρe = static_energy_density(model) ρE = Field(Average(ρe, dims=1)) ρe′ = Field(ρe - ρE) @@ -64,7 +64,7 @@ simulation = Simulation(model; Δt=2, stop_time=25minutes) conjure_time_step_wizard!(simulation, cfl=0.7) function progress(sim) - ρe = energy_density(sim.model) + ρe = static_energy_density(sim.model) u, v, w = sim.model.velocities msg = @sprintf("Iter: %d, t: %s, Δt: %s, extrema(ρe): (%.2f, %.2f) J/kg, max|u|: %.2f m/s, max|w|: %.2f m/s", diff --git a/examples/time_step_atmos_model.jl b/examples/time_step_atmos_model.jl index 78a39b67..e792e749 100644 --- a/examples/time_step_atmos_model.jl +++ b/examples/time_step_atmos_model.jl @@ -80,7 +80,7 @@ set!(θ_bg_field, z -> Tₛ + dθdz * z) T = model.temperature ρʳ = model.formulation.reference_state.density cᵖᵈ = model.thermodynamic_constants.dry_air.heat_capacity -ρe = energy_density(model) +ρe = static_energy_density(model) θ = ρe / (ρʳ * cᵖᵈ) θ′ = θ - θ_bg_field diff --git a/src/AtmosphereModels/AtmosphereModels.jl b/src/AtmosphereModels/AtmosphereModels.jl index 0b747b3a..1d828659 100644 --- a/src/AtmosphereModels/AtmosphereModels.jl +++ b/src/AtmosphereModels/AtmosphereModels.jl @@ -5,10 +5,12 @@ export AtmosphereModelBuoyancy, AnelasticFormulation, StaticEnergyThermodynamics, - PotentialTemperature, - PotentialTemperatureField, - energy_density, - specific_energy + LiquidIcePotentialTemperatureThermodynamics, + static_energy_density, + static_energy, + total_energy, + potential_temperature_density, + liquid_ice_potential_temperature using DocStringExtensions: TYPEDSIGNATURES using Adapt: Adapt, adapt @@ -16,12 +18,14 @@ using Adapt: Adapt, adapt include("atmosphere_model.jl") include("diagnostic_fields.jl") +include("set_atmosphere_model.jl") include("anelastic_formulation.jl") +include("static_energy_thermodynamics.jl") +include("potential_temperature_thermodynamics.jl") include("atmosphere_model_buoyancy.jl") include("microphysics_interface.jl") include("dynamics_kernel_functions.jl") include("update_atmosphere_model_state.jl") -include("set_atmosphere_model.jl") include("compute_hydrostatic_pressure.jl") end diff --git a/src/AtmosphereModels/anelastic_formulation.jl b/src/AtmosphereModels/anelastic_formulation.jl index 6b052d80..7b9f9598 100644 --- a/src/AtmosphereModels/anelastic_formulation.jl +++ b/src/AtmosphereModels/anelastic_formulation.jl @@ -1,6 +1,7 @@ using ..Thermodynamics: MoistureMassFractions, StaticEnergyState, + LiquidIcePotentialTemperatureState, ThermodynamicConstants, ReferenceState, mixture_gas_constant, @@ -16,6 +17,7 @@ using Oceananigans.Solvers: solve!, AbstractHomogeneousNeumannFormulation using KernelAbstractions: @kernel, @index using Adapt: Adapt, adapt +import Oceananigans.BoundaryConditions: fill_halo_regions! import Oceananigans.Solvers: tridiagonal_direction, compute_main_diagonal!, compute_lower_diagonal! import Oceananigans.TimeSteppers: compute_pressure_correction!, make_pressure_correction! @@ -37,23 +39,17 @@ struct AnelasticFormulation{T, R, P} pressure_anomaly :: P end -struct StaticEnergyThermodynamics{E, S} - energy_density :: E - specific_energy :: S -end +const valid_thermodynamics_types = (:LiquidIcePotentialTemperature, :StaticEnergy) """ $(TYPEDSIGNATURES) -Construct a "stub" `AnelasticFormulation` with just the `reference_state`. +Construct an un-materialized "stub" `AnelasticFormulation` with `reference_state` and `thermodynamics`. The thermodynamics and pressure fields are materialized later in the model constructor. """ -AnelasticFormulation(reference_state) = - AnelasticFormulation(nothing, reference_state, nothing) - -Adapt.adapt_structure(to, thermo::StaticEnergyThermodynamics) = - StaticEnergyThermodynamics(adapt(to, thermo.energy_density), - adapt(to, thermo.specific_energy)) +function AnelasticFormulation(reference_state; thermodynamics=:StaticEnergy) + return AnelasticFormulation(thermodynamics, reference_state, nothing) +end Adapt.adapt_structure(to, formulation::AnelasticFormulation) = AnelasticFormulation(adapt(to, formulation.thermodynamics), @@ -62,6 +58,33 @@ Adapt.adapt_structure(to, formulation::AnelasticFormulation) = const AnelasticModel = AtmosphereModel{<:AnelasticFormulation} +# Type aliases for convenience + +function prognostic_field_names(formulation::AnelasticFormulation{<:Symbol}) + if formulation.thermodynamics == :StaticEnergy + return tuple(:ρe) + elseif formulation.thermodynamics == :LiquidIcePotentialTemperature + return tuple(:ρθ) + else + throw(ArgumentError("Got $(formulation.thermodynamics) thermodynamics, which is not one of \ + the valid types $valid_thermodynamics_types.")) + end +end + +function additional_field_names(formulation::AnelasticFormulation{<:Symbol}) + if formulation.thermodynamics == :StaticEnergy + return tuple(:e) + elseif formulation.thermodynamics == :LiquidIcePotentialTemperature + return tuple(:θ) + end +end + +""" + $(TYPEDSIGNATURES) + +Construct a "stub" `AnelasticFormulation` with just the `reference_state`. +The thermodynamics and pressure fields are materialized later in the model constructor. +""" function default_formulation(grid, constants) reference_state = ReferenceState(grid, constants) return AnelasticFormulation(reference_state) @@ -71,65 +94,29 @@ end $(TYPEDSIGNATURES) Materialize a stub `AnelasticFormulation` into a full formulation with thermodynamic fields -(`energy_density`, `specific_energy`) and the pressure anomaly field. +and the pressure anomaly field. The thermodynamic fields depend on the type of thermodynamics +specified in the stub (`:static_energy` or `:potential_temperature`). """ function materialize_formulation(stub::AnelasticFormulation, grid, boundary_conditions) - energy_density = CenterField(grid, boundary_conditions=boundary_conditions.ρe) - specific_energy = CenterField(grid) # e = ρe / ρᵣ (diagnostic per-mass energy) - thermodynamics = StaticEnergyThermodynamics(energy_density, specific_energy) + thermo_type = stub.thermodynamics pressure_anomaly = CenterField(grid) + thermodynamics = materialize_thermodynamics(Val(thermo_type), grid, boundary_conditions) return AnelasticFormulation(thermodynamics, stub.reference_state, pressure_anomaly) end -function Base.summary(formulation::AnelasticFormulation) - p₀ = formulation.reference_state.base_pressure - θ₀ = formulation.reference_state.potential_temperature - return string("AnelasticFormulation(p₀=", prettysummary(p₀), - ", θ₀=", prettysummary(θ₀), ")") +function materialize_thermodynamics(::Val{T}, grid, boundary_conditions) where T + throw(ArgumentError("Got $T thermodynamics, which is not one of \ + the valid types $valid_thermodynamics_types.")) + return nothing end -Base.show(io::IO, formulation::AnelasticFormulation) = print(io, "AnelasticFormulation") - -##### -##### Thermodynamic state -##### - -""" - $(TYPEDSIGNATURES) - -Return `StaticEnergyState` computed from the prognostic state including -energy density, moisture density, and microphysical fields. -""" -function diagnose_thermodynamic_state(i, j, k, grid, formulation::AnelasticFormulation, - microphysics, - microphysical_fields, - constants, - specific_energy, - specific_moisture) - @inbounds begin - e = specific_energy[i, j, k] - pᵣ = formulation.reference_state.pressure[i, j, k] - ρᵣ = formulation.reference_state.density[i, j, k] - qᵗ = specific_moisture[i, j, k] - end - - q = compute_moisture_fractions(i, j, k, grid, microphysics, ρᵣ, qᵗ, microphysical_fields) - z = znode(i, j, k, grid, c, c, c) - - return StaticEnergyState(e, q, z, pᵣ) +function Base.summary(formulation::AnelasticFormulation) + p₀_str = prettysummary(formulation.reference_state.base_pressure) + θ₀_str = prettysummary(formulation.reference_state.potential_temperature) + return string("AnelasticFormulation(p₀=", p₀_str, ", θ₀=", θ₀_str, ")") end - -function collect_prognostic_fields(::AnelasticFormulation, - momentum, - energy_density, - moisture_density, - microphysical_fields, - tracers) - - thermodynamic_variables = (ρe=energy_density, ρqᵗ=moisture_density) - return merge(momentum, thermodynamic_variables, microphysical_fields, tracers) -end +Base.show(io::IO, formulation::AnelasticFormulation) = print(io, "AnelasticFormulation") function materialize_momentum_and_velocities(formulation::AnelasticFormulation, grid, boundary_conditions) ρu = XFaceField(grid, boundary_conditions=boundary_conditions.ρu) diff --git a/src/AtmosphereModels/atmosphere_model.jl b/src/AtmosphereModels/atmosphere_model.jl index 993b3f7e..92c36e92 100644 --- a/src/AtmosphereModels/atmosphere_model.jl +++ b/src/AtmosphereModels/atmosphere_model.jl @@ -3,11 +3,11 @@ using ..Thermodynamics: Thermodynamics, ThermodynamicConstants, ReferenceState using Oceananigans: AbstractModel, Center, CenterField, Clock, Field using Oceananigans: Centered, XFaceField, YFaceField, ZFaceField using Oceananigans.Advection: adapt_advection_order +using Oceananigans.AbstractOperations: @at using Oceananigans.Forcings: materialize_forcing using Oceananigans.BoundaryConditions: FieldBoundaryConditions, regularize_field_boundary_conditions using Oceananigans.Grids: ZDirection using Oceananigans.Models: validate_model_halo, validate_tracer_advection -import Oceananigans.Models.HydrostaticFreeSurfaceModels: validate_momentum_advection using Oceananigans.Solvers: FourierTridiagonalPoissonSolver using Oceananigans.TimeSteppers: TimeStepper using Oceananigans.TurbulenceClosures: implicit_diffusion_solver, time_discretization, build_closure_fields @@ -15,6 +15,7 @@ using Oceananigans.Utils: launch!, prettytime, prettykeys, with_tracers import Oceananigans: fields, prognostic_fields import Oceananigans.Advection: cell_advection_timescale +import Oceananigans.Models.HydrostaticFreeSurfaceModels: validate_momentum_advection struct DefaultValue end @@ -132,7 +133,6 @@ function AtmosphereModel(grid; # Next, we form a list of default boundary conditions: prognostic_names = prognostic_field_names(formulation, microphysics, tracers) - FT = eltype(grid) default_boundary_conditions = NamedTuple{prognostic_names}(FieldBoundaryConditions() for _ in prognostic_names) boundary_conditions = merge(default_boundary_conditions, boundary_conditions) all_names = field_names(formulation, microphysics, tracers) @@ -151,9 +151,6 @@ function AtmosphereModel(grid; moisture_density = CenterField(grid, boundary_conditions=boundary_conditions.ρqᵗ) end - # Access energy fields from the formulation - energy_density = formulation.thermodynamics.energy_density - # Diagnostic fields specific_moisture = CenterField(grid) temperature = CenterField(grid) @@ -162,7 +159,6 @@ function AtmosphereModel(grid; prognostic_microphysical_fields = NamedTuple(microphysical_fields[name] for name in prognostic_field_names(microphysics)) prognostic_fields = collect_prognostic_fields(formulation, momentum, - energy_density, moisture_density, prognostic_microphysical_fields, tracers) @@ -174,8 +170,9 @@ function AtmosphereModel(grid; model_fields = merge(prognostic_fields, velocities, (; T=temperature, qᵗ=specific_moisture)) forcing = atmosphere_model_forcing(forcing, prognostic_fields, model_fields) - # Include ρe, ρqᵗ plus user tracers for closure field construction - scalar_names = tuple(:ρe, :ρqᵗ, tracer_names...) + # Include thermodynamic density (ρe or ρθ), ρqᵗ plus user tracers for closure field construction + closure_thermo_name = thermodynamic_density_name(formulation) + scalar_names = tuple(closure_thermo_name, :ρqᵗ, tracer_names...) closure = Oceananigans.Utils.with_tracers(scalar_names, closure) closure_fields = build_closure_fields(nothing, grid, clock, scalar_names, boundary_conditions, closure) @@ -250,16 +247,19 @@ end cell_advection_timescale(model::AtmosphereModel) = cell_advection_timescale(model.grid, model.velocities) +# Default prognostic field names - overloaded by formulation-specific files function prognostic_field_names(formulation, microphysics, tracer_names) - default_names = (:ρu, :ρv, :ρw, :ρe, :ρqᵗ) + default_names = (:ρu, :ρv, :ρw, :ρqᵗ) + formulation_names = prognostic_field_names(formulation) microphysical_names = prognostic_field_names(microphysics) - return tuple(default_names..., microphysical_names..., tracer_names...) + return tuple(default_names..., formulation_names..., microphysical_names..., tracer_names...) end function field_names(formulation, microphysics, tracer_names) - prognostic_names = prognostic_field_names(formulation, microphysics, tracer_names) - additional_names = (:u, :v, :w, :e, :T, :qᵗ) - return tuple(prognostic_names..., additional_names...) + prog_names = prognostic_field_names(formulation, microphysics, tracer_names) + default_additional_names = (:u, :v, :w, :T, :qᵗ) + formulation_additional_names = additional_field_names(formulation) + return tuple(prog_names..., default_additional_names..., formulation_additional_names...) end function atmosphere_model_forcing(user_forcings, prognostic_fields, model_fields) @@ -300,39 +300,60 @@ function atmosphere_model_forcing(user_forcings::NamedTuple, prognostic_fields, end function fields(model::AtmosphereModel) - specific_energy = model.formulation.thermodynamics.specific_energy - auxiliary_thermodynamic_fields = (e=specific_energy, T=model.temperature, qᵗ=model.specific_moisture) - return merge(prognostic_fields(model), - model.velocities, - auxiliary_thermodynamic_fields) + formulation_fields = fields(model.formulation) + auxiliary = (; T=model.temperature, qᵗ=model.specific_moisture) + return merge(prognostic_fields(model), formulation_fields, model.velocities, auxiliary) end function prognostic_fields(model::AtmosphereModel) - energy_density = model.formulation.thermodynamics.energy_density - thermodynamic_fields = (ρe=energy_density, ρqᵗ=model.moisture_density) - microphysical_names = prognostic_field_names(model.microphysics) - prognostic_microphysical_fields = NamedTuple{microphysical_names}( - model.microphysical_fields[name] for name in microphysical_names) - - return merge(model.momentum, thermodynamic_fields, prognostic_microphysical_fields, model.tracers) + prognostic_formulation_fields = prognostic_fields(model.formulation) + thermodynamic_fields = merge(prognostic_formulation_fields, (; ρqᵗ=model.moisture_density)) + μ_names = prognostic_field_names(model.microphysics) + μ_fields= NamedTuple{μ_names}(model.microphysical_fields[name] for name in μ_names) + return merge(model.momentum, thermodynamic_fields, μ_fields, model.tracers) end ##### ##### Helper functions for accessing thermodynamic fields ##### +# Stub function - implementation in anelastic_formulation.jl +function thermodynamic_density_name end + +""" + static_energy_density(model::AtmosphereModel) + +Return an `AbstractField` representing static energy density for `model`. """ - energy_density(model::AtmosphereModel) +function static_energy_density end -Return the energy density field `ρe` from the model formulation. -This is a convenience function for accessing `model.formulation.thermodynamics.energy_density`. """ -energy_density(model::AtmosphereModel) = model.formulation.thermodynamics.energy_density + static_energy(model::AtmosphereModel) +Return an `AbstractField` representing the (specific) static energy +for `model`. """ - specific_energy(model::AtmosphereModel) +function static_energy end -Return the specific energy field `e` from the model formulation. -This is a convenience function for accessing `model.formulation.thermodynamics.specific_energy`. """ -specific_energy(model::AtmosphereModel) = model.formulation.thermodynamics.specific_energy + potential_temperature_density(model::AtmosphereModel) + +Return an `AbstractField` representing potential temperature density +for `model`. +""" +function potential_temperature_density end + +""" + liquid_ice_potential_temperature(model::AtmosphereModel) + +Return an `AbstractField` representing potential temperature `θ` +for `model`. +""" +function liquid_ice_potential_temperature end + +function total_energy(model) + u, v, w = model.velocities + k = @at (Center, Center, Center) (u^2 + v^2 + w^2) / 2 |> Field + e = static_energy(model) |> Field + return k + e +end diff --git a/src/AtmosphereModels/atmosphere_model_buoyancy.jl b/src/AtmosphereModels/atmosphere_model_buoyancy.jl index 7386f1fc..88b2519f 100644 --- a/src/AtmosphereModels/atmosphere_model_buoyancy.jl +++ b/src/AtmosphereModels/atmosphere_model_buoyancy.jl @@ -34,8 +34,7 @@ function OceanTurbulenceClosures.buoyancy_tracers(model::AtmosphereModel) # Diagnostic fields for buoyancy gradient calculation buoyancy_tracers = (; T = model.temperature, qᵗ = model.specific_moisture) # Prognostic tracer fields for diffusivity computation - energy_density = model.formulation.thermodynamics.energy_density - prognostic_tracers = (; ρe = energy_density, ρqᵗ = model.moisture_density) + prognostic_tracers = merge(prognostic_fields(model.formulation), (; ρqᵗ = model.moisture_density)) # Merge with user tracers all_prognostic = merge(prognostic_tracers, model.tracers) # Final merge - buoyancy tracers at end for named access in ∂z_b diff --git a/src/AtmosphereModels/diagnostic_fields.jl b/src/AtmosphereModels/diagnostic_fields.jl index 0c3ff54f..357fa1e4 100644 --- a/src/AtmosphereModels/diagnostic_fields.jl +++ b/src/AtmosphereModels/diagnostic_fields.jl @@ -3,73 +3,99 @@ using Oceananigans: Center, Field using Oceananigans.AbstractOperations: KernelFunctionOperation using Adapt: Adapt, adapt -""" - SaturationSpecificHumidityKernel(temperature) +struct Specific end +struct Density end -A callable object for diagnosing the saturation specific humidity field, -given a temperature field, designed for use as a `KernelFunctionOperation` -in Oceananigans. Parameters are captured within the kernel struct for GPU friendliness. -""" -struct SaturationSpecificHumidityKernelFunction{R, T, μ, M, MF, TH} +struct LiquidIcePotentialTemperatureKernelFunction{F, R, μ, M, MF, TMP, TH} + flavor :: F reference_state :: R - temperature :: T microphysics :: μ microphysical_fields :: M specific_moisture :: MF + temperature :: TMP thermodynamic_constants :: TH end -Adapt.adapt_structure(to, k::SaturationSpecificHumidityKernelFunction) = - SaturationSpecificHumidityKernelFunction(adapt(to, k.reference_state), - adapt(to, k.temperature), - adapt(to, k.microphysics), - adapt(to, k.microphysical_fields), - adapt(to, k.specific_moisture), - adapt(to, k.thermodynamic_constants)) +Adapt.adapt_structure(to, k::LiquidIcePotentialTemperatureKernelFunction) = + LiquidIcePotentialTemperatureKernelFunction(adapt(to, k.flavor), + adapt(to, k.reference_state), + adapt(to, k.microphysics), + adapt(to, k.microphysical_fields), + adapt(to, k.specific_moisture), + adapt(to, k.temperature), + adapt(to, k.thermodynamic_constants)) -const C = Center - -const SaturationSpecificHumidityOperation = KernelFunctionOperation{C, C, C, <:Any, <:Any, <:SaturationSpecificHumidityKernelFunction} -const SaturationSpecificHumidityField = Field{C, C, C, <:SaturationSpecificHumidityOperation} +const LiquidIcePotentialTemperature = KernelFunctionOperation{Center, Center, Center, <:Any, <:Any, <:LiquidIcePotentialTemperatureKernelFunction} +const LiquidIcePotentialTemperatureField = Field{Center, Center, Center, <:LiquidIcePotentialTemperature} """ -$(TYPEDSIGNATURES) + $(TYPEDSIGNATURES) -Return a field for the saturation specific humidity. +Return a `KernelFunctionOperation` representing liquid-ice potential temperature. """ -function SaturationSpecificHumidityField(model) - func = SaturationSpecificHumidityKernelFunction(model.formulation.reference_state, - model.temperature, - model.microphysics, - model.microphysical_fields, - model.specific_moisture, - model.thermodynamic_constants) +function LiquidIcePotentialTemperature(model::AtmosphereModel, flavor_symbol=:specific) + + flavor = if flavor_symbol === :specific + Specific() + elseif flavor_symbol === :density + Density() + else + error("Unknown $flavor_symbol") + end - op = KernelFunctionOperation{Center, Center, Center}(func, model.grid) + func = LiquidIcePotentialTemperatureKernelFunction(flavor, + model.formulation.reference_state, + model.microphysics, + model.microphysical_fields, + model.specific_moisture, + model.temperature, + model.thermodynamic_constants) - return Field(op) + return KernelFunctionOperation{Center, Center, Center}(func, model.grid) end -function (d::SaturationSpecificHumidityKernelFunction)(i, j, k, grid) +""" + $(TYPEDSIGNATURES) + +Return a `Field` representing potential temperature. +""" +LiquidIcePotentialTemperatureField(model, flavor_symbol=:specific) = + LiquidIcePotentialTemperature(model, flavor_symbol) |> Field + +function (d::LiquidIcePotentialTemperatureKernelFunction)(i, j, k, grid) @inbounds begin pᵣ = d.reference_state.pressure[i, j, k] - T = d.temperature[i, j, k] - qᵗ = d.specific_moisture[i, j, k] ρᵣ = d.reference_state.density[i, j, k] + qᵗ = d.specific_moisture[i, j, k] + p₀ = d.reference_state.base_pressure + T = d.temperature[i, j, k] end + q = compute_moisture_fractions(i, j, k, grid, d.microphysics, ρᵣ, qᵗ, d.microphysical_fields) - ρ = Thermodynamics.density(pᵣ, T, q, d.thermodynamic_constants) - return Thermodynamics.saturation_specific_humidity(T, ρ, d.thermodynamic_constants, d.thermodynamic_constants.liquid) + Rᵐ = Thermodynamics.mixture_gas_constant(q, d.thermodynamic_constants) + cᵖᵐ = Thermodynamics.mixture_heat_capacity(q, d.thermodynamic_constants) + Π = (pᵣ / p₀)^(Rᵐ / cᵖᵐ) + + ℒˡᵣ = d.thermodynamic_constants.liquid.reference_latent_heat + ℒⁱᵣ = d.thermodynamic_constants.ice.reference_latent_heat + qˡ = q.liquid + qⁱ = q.ice + + θ = (T - (ℒˡᵣ * qˡ + ℒⁱᵣ * qⁱ) / cᵖᵐ) / Π + + if d.flavor isa Specific + return θ + elseif d.flavor isa Density + return ρᵣ * θ + end end -""" - PotentialTemperatureKernel(temperature) +##### +##### Static energy +##### -A callable object for diagnosing the potential temperature field, -given a temperature field, designed for use as a `KernelFunctionOperation` -in Oceananigans. Follows the pattern used for saturation diagnostics. -""" -struct PotentialTemperatureKernelFunction{R, μ, M, MF, TMP, TH} +struct StaticEnergyKernelFunction{F, R, μ, M, MF, TMP, TH} + flavor :: F reference_state :: R microphysics :: μ microphysical_fields :: M @@ -77,32 +103,43 @@ struct PotentialTemperatureKernelFunction{R, μ, M, MF, TMP, TH} temperature :: TMP thermodynamic_constants :: TH end - -Adapt.adapt_structure(to, k::PotentialTemperatureKernelFunction) = - PotentialTemperatureKernelFunction(adapt(to, k.reference_state), - adapt(to, k.microphysics), - adapt(to, k.microphysical_fields), - adapt(to, k.specific_moisture), - adapt(to, k.temperature), - adapt(to, k.thermodynamic_constants)) - -const PotentialTemperature = KernelFunctionOperation{Center, Center, Center, <:Any, <:Any, <:PotentialTemperatureKernelFunction} -const PotentialTemperatureField = Field{Center, Center, Center, <:PotentialTemperature} + +Adapt.adapt_structure(to, k::StaticEnergyKernelFunction) = + StaticEnergyKernelFunction(adapt(to, k.flavor), + adapt(to, k.reference_state), + adapt(to, k.microphysics), + adapt(to, k.microphysical_fields), + adapt(to, k.specific_moisture), + adapt(to, k.temperature), + adapt(to, k.thermodynamic_constants)) + +const StaticEnergy = KernelFunctionOperation{Center, Center, Center, <:Any, <:Any, <:StaticEnergyKernelFunction} +const StaticEnergyField = Field{Center, Center, Center, <:StaticEnergy} """ $(TYPEDSIGNATURES) Return a `KernelFunctionOperation` representing potential temperature. """ -function PotentialTemperature(model) - grid = model.grid - func = PotentialTemperatureKernelFunction(model.formulation.reference_state, - model.microphysics, - model.microphysical_fields, - model.specific_moisture, - model.temperature, - model.thermodynamic_constants) - return KernelFunctionOperation{Center, Center, Center}(func, grid) +function StaticEnergy(model, flavor_symbol=:specific) + + flavor = if flavor_symbol === :specific + Specific() + elseif flavor_symbol === :density + Density() + else + error("Unknown $flavor_symbol") + end + + func = StaticEnergyKernelFunction(flavor, + model.formulation.reference_state, + model.microphysics, + model.microphysical_fields, + model.specific_moisture, + model.temperature, + model.thermodynamic_constants) + + return KernelFunctionOperation{Center, Center, Center}(func, model.grid) end """ @@ -110,11 +147,11 @@ end Return a `Field` representing potential temperature. """ -PotentialTemperatureField(model) = Field(PotentialTemperature(model)) +StaticEnergyField(model, flavor_symbol=:specific) = + StaticEnergy(model, flavor_symbol) |> Field -function (d::PotentialTemperatureKernelFunction)(i, j, k, grid) +function (d::StaticEnergyKernelFunction)(i, j, k, grid) @inbounds begin - pᵣ = d.reference_state.pressure[i, j, k] ρᵣ = d.reference_state.density[i, j, k] qᵗ = d.specific_moisture[i, j, k] p₀ = d.reference_state.base_pressure @@ -122,14 +159,22 @@ function (d::PotentialTemperatureKernelFunction)(i, j, k, grid) end q = compute_moisture_fractions(i, j, k, grid, d.microphysics, ρᵣ, qᵗ, d.microphysical_fields) - Rᵐ = Thermodynamics.mixture_gas_constant(q, d.thermodynamic_constants) cᵖᵐ = Thermodynamics.mixture_heat_capacity(q, d.thermodynamic_constants) - Π = (pᵣ / p₀)^(Rᵐ / cᵖᵐ) + + g = d.thermodynamic_constants.gravitational_acceleration + z = znode(i, j, k, grid, c, c, c) ℒˡᵣ = d.thermodynamic_constants.liquid.reference_latent_heat ℒⁱᵣ = d.thermodynamic_constants.ice.reference_latent_heat qˡ = q.liquid qⁱ = q.ice - return (T - (ℒˡᵣ * qˡ + ℒⁱᵣ * qⁱ) / cᵖᵐ) / Π + # Moist static energy + e = cᵖᵐ * T + g * z - ℒˡᵣ * qˡ + ℒⁱᵣ * qⁱ + + if d.flavor isa Specific + return e + elseif d.flavor isa Density + return ρᵣ * e + end end diff --git a/src/AtmosphereModels/dynamics_kernel_functions.jl b/src/AtmosphereModels/dynamics_kernel_functions.jl index 5000fec4..7b020d76 100644 --- a/src/AtmosphereModels/dynamics_kernel_functions.jl +++ b/src/AtmosphereModels/dynamics_kernel_functions.jl @@ -126,7 +126,6 @@ end advection, formulation, constants, - specific_energy, specific_moisture, velocities, microphysics, @@ -147,7 +146,6 @@ end microphysics, microphysical_fields, constants, - specific_energy, specific_moisture) return ( - div_ρUc(i, j, k, grid, advection, ρ, Uᵗ, c) @@ -155,44 +153,3 @@ end + microphysical_tendency(i, j, k, grid, microphysics, name, microphysical_fields, 𝒰, constants) + c_forcing(i, j, k, grid, clock, model_fields)) end - -@inline function static_energy_tendency(i, j, k, grid, - id, - ρe_forcing, - advection, - formulation, - constants, - specific_energy, - specific_moisture, - velocities, - microphysics, - microphysical_fields, - closure, - closure_fields, - clock, - model_fields, - temperature) - - 𝒰 = diagnose_thermodynamic_state(i, j, k, grid, - formulation, - microphysics, - microphysical_fields, - constants, - specific_energy, - specific_moisture) - - ρ = formulation.reference_state.density - - # Compute the buoyancy flux term, ρᵣ w b - buoyancy_flux = ℑzᵃᵃᶜ(i, j, k, grid, ρ_w_bᶜᶜᶠ, - velocities.w, formulation, ρ, temperature, specific_moisture, - microphysics, microphysical_fields, constants) - - closure_buoyancy = AtmosphereModelBuoyancy(formulation, constants) - - return ( - div_ρUc(i, j, k, grid, advection, ρ, velocities, specific_energy) - + buoyancy_flux - - ∇_dot_Jᶜ(i, j, k, grid, ρ, closure, closure_fields, id, specific_energy, clock, model_fields, closure_buoyancy) - + microphysical_tendency(i, j, k, grid, microphysics, Val(:ρe), microphysical_fields, 𝒰, constants) - + ρe_forcing(i, j, k, grid, clock, model_fields)) -end diff --git a/src/AtmosphereModels/potential_temperature_thermodynamics.jl b/src/AtmosphereModels/potential_temperature_thermodynamics.jl new file mode 100644 index 00000000..08a26269 --- /dev/null +++ b/src/AtmosphereModels/potential_temperature_thermodynamics.jl @@ -0,0 +1,202 @@ +using Breeze.Thermodynamics: LiquidIcePotentialTemperatureState, with_temperature + +struct LiquidIcePotentialTemperatureThermodynamics{F, T} + potential_temperature_density :: F # ρθ (prognostic) + potential_temperature :: T # θ = ρθ / ρᵣ (diagnostic) +end + +Adapt.adapt_structure(to, thermo::LiquidIcePotentialTemperatureThermodynamics) = + LiquidIcePotentialTemperatureThermodynamics(adapt(to, thermo.potential_temperature_density), + adapt(to, thermo.potential_temperature)) + +function fill_halo_regions!(thermo::LiquidIcePotentialTemperatureThermodynamics) + fill_halo_regions!(thermo.potential_temperature_density) + fill_halo_regions!(thermo.potential_temperature) + return nothing +end + +const APTF = AnelasticFormulation{<:LiquidIcePotentialTemperatureThermodynamics} + +prognostic_field_names(formulation::APTF) = tuple(:ρθ) +additional_field_names(formulation::APTF) = tuple(:θ) +thermodynamic_density_name(::APTF) = :ρθ +thermodynamic_density(formulation::APTF) = formulation.thermodynamics.potential_temperature_density +fields(formulation::APTF) = (; θ=formulation.thermodynamics.potential_temperature) +prognostic_fields(formulation::APTF) = (; ρθ=formulation.thermodynamics.potential_temperature_density) + +function materialize_thermodynamics(::Val{:LiquidIcePotentialTemperature}, grid, boundary_conditions) + potential_temperature_density = CenterField(grid, boundary_conditions=boundary_conditions.ρθ) + potential_temperature = CenterField(grid) # θ = ρθ / ρᵣ (diagnostic) + return LiquidIcePotentialTemperatureThermodynamics(potential_temperature_density, potential_temperature) +end + +function compute_auxiliary_thermodynamic_variables!(formulation::APTF, i, j, k, grid) + @inbounds begin + ρᵣ = formulation.reference_state.density[i, j, k] + ρθ = formulation.thermodynamics.potential_temperature_density[i, j, k] + formulation.thermodynamics.potential_temperature[i, j, k] = ρθ / ρᵣ + end + return nothing +end + +function diagnose_thermodynamic_state(i, j, k, grid, formulation::APTF, + microphysics, + microphysical_fields, + constants, + specific_moisture) + + θ = @inbounds formulation.thermodynamics.potential_temperature[i, j, k] + pᵣ = @inbounds formulation.reference_state.pressure[i, j, k] + ρᵣ = @inbounds formulation.reference_state.density[i, j, k] + p₀ = formulation.reference_state.base_pressure + qᵗ = @inbounds specific_moisture[i, j, k] + + q = compute_moisture_fractions(i, j, k, grid, microphysics, ρᵣ, qᵗ, microphysical_fields) + + return LiquidIcePotentialTemperatureState(θ, q, p₀, pᵣ) +end + +function collect_prognostic_fields(formulation::APTF, + momentum, + moisture_density, + microphysical_fields, + tracers) + + ρθ = formulation.thermodynamics.potential_temperature_density + thermodynamic_variables = (ρθ=ρθ, ρqᵗ=moisture_density) + return merge(momentum, thermodynamic_variables, microphysical_fields, tracers) +end + +const LiquidIcePotentialTemperatureAnelasticModel = AtmosphereModel{<:APTF} +const LIPTAM = LiquidIcePotentialTemperatureAnelasticModel + +liquid_ice_potential_temperature_density(model::LIPTAM) = model.formulation.thermodynamics.potential_temperature_density +liquid_ice_potential_temperature(model::LIPTAM) = model.formulation.thermodynamics.potential_temperature +static_energy(model::LIPTAM) = StaticEnergyField(model, :specific) +static_energy_density(model::LIPTAM) = StaticEnergyField(model, :density) + +function compute_thermodynamic_tendency!(model::LiquidIcePotentialTemperatureAnelasticModel, common_args) + grid = model.grid + arch = grid.architecture + + ρθ_args = ( + Val(1), + model.forcing.ρθ, + model.advection.ρθ, + common_args..., + model.temperature) + + Gρθ = model.timestepper.Gⁿ.ρθ + launch!(arch, grid, :xyz, compute_potential_temperature_tendency!, Gρθ, grid, ρθ_args) + return nothing +end + +@inline function potential_temperature_tendency(i, j, k, grid, + id, + ρθ_forcing, + advection, + formulation, + constants, + specific_moisture, + velocities, + microphysics, + microphysical_fields, + closure, + closure_fields, + clock, + model_fields, + temperature) + + potential_temperature = formulation.thermodynamics.potential_temperature + ρᵣ = formulation.reference_state.density + + 𝒰 = diagnose_thermodynamic_state(i, j, k, grid, + formulation, + microphysics, + microphysical_fields, + constants, + specific_moisture) + + closure_buoyancy = AtmosphereModelBuoyancy(formulation, constants) + + return ( - div_ρUc(i, j, k, grid, advection, ρᵣ, velocities, potential_temperature) + - ∇_dot_Jᶜ(i, j, k, grid, ρᵣ, closure, closure_fields, id, potential_temperature, clock, model_fields, closure_buoyancy) + + microphysical_tendency(i, j, k, grid, microphysics, Val(:ρθ), microphysical_fields, 𝒰, constants) + + ρθ_forcing(i, j, k, grid, clock, model_fields)) +end + +##### +##### Set +##### + +set_thermodynamic_variable!(model::LiquidIcePotentialTemperatureAnelasticModel, ::Val{:ρθ}, value) = + set!(model.formulation.thermodynamics.potential_temperature_density, value) + +function set_thermodynamic_variable!(model::LiquidIcePotentialTemperatureAnelasticModel, ::Val{:θ}, value) + set!(model.formulation.thermodynamics.potential_temperature, value) + ρᵣ = model.formulation.reference_state.density + θ = model.formulation.thermodynamics.potential_temperature + set!(model.formulation.thermodynamics.potential_temperature_density, ρᵣ * θ) + return nothing +end + +# Setting :θ (potential temperature) +function set_thermodynamic_variable!(model::LiquidIcePotentialTemperatureAnelasticModel, ::Val{:e}, value) + thermo = model.formulation.thermodynamics + e = model.temperature # scratch space + set!(e, value) + + grid = model.grid + arch = grid.architecture + launch!(arch, grid, :xyz, + _potential_temperature_from_energy!, + thermo.potential_temperature_density, + thermo.potential_temperature, + grid, + e, + model.specific_moisture, + model.formulation, + model.microphysics, + model.microphysical_fields, + model.thermodynamic_constants) + + return nothing +end + +function set_thermodynamic_variable!(model::LiquidIcePotentialTemperatureAnelasticModel, ::Val{:ρe}, value) + ρe = model.temperature # scratch space + set!(ρe, value) + ρᵣ = model.formulation.reference_state.density + return set_thermodynamic_variable!(model, Val(:e), ρe / ρᵣ) +end + +@kernel function _potential_temperature_from_energy!(potential_temperature_density, + potential_temperature, + grid, + specific_energy, + specific_moisture, + formulation, + microphysics, + microphysical_fields, + constants) + i, j, k = @index(Global, NTuple) + + @inbounds begin + pᵣ = formulation.reference_state.pressure[i, j, k] + ρᵣ = formulation.reference_state.density[i, j, k] + qᵗ = specific_moisture[i, j, k] + e = specific_energy[i, j, k] + end + + z = znode(i, j, k, grid, c, c, c) + q = compute_moisture_fractions(i, j, k, grid, microphysics, ρᵣ, qᵗ, microphysical_fields) + 𝒰e₀ = StaticEnergyState(e, q, z, pᵣ) + 𝒰e₁ = maybe_adjust_thermodynamic_state(𝒰e₀, microphysics, microphysical_fields, qᵗ, constants) + T = temperature(𝒰e₁, constants) + + p₀ = formulation.reference_state.base_pressure + q₁ = 𝒰e₁.moisture_mass_fractions + 𝒰θ = LiquidIcePotentialTemperatureState(zero(T), q₁, p₀, pᵣ) + @inbounds potential_temperature[i, j, k] = with_temperature(𝒰θ, T, constants).potential_temperature + @inbounds potential_temperature_density[i, j, k] = ρᵣ * with_temperature(𝒰θ, T, constants).potential_temperature +end diff --git a/src/AtmosphereModels/set_atmosphere_model.jl b/src/AtmosphereModels/set_atmosphere_model.jl index c5af3750..f34f9c10 100644 --- a/src/AtmosphereModels/set_atmosphere_model.jl +++ b/src/AtmosphereModels/set_atmosphere_model.jl @@ -4,9 +4,10 @@ using Oceananigans.BoundaryConditions: fill_halo_regions! using Oceananigans.TimeSteppers: compute_pressure_correction!, make_pressure_correction!, update_state! using ..Thermodynamics: - PotentialTemperatureState, + LiquidIcePotentialTemperatureState, MoistureMassFractions, mixture_heat_capacity, + mixture_gas_constant, temperature import Oceananigans.Fields: set! @@ -25,6 +26,8 @@ function prioritize_names(names) return names end +function set_thermodynamic_variable! end + function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) names = collect(keys(kw)) prioritized = prioritize_names(names) @@ -42,8 +45,10 @@ function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) set!(c, value) elseif name == :ρe - energy_density = model.formulation.thermodynamics.energy_density - set!(energy_density, value) + set_thermodynamic_variable!(model, Val(:ρe), value) + + elseif name == :ρθ + set_thermodynamic_variable!(model, Val(:ρθ), value) elseif name == :ρqᵗ set!(model.moisture_density, value) @@ -72,33 +77,10 @@ function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) set!(ϕ, value) elseif name == :e - # Set specific energy directly - specific_energy = model.formulation.thermodynamics.specific_energy - energy_density = model.formulation.thermodynamics.energy_density - set!(specific_energy, value) - ρᵣ = model.formulation.reference_state.density - set!(energy_density, ρᵣ * specific_energy) + set_thermodynamic_variable!(model, Val(:e), value) elseif name == :θ - θ = model.temperature # use scratch - set!(θ, value) - - grid = model.grid - arch = grid.architecture - energy_density = model.formulation.thermodynamics.energy_density - specific_energy = model.formulation.thermodynamics.specific_energy - - launch!(arch, grid, :xyz, - _energy_density_from_potential_temperature!, - energy_density, - specific_energy, - grid, - θ, - model.specific_moisture, - model.formulation, - model.microphysics, - model.microphysical_fields, - model.thermodynamic_constants) + set_thermodynamic_variable!(model, Val(:θ), value) else prognostic_names = keys(prognostic_fields(model)) @@ -125,48 +107,5 @@ function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) update_state!(model, compute_tendencies=false) end - energy_density = model.formulation.thermodynamics.energy_density - fill_halo_regions!(energy_density) - return nothing end - -@kernel function _energy_density_from_potential_temperature!(energy_density, - specific_energy, - grid, - potential_temperature, - specific_moisture, - formulation::AnelasticFormulation, - microphysics, - microphysical_fields, - constants) - i, j, k = @index(Global, NTuple) - - @inbounds begin - pᵣ = formulation.reference_state.pressure[i, j, k] - ρᵣ = formulation.reference_state.density[i, j, k] - qᵗ = specific_moisture[i, j, k] - θ = potential_temperature[i, j, k] - end - - g = constants.gravitational_acceleration - z = znode(i, j, k, grid, c, c, c) - p₀ = formulation.reference_state.base_pressure - - q = compute_moisture_fractions(i, j, k, grid, microphysics, ρᵣ, qᵗ, microphysical_fields) - 𝒰₀ = PotentialTemperatureState(θ, q, p₀, pᵣ) - 𝒰 = maybe_adjust_thermodynamic_state(𝒰₀, microphysics, microphysical_fields, qᵗ, constants) - - T = temperature(𝒰, constants) - q = 𝒰.moisture_mass_fractions - cᵖᵐ = mixture_heat_capacity(q, constants) - - ℒˡᵣ = constants.liquid.reference_latent_heat - ℒⁱᵣ = constants.ice.reference_latent_heat - qˡ = q.liquid - qⁱ = q.ice - - e = cᵖᵐ * T + g * z - ℒˡᵣ * qˡ - ℒⁱᵣ * qⁱ - @inbounds specific_energy[i, j, k] = e - @inbounds energy_density[i, j, k] = ρᵣ * e -end diff --git a/src/AtmosphereModels/static_energy_thermodynamics.jl b/src/AtmosphereModels/static_energy_thermodynamics.jl new file mode 100644 index 00000000..e3d13fc5 --- /dev/null +++ b/src/AtmosphereModels/static_energy_thermodynamics.jl @@ -0,0 +1,205 @@ +using Breeze.Thermodynamics: StaticEnergyState, with_temperature + +struct StaticEnergyThermodynamics{E, S} + energy_density :: E + specific_energy :: S +end + +Adapt.adapt_structure(to, thermo::StaticEnergyThermodynamics) = + StaticEnergyThermodynamics(adapt(to, thermo.energy_density), + adapt(to, thermo.specific_energy)) + +function fill_halo_regions!(thermo::StaticEnergyThermodynamics) + fill_halo_regions!(thermo.energy_density) + fill_halo_regions!(thermo.specific_energy) + return nothing +end + +const ASEF = AnelasticFormulation{<:StaticEnergyThermodynamics} + +prognostic_field_names(formulation::ASEF) = tuple(:ρe) +additional_field_names(formulation::ASEF) = tuple(:e) +thermodynamic_density_name(::ASEF) = :ρe +thermodynamic_density(formulation::ASEF) = formulation.thermodynamics.energy_density +fields(formulation::ASEF) = (; e=formulation.thermodynamics.specific_energy) +prognostic_fields(formulation::ASEF) = (; ρe=formulation.thermodynamics.energy_density) + +function materialize_thermodynamics(::Val{:StaticEnergy}, grid, boundary_conditions) + energy_density = CenterField(grid, boundary_conditions=boundary_conditions.ρe) + specific_energy = CenterField(grid) # e = ρe / ρᵣ (diagnostic per-mass energy) + return StaticEnergyThermodynamics(energy_density, specific_energy) +end + +function compute_auxiliary_thermodynamic_variables!(formulation::ASEF, i, j, k, grid) + @inbounds begin + ρᵣ = formulation.reference_state.density[i, j, k] + ρe = formulation.thermodynamics.energy_density[i, j, k] + formulation.thermodynamics.specific_energy[i, j, k] = ρe / ρᵣ + end + return nothing +end + +function diagnose_thermodynamic_state(i, j, k, grid, formulation::ASEF, + microphysics, + microphysical_fields, + constants, + specific_moisture) + + e = @inbounds formulation.thermodynamics.specific_energy[i, j, k] + pᵣ = @inbounds formulation.reference_state.pressure[i, j, k] + ρᵣ = @inbounds formulation.reference_state.density[i, j, k] + qᵗ = @inbounds specific_moisture[i, j, k] + + q = compute_moisture_fractions(i, j, k, grid, microphysics, ρᵣ, qᵗ, microphysical_fields) + z = znode(i, j, k, grid, c, c, c) + + return StaticEnergyState(e, q, z, pᵣ) +end + +function collect_prognostic_fields(formulation::ASEF, + momentum, + moisture_density, + microphysical_fields, + tracers) + ρe = formulation.thermodynamics.energy_density + thermodynamic_variables = (ρe=ρe, ρqᵗ=moisture_density) + return merge(momentum, thermodynamic_variables, microphysical_fields, tracers) +end + +const StaticEnergyAnelasticModel = AtmosphereModel{<:ASEF} +const SEAM = StaticEnergyAnelasticModel + +liquid_ice_potential_temperature(model::SEAM) = LiquidIcePotentialTemperature(model, :specific) +potential_temperature_density(model::SEAM) = LiquidIcePotentialTemperature(model, :density) +static_energy(model::SEAM) = model.formulation.thermodynamics.specific_energy +static_energy_density(model::SEAM) = model.formulation.thermodynamics.energy_density + +function compute_thermodynamic_tendency!(model::StaticEnergyAnelasticModel, common_args) + grid = model.grid + arch = grid.architecture + + ρe_args = ( + Val(1), + model.forcing.ρe, + model.advection.ρe, + common_args..., + model.temperature) + + Gρe = model.timestepper.Gⁿ.ρe + launch!(arch, grid, :xyz, compute_static_energy_tendency!, Gρe, grid, ρe_args) + return nothing +end + +@inline function static_energy_tendency(i, j, k, grid, + id, + ρe_forcing, + advection, + formulation, + constants, + specific_moisture, + velocities, + microphysics, + microphysical_fields, + closure, + closure_fields, + clock, + model_fields, + temperature) + + specific_energy = formulation.thermodynamics.specific_energy + + 𝒰 = diagnose_thermodynamic_state(i, j, k, grid, + formulation, + microphysics, + microphysical_fields, + constants, + specific_moisture) + + ρ = formulation.reference_state.density + + # Compute the buoyancy flux term, ρᵣ w b + buoyancy_flux = ℑzᵃᵃᶜ(i, j, k, grid, ρ_w_bᶜᶜᶠ, + velocities.w, formulation, ρ, temperature, specific_moisture, + microphysics, microphysical_fields, constants) + + closure_buoyancy = AtmosphereModelBuoyancy(formulation, constants) + + return ( - div_ρUc(i, j, k, grid, advection, ρ, velocities, specific_energy) + + buoyancy_flux + - ∇_dot_Jᶜ(i, j, k, grid, ρ, closure, closure_fields, id, specific_energy, clock, model_fields, closure_buoyancy) + + microphysical_tendency(i, j, k, grid, microphysics, Val(:ρe), microphysical_fields, 𝒰, constants) + + ρe_forcing(i, j, k, grid, clock, model_fields)) +end + +##### +##### Dispatch for setting thermodynamic variables +##### + +# StaticEnergyThermodynamics: :ρe sets energy density directly +set_thermodynamic_variable!(model::StaticEnergyAnelasticModel, ::Val{:ρe}, value) = + set!(model.formulation.thermodynamics.energy_density, value) + +function set_thermodynamic_variable!(model::StaticEnergyAnelasticModel, ::Val{:e}, value) + set!(model.formulation.thermodynamics.specific_energy, value) + ρᵣ = model.formulation.reference_state.density + e = model.formulation.thermodynamics.specific_energy + set!(model.formulation.thermodynamics.energy_density, ρᵣ * e) + return nothing +end + +# Setting :θ (potential temperature) +function set_thermodynamic_variable!(model::StaticEnergyAnelasticModel, ::Val{:θ}, value) + thermo = model.formulation.thermodynamics + θ = model.temperature # scratch space + set!(θ, value) + + grid = model.grid + arch = grid.architecture + launch!(arch, grid, :xyz, + _energy_density_from_potential_temperature!, + thermo.energy_density, + thermo.specific_energy, + grid, + θ, + model.specific_moisture, + model.formulation, + model.microphysics, + model.microphysical_fields, + model.thermodynamic_constants) + + return nothing +end + +@kernel function _energy_density_from_potential_temperature!(energy_density, + specific_energy, + grid, + potential_temperature, + specific_moisture, + formulation::AnelasticFormulation, + microphysics, + microphysical_fields, + constants) + i, j, k = @index(Global, NTuple) + + @inbounds begin + pᵣ = formulation.reference_state.pressure[i, j, k] + ρᵣ = formulation.reference_state.density[i, j, k] + qᵗ = specific_moisture[i, j, k] + θ = potential_temperature[i, j, k] + end + + p₀ = formulation.reference_state.base_pressure + q = compute_moisture_fractions(i, j, k, grid, microphysics, ρᵣ, qᵗ, microphysical_fields) + 𝒰θ₀ = LiquidIcePotentialTemperatureState(θ, q, p₀, pᵣ) + 𝒰θ₁ = maybe_adjust_thermodynamic_state(𝒰θ₀, microphysics, microphysical_fields, qᵗ, constants) + T = temperature(𝒰θ₁, constants) + + z = znode(i, j, k, grid, c, c, c) + q₁ = 𝒰θ₁.moisture_mass_fractions + 𝒰e₀ = StaticEnergyState(zero(T), q₁, z, pᵣ) + 𝒰e₁ = with_temperature(𝒰e₀, T, constants) + e = 𝒰e₁.static_energy + + @inbounds specific_energy[i, j, k] = e + @inbounds energy_density[i, j, k] = ρᵣ * e +end diff --git a/src/AtmosphereModels/update_atmosphere_model_state.jl b/src/AtmosphereModels/update_atmosphere_model_state.jl index b1e3e698..35a91108 100644 --- a/src/AtmosphereModels/update_atmosphere_model_state.jl +++ b/src/AtmosphereModels/update_atmosphere_model_state.jl @@ -50,7 +50,7 @@ Compute auxiliary model variables: - thermodynamic variables from the prognostic thermodynamic state, * temperature ``T``, possibly involving saturation adjustment - * moist static energy ``e = ρe / ρ`` + * specific thermodynamic variable (``e = ρe / ρ`` or ``θ = ρθ / ρ``) * moisture mass fraction ``qᵗ = ρqᵗ / ρ`` @@ -69,33 +69,36 @@ function compute_auxiliary_variables!(model) fill_halo_regions!(model.velocities) foreach(mask_immersed_field!, model.velocities) - energy_density = model.formulation.thermodynamics.energy_density - specific_energy = model.formulation.thermodynamics.specific_energy + # Dispatch on thermodynamics type + compute_auxiliary_thermodynamic_variables!(model) + + # Compute diffusivities + compute_diffusivities!(model.closure_fields, model.closure, model) + + # TODO: should we mask the auxiliary variables? They can also be masked in the kernel + + return nothing +end + +function compute_auxiliary_thermodynamic_variables!(model::AtmosphereModel) + grid = model.grid + arch = grid.architecture launch!(arch, grid, :xyz, _compute_auxiliary_thermodynamic_variables!, model.temperature, - specific_energy, model.specific_moisture, + model.formulation, grid, model.thermodynamic_constants, - model.formulation, model.microphysics, model.microphysical_fields, - energy_density, model.moisture_density) - # TODO: Can we compute the thermodynamic variable within halos as well, and avoid - # halo filling later on? fill_halo_regions!(model.temperature) - fill_halo_regions!(specific_energy) fill_halo_regions!(model.specific_moisture) fill_halo_regions!(model.microphysical_fields) - - # Compute diffusivities - compute_diffusivities!(model.closure_fields, model.closure, model) - - # TODO: should we mask the auxiliary variables? They can also be masked in the kernel + fill_halo_regions!(model.formulation.thermodynamics) return nothing end @@ -118,25 +121,63 @@ end end @kernel function _compute_auxiliary_thermodynamic_variables!(temperature, - specific_energy, specific_moisture, + formulation, grid, constants, - formulation, microphysics, microphysical_fields, - energy_density, moisture_density) i, j, k = @index(Global, NTuple) + compute_auxiliary_thermodynamic_variables!(formulation, i, j, k, grid) + + @inbounds begin + ρ = formulation.reference_state.density[i, j, k] + ρqᵗ = moisture_density[i, j, k] + qᵗ = ρqᵗ / ρ + specific_moisture[i, j, k] = qᵗ + end + + 𝒰₀ = diagnose_thermodynamic_state(i, j, k, grid, + formulation, + microphysics, + microphysical_fields, + constants, + specific_moisture) + + # Adjust the thermodynamic state if using a microphysics scheme + # that invokes saturation adjustment + 𝒰₁ = maybe_adjust_thermodynamic_state(𝒰₀, microphysics, microphysical_fields, qᵗ, constants) + + update_microphysical_fields!(microphysical_fields, microphysics, + i, j, k, grid, + ρ, 𝒰₁, constants) + + T = Thermodynamics.temperature(𝒰₁, constants) + @inbounds temperature[i, j, k] = T +end + +@kernel function _compute_potential_temperature_auxiliary_variables!(temperature, + potential_temperature, + specific_moisture, + grid, + constants, + formulation, + microphysics, + microphysical_fields, + potential_temperature_density, + moisture_density) + i, j, k = @index(Global, NTuple) + @inbounds begin - ρe = energy_density[i, j, k] + ρθ = potential_temperature_density[i, j, k] ρqᵗ = moisture_density[i, j, k] ρ = formulation.reference_state.density[i, j, k] - e = ρe / ρ + θ = ρθ / ρ qᵗ = ρqᵗ / ρ - specific_energy[i, j, k] = e + potential_temperature[i, j, k] = θ specific_moisture[i, j, k] = qᵗ end @@ -145,7 +186,6 @@ end microphysics, microphysical_fields, constants, - specific_energy, specific_moisture) # Adjust the thermodynamic state if using a microphysics scheme @@ -201,13 +241,10 @@ function compute_tendencies!(model::AnelasticModel) launch!(arch, grid, :xyz, compute_y_momentum_tendency!, Gρv, grid, v_args) launch!(arch, grid, :xyz, compute_z_momentum_tendency!, Gρw, grid, w_args) - specific_energy = model.formulation.thermodynamics.specific_energy - # Arguments common to energy density, moisture density, and tracer density tendencies: common_args = ( model.formulation, model.thermodynamic_constants, - specific_energy, model.specific_moisture, model.velocities, model.microphysics, @@ -218,18 +255,10 @@ function compute_tendencies!(model::AnelasticModel) model_fields) ##### - ##### Energy density tendency + ##### Thermodynamic density tendency (dispatches on thermodynamics type) ##### - ρe_args = ( - Val(1), - model.forcing.ρe, - model.advection.ρe, - common_args..., - model.temperature) - - Gρe = model.timestepper.Gⁿ.ρe - launch!(arch, grid, :xyz, compute_static_energy_tendency!, Gρe, grid, ρe_args) + compute_thermodynamic_tendency!(model, common_args) ##### ##### Moisture density tendency @@ -279,6 +308,11 @@ end @inbounds Gρe[i, j, k] = static_energy_tendency(i, j, k, grid, args...) end +@kernel function compute_potential_temperature_tendency!(Gρθ, grid, args) + i, j, k = @index(Global, NTuple) + @inbounds Gρθ[i, j, k] = potential_temperature_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...) diff --git a/src/Breeze.jl b/src/Breeze.jl index f9af842a..606d9e93 100644 --- a/src/Breeze.jl +++ b/src/Breeze.jl @@ -11,8 +11,8 @@ export ReferenceState, AnelasticFormulation, AtmosphereModel, - PotentialTemperature, - PotentialTemperatureField, + StaticEnergyThermodynamics, + LiquidIcePotentialTemperatureThermodynamics, TemperatureField, IdealGas, CondensedPhase, @@ -21,9 +21,14 @@ export SaturationAdjustment, MixedPhaseEquilibrium, WarmPhaseEquilibrium, + SaturationSpecificHumidity, + SaturationSpecificHumidityField, BulkMicrophysics, - energy_density, - specific_energy + static_energy_density, + static_energy, + total_energy, + potential_temperature_density, + liquid_ice_potential_temperature using Oceananigans: Oceananigans, @at, AnisotropicMinimumDissipation, Average, AveragedTimeInterval, BackgroundField, BetaPlane, Bounded, diff --git a/src/Microphysics/Microphysics.jl b/src/Microphysics/Microphysics.jl index 79350fd6..09cb5a9d 100644 --- a/src/Microphysics/Microphysics.jl +++ b/src/Microphysics/Microphysics.jl @@ -7,7 +7,9 @@ export MixedPhaseEquilibrium, WarmPhaseEquilibrium, BulkMicrophysics, - FourCategories + FourCategories, + SaturationSpecificHumidity, + SaturationSpecificHumidityField import ..AtmosphereModels: maybe_adjust_thermodynamic_state, @@ -20,5 +22,6 @@ import ..AtmosphereModels: include("saturation_adjustment.jl") include("bulk_microphysics.jl") +include("microphysics_diagnostics.jl") end # module Microphysics diff --git a/src/Microphysics/microphysics_diagnostics.jl b/src/Microphysics/microphysics_diagnostics.jl new file mode 100644 index 00000000..f74a22ef --- /dev/null +++ b/src/Microphysics/microphysics_diagnostics.jl @@ -0,0 +1,194 @@ +using Adapt: Adapt, adapt +using Oceananigans.AbstractOperations: KernelFunctionOperation +using Oceananigans.Fields: Field, Center + +using Breeze.Thermodynamics: + saturation_specific_humidity, + dry_air_gas_constant, + vapor_gas_constant, + density, + saturation_vapor_pressure + +import Oceananigans.Grids: prettysummary + +struct SaturationSpecificHumidityKernelFunction{μ, FL, M, MF, T, R, TH} + flavor :: FL + microphysics :: μ + microphysical_fields :: M + specific_moisture :: MF + temperature :: T + reference_state :: R + thermodynamic_constants :: TH +end + +prettysummary(kf::SaturationSpecificHumidityKernelFunction) = "$(kf.flavor) SaturationSpecificHumidityKernelFunction" + +Adapt.adapt_structure(to, k::SaturationSpecificHumidityKernelFunction) = + SaturationSpecificHumidityKernelFunction(adapt(to, k.flavor), + adapt(to, k.microphysics), + adapt(to, k.microphysical_fields), + adapt(to, k.specific_moisture), + adapt(to, k.temperature), + adapt(to, k.reference_state), + adapt(to, k.thermodynamic_constants)) + +const C = Center +const SaturationSpecificHumidity = KernelFunctionOperation{C, C, C, <:Any, <:Any, <:SaturationSpecificHumidityKernelFunction} + +struct Prognostic end +struct Equilibrium end +struct TotalMoisture end + +""" +$(TYPEDSIGNATURES) + +Return a `KernelFunctionOperation` representing the specified flavor +of *saturation specific humidity* ``qᵛ⁺`` which correpsonds to `model.microphysics`. +If `model.microphysics` is not a saturation adjustment scheme, then +a warm phase scheme is assumed which computes the saturation specific humidity +over a planar liquid surface. + +## Flavor options + +### `:prognostic` + +Returns the *saturation specific humidity* corresponding to the `model`'s prognostic state. +This is the same as the equilibrium saturation specific humidity for saturated conditions +and a model that uses saturation adjustment microphysics. + +### `:equilibrium` + +Returns the *saturation specific humidity* in saturated conditions, using the +`model.specific_moisture`. This is equivalent to the `:total_moisture` flavor +under saturated conditions with no condensate; or in other words, if `model.specific_moisture` happens +to be equal to the saturation specific humidity. + +### `:total_moisture` + +Returns *saturation specific humidity* in the case that the total specific moisture is +equal to the saturation specific humidity and there is no condensate. +This is useful for manufacturing perfectly saturated initial conditions. + +## Examples + +```jldoctestssh +using Breeze +grid = RectilinearGrid(size=(1, 1, 128), extent=(1e3, 1e3, 1e3)) +microphysics = SaturationAdjustment() +model = AtmosphereModel(grid; microphysics) +set!(model, θ=300) +qᵛ⁺ = SaturationSpecificHumidity(model, :prognostic) + +# output +KernelFunctionOperation at (Center, Center, Center) +├── grid: 1×1×128 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo +├── kernel_function: Breeze.Microphysics.Prognostic() SaturationSpecificHumidityKernelFunction +└── arguments: () +``` + +As `SaturationSpecificHumidity` it may be wrapped in `Field` to store the result +of its computation. For example, a `Field` representing the equilibrium saturation specific +humidity may be formed via, + +```jldoctest ssh +qᵛ = SaturationSpecificHumidity(model, :equilibrium) |> Field + +# output +1×1×128 Field{Center, Center, Center} on RectilinearGrid on CPU +├── grid: 1×1×128 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo +├── boundary conditions: FieldBoundaryConditions +│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing +├── operand: KernelFunctionOperation at (Center, Center, Center) +├── status: time=0.0 +└── data: 3×3×134 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:131) with eltype Float64 with indices 0:2×0:2×-2:131 + └── max=0.0361828, min=0.0224965, mean=0.028878 +``` + +We also provide a constructor and type alias for the `Field` itself. +For example, to build a `Field` representing the saturation specific humidity +in the case that the total specific moisture is exactly at saturation, + +```@example ssh +qᵗ = SaturationSpecificHumidityField(model, :total_moisture) + +# output +1×1×128 Field{Center, Center, Center} on RectilinearGrid on CPU +├── grid: 1×1×128 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo +├── boundary conditions: FieldBoundaryConditions +│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing +├── operand: KernelFunctionOperation at (Center, Center, Center) +├── status: time=0.0 +└── data: 3×3×134 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:131) with eltype Float64 with indices 0:2×0:2×-2:131 + └── max=0.0561539, min=0.0353807, mean=0.0451121 +``` +""" +function SaturationSpecificHumidity(model, flavor_symbol=:prognostic) + + flavor = if flavor_symbol == :prognostic + Prognostic() + elseif flavor_symbol == :equilibrium + Equilibrium() + elseif flavor_symbol == :total_moisture + TotalMoisture() + else + valid_flavors = (:prognostic, :equilibrium, :total_moisture) + throw(ArgumentError("Flavor $flavor_symbol is not one of the valid flavors $valid_flavors")) + end + + microphysics = if model.microphysics isa SaturationAdjustment + model.microphysics + else + SaturationAdjustment(equilibrium=WarmPhaseEquilibrium()) + end + + func = SaturationSpecificHumidityKernelFunction(flavor, + microphysics, + model.microphysical_fields, + model.specific_moisture, + model.temperature, + model.formulation.reference_state, + model.thermodynamic_constants) + + return KernelFunctionOperation{Center, Center, Center}(func, model.grid) +end + +@inline function saturation_total_specific_moisture(T, pᵣ, constants, equil) + surface = equilibrated_surface(equil, T) + pᵛ⁺ = saturation_vapor_pressure(T, constants, surface) + Rᵈ = dry_air_gas_constant(constants) + Rᵛ = vapor_gas_constant(constants) + δᵈᵛ = Rᵈ / Rᵛ - 1 + return pᵛ⁺ / (pᵣ + δᵈᵛ * pᵛ⁺) +end + +const AdjustmentSH = SaturationSpecificHumidityKernelFunction{<:SaturationAdjustment} + +function (d::AdjustmentSH)(i, j, k, grid) + @inbounds begin + pᵣ = d.reference_state.pressure[i, j, k] + ρᵣ = d.reference_state.density[i, j, k] + T = d.temperature[i, j, k] + end + + constants = d.thermodynamic_constants + equil = d.microphysics.equilibrium + + if d.flavor isa Prognostic + qᵗ = @inbounds d.specific_moisture[i, j, k] + q = compute_moisture_fractions(i, j, k, grid, d.microphysics, ρᵣ, qᵗ, d.microphysical_fields) + ρ = density(pᵣ, T, q, constants) + surface = equilibrated_surface(equil, T) + return saturation_specific_humidity(T, ρ, constants, surface) + + elseif d.flavor isa Equilibrium + qᵗ = @inbounds d.specific_moisture[i, j, k] + return equilibrium_saturation_specific_humidity(T, pᵣ, qᵗ, constants, equil) + + elseif d.flavor isa TotalMoisture + return saturation_total_specific_moisture(T, pᵣ, constants, equil) + + end +end + +const SaturationSpecificHumidityField = Field{C, C, C, <:SaturationSpecificHumidity} +SaturationSpecificHumidityField(model, flavor_symbol=:prognostic) = Field(SaturationSpecificHumidity(model, flavor_symbol)) diff --git a/src/Microphysics/saturation_adjustment.jl b/src/Microphysics/saturation_adjustment.jl index 4ef52f86..43a1defc 100644 --- a/src/Microphysics/saturation_adjustment.jl +++ b/src/Microphysics/saturation_adjustment.jl @@ -176,7 +176,7 @@ end return saturation_specific_humidity(T, ρ, constants, surface) end -@inline function adjustment_saturation_specific_humidity(T, pᵣ, qᵗ, constants, equil) +@inline function equilibrium_saturation_specific_humidity(T, pᵣ, qᵗ, constants, equil) surface = equilibrated_surface(equil, T) pᵛ⁺ = saturation_vapor_pressure(T, constants, surface) Rᵈ = dry_air_gas_constant(constants) @@ -188,7 +188,7 @@ end @inline function adjust_state(𝒰₀, T, constants, equilibrium) pᵣ = 𝒰₀.reference_pressure qᵗ = total_specific_moisture(𝒰₀) - qᵛ⁺ = adjustment_saturation_specific_humidity(T, pᵣ, qᵗ, constants, equilibrium) + qᵛ⁺ = equilibrium_saturation_specific_humidity(T, pᵣ, qᵗ, constants, equilibrium) q₁ = equilibrated_moisture_mass_fractions(T, qᵗ, qᵛ⁺, equilibrium) return with_moisture(𝒰₀, q₁) end diff --git a/src/MoistAirBuoyancies.jl b/src/MoistAirBuoyancies.jl index 57d59b68..702bdbc2 100644 --- a/src/MoistAirBuoyancies.jl +++ b/src/MoistAirBuoyancies.jl @@ -7,7 +7,7 @@ export SaturationField using ..Thermodynamics: - PotentialTemperatureState, + LiquidIcePotentialTemperatureState, MoistureMassFractions, total_specific_moisture, dry_air_gas_constant, @@ -124,7 +124,7 @@ const c = Center() z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) p₀ = mb.reference_state.base_pressure q = MoistureMassFractions(qᵗ) - 𝒰 = PotentialTemperatureState(θ, q, p₀, pᵣ) + 𝒰 = LiquidIcePotentialTemperatureState(θ, q, p₀, pᵣ) # Perform saturation adjustment T = compute_boussinesq_adjustment_temperature(𝒰, mb.thermodynamic_constants) @@ -177,7 +177,7 @@ r(T) ≡ T - θ Π - ℒˡᵣ qˡ / cᵖᵐ . Solution of ``r(T) = 0`` is found via the [secant method](https://en.wikipedia.org/wiki/Secant_method). """ -@inline function compute_boussinesq_adjustment_temperature(𝒰₀::PotentialTemperatureState{FT}, constants) where FT +@inline function compute_boussinesq_adjustment_temperature(𝒰₀::LiquidIcePotentialTemperatureState{FT}, constants) where FT θ = 𝒰₀.potential_temperature θ == 0 && return zero(FT) @@ -198,7 +198,7 @@ Solution of ``r(T) = 0`` is found via the [secant method](https://en.wikipedia.o # has to be modified to consistently include the liquid mass fraction. # Subsequent computations will assume that the specific humidity # is given by the saturation specific humidity, eg ``qᵛ = qᵛ⁺``. - qᵛ⁺₁ = adjustment_saturation_specific_humidity(T₁, 𝒰₁, constants) + qᵛ⁺₁ = equilibrium_saturation_specific_humidity(T₁, 𝒰₁, constants) qˡ₁ = qᵗ - qᵛ⁺₁ q₁ = MoistureMassFractions(qᵛ⁺₁, qˡ₁) 𝒰₁ = with_moisture(𝒰₀, q₁) @@ -246,7 +246,7 @@ end # This consideration culminates in a new expression for saturation specific humidity # used below, and also written in Pressel et al 2015, equation 37. # (There is an error in the description below it, but the equation 37 is correct.) -@inline function adjustment_saturation_specific_humidity(T, 𝒰, constants) +@inline function equilibrium_saturation_specific_humidity(T, 𝒰, constants) pᵛ⁺ = saturation_vapor_pressure(T, constants, constants.liquid) pᵣ = 𝒰.reference_pressure qᵗ = total_specific_moisture(𝒰) @@ -257,7 +257,7 @@ end end @inline function adjust_state(𝒰₀, T, constants) - qᵛ⁺ = adjustment_saturation_specific_humidity(T, 𝒰₀, constants) + qᵛ⁺ = equilibrium_saturation_specific_humidity(T, 𝒰₀, constants) qᵗ = total_specific_moisture(𝒰₀) qˡ = max(0, qᵗ - qᵛ⁺) qᵛ = qᵗ - qˡ @@ -293,7 +293,7 @@ const c = Center() z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) p₀ = mb.reference_state.base_pressure q = MoistureMassFractions(qᵗᵢ) - 𝒰 = PotentialTemperatureState(θᵢ, q, p₀, pᵣ) + 𝒰 = LiquidIcePotentialTemperatureState(θᵢ, q, p₀, pᵣ) return compute_boussinesq_adjustment_temperature(𝒰, mb.thermodynamic_constants) end @@ -370,12 +370,12 @@ Adapt.adapt_structure(to, ck::CondensateKernel) = CondensateKernel(adapt(to, ck. z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) p₀ = mb.reference_state.base_pressure q = MoistureMassFractions(qᵗᵢ) - 𝒰 = PotentialTemperatureState(Tᵢ, q, p₀, pᵣ) + 𝒰 = LiquidIcePotentialTemperatureState(Tᵢ, q, p₀, pᵣ) Π = exner_function(𝒰, mb.thermodynamic_constants) Tᵢ <= Π * θᵢ + 10 * eps(Tᵢ) && return zero(qᵗᵢ) # Next assume a saturation value - qᵛ⁺ = adjustment_saturation_specific_humidity(Tᵢ, 𝒰, mb.thermodynamic_constants) + qᵛ⁺ = equilibrium_saturation_specific_humidity(Tᵢ, 𝒰, mb.thermodynamic_constants) return max(0, qᵗᵢ - qᵛ⁺) end diff --git a/src/Thermodynamics/dynamic_states.jl b/src/Thermodynamics/dynamic_states.jl index f71858e7..e5eb982c 100644 --- a/src/Thermodynamics/dynamic_states.jl +++ b/src/Thermodynamics/dynamic_states.jl @@ -2,49 +2,80 @@ abstract type AbstractThermodynamicState{FT} end @inline Base.eltype(::AbstractThermodynamicState{FT}) where FT = FT -struct PotentialTemperatureState{FT} <: AbstractThermodynamicState{FT} +@inline function density(𝒰::AbstractThermodynamicState, constants) + pᵣ = 𝒰.reference_pressure + T = temperature(𝒰, constants) + q = 𝒰.moisture_mass_fractions + return density(pᵣ, T, q, constants) +end + +@inline function saturation_specific_humidity(𝒰::AbstractThermodynamicState, constants, equil) + T = temperature(𝒰, constants) + ρ = density(𝒰, constants) + return saturation_specific_humidity(T, ρ, constants, equil) +end + +##### +##### Liquid-ice potential temperature state +##### + +struct LiquidIcePotentialTemperatureState{FT} <: AbstractThermodynamicState{FT} potential_temperature :: FT moisture_mass_fractions :: MoistureMassFractions{FT} base_pressure :: FT reference_pressure :: FT end -@inline is_absolute_zero(𝒰::PotentialTemperatureState) = 𝒰.potential_temperature == 0 +@inline is_absolute_zero(𝒰::LiquidIcePotentialTemperatureState) = 𝒰.potential_temperature == 0 -@inline function exner_function(𝒰::PotentialTemperatureState, thermo::ThermodynamicConstants) +@inline function exner_function(𝒰::LiquidIcePotentialTemperatureState, constants::ThermodynamicConstants) q = 𝒰.moisture_mass_fractions - Rᵐ = mixture_gas_constant(q, thermo) - cᵖᵐ = mixture_heat_capacity(q, thermo) + Rᵐ = mixture_gas_constant(q, constants) + cᵖᵐ = mixture_heat_capacity(q, constants) pᵣ = 𝒰.reference_pressure p₀ = 𝒰.base_pressure return (pᵣ / p₀)^(Rᵐ / cᵖᵐ) end -@inline total_specific_moisture(state::PotentialTemperatureState) = +@inline total_specific_moisture(state::LiquidIcePotentialTemperatureState) = total_specific_moisture(state.moisture_mass_fractions) -@inline with_moisture(𝒰::PotentialTemperatureState{FT}, q::MoistureMassFractions{FT}) where FT = - PotentialTemperatureState{FT}(𝒰.potential_temperature, q, 𝒰.base_pressure, 𝒰.reference_pressure) +@inline with_moisture(𝒰::LiquidIcePotentialTemperatureState{FT}, q::MoistureMassFractions{FT}) where FT = + LiquidIcePotentialTemperatureState{FT}(𝒰.potential_temperature, q, 𝒰.base_pressure, 𝒰.reference_pressure) -@inline function temperature(𝒰::PotentialTemperatureState, thermo::ThermodynamicConstants) +@inline function temperature(𝒰::LiquidIcePotentialTemperatureState, constants::ThermodynamicConstants) θ = 𝒰.potential_temperature - Π = exner_function(𝒰, thermo) + Π = exner_function(𝒰, constants) q = 𝒰.moisture_mass_fractions - cᵖᵐ = mixture_heat_capacity(q, thermo) - ℒˡᵣ = thermo.liquid.reference_latent_heat - ℒⁱᵣ = thermo.ice.reference_latent_heat + cᵖᵐ = mixture_heat_capacity(q, constants) + ℒˡᵣ = constants.liquid.reference_latent_heat + ℒⁱᵣ = constants.ice.reference_latent_heat qˡ = q.liquid qⁱ = q.ice - return Π*θ + (ℒˡᵣ*qˡ + ℒⁱᵣ*qⁱ) / cᵖᵐ + return Π * θ + (ℒˡᵣ * qˡ + ℒⁱᵣ * qⁱ) / cᵖᵐ end -@inline function density(𝒰::PotentialTemperatureState, thermo) +@inline function with_temperature(𝒰::LiquidIcePotentialTemperatureState, T, constants) + Π = exner_function(𝒰, constants) + q = 𝒰.moisture_mass_fractions + cᵖᵐ = mixture_heat_capacity(q, constants) + ℒˡᵣ = constants.liquid.reference_latent_heat + ℒⁱᵣ = constants.ice.reference_latent_heat + qˡ = q.liquid + qⁱ = q.ice + + θ = (T - (ℒˡᵣ * qˡ + ℒⁱᵣ * qⁱ) / cᵖᵐ) / Π + + return LiquidIcePotentialTemperatureState(θ, q, 𝒰.base_pressure, 𝒰.reference_pressure) +end + +@inline function density(𝒰::LiquidIcePotentialTemperatureState, constants) pᵣ = 𝒰.reference_pressure - T = temperature(𝒰, thermo) + T = temperature(𝒰, constants) q = 𝒰.moisture_mass_fractions - return density(pᵣ, T, q, thermo) + return density(pᵣ, T, q, constants) end ##### @@ -64,32 +95,34 @@ end @inline with_moisture(𝒰::StaticEnergyState{FT}, q::MoistureMassFractions{FT}) where FT = StaticEnergyState{FT}(𝒰.static_energy, q, 𝒰.height, 𝒰.reference_pressure) -@inline function temperature(𝒰::StaticEnergyState, thermo::ThermodynamicConstants) +@inline function temperature(𝒰::StaticEnergyState, constants::ThermodynamicConstants) e = 𝒰.static_energy q = 𝒰.moisture_mass_fractions - cᵖᵐ = mixture_heat_capacity(q, thermo) + cᵖᵐ = mixture_heat_capacity(q, constants) - g = thermo.gravitational_acceleration + g = constants.gravitational_acceleration z = 𝒰.height - ℒˡᵣ = thermo.liquid.reference_latent_heat - ℒⁱᵣ = thermo.ice.reference_latent_heat + ℒˡᵣ = constants.liquid.reference_latent_heat + ℒⁱᵣ = constants.ice.reference_latent_heat qˡ = q.liquid qⁱ = q.ice # e = cᵖᵐ * T + g * z - ℒˡᵣ * qˡ - ℒⁱᵣ * qⁱ - return (e - g*z + ℒˡᵣ*qˡ + ℒⁱᵣ*qⁱ) / cᵖᵐ + return (e - g * z + ℒˡᵣ * qˡ + ℒⁱᵣ * qⁱ) / cᵖᵐ end -@inline function density(𝒰::AbstractThermodynamicState, thermo) - pᵣ = 𝒰.reference_pressure - T = temperature(𝒰, thermo) +@inline function with_temperature(𝒰::StaticEnergyState, T, constants) q = 𝒰.moisture_mass_fractions - return density(pᵣ, T, q, thermo) -end + cᵖᵐ = mixture_heat_capacity(q, constants) + g = constants.gravitational_acceleration + z = 𝒰.height + ℒˡᵣ = constants.liquid.reference_latent_heat + ℒⁱᵣ = constants.ice.reference_latent_heat + qˡ = q.liquid + qⁱ = q.ice -@inline function saturation_specific_humidity(𝒰::AbstractThermodynamicState, thermo, equil) - T = temperature(𝒰, thermo) - ρ = density(𝒰, thermo) - return saturation_specific_humidity(T, ρ, thermo, equil) -end + e = cᵖᵐ * T + g * z - ℒˡᵣ * qˡ - ℒⁱᵣ * qⁱ + + return StaticEnergyState(e, q, z, 𝒰.reference_pressure) +end \ No newline at end of file diff --git a/src/Thermodynamics/thermodynamics_constants.jl b/src/Thermodynamics/thermodynamics_constants.jl index a2926f62..1b5b7d92 100644 --- a/src/Thermodynamics/thermodynamics_constants.jl +++ b/src/Thermodynamics/thermodynamics_constants.jl @@ -12,6 +12,7 @@ A struct representing an ideal gas with molar mass and specific heat capacity. # Examples ```jldoctest +using Breeze dry_air = IdealGas(molar_mass=0.02897, heat_capacity=1005) # output @@ -67,7 +68,7 @@ Adapt.adapt_structure(to, pt::CondensedPhase) = """ $(TYPEDSIGNATURES) -Returns `CondensedPhase` with specified parameters converted to `FT`. +Return `CondensedPhase` with specified parameters converted to `FT`. Two examples of `CondensedPhase` are liquid and ice. When matter is converted from vapor to liquid, water molecules in the diff --git a/test/atmosphere_model_unit.jl b/test/atmosphere_model_unit.jl index 26a80037..3891ed04 100644 --- a/test/atmosphere_model_unit.jl +++ b/test/atmosphere_model_unit.jl @@ -7,10 +7,10 @@ using Test grid = RectilinearGrid(default_arch, FT; size=(8, 8, 8), x=(0, 1_000), y=(0, 1_000), z=(0, 1_000)) constants = ThermodynamicConstants(FT) - for p₀ in (101325, 100000), θ₀ in (288, 300) - @testset let p₀ = p₀, θ₀ = θ₀ + for p₀ in (101325, 100000), θ₀ in (288, 300), thermodynamics in (:LiquidIcePotentialTemperature, :StaticEnergy) + @testset let p₀ = p₀, θ₀ = θ₀, thermodynamics = thermodynamics reference_state = ReferenceState(grid, constants, base_pressure=p₀, potential_temperature=θ₀) - formulation = AnelasticFormulation(reference_state) + formulation = AnelasticFormulation(reference_state; thermodynamics) model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation) # test set! @@ -19,15 +19,17 @@ using Test ρeᵢ = ρᵣ * cᵖᵈ * θ₀ set!(model; θ = θ₀) - ρe₁ = deepcopy(energy_density(model)) + ρe₁ = deepcopy(static_energy_density(model)) + θ₁ = deepcopy(liquid_ice_potential_temperature(model)) set!(model; ρe = ρeᵢ) - @test energy_density(model) ≈ ρe₁ + @test static_energy_density(model) ≈ ρe₁ + @test liquid_ice_potential_temperature(model) ≈ θ₁ end end end -@testset "PotentialTemperatureField (no microphysics) [$(FT)]" for FT in (Float32, Float64) +@testset "liquid_ice_potential_temperature no microphysics) [$(FT)]" for FT in (Float32, Float64), thermodynamics in (:LiquidIcePotentialTemperature, :StaticEnergy) Oceananigans.defaults.FloatType = FT grid = RectilinearGrid(default_arch; size=(8, 8, 8), x=(0, 1_000), y=(0, 1_000), z=(0, 1_000)) constants = ThermodynamicConstants() @@ -35,7 +37,7 @@ end p₀ = FT(101325) θ₀ = FT(300) reference_state = ReferenceState(grid, constants, base_pressure=p₀, potential_temperature=θ₀) - formulation = AnelasticFormulation(reference_state) + formulation = AnelasticFormulation(reference_state; thermodynamics) model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation) # Initialize with potential temperature and dry air @@ -43,11 +45,11 @@ end set!(θᵢ, (x, y, z) -> θ₀ + rand()) set!(model; θ=θᵢ) - θ_model = Breeze.AtmosphereModels.PotentialTemperatureField(model) + θ_model = liquid_ice_potential_temperature(model) |> Field @test θ_model ≈ θᵢ end -@testset "Saturation and PotentialTemperatureField (WarmPhase) [$(FT)]" for FT in (Float32, Float64) +@testset "Saturation and LiquidIcePotentialTemperatureField (WarmPhase) [$(FT)]" for FT in (Float32, Float64), thermodynamics in (:LiquidIcePotentialTemperature, :StaticEnergy) Oceananigans.defaults.FloatType = FT grid = RectilinearGrid(default_arch; size=(8, 8, 8), x=(0, 1_000), y=(0, 1_000), z=(0, 1_000)) constants = ThermodynamicConstants() @@ -55,7 +57,7 @@ end p₀ = FT(101325) θ₀ = FT(300) reference_state = ReferenceState(grid, constants, base_pressure=p₀, potential_temperature=θ₀) - formulation = AnelasticFormulation(reference_state) + formulation = AnelasticFormulation(reference_state; thermodynamics) microphysics = SaturationAdjustment() model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation, microphysics) @@ -63,7 +65,7 @@ end set!(model; θ=θ₀) # Check SaturationSpecificHumidityField matches direct thermodynamics - qᵛ⁺ = Breeze.AtmosphereModels.SaturationSpecificHumidityField(model) + qᵛ⁺ = SaturationSpecificHumidityField(model) # Sample mid-level cell _, _, Nz = size(grid) @@ -88,72 +90,82 @@ end p₀ = FT(101325) θ₀ = FT(300) reference_state = ReferenceState(grid, constants, base_pressure=p₀, potential_temperature=θ₀) - formulation = AnelasticFormulation(reference_state) + potential_temperature_formulation = AnelasticFormulation(reference_state; thermodynamics=:LiquidIcePotentialTemperature) + static_energy_formulation = AnelasticFormulation(reference_state; thermodynamics=:StaticEnergy) @testset "Default advection schemes" begin - model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation) - @test model.advection.momentum isa Centered - @test model.advection.ρe isa Centered - @test model.advection.ρqᵗ isa Centered + static_energy_model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation=static_energy_formulation) + potential_temperature_model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation=potential_temperature_formulation) - time_step!(model, 1) - @test true + @test static_energy_model.advection.ρe isa Centered + @test potential_temperature_model.advection.ρθ isa Centered + + for model in (static_energy_model, potential_temperature_model) + @test model.advection.momentum isa Centered + @test model.advection.ρqᵗ isa Centered + time_step!(model, 1) + end end @testset "Unified advection parameter" begin - model_weno = AtmosphereModel(grid; thermodynamic_constants=constants, formulation, advection=WENO()) - @test model_weno.advection.momentum isa WENO - @test model_weno.advection.ρe isa WENO - @test model_weno.advection.ρqᵗ isa WENO - time_step!(model_weno, 1) - @test true - - model_centered = AtmosphereModel(grid; thermodynamic_constants=constants, formulation, advection=Centered(order=4)) - @test model_centered.advection.momentum isa Centered - @test model_centered.advection.ρe isa Centered - time_step!(model_centered, 1) - @test true + static_energy_model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation=static_energy_formulation, advection=WENO()) + potential_temperature_model= AtmosphereModel(grid; thermodynamic_constants=constants, formulation=potential_temperature_formulation, advection=WENO()) + + @test static_energy_model.advection.ρe isa WENO + @test potential_temperature_model.advection.ρθ isa WENO + + for model in (static_energy_model, potential_temperature_model) + @test model.advection.momentum isa WENO + @test model.advection.ρqᵗ isa WENO + time_step!(model, 1) + end end @testset "Separate momentum and tracer advection" begin - model = AtmosphereModel(grid; - thermodynamic_constants = constants, - formulation, - momentum_advection = WENO(), - scalar_advection = Centered(order=2)) - @test model.advection.momentum isa WENO - @test model.advection.ρe isa Centered - @test model.advection.ρqᵗ isa Centered - time_step!(model, 1) - @test true + kw = (thermodynamic_constants=constants, momentum_advection = WENO(), scalar_advection = Centered(order=2)) + static_energy_model = AtmosphereModel(grid; formulation=static_energy_formulation, kw...) + potential_temperature_model = AtmosphereModel(grid; formulation=potential_temperature_formulation, kw...) + + @test static_energy_model.advection.ρe isa Centered + @test potential_temperature_model.advection.ρθ isa Centered + + for model in (static_energy_model, potential_temperature_model) + @test model.advection.momentum isa WENO + @test model.advection.ρqᵗ isa Centered + time_step!(model, 1) + end end @testset "Tracer advection with user tracers" begin - model = AtmosphereModel(grid; - thermodynamic_constants = constants, - formulation, - tracers = :c, - scalar_advection = UpwindBiased(order=1)) - @test model.advection.momentum isa Centered - @test model.advection.ρe isa UpwindBiased - @test model.advection.ρqᵗ isa UpwindBiased - @test model.advection.c isa UpwindBiased - time_step!(model, 1) - @test true + kw = (thermodynamic_constants=constants, tracers = :c, scalar_advection = UpwindBiased(order=1)) + static_energy_model = AtmosphereModel(grid; formulation=static_energy_formulation, kw...) + potential_temperature_model = AtmosphereModel(grid; formulation=potential_temperature_formulation, kw...) + + @test static_energy_model.advection.ρe isa UpwindBiased + @test potential_temperature_model.advection.ρθ isa UpwindBiased + + for model in (static_energy_model, potential_temperature_model) + @test model.advection.momentum isa Centered + @test model.advection.ρqᵗ isa UpwindBiased + @test model.advection.c isa UpwindBiased + time_step!(model, 1) + end end @testset "Mixed configuration with tracers" begin - model = AtmosphereModel(grid; - thermodynamic_constants = constants, - formulation, - tracers = :c, - momentum_advection = WENO(), - scalar_advection = Centered(order=2)) - @test model.advection.momentum isa WENO - @test model.advection.ρe isa Centered - @test model.advection.ρqᵗ isa Centered - @test model.advection.c isa Centered - time_step!(model, 1) - @test true + scalar_advection = (; c=Centered(order=2), ρqᵗ=WENO()) + kw = (thermodynamic_constants=constants, tracers = :c, momentum_advection = WENO(), scalar_advection) + static_energy_model = AtmosphereModel(grid; formulation=static_energy_formulation, kw...) + potential_temperature_model = AtmosphereModel(grid; formulation=potential_temperature_formulation, kw...) + + @test static_energy_model.advection.ρe isa Centered + @test potential_temperature_model.advection.ρθ isa Centered + + for model in (static_energy_model, potential_temperature_model) + @test model.advection.momentum isa WENO + @test model.advection.ρqᵗ isa WENO + @test model.advection.c isa Centered + time_step!(model, 1) + end end end diff --git a/test/forcing.jl b/test/forcing.jl index d1ce4c19..ffce7310 100644 --- a/test/forcing.jl +++ b/test/forcing.jl @@ -38,9 +38,9 @@ increment_tolerance(::Type{Float64}) = 1e-10 e_forcing = (; ρe=forcing) model = setup_forcing_model(grid, e_forcing) - ρe_before = deepcopy(energy_density(model)) + ρe_before = deepcopy(static_energy_density(model)) time_step!(model, Δt) - @test maximum(energy_density(model)) ≈ maximum(ρe_before) + Δt + @test maximum(static_energy_density(model)) ≈ maximum(ρe_before) + Δt q_forcing = (; ρqᵗ=forcing) model = setup_forcing_model(grid, q_forcing) diff --git a/test/saturation_adjustment.jl b/test/saturation_adjustment.jl index b71af94a..07796d77 100644 --- a/test/saturation_adjustment.jl +++ b/test/saturation_adjustment.jl @@ -5,7 +5,7 @@ using Test using Breeze.Thermodynamics: MoistureMassFractions, - PotentialTemperatureState, + LiquidIcePotentialTemperatureState, StaticEnergyState, exner_function, density, @@ -18,14 +18,16 @@ using Breeze.MoistAirBuoyancies: compute_boussinesq_adjustment_temperature using Breeze.Microphysics: compute_temperature using Breeze.Microphysics: - adjustment_saturation_specific_humidity + equilibrium_saturation_specific_humidity solver_tol(::Type{Float64}) = 1e-6 solver_tol(::Type{Float32}) = 1e-3 test_tol(FT::Type{Float64}) = 10 * sqrt(solver_tol(FT)) test_tol(FT::Type{Float32}) = sqrt(solver_tol(FT)) -@testset "Warm-phase saturation adjustment (AtmosphereModel) [$(FT)]" for FT in (Float32, Float64) +test_thermodynamics = (:StaticEnergy, :LiquidIcePotentialTemperature) + +@testset "Warm-phase saturation adjustment [$(FT)]" for FT in (Float32, Float64) grid = RectilinearGrid(default_arch, FT; size=(1, 1, 1), x=(0, 1), y=(0, 1), z=(0, 1)) constants = ThermodynamicConstants(FT) reference_state = ReferenceState(grid, constants; base_pressure=101325, potential_temperature=288) @@ -57,38 +59,40 @@ test_tol(FT::Type{Float32}) = sqrt(solver_tol(FT)) @test compute_temperature(𝒰₁, microphysics, constants) ≈ T₁ atol=atol @test compute_temperature(𝒰₁, nothing, constants) ≈ T₁ atol=atol - formulation = AnelasticFormulation(reference_state) - model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation, microphysics) - ρᵣ = @allowscalar first(reference_state.density) - - # Many more tests that touch saturated conditions - for T₂ in 270:4:320, qᵗ₂ in 1e-2:2e-3:5e-2 - @testset let T₂=T₂, qᵗ₂=qᵗ₂ - T₂ = convert(FT, T₂) - qᵗ₂ = convert(FT, qᵗ₂) - qᵛ⁺₂ = adjustment_saturation_specific_humidity(T₂, pᵣ, qᵗ₂, constants, microphysics.equilibrium) - @test qᵛ⁺₂ isa FT - - if qᵗ₂ > qᵛ⁺₂ # saturated conditions - qˡ₂ = qᵗ₂ - qᵛ⁺₂ - q₂ = MoistureMassFractions(qᵛ⁺₂, qˡ₂) - cᵖᵐ = mixture_heat_capacity(q₂, constants) - ℒˡᵣ = constants.liquid.reference_latent_heat - e₂ = cᵖᵐ * T₂ + g * z - ℒˡᵣ * qˡ₂ - - 𝒰₂ = StaticEnergyState(e₂, q₂, z, pᵣ) - T★ = compute_temperature(𝒰₂, microphysics, constants) - @test T★ ≈ T₂ atol=atol - - # Parcel test for AtmosphereModel - set!(model, ρe = ρᵣ * e₂, qᵗ = qᵗ₂) - T★ = @allowscalar first(model.temperature) - qᵛ = @allowscalar first(model.microphysical_fields.qᵛ) - qˡ = @allowscalar first(model.microphysical_fields.qˡ) + @testset "AtmosphereModel with $thermodynamics thermodynamics [$FT]" for thermodynamics in test_thermodynamics + formulation = AnelasticFormulation(reference_state; thermodynamics) + model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation, microphysics) + ρᵣ = @allowscalar first(reference_state.density) + + # Many more tests that touch saturated conditions + for T₂ in 270:4:320, qᵗ₂ in 1e-2:2e-3:5e-2 + @testset let T₂=T₂, qᵗ₂=qᵗ₂ + T₂ = convert(FT, T₂) + qᵗ₂ = convert(FT, qᵗ₂) + qᵛ⁺₂ = equilibrium_saturation_specific_humidity(T₂, pᵣ, qᵗ₂, constants, microphysics.equilibrium) + @test qᵛ⁺₂ isa FT + + if qᵗ₂ > qᵛ⁺₂ # saturated conditions + qˡ₂ = qᵗ₂ - qᵛ⁺₂ + q₂ = MoistureMassFractions(qᵛ⁺₂, qˡ₂) + cᵖᵐ = mixture_heat_capacity(q₂, constants) + ℒˡᵣ = constants.liquid.reference_latent_heat + e₂ = cᵖᵐ * T₂ + g * z - ℒˡᵣ * qˡ₂ + + 𝒰₂ = StaticEnergyState(e₂, q₂, z, pᵣ) + T★ = compute_temperature(𝒰₂, microphysics, constants) + @test T★ ≈ T₂ atol=atol + + # Parcel test for AtmosphereModel + set!(model, ρe = ρᵣ * e₂, qᵗ = qᵗ₂) + T★ = @allowscalar first(model.temperature) + qᵛ = @allowscalar first(model.microphysical_fields.qᵛ) + qˡ = @allowscalar first(model.microphysical_fields.qˡ) - @test T★ ≈ T₂ atol=atol - @test qᵛ ≈ qᵛ⁺₂ atol=atol - @test qˡ ≈ qˡ₂ atol=atol + @test T★ ≈ T₂ atol=atol + @test qᵛ ≈ qᵛ⁺₂ atol=atol + @test qˡ ≈ qˡ₂ atol=atol + end end end end @@ -101,28 +105,29 @@ end @testset "Mixed-phase saturation adjustment (AtmosphereModel) [$(FT)]" for FT in (Float32, Float64) grid = RectilinearGrid(default_arch, FT; size=(1, 1, 1), x=(0, 1), y=(0, 1), z=(0, 1)) + constants = ThermodynamicConstants(FT) + ℒˡᵣ = constants.liquid.reference_latent_heat + ℒⁱᵣ = constants.ice.reference_latent_heat + g = constants.gravitational_acceleration + z = zero(FT) + reference_state = ReferenceState(grid, constants; base_pressure=101325, potential_temperature=288) - atol = test_tol(FT) + pᵣ = @allowscalar first(reference_state.pressure) + ρᵣ = @allowscalar first(reference_state.density) - Tᶠ = FT(273.15) # Freezing temperature + atol = test_tol(FT) Tʰ = FT(233.15) # Homogeneous ice nucleation temperature + Tᶠ = FT(273.15) # Freezing temperature equilibrium = MixedPhaseEquilibrium(FT; freezing_temperature=Tᶠ, homogeneous_ice_nucleation_temperature=Tʰ) - microphysics = SaturationAdjustment(FT; tolerance=solver_tol(FT), equilibrium=equilibrium) - formulation = AnelasticFormulation(reference_state) - model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation, microphysics) - ρᵣ = @allowscalar first(reference_state.density) + microphysics = SaturationAdjustment(FT; tolerance=solver_tol(FT), equilibrium) - # Sample a single cell - pᵣ = @allowscalar first(reference_state.pressure) - g = constants.gravitational_acceleration - z = zero(FT) - ℒˡᵣ = constants.liquid.reference_latent_heat - ℒⁱᵣ = constants.ice.reference_latent_heat + @testset "AtmosphereModel with $thermodynamics thermodynamics [$FT]" for thermodynamics in test_thermodynamics + formulation = AnelasticFormulation(reference_state; thermodynamics) + model = AtmosphereModel(grid; thermodynamic_constants=constants, formulation, microphysics) - # Test 1: Constructor and equilibrated_surface utility - @testset "Constructor and equilibrated_surface" begin + # Test 1: Constructor and equilibrated_surface utility @test microphysics isa SaturationAdjustment @test microphysics.equilibrium isa MixedPhaseEquilibrium{FT} @test microphysics.equilibrium.freezing_temperature == Tᶠ @@ -133,129 +138,132 @@ end @test model.microphysics.equilibrium.freezing_temperature == Tᶠ @test model.microphysics.equilibrium.homogeneous_ice_nucleation_temperature == Tʰ - # Test equilibrated_surface at different temperatures - surface_above_freezing = Breeze.Microphysics.equilibrated_surface(equilibrium, FT(300)) - @test surface_above_freezing isa PlanarMixedPhaseSurface{FT} - @test surface_above_freezing.liquid_fraction == 1 # Above freezing, all liquid + @testset "equilibrated_surface" begin + # Test equilibrated_surface at different temperatures + surface_above_freezing = Breeze.Microphysics.equilibrated_surface(equilibrium, FT(300)) + @test surface_above_freezing isa PlanarMixedPhaseSurface{FT} + @test surface_above_freezing.liquid_fraction == 1 # Above freezing, all liquid + + surface_below_homogeneous_ice_nucleation = Breeze.Microphysics.equilibrated_surface(equilibrium, FT(200)) + @test surface_below_homogeneous_ice_nucleation isa PlanarMixedPhaseSurface{FT} + @test surface_below_homogeneous_ice_nucleation.liquid_fraction == 0 # Below homogeneous nucleation, all ice + + T_mid = FT(253.15) # Midway between Tᶠ and Tʰ + surface_midway = Breeze.Microphysics.equilibrated_surface(equilibrium, T_mid) + @test surface_midway isa PlanarMixedPhaseSurface{FT} + λ_expected = test_liquid_fraction(T_mid, Tᶠ, Tʰ) + @test surface_midway.liquid_fraction ≈ λ_expected + end - surface_below_homogeneous_ice_nucleation = Breeze.Microphysics.equilibrated_surface(equilibrium, FT(200)) - @test surface_below_homogeneous_ice_nucleation isa PlanarMixedPhaseSurface{FT} - @test surface_below_homogeneous_ice_nucleation.liquid_fraction == 0 # Below homogeneous nucleation, all ice + # Test 2: Temperatures above freezing - should match warm phase behavior + @testset "Temperatures above freezing (warm phase equivalence)" begin + T_warm = FT(300) + qᵗ = FT(0.02) + qᵛ⁺ = equilibrium_saturation_specific_humidity(T_warm, pᵣ, qᵗ, constants, equilibrium) + atol = test_tol(FT) - T_mid = FT(253.15) # Midway between Tᶠ and Tʰ - surface_midway = Breeze.Microphysics.equilibrated_surface(equilibrium, T_mid) - @test surface_midway isa PlanarMixedPhaseSurface{FT} - λ_expected = test_liquid_fraction(T_mid, Tᶠ, Tʰ) - @test surface_midway.liquid_fraction ≈ λ_expected - end + if qᵗ > qᵛ⁺ # saturated conditions + # For warm temperatures, all condensate should be liquid + qˡ = qᵗ - qᵛ⁺ + q = MoistureMassFractions(qᵛ⁺, qˡ) + cᵖᵐ = mixture_heat_capacity(q, constants) + e = cᵖᵐ * T_warm + g * z - ℒˡᵣ * qˡ - # Test 2: Temperatures above freezing - should match warm phase behavior - @testset "Temperatures above freezing (warm phase equivalence)" begin - T_warm = FT(300) - qᵗ = FT(0.02) - qᵛ⁺ = adjustment_saturation_specific_humidity(T_warm, pᵣ, qᵗ, constants, equilibrium) - atol = test_tol(FT) + 𝒰 = StaticEnergyState(e, q, z, pᵣ) + T★ = compute_temperature(𝒰, microphysics, constants) + @test T★ ≈ T_warm atol=atol - if qᵗ > qᵛ⁺ # saturated conditions - # For warm temperatures, all condensate should be liquid - qˡ = qᵗ - qᵛ⁺ - q = MoistureMassFractions(qᵛ⁺, qˡ) - cᵖᵐ = mixture_heat_capacity(q, constants) - e = cᵖᵐ * T_warm + g * z - ℒˡᵣ * qˡ - - 𝒰 = StaticEnergyState(e, q, z, pᵣ) - T★ = compute_temperature(𝒰, microphysics, constants) - @test T★ ≈ T_warm atol=atol - - # Parcel test for AtmosphereModel - set!(model, ρe = ρᵣ * e, qᵗ = qᵗ) - T★ = @allowscalar first(model.temperature) - qᵛm = @allowscalar first(model.microphysical_fields.qᵛ) - qˡm = @allowscalar first(model.microphysical_fields.qˡ) - qⁱm = @allowscalar first(model.microphysical_fields.qⁱ) - - @test T★ ≈ T_warm atol=atol - @test qᵛm ≈ qᵛ⁺ atol=atol - @test qˡm ≈ qˡ atol=atol - @test qⁱm ≈ zero(FT) atol=atol + # Parcel test for AtmosphereModel + set!(model, ρe = ρᵣ * e, qᵗ = qᵗ) + T★ = @allowscalar first(model.temperature) + qᵛm = @allowscalar first(model.microphysical_fields.qᵛ) + qˡm = @allowscalar first(model.microphysical_fields.qˡ) + qⁱm = @allowscalar first(model.microphysical_fields.qⁱ) + + @test T★ ≈ T_warm atol=atol + @test qᵛm ≈ qᵛ⁺ atol=atol + @test qˡm ≈ qˡ atol=atol + @test qⁱm ≈ zero(FT) atol=atol + end end - end - - # Test 3: Temperatures below homogeneous ice nucleation - all ice - @testset "Temperatures below homogeneous ice nucleation (all ice)" begin - T_cold = FT(220) # Below Tʰ - qᵗ = FT(0.01) - qᵛ⁺ = adjustment_saturation_specific_humidity(T_cold, pᵣ, qᵗ, constants, equilibrium) - atol = test_tol(FT) - if qᵗ > qᵛ⁺ # saturated conditions - # All condensate should be ice - qⁱ = qᵗ - qᵛ⁺ - q = MoistureMassFractions(qᵛ⁺, zero(FT), qⁱ) - cᵖᵐ = mixture_heat_capacity(q, constants) - e = cᵖᵐ * T_cold + g * z - ℒⁱᵣ * qⁱ - - 𝒰 = StaticEnergyState(e, q, z, pᵣ) - T★ = compute_temperature(𝒰, microphysics, constants) - @test T★ ≈ T_cold atol=atol - - set!(model, ρe = ρᵣ * e, qᵗ = qᵗ) - T★ = @allowscalar first(model.temperature) - qᵛm = @allowscalar first(model.microphysical_fields.qᵛ) - qˡm = @allowscalar first(model.microphysical_fields.qˡ) - qⁱm = @allowscalar first(model.microphysical_fields.qⁱ) - - @test T★ ≈ T_cold atol=atol - @test qᵛm ≈ qᵛ⁺ atol=atol - @test qˡm ≈ zero(FT) atol=atol - @test qⁱm ≈ qⁱ atol=atol + # Test 3: Temperatures below homogeneous ice nucleation - all ice + @testset "Temperatures below homogeneous ice nucleation (all ice)" begin + T_cold = FT(220) # Below Tʰ + qᵗ = FT(0.01) + qᵛ⁺ = equilibrium_saturation_specific_humidity(T_cold, pᵣ, qᵗ, constants, equilibrium) + atol = test_tol(FT) + + if qᵗ > qᵛ⁺ # saturated conditions + # All condensate should be ice + qⁱ = qᵗ - qᵛ⁺ + q = MoistureMassFractions(qᵛ⁺, zero(FT), qⁱ) + cᵖᵐ = mixture_heat_capacity(q, constants) + e = cᵖᵐ * T_cold + g * z - ℒⁱᵣ * qⁱ + + 𝒰 = StaticEnergyState(e, q, z, pᵣ) + T★ = compute_temperature(𝒰, microphysics, constants) + @test T★ ≈ T_cold atol=atol + + set!(model, ρe = ρᵣ * e, qᵗ = qᵗ) + T★ = @allowscalar first(model.temperature) + qᵛm = @allowscalar first(model.microphysical_fields.qᵛ) + qˡm = @allowscalar first(model.microphysical_fields.qˡ) + qⁱm = @allowscalar first(model.microphysical_fields.qⁱ) + + @test T★ ≈ T_cold atol=atol + @test qᵛm ≈ qᵛ⁺ atol=atol + @test qˡm ≈ zero(FT) atol=atol + @test qⁱm ≈ qⁱ atol=atol + end end - end - - # Test 4: Mixed-phase range temperatures with moist static energy verification - @testset "Mixed-phase range temperatures with moist static energy" begin - atol = test_tol(FT) - for T in 240:10:270 - @testset let T=T - T = convert(FT, T) - λ = test_liquid_fraction(T, Tᶠ, Tʰ) - qᵗ = FT(0.015) - qᵛ⁺ = adjustment_saturation_specific_humidity(T, pᵣ, qᵗ, constants, equilibrium) - - if qᵗ > qᵛ⁺ # saturated conditions - # Partition condensate between liquid and ice based on λ - q_condensate = qᵗ - qᵛ⁺ - qˡ = λ * q_condensate - qⁱ = (1 - λ) * q_condensate - q = MoistureMassFractions(qᵛ⁺, qˡ, qⁱ) - - # Verify partitioning sums correctly - @test q.vapor + q.liquid + q.ice ≈ qᵗ - - # Compute moist static energy: e = cᵖᵐ*T + g*z - ℒˡᵣ*qˡ - ℒⁱᵣ*qⁱ - cᵖᵐ = mixture_heat_capacity(q, constants) - e = cᵖᵐ * T + g * z - ℒˡᵣ * qˡ - ℒⁱᵣ * qⁱ - # Verify moist static energy can recover temperature - 𝒰 = StaticEnergyState(e, q, z, pᵣ) - T_recovered = (e - g * z + ℒˡᵣ * q.liquid + ℒⁱᵣ * q.ice) / mixture_heat_capacity(q, constants) - @test T_recovered ≈ T - - # Test saturation adjustment recovers temperature - 𝒰_unadjusted = StaticEnergyState(e, MoistureMassFractions(qᵗ), z, pᵣ) - T★ = compute_temperature(𝒰_unadjusted, microphysics, constants) - @test T★ ≈ T atol=atol - - set!(model, ρe = ρᵣ * e, qᵗ = qᵗ) - T★ = @allowscalar first(model.temperature) - qᵛm = @allowscalar first(model.microphysical_fields.qᵛ) - qˡm = @allowscalar first(model.microphysical_fields.qˡ) - qⁱm = @allowscalar first(model.microphysical_fields.qⁱ) - - @test T★ ≈ T atol=atol - @test qᵛm ≈ qᵛ⁺ atol=atol - @test qˡm ≈ qˡ atol=atol - @test qⁱm ≈ qⁱ atol=atol + # Test 4: Mixed-phase range temperatures with moist static energy verification + @testset "Mixed-phase range temperatures with moist static energy" begin + atol = test_tol(FT) + + for T in 240:10:270 + @testset let T=T + T = convert(FT, T) + λ = test_liquid_fraction(T, Tᶠ, Tʰ) + qᵗ = FT(0.015) + qᵛ⁺ = equilibrium_saturation_specific_humidity(T, pᵣ, qᵗ, constants, equilibrium) + + if qᵗ > qᵛ⁺ # saturated conditions + # Partition condensate between liquid and ice based on λ + q_condensate = qᵗ - qᵛ⁺ + qˡ = λ * q_condensate + qⁱ = (1 - λ) * q_condensate + q = MoistureMassFractions(qᵛ⁺, qˡ, qⁱ) + + # Verify partitioning sums correctly + @test q.vapor + q.liquid + q.ice ≈ qᵗ + + # Compute moist static energy: e = cᵖᵐ*T + g*z - ℒˡᵣ*qˡ - ℒⁱᵣ*qⁱ + cᵖᵐ = mixture_heat_capacity(q, constants) + e = cᵖᵐ * T + g * z - ℒˡᵣ * qˡ - ℒⁱᵣ * qⁱ + + # Verify moist static energy can recover temperature + 𝒰 = StaticEnergyState(e, q, z, pᵣ) + T_recovered = (e - g * z + ℒˡᵣ * q.liquid + ℒⁱᵣ * q.ice) / mixture_heat_capacity(q, constants) + @test T_recovered ≈ T + + # Test saturation adjustment recovers temperature + 𝒰_unadjusted = StaticEnergyState(e, MoistureMassFractions(qᵗ), z, pᵣ) + T★ = compute_temperature(𝒰_unadjusted, microphysics, constants) + @test T★ ≈ T atol=atol + + set!(model, ρe = ρᵣ * e, qᵗ = qᵗ) + T★ = @allowscalar first(model.temperature) + qᵛm = @allowscalar first(model.microphysical_fields.qᵛ) + qˡm = @allowscalar first(model.microphysical_fields.qˡ) + qⁱm = @allowscalar first(model.microphysical_fields.qⁱ) + + @test T★ ≈ T atol=atol + @test qᵛm ≈ qᵛ⁺ atol=atol + @test qˡm ≈ qˡ atol=atol + @test qⁱm ≈ qⁱ atol=atol + end end end end @@ -269,7 +277,7 @@ end for qᵗ in FT.(5e-3:5e-3:3e-2) @testset let qᵗ=qᵗ - qᵛ⁺ = adjustment_saturation_specific_humidity(T, pᵣ, qᵗ, constants, equilibrium) + qᵛ⁺ = equilibrium_saturation_specific_humidity(T, pᵣ, qᵗ, constants, equilibrium) if qᵗ > qᵛ⁺ # saturated conditions qᶜ = qᵗ - qᵛ⁺ @@ -299,7 +307,7 @@ end λ_expected = test_liquid_fraction(T_partition, Tᶠ, Tʰ) qᵗ = FT(0.02) - qᵛ⁺ = adjustment_saturation_specific_humidity(T_partition, pᵣ, qᵗ, constants, equilibrium) + qᵛ⁺ = equilibrium_saturation_specific_humidity(T_partition, pᵣ, qᵗ, constants, equilibrium) if qᵗ > qᵛ⁺ # saturated conditions q_condensate = qᵗ - qᵛ⁺ @@ -344,7 +352,7 @@ end # Case 0: Absolute zero potential temperature returns zero temperature θ₀ = zero(FT) q₀ = MoistureMassFractions{FT} |> zero - 𝒰₀ = PotentialTemperatureState(θ₀, q₀, p₀, pᵣ) + 𝒰₀ = LiquidIcePotentialTemperatureState(θ₀, q₀, p₀, pᵣ) T₀ = compute_boussinesq_adjustment_temperature(𝒰₀, constants) @test T₀ == 0 @@ -352,7 +360,7 @@ end θ₁ = FT(300) qᵗ₁ = zero(FT) q₁ = MoistureMassFractions(qᵗ₁) - 𝒰₁ = PotentialTemperatureState(θ₁, q₁, p₀, pᵣ) + 𝒰₁ = LiquidIcePotentialTemperatureState(θ₁, q₁, p₀, pᵣ) Π₁ = exner_function(𝒰₁, constants) T_dry₁ = Π₁ * θ₁ @@ -362,7 +370,7 @@ end # Case 2: Unsaturated, humid but below saturation at dry temperature θ₂ = FT(300) q₂ = MoistureMassFractions{FT} |> zero - 𝒰₂ = PotentialTemperatureState(θ₂, q₂, p₀, pᵣ) + 𝒰₂ = LiquidIcePotentialTemperatureState(θ₂, q₂, p₀, pᵣ) Π₂ = exner_function(𝒰₂, constants) T_dry₂ = Π₂ * θ₂ @@ -382,8 +390,8 @@ end T₃ = θ̃ = FT(300) qᵗ = FT(0.025) q̃ = MoistureMassFractions(qᵗ) - 𝒰 = PotentialTemperatureState(θ̃, q̃, p₀, pᵣ) - qᵛ⁺ = Breeze.MoistAirBuoyancies.adjustment_saturation_specific_humidity(T₃, 𝒰, constants) + 𝒰 = LiquidIcePotentialTemperatureState(θ̃, q̃, p₀, pᵣ) + qᵛ⁺ = Breeze.MoistAirBuoyancies.equilibrium_saturation_specific_humidity(T₃, 𝒰, constants) @test qᵗ > qᵛ⁺ # otherwise the test is wrong qˡ = qᵗ - qᵛ⁺ @@ -393,7 +401,7 @@ end cᵖᵐ = mixture_heat_capacity(q₃, constants) ℒˡᵣ = constants.liquid.reference_latent_heat θ₃ = (T₃ - ℒˡᵣ / cᵖᵐ * qˡ) / Π₃ - 𝒰₃ = PotentialTemperatureState(θ₃, q₃, p₀, pᵣ) + 𝒰₃ = LiquidIcePotentialTemperatureState(θ₃, q₃, p₀, pᵣ) T₃_solve = compute_boussinesq_adjustment_temperature(𝒰₃, constants) @test isapprox(T₃_solve, T₃; atol=atol) diff --git a/test/turbulence_closures.jl b/test/turbulence_closures.jl index 2b38a873..4c4f242b 100644 --- a/test/turbulence_closures.jl +++ b/test/turbulence_closures.jl @@ -2,6 +2,8 @@ using Breeze using Oceananigans using Test +test_thermodynamics = (:StaticEnergy, :LiquidIcePotentialTemperature) + @testset "Time stepping with TurbulenceClosures [$(FT)]" for FT in (Float32, Float64) Oceananigans.defaults.FloatType = FT grid = RectilinearGrid(default_arch; size=(8, 8, 8), x=(0, 100), y=(0, 100), z=(0, 100)) @@ -22,77 +24,84 @@ using Test end end + constants = ThermodynamicConstants() + reference_state = ReferenceState(grid, constants) etd = Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization() discretizations = (vitd, etd) - @testset "Implicit diffusion solver with ScalarDiffusivity [$(FT), $(typeof(disc))]" for disc in discretizations - closure = ScalarDiffusivity(disc, ν=1, κ=1) - model = AtmosphereModel(grid; closure, tracers=:ρc) - # Set uniform specific energy for no diffusion - θ₀ = model.formulation.reference_state.potential_temperature - cᵖᵈ = model.thermodynamic_constants.dry_air.heat_capacity - e₀ = cᵖᵈ * θ₀ - set!(model; e=e₀) - ρe₀ = deepcopy(energy_density(model)) - time_step!(model, 1) - # Use rtol for implicit solver which may have small numerical effects - @test isapprox(energy_density(model), ρe₀, rtol=1e-5) - end - @testset "Closure flux affects momentum tendency [$(FT)]" begin - closure = ScalarDiffusivity(ν=1e4) - model = AtmosphereModel(grid; advection=nothing, closure) - set!(model; ρu = (x, y, z) -> exp((z - 50)^2 / (2 * 20^2))) - Breeze.AtmosphereModels.compute_tendencies!(model) - Gρu = model.timestepper.Gⁿ.ρu - @test maximum(abs, Gρu) > 0 - end + @testset "AtmosphereModel with $thermodynamics thermodynamics [$FT]" for thermodynamics in test_thermodynamics + formulation = AnelasticFormulation(reference_state; thermodynamics) - @testset "SmagorinskyLilly with velocity gradients [$(FT)]" begin - model = AtmosphereModel(grid; closure=SmagorinskyLilly()) - θ₀ = model.formulation.reference_state.potential_temperature - set!(model; θ=θ₀, ρu = (x, y, z) -> z / 100) - Breeze.AtmosphereModels.update_state!(model) - @test maximum(abs, model.closure_fields.νₑ) > 0 - end + @testset "Implicit diffusion solver with ScalarDiffusivity [$thermodynamics, $(FT), $(typeof(disc))]" for disc in discretizations + closure = ScalarDiffusivity(disc, ν=1, κ=1) + model = AtmosphereModel(grid; closure, tracers=:ρc) + # Set uniform specific energy for no diffusion + θ₀ = model.formulation.reference_state.potential_temperature + cᵖᵈ = model.thermodynamic_constants.dry_air.heat_capacity + e₀ = cᵖᵈ * θ₀ + set!(model; e=e₀) + ρe₀ = deepcopy(static_energy_density(model)) + time_step!(model, 1) + # Use rtol for implicit solver which may have small numerical effects + @test isapprox(static_energy_density(model), ρe₀, rtol=1e-5) + end - @testset "AnisotropicMinimumDissipation with velocity gradients [$(FT)]" begin - model = AtmosphereModel(grid; closure=AnisotropicMinimumDissipation()) - set!(model; ρu = (x, y, z) -> z / 100) - Breeze.AtmosphereModels.update_state!(model) - @test haskey(model.closure_fields, :νₑ) || haskey(model.closure_fields, :κₑ) - end + @testset "Closure flux affects momentum tendency [$thermodynamics, $(FT)]" begin + closure = ScalarDiffusivity(ν=1e4) + model = AtmosphereModel(grid; advection=nothing, closure) + set!(model; ρu = (x, y, z) -> exp((z - 50)^2 / (2 * 20^2))) + Breeze.AtmosphereModels.compute_tendencies!(model) + Gρu = model.timestepper.Gⁿ.ρu + @test maximum(abs, Gρu) > 0 + end + + @testset "SmagorinskyLilly with velocity gradients [$thermodynamics, $(FT)]" begin + model = AtmosphereModel(grid; closure=SmagorinskyLilly()) + θ₀ = model.formulation.reference_state.potential_temperature + set!(model; θ=θ₀, ρu = (x, y, z) -> z / 100) + Breeze.AtmosphereModels.update_state!(model) + @test maximum(abs, model.closure_fields.νₑ) > 0 + end - # Test LES scalar diffusion with advection=nothing - # This isolates the effect of the closure on scalar fields - les_closures = (SmagorinskyLilly(), AnisotropicMinimumDissipation()) - - @testset "LES scalar diffusion without advection [$(FT), $(nameof(typeof(closure)))]" for closure in les_closures - model = AtmosphereModel(grid; closure, advection=nothing, tracers=:ρc) - - # Set random velocity field to trigger non-zero eddy diffusivity - Ξ(x, y, z) = randn() - - # Set scalar gradients for energy, moisture, and passive tracer - θ₀ = model.formulation.reference_state.potential_temperature - z₀, dz = 50, 10 - gaussian(z) = exp(- (z - z₀)^2 / 2dz^2) - θᵢ(x, y, z) = θ₀ + 10 * gaussian(z) - qᵗᵢ(x, y, z) = 0.01 + 1e-3 * gaussian(z) - ρcᵢ(x, y, z) = gaussian(z) - set!(model; θ = θᵢ, ρqᵗ = qᵗᵢ, ρc = ρcᵢ, ρu = Ξ, ρv = Ξ, ρw = Ξ) - - # Store initial scalar fields (using copy of data to avoid reference issues) - ρe₀ = copy(interior(energy_density(model))) - ρqᵗ₀ = copy(interior(model.moisture_density)) - ρc₀ = copy(interior(model.tracers.ρc)) - - # Take a time step - time_step!(model, 1) - - # Scalars should change due to diffusion (not advection since advection=nothing) - # Use explicit maximum difference check instead of ≈ to handle Float32 - @test maximum(abs, interior(energy_density(model)) .- ρe₀) > 0 - @test maximum(abs, interior(model.moisture_density) .- ρqᵗ₀) > 0 - @test maximum(abs, interior(model.tracers.ρc) .- ρc₀) > 0 + @testset "AnisotropicMinimumDissipation with velocity gradients [$thermodynamics, $(FT)]" begin + model = AtmosphereModel(grid; closure=AnisotropicMinimumDissipation()) + set!(model; ρu = (x, y, z) -> z / 100) + Breeze.AtmosphereModels.update_state!(model) + @test haskey(model.closure_fields, :νₑ) || haskey(model.closure_fields, :κₑ) + end + + # Test LES scalar diffusion with advection=nothing + # This isolates the effect of the closure on scalar fields + les_closures = (SmagorinskyLilly(), AnisotropicMinimumDissipation()) + + @testset "LES scalar diffusion without advection [$thermodynamics, $(FT), $(nameof(typeof(closure)))]" for closure in les_closures + model = AtmosphereModel(grid; closure, advection=nothing, tracers=:ρc) + + # Set random velocity field to trigger non-zero eddy diffusivity + Ξ(x, y, z) = randn() + + # Set scalar gradients for energy, moisture, and passive tracer + θ₀ = model.formulation.reference_state.potential_temperature + z₀, dz = 50, 10 + gaussian(z) = exp(- (z - z₀)^2 / 2dz^2) + θᵢ(x, y, z) = θ₀ + 10 * gaussian(z) + qᵗᵢ(x, y, z) = 0.01 + 1e-3 * gaussian(z) + ρcᵢ(x, y, z) = gaussian(z) + set!(model; θ = θᵢ, ρqᵗ = qᵗᵢ, ρc = ρcᵢ, ρu = Ξ, ρv = Ξ, ρw = Ξ) + + # Store initial scalar fields (using copy of data to avoid reference issues) + ρe₀ = copy(interior(static_energy_density(model))) + ρqᵗ₀ = copy(interior(model.moisture_density)) + ρc₀ = copy(interior(model.tracers.ρc)) + + # Take a time step + time_step!(model, 1) + + # Scalars should change due to diffusion (not advection since advection=nothing) + # Use explicit maximum difference check instead of ≈ to handle Float32 + @test maximum(abs, interior(static_energy_density(model)) .- ρe₀) > 0 + @test maximum(abs, interior(model.moisture_density) .- ρqᵗ₀) > 0 + @test maximum(abs, interior(model.tracers.ρc) .- ρc₀) > 0 + end end end