[Dycore] Implement tendency term for moist static energy following Pauluis (2008)

docs/make.jl

@@ -20,6 +20,10 @@ makedocs(sitename="Breeze",
             "Overview" => "microphysics/microphysics_overview.md",
             "Warm phase saturation adjustment" => "microphysics/saturation_adjustment.md",
         ],
+        "Dycore equations and algorithms" => "dynamics.md",
+        "Appendix" => Any[
+            "Notation" => "appendix/notation.md",
+        ],
         "References" => "references.md",
         "API" => "api.md",
     ]

docs/src/appendix/notation.md

@@ -0,0 +1,48 @@
+# Notation
+
+This appendix establishes a common notation across the documentation and source.
+Each entry lists a mathematical symbol, the Unicode form commonly used in
+the codebase, and a brief description. Mathematical symbols are shown with
+inline math, while the Unicode column shows the exact glyphs used in code.
+
+| mathematical symbol | unicode | description |
+| --- | --- | --- |
+| ``p`` | `p` | Pressure |
+| ``p_r(z)`` | `pᵣ` | Hydrostatic reference pressure at height ``z`` |
+| ``p_0`` | `p₀` | Base (surface) reference pressure |
+| ``p_h'`` | `pₕ′` | Hydrostatic pressure anomaly, ``∂_z p_h' = - ρ_r b`` |
+| ``p_n`` | `pₙ` | Nonhydrostatic pressure (projection/correction potential) |
+| ``\rho`` | `ρ` | Density |
+| ``\rho_r(z)`` | `ρᵣ` | Reference density at height ``z`` |
+| ``\boldsymbol{u} = (u,v,w)`` | `u, v, w` | Velocity components (x, y, z) |
+| ``\rho_r\,u,\; \rho_r\,v,\; \rho_r\,w`` | `ρu, ρv, ρw` | Momentum components (ρᵣ-weighted velocities) |
+| ``g`` | `g` | Gravitational acceleration |
+| ``T`` | `T` | Temperature |
+| ``\theta`` | `θ` | Potential temperature |
+| ``\theta_r`` | `θᵣ` | Reference potential temperature (constant) |
+| ``\alpha`` | `α` | Specific volume, ``α = 1/ρ`` or ``α = R^{m} T / p_r`` in anelastic |
+| ``\alpha_{r}`` | `αᵣ` | Reference specific volume, ``α_{r} = R^{d} θ_r / p_r`` |
+| ``\Pi`` | `Π` | Exner function, ``\Pi = (p_r/p_0)^{R^{m}/c^{pm}}`` |
+| ``R^{d}`` | `Rᵈ` | Dry air gas constant |
+| ``R^{v}`` | `Rᵛ` | Water vapor gas constant |
+| ``R^{m}`` | `Rᵐ` | Mixture gas constant, function of ``q`` |
+| ``c^{pd}`` | `cᵖᵈ` | Heat capacity of dry air at constant pressure |
+| ``c^{pv}`` | `cᵖᵛ` | Heat capacity of vapor at constant pressure |
+| ``c^{pm}`` | `cᵖᵐ` | Mixture heat capacity at constant pressure |
+| ``q`` | `q` | Specific humidity (vapor mass fraction) |
+| ``q^{t}`` | `qᵗ` | Total specific humidity (vapor + condensates) |
+| ``q^{v}`` | `qᵛ` | Water vapor specific humidity |
+| ``q^{\ell}`` | `qˡ` | Liquid water specific humidity |
+| ``q^{i}`` | `qi` | Ice specific humidity |
+| ``q^{v+}`` | `qᵛ⁺` | Saturation specific humidity over liquid/ice (context-dependent) |
+| ``\mathcal{L}^{v}`` | `ℒᵛ` | Latent heat of vaporization |
+| ``\mathcal{R}`` | `ℛ` | Universal (molar) gas constant |
+| ``b`` | `b` | Buoyancy, ``b = g (α - α_{r}) / α_{r}`` |
+| ``\Delta t`` | `Δt` | Time step |
+| ``\Delta z`` | `Δz` | Vertical grid spacing |
+
+
+Notes:
+- Reference-state quantities use a subscript ``r`` (e.g., ``p_r``, ``\rho_r``), following the Thermodynamics docs and code.
+- Phase or mixture identifiers (``d``, ``v``, ``m``) appear as superscripts (e.g., ``R^{d}``, ``c^{pm}``), matching usage in the codebase (e.g., `Rᵈ`, `cᵖᵐ`).
+- Conservative variables are typically stored in ρᵣ-weighted form in the code (e.g., `ρu`, `ρv`, `ρw`, `ρe`, `ρqᵗ`).

docs/src/dynamics.md

@@ -0,0 +1,128 @@
+# Dynamics
+
+This page summarizes the governing equations behind Breeze’s atmospheric dynamics and the anelastic formulation used by `AtmosphereModel`, following the thermodynamically consistent framework of Pauluis (2008) [Pauluis2008](@citet).
+
+We begin with the compressible Navier–Stokes momentum equations and reduce them to an anelastic, conservative form. We then introduce the moist static energy equation and outline the time-discretized pressure correction used to enforce the anelastic constraint.
+
+## Compressible momentum equations
+
+Let ``ρ`` denote density, ``\boldsymbol{u}`` velocity, ``p`` pressure, ``\boldsymbol{f}`` non-pressure body forces (e.g., Coriolis), and ``\boldsymbol{\tau}`` subgrid/viscous stresses. With gravity ``g \hat{\boldsymbol{z}}``, the inviscid compressible equations in flux form are
+
+```math
+\begin{aligned}
+&\text{Mass:} && \partial_t \rho + \nabla\!\cdot(\rho\,\boldsymbol{u}) = 0. \\
+&\text{Momentum:} && \partial_t(\rho\,\boldsymbol{u}) + \nabla\!\cdot(\rho\,\boldsymbol{u}\,\boldsymbol{u}) + \nabla p 
+\;=\; \rho\,g\,\hat{\boldsymbol{z}} + \rho\,\boldsymbol{f} + \nabla\!\cdot\boldsymbol{\tau}. 
+\end{aligned}
+```
+
+For moist flows we also track total water (vapor + condensates) via
+
+```math
+\partial_t(\rho q^t) + \nabla\!\cdot(\rho q^t\,\boldsymbol{u}) = S_q ,
+```
+
+where ``q^t`` is total specific humidity and ``S_q`` accounts for sources/sinks from microphysics and boundary fluxes.
+
+Thermodynamic relations (mixture gas constant ``R^{m}``, heat capacity ``c^{pm}``, Exner function, etc.) are summarized on the Thermodynamics page.
+
+## Anelastic approximation
+
+To filter acoustic waves while retaining compressibility effects in buoyancy and thermodynamics, we linearize about a hydrostatic, horizontally uniform reference state ``(p_r(z), \rho_r(z))`` with constant reference potential temperature ``\theta_r``. The key assumptions are
+
+- Small Mach number and small relative density perturbations except in buoyancy.
+- Hydrostatic reference balance: ``\partial_z p_r = -\rho_r g``.
+- Mass flux divergence constraint: ``\nabla\!\cdot(\rho_r\,\boldsymbol{u}) = 0``.
+
+Define the specific volume of moist air and its reference value as
+
+```math
+\alpha = \frac{R^{m} T}{p_r} , \qquad \alpha_{r} = \frac{R^{d}\,\theta_r}{p_r} ,
+```
+
+where ``R^{m}`` is the mixture gas constant and ``R^{d}`` is the dry-air gas constant. The buoyancy appearing in the vertical momentum is
+
+```math
+b \equiv g\,\frac{\alpha - \alpha_{r}}{\alpha_{r}} .
+```
+
+## Conservative anelastic system
+
+With ``\rho_r(z)`` fixed by the reference state, the prognostic equations advanced in Breeze are written in conservative form for the ``\rho_r``-weighted fields:
+
+- Continuity (constraint):
+
+```math
+\nabla\!\cdot(\rho_r\,\boldsymbol{u}) = 0 .
+```
+
+- Momentum:
+
+```math
+\partial_t(\rho_r\,\boldsymbol{u}) + \nabla\!\cdot(\rho_r\,\boldsymbol{u}\,\boldsymbol{u}) 
+\;=\; -\,\rho_r\,\nabla \phi + \rho_r\,b\,\hat{\boldsymbol{z}} + \rho_r\,\boldsymbol{f} + \nabla\!\cdot\boldsymbol{\tau} ,
+```
+
+where ``\phi`` is a nonhydrostatic pressure correction potential defined by the projection step (see below). Pressure is decomposed as ``p = p_r(z) + p_h'(x,y,z,t) + p_n``, where ``p_h'`` is a hydrostatic anomaly (obeying ``\partial_z p_h' = -\rho_r b``) and ``p_n`` is the nonhydrostatic component responsible for enforcing the anelastic constraint. In the discrete formulation used here, ``\phi`` coincides with the pressure correction variable.
+
+- Total water:
+
+```math
+\partial_t(\rho_r q^t) + \nabla\!\cdot(\rho_r q^t\,\boldsymbol{u}) = S_q .
+```
+
+### Moist static energy (Pauluis, 2008)
+
+Following [Pauluis2008](@citet), Breeze advances a conservative moist static energy density
+
+```math
+\rho_r e \equiv \rho_r\,c^{pd}\,\theta ,
+```
+
+where ``c^{p^{d}}`` is the dry-air heat capacity at constant pressure and ``\theta`` is the (moist) potential temperature. The prognostic equation reads
+
+```math
+\partial_t(\rho_r e) + \nabla\!\cdot(\rho_r e\,\boldsymbol{u}) 
+\;=\; \rho_r\,w\,b \; + \; S_e ,
+```
+
+with vertical velocity ``w``, buoyancy ``b`` as above, and ``S_e`` including microphysical and external energy sources/sinks. The ``\rho_r w b`` term is the buoyancy flux that links the energy and momentum budgets in the anelastic limit.
+
+Thermodynamic closures needed for ``R^{m}``, ``c^{pm}`` and the Exner function ``\Pi = (p_r/p_0)^{R^{m}/c^{pm}}`` are given in Thermodynamics.
+
+## Time discretization and pressure correction
+
+Breeze uses an explicit multi-stage time integrator for advection, Coriolis, buoyancy, forcing, and tracer terms, coupled with a projection step to enforce the anelastic constraint at each substep. Denote the predicted momentum by ``\widetilde{(\rho_r\,\boldsymbol{u})}``. The projection is
+
+1. Solve the variable-coefficient Poisson problem for the pressure correction potential ``\phi``:
+
+```math
+\nabla\!\cdot\big( \rho_r \, \nabla \phi \big) 
+\;=\; \frac{1}{\Delta t} \, \nabla\!\cdot\,\widetilde{(\rho_r\,\boldsymbol{u})} ,
+```
+
+   with periodic lateral boundaries and homogeneous Neumann boundary conditions in ``z``.
+
+2. Update momentum to enforce ``\nabla\!\cdot(\rho_r\,\boldsymbol{u}^{\,n+1})=0``:
+
+```math
+\rho_r\,\boldsymbol{u}^{\,n+1} \;=\; \widetilde{(\rho_r\,\boldsymbol{u})} \; - \; \Delta t\, \rho_r \, \nabla \phi .
+```
+
+In Breeze this projection is implemented as a Fourier–tridiagonal solve in the vertical with variable ``\rho_r(z)``, aligning with the hydrostatic reference state. The hydrostatic pressure anomaly ``p_h'`` can be obtained diagnostically by vertical integration of buoyancy and used when desired to separate hydrostatic and nonhydrostatic pressure effects.
+
+## Symbols and closures used here
+
+- ``\rho_r(z)``, ``p_r(z)``: Reference density and pressure satisfying hydrostatic balance for a constant ``\theta_r``.
+- ``\alpha = R^{m} T/p_r``, ``\alpha_{r} = R^{d}\,\theta_r/p_r``: Specific volume and its reference value.
+- ``b = g (\alpha - \alpha_{r})/\alpha_{r}``: Buoyancy.
+- ``e = c^{pd}\,\theta``: Energy variable used for moist static energy in the conservative equation.
+- ``q^t``: Total specific humidity (vapor + condensates).
+- ``\phi``: Nonhydrostatic pressure correction potential used by the projection.
+
+See Thermodynamics for definitions of ``R^{m}(q)``, ``c^{pm}(q)``, and ``\Pi``.
+
+## References
+
+```@bibliography
+```

docs/src/index.md

@@ -26,7 +26,7 @@ Pkg.add("https://github.com/NumericalEarth/Breeze.jl.git")
 
 A basic free convection simulation:
 
-```@example intro
+```julia
 using Oceananigans
 using Oceananigans.Units
 using CairoMakie

examples/atmos_convection.jl

@@ -0,0 +1,75 @@
+using Breeze
+using Oceananigans.Units
+using Printf
+using GLMakie
+
+Nx = Nz = 64
+Δx, Δz = 100, 50
+Lx, Lz = Δx * Nx, Δz * Nz
+grid = RectilinearGrid(size=(Nx, Nz), x=(0, Lx), z=(0, Lz), topology=(Periodic, Flat, Bounded))
+
+ρe_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1000))
+model = AtmosphereModel(grid, advection=WENO(), boundary_conditions=(; ρe=ρe_bcs))
+
+Tₛ = model.formulation.constants.reference_potential_temperature
+g = model.thermodynamics.gravitational_acceleration
+N² = 1e-6
+dθdz = N² * Tₛ / g  # Background potential temperature gradient
+θᵢ(x, z) = Tₛ + dθdz * z # background stratification
+Ξᵢ(x, z) = 1e-2 * randn()
+set!(model, θ=θᵢ, u=Ξᵢ, v=Ξᵢ)
+
+simulation = Simulation(model, Δt=10, stop_time=20minutes)
+conjure_time_step_wizard!(simulation, cfl=0.7)
+∫ρe = Field(Integral(model.energy))
+
+function progress(sim)
+    compute!(∫ρe)
+    max_w = maximum(abs, model.velocities.w)
+
+    @info @sprintf("%s: t = %s, ∫ρe dV = %.9e J/kg, max|w| = %.2e m/s",
+                   iteration(sim), prettytime(sim), first(∫ρe), max_w)
+
+    return nothing
+end
+
+add_callback!(simulation, progress, IterationInterval(5))
+
+∫ρe_series = Float64[]
+w_series = []
+t = []
+
+function record_data!(sim)
+    compute!(∫ρe)
+    push!(∫ρe_series, first(∫ρe))
+    push!(w_series, deepcopy(model.velocities.w))
+    push!(t, sim.model.clock.time)
+    return nothing
+end
+add_callback!(simulation, record_data!)
+
+run!(simulation)
+
+fig = Figure(size=(800, 550), fontsize=12)
+
+Nt = length(∫ρe_series)
+axw = Axis(fig[1, 1], xlabel="x (m)", ylabel="z (m)", aspect=DataAspect())
+axe = Axis(fig[2, 1], xlabel="Time (s)", ylabel="∫ρe dV / ∫ρe dV |t=0", height=Relative(0.2))
+slider = Slider(fig[3, 1], range=1:Nt, startvalue=1)
+n = slider.value
+
+∫ρe_n = @lift ∫ρe_series[$n] / ∫ρe_series[1]
+t_n = @lift t[$n]
+w_n = @lift w_series[$n]
+w_lim = maximum(abs, model.velocities.w)
+
+lines!(axe, t, ∫ρe_series / ∫ρe_series[1])
+scatter!(axe, t_n, ∫ρe_n, markersize=10)
+
+hm = heatmap!(axw, w_n, colorrange=(-w_lim, w_lim), colormap=:balance)
+Colorbar(fig[1, 2], hm, label="w (m/s)")
+
+rowgap!(fig.layout, 1, Relative(-0.2))
+rowgap!(fig.layout, 2, Relative(-0.2))
+
+display(fig)

examples/mountain_wave.jl

@@ -0,0 +1,73 @@
+using Breeze
+using Oceananigans.Grids: ImmersedBoundaryGrid, PartialCellBottom
+using Oceananigans.Units
+using Printf
+
+Nx, Nz = 512, 512
+H, L = 20kilometers, 200kilometers
+
+underlying_grid = RectilinearGrid(size = (Nx, Nz), halo = (4, 4),
+                                  x = (-L, L), z = (0, H),
+                                  topology = (Periodic, Flat, Bounded))
+
+h₀ = 250meters
+a = 5kilometers
+λ = 4kilometers
+hill(x) = h₀ * exp(-(x / a)^2) * cos(π * x / λ)^2
+grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(hill))
+
+model = AtmosphereModel(grid, advection = WENO())
+
+# Initial conditions
+θ₀ = model.formulation.constants.reference_potential_temperature
+g = model.thermodynamics.gravitational_acceleration
+N² = 1e-4           # Brunt-Väisälä frequency squared (s⁻²)
+dθdz = N² * θ₀ / g  # Background potential temperature gradient
+θᵢ(x, z) = θ₀ + dθdz * z # background stratification
+Uᵢ = 10
+set!(model, θ=θᵢ, u=Uᵢ)
+
+Δt = 1 # seconds
+stop_iteration = 1000
+simulation = Simulation(model; Δt, stop_iteration)
+conjure_time_step_wizard!(simulation, cfl=0.7)
+
+wall_clock = Ref(time_ns())
+
+ρu, ρv, ρw = model.momentum
+δ = Field(∂x(ρu) + ∂y(ρv) + ∂z(ρw))
+
+function progress(sim)
+    compute!(δ)
+    elapsed = 1e-9 * (time_ns() - wall_clock[])
+
+    msg = @sprintf("Iter: %d, time: %s, wall time: %s, max|w|: %6.3e, m s⁻¹, max|δ|: %6.3e s⁻¹\n",
+                   iteration(sim), prettytime(sim), prettytime(elapsed),
+                   maximum(abs, sim.model.velocities.w), maximum(abs, δ))
+
+    wall_clock[] = time_ns()
+
+    @info msg
+
+    return nothing
+end
+
+add_callback!(simulation, progress, name=:progress, IterationInterval(200))
+
+filename = "mountain_waves"
+outputs = merge(model.velocities, (; δ))
+simulation.output_writers[:fields] = JLD2Writer(model, outputs; filename,
+                                                schedule = IterationInterval(10),
+                                                overwrite_existing = true)
+
+run!(simulation)
+
+using GLMakie
+
+fig = Figure()
+axw = Axis(fig[1, 1], title="Vertical Velocity w (m/s)")
+axδ = Axis(fig[1, 2], title="Divergence δ (s⁻¹)")
+
+heatmap!(axw, model.velocities.w)
+heatmap!(axδ, δ)
+fig