More general set! plus MoistureMassFractions abstraction (#72)

* More general set! plus SpecificHumidities abstraction

* try to get SpecificHumidities working

* reduce thermodynamic states

* get SpecificHumidities into docs

* SpecificHumidities -> MassRatios

* humidities -> mass_ratios

* mass ratio -> specific moisture content

* lets try "mass fractions"

* change "latent heat" to "reference latent heat"

* reorganize stuff

* cleanup

* fixing trhe reference state

* bomex running

* thermal bubble runss

* fix thermal bubble

* fix tests

* fix

* use approx

* more test fixes

* rm unused import

* fix docs

* fix

* fix saturation adj bug

* fix saturation adjustment solver for boussinesq model

* fix broken imports

* another import fix

* clean up

* fix docs example

* fix another saturation adj bug

* add comments

* dont show iter

* fix diagnostics

* update saturation adj docs

* use rand not randn

* add tests and change algorithm

* rm show methods

* change ref from chammas to pressel

* fix quick start

* another fix

Author: glwagner

docs/src/breeze.bib

@@ -28,3 +28,14 @@ @book{Straka2009
     year = {2009},
     doi = {10.1017/CBO9780511581168}
 }
+
+@article{Pressel2015,
+  title={Large-eddy simulation in an anelastic framework with closed water and entropy balances},
+  author={Pressel, Kyle G and Kaul, Colleen M and Schneider, Tapio and Tan, Zhihong and Mishra, Siddhartha},
+  journal={Journal of Advances in Modeling Earth Systems},
+  volume={7},
+  number={3},
+  pages={1425--1456},
+  year={2015},
+  publisher={Wiley Online Library}
+}

docs/src/index.md

@@ -39,28 +39,27 @@ seed!(42)
 
 Nx = Nz = 64
 Lz = 4 * 1024
-grid = RectilinearGrid(CPU(), size=(Nx, Nz), x=(0, 2Lz), z=(0, Lz), topology=(Periodic, Flat, Bounded))
+grid = RectilinearGrid(size=(Nx, Nz), x=(0, 2Lz), z=(0, Lz), topology=(Periodic, Flat, Bounded))
 
-reference_constants = Breeze.Thermodynamics.ReferenceStateConstants(base_pressure=1e5, potential_temperature=288)
-buoyancy = Breeze.MoistAirBuoyancy(; reference_constants)
+p₀, θᵣ = 1e5, 288 # reference state parameters
+buoyancy = Breeze.MoistAirBuoyancy(grid, base_pressure=p₀, reference_potential_temperature=θᵣ)
 
-Q₀ = 1000 # heat flux in W / m²
-ρ₀ = Breeze.MoistAirBuoyancies.base_density(buoyancy) # air density at z=0
+thermo = buoyancy.thermodynamics
+ρ₀ = Breeze.Thermodynamics.base_density(p₀, θᵣ, thermo)
 cₚ = buoyancy.thermodynamics.dry_air.heat_capacity
+Q₀ = 1000 # heat flux in W / m²
 θ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(Q₀ / (ρ₀ * cₚ)))
-q_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1e-2))
+qᵗ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1e-2))
 
 advection = WENO()
-tracers = (:θ, :q)
 model = NonhydrostaticModel(; grid, advection, buoyancy,
-                            tracers = (:θ, :q),
-                            boundary_conditions = (θ=θ_bcs, q=q_bcs))
+                            tracers = (:θ, :qᵗ),
+                            boundary_conditions = (θ=θ_bcs, qᵗ=qᵗ_bcs))
 
 Δθ = 2 # ᵒK
-Tₛ = reference_constants.reference_potential_temperature # K
+Tₛ = buoyancy.reference_state.potential_temperature # K
 θᵢ(x, z) = Tₛ + Δθ * z / grid.Lz + 2e-2 * Δθ * (rand() - 0.5)
-qᵢ(x, z) = 0 # 1e-2 + 1e-5 * rand()
-set!(model, θ=θᵢ, q=qᵢ)
+set!(model, θ=θᵢ)
 
 simulation = Simulation(model, Δt=10, stop_time=2hours)
 conjure_time_step_wizard!(simulation, cfl=0.7)

docs/src/microphysics/saturation_adjustment.md

@@ -1,11 +1,11 @@
 # Warm-phase saturation adjustment
 
 Warm-phase saturation adjustment is a model for water droplet nucleation that assumes that water vapor in excess of the saturation specific humidity is instantaneously converted to liquid water.
-Mixed-phase saturation adjustment is described by [Chammas2023](@citet).
+Mixed-phase saturation adjustment is described by [Pressel2015](@citet).
 Saturation adjustment may be formulated as a nonlinear algebraic equation that relates temperature, potential temperature, and total specific humidity, derived from the definition of liquid potential temperature,
 
 ```math
-θ = \frac{T}{Π} \left \{1 - \frac{ℒᵥ₀}{cᵖᵐ T} \max \left [0, qᵗ - qᵛ⁺(T) \right ] \right \} ,
+θ = \frac{1}{Π} \left \{T - \frac{ℒᵥ₀}{cᵖᵐ} \max \left [0, qᵗ - qᵛ⁺(T) \right ] \right \} ,
 ```
 
 where ``Π`` is the Exner function, ``θ`` is potential temperature, ``T`` is temperature,
@@ -32,14 +32,15 @@ The saturation specific humidity is then
 
 ```@example microphysics
 using Breeze
-using Breeze.MoistAirBuoyancies: saturation_specific_humidity, HeightReferenceThermodynamicState
+using Breeze.Thermodynamics: saturation_specific_humidity
 
 thermo = ThermodynamicConstants()
-ref = ReferenceStateConstants(base_pressure=101325, potential_temperature=288)
 
-z = 0.0    # [m] height
-θ = 290.0  # [ᵒK] potential temperature
-qᵛ⁺₀ = saturation_specific_humidity(θ, z, ref, thermo, thermo.liquid)
+p₀ = 101325.0
+θ₀ = 288.0
+Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo)
+ρ₀ = p₀ / (Rᵈ * θ₀)
+qᵛ⁺₀ = saturation_specific_humidity(θ₀, ρ₀, thermo, thermo.liquid)
 ```
 
 Recall that the specific humidity is unitless, or has units "``kg / kg``": kg of water vapor
@@ -49,17 +50,20 @@ given a total specific humidity slightly above saturation:
 
 ```@example microphysics
 using Breeze.MoistAirBuoyancies: temperature
+using Breeze.Thermodynamics: PotentialTemperatureState, MoistureMassFractions
 
+z = 0.0
 qᵗ = 0.012   # [kg kg⁻¹] total specific humidity
-U = HeightReferenceThermodynamicState(θ, qᵗ, z)
-T = temperature(U, ref, thermo)
+q = MoistureMassFractions(qᵗ, zero(qᵗ), zero(qᵗ))
+U = PotentialTemperatureState(θ₀, q, z, p₀, p₀, ρ₀)
+T = temperature(U, thermo)
 ```
 
 Finally, we recover the amount of liquid condensate by subtracting the saturation
 specific humidity from the total:
 
 ```@example microphysics
-qᵛ⁺ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid)
+qᵛ⁺ = saturation_specific_humidity(T, ρ₀, thermo, thermo.liquid)
 qˡ = qᵗ - qᵛ⁺
 ```
 
@@ -69,14 +73,15 @@ As a second example, we examine the dependence of temperature on total specific
 when the potential temperature is constant:
 
 ```@example microphysics
-qᵗ = 0:1e-4:0.04 # [kg kg⁻¹] total specific humidity
-U = [HeightReferenceThermodynamicState(θ, qᵗⁱ, z) for qᵗⁱ in qᵗ]
-T = [temperature(Uⁱ, ref, thermo) for Uⁱ in U]
+qᵗ = 0:1e-4:0.035 # [kg kg⁻¹] total specific humidity
+q = [MoistureMassFractions(qᵗⁱ, 0.0, 0.0) for qᵗⁱ in qᵗ]
+U = [PotentialTemperatureState(θ₀, qⁱ, z, p₀, p₀, ρ₀) for qⁱ in q]
+T = [temperature(Uⁱ, thermo) for Uⁱ in U]
 
 ## Compare with a simple piecewise linear model
-ℒᵥ₀ = thermo.liquid.latent_heat
+ℒᵥ₀ = thermo.liquid.reference_latent_heat
 cᵖᵈ = thermo.dry_air.heat_capacity
-T̃ = [290 + ℒᵥ₀ / cᵖᵈ * max(0, qᵗⁱ - qᵛ⁺₀) for qᵗⁱ in qᵗ]
+T̃ = [288 + ℒᵥ₀ / cᵖᵈ * max(0, qᵗⁱ - qᵛ⁺₀) for qᵗⁱ in qᵗ]
 
 using CairoMakie
 
@@ -94,13 +99,32 @@ For a third example, we consider a state with constant potential temperature and
 but at varying heights:
 
 ```@example microphysics
+grid = RectilinearGrid(size=100, z=(0, 1e4), topology=(Flat, Flat, Bounded))
+thermo = ThermodynamicConstants()
+reference_state = ReferenceState(grid, thermo)
+
+θᵣ = reference_state.potential_temperature
+p₀ = reference_state.base_pressure
 qᵗ = 0.005
-z = 0:100:10e3
+q = MoistureMassFractions(qᵗ, 0.0, 0.0)
+
+z = znodes(grid, Center())
+T = zeros(grid.Nz)
+qᵛ⁺ = zeros(grid.Nz)
+qˡ = zeros(grid.Nz)
+rh = zeros(grid.Nz)
+
+for k = 1:grid.Nz
+    ρᵣ = reference_state.density[1, 1, k]
+    pᵣ = reference_state.pressure[1, 1, k]
+
+    U = PotentialTemperatureState(θᵣ, q, z[k], p₀, pᵣ, ρᵣ)
+    T[k] = temperature(U, thermo)
 
-T = [temperature(HeightReferenceThermodynamicState(θ, qᵗ, zᵏ), ref, thermo) for zᵏ in z]
-qᵛ⁺ = [saturation_specific_humidity(T[k], z[k], ref, thermo, thermo.liquid) for k = 1:length(z)]
-qˡ = [max(0, qᵗ - qᵛ⁺ᵏ) for qᵛ⁺ᵏ in qᵛ⁺]
-rh = [100 * min(qᵗ, qᵛ⁺ᵏ) / qᵛ⁺ᵏ for qᵛ⁺ᵏ in qᵛ⁺]
+    qᵛ⁺[k] = saturation_specific_humidity(T[k], ρᵣ, thermo, thermo.liquid)
+    qˡ[k] = max(0, qᵗ - qᵛ⁺[k])
+    rh[k] = 100 * min(qᵗ, qᵛ⁺[k]) / qᵛ⁺[k]
+end
 
 cᵖᵈ = thermo.dry_air.heat_capacity
 g = thermo.gravitational_acceleration

docs/src/thermodynamics.md

@@ -187,23 +187,19 @@ with ``Rᵈ = 286.71 \; \mathrm{J} \, \mathrm{K}^{-1}``:
 
 ```@example reference_state
 using Breeze
-using Breeze.Thermodynamics: reference_pressure, reference_density
 using CairoMakie
 
-thermo = ThermodynamicConstants()
-constants = ReferenceStateConstants(base_pressure=101325, potential_temperature=288)
 grid = RectilinearGrid(size=160, z=(0, 12_000), topology=(Flat, Flat, Bounded))
+thermo = ThermodynamicConstants()
+reference_state = ReferenceState(grid, thermo, base_pressure=101325, potential_temperature=288)
 
-pᵣ = CenterField(grid)
-ρᵣ = CenterField(grid)
-
-set!(pᵣ, z -> reference_pressure(z, constants, thermo))
-set!(ρᵣ, z -> reference_density(z, constants, thermo))
+pᵣ = reference_state.pressure
+ρᵣ = reference_state.density
 
 Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo)
 cᵖᵈ = thermo.dry_air.heat_capacity
-p₀ = constants.base_pressure
-θ₀ = constants.reference_potential_temperature
+p₀ = reference_state.base_pressure
+θ₀ = reference_state.potential_temperature
 g = thermo.gravitational_acceleration
 
 # Verify that Tᵣ = θ₀ (1 - g z / (cᵖᵈ θ₀))
@@ -280,17 +276,18 @@ More generally, ``qᵈ = 1 - qᵛ - qᶜ``, where ``qᶜ`` is the total mass
 ratio of condensed species. In most situations on Earth, ``qᶜ ≪ qᵛ``.
 
 ```@example thermo
+using Breeze.Thermodynamics: MoistureMassFractions
+
 # Compute mixture properties for air with 0.01 specific humidity
-qᵛ = 0.01 # 1% water vapor by mass
-Rᵐ = mixture_gas_constant(qᵛ, thermo)
+qᵗ = 0.01 # 1% water vapor by mass
+q = MoistureMassFractions(qᵗ, zero(qᵗ), zero(qᵗ))
+Rᵐ = mixture_gas_constant(q, thermo)
 ```
 
 We likewise define a mixture heat capacity via ``cᵖᵐ = qᵈ cᵖᵈ + qᵛ cᵖᵛ``,
 
-
 ```@example thermo
-q = 0.01 # 1% water vapor by mass
-cᵖᵐ = mixture_heat_capacity(qᵛ, thermo)
+cᵖᵐ = mixture_heat_capacity(q, thermo)
 ```
 
 ## Liquid-ice potential temperature
@@ -313,7 +310,7 @@ heats, the latent heat of a phase transition is linear in temperature.
 For example, for phase change from vapor to liquid,
 
 ```math
-ℒˡ(T) = ℒˡ(T=0) + \big ( \underbrace{cᵖᵛ - cᵖˡ}_{≡Δcˡ} \big ) T ,
+ℒˡ(T) = ℒˡ(T=0) + \big ( \underbrace{cᵖᵛ - cˡ}_{≡Δcˡ} \big ) T ,
 ```
 
 where ``ℒˡ(T=0)`` is the latent heat at absolute zero, ``T = 0 \; \mathrm{K}``.
@@ -327,13 +324,13 @@ Consider parameters for liquid water,
 
 ```@example thermo
 using Breeze.Thermodynamics: CondensedPhase
-liquid_water = CondensedPhase(latent_heat=2500800, heat_capacity=4181)
+liquid_water = CondensedPhase(reference_latent_heat=2500800, heat_capacity=4181)
 ```
 
 or water ice,
 
 ```@example thermo
-water_ice = CondensedPhase(latent_heat=2834000, heat_capacity=2108)
+water_ice = CondensedPhase(reference_latent_heat=2834000, heat_capacity=2108)
 ```
 
 The saturation vapor pressure is
@@ -369,14 +366,19 @@ and this is what it looks like:
 
 ```@example
 using Breeze
-using Breeze.MoistAirBuoyancies: saturation_specific_humidity
+using Breeze.Thermodynamics: saturation_specific_humidity
 
 thermo = ThermodynamicConstants()
-ref = ReferenceStateConstants(base_pressure=101325, potential_temperature=288)
 
-z = 0
+p₀ = 101325
+Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo)
 T = collect(273.2:0.1:313.2)
-qᵛ⁺ = [saturation_specific_humidity(Tⁱ, z, ref, thermo, thermo.liquid) for Tⁱ in T]
+qᵛ⁺ = zeros(length(T))
+
+for i = 1:length(T)
+    ρ = p₀ / (Rᵈ * T[i])
+    qᵛ⁺[i] = saturation_specific_humidity(T[i], ρ, thermo, thermo.liquid)
+end
 
 using CairoMakie
 

examples/JRA55_saturation_specific_humidity.jl

@@ -24,7 +24,7 @@ ax4 = Axis(fig[2, 2], title="Liquid specific humidity")
 # Compute cloudiness for instantaneous drop
 θ⁻ = T .- 10
 Ψ = Breeze.ThermodynamicState{Float64}.(θ⁻, q, 0)
-ℛ = Breeze.ReferenceStateConstants{Float64}(101325, 20)
+ℛ = Breeze.ReferenceState{Float64}(101325, 20)
 T⁻ = Breeze.temperature.(Ψ, Ref(ℛ), Ref(thermo))
 qᵛ★ = saturation_specific_humidity.(T⁻, 1.2, Ref(thermo))
 qˡ = @. max(0, q - qᵛ★)

examples/atmos_convection.jl

@@ -11,7 +11,7 @@ grid = RectilinearGrid(size=(Nx, Nz), x=(0, Lx), z=(0, Lz), topology=(Periodic,
 ρe_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1000))
 model = AtmosphereModel(grid, advection=WENO(), boundary_conditions=(; ρe=ρe_bcs))
 
-Tₛ = model.formulation.constants.reference_potential_temperature
+Tₛ = model.formulation.reference_state.potential_temperature
 g = model.thermodynamics.gravitational_acceleration
 N² = 1e-6
 dθdz = N² * Tₛ / g  # Background potential temperature gradient