diff --git a/docs/make.jl b/docs/make.jl index 26704464..6c49eb4f 100644 --- a/docs/make.jl +++ b/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", ] diff --git a/docs/src/appendix/notation.md b/docs/src/appendix/notation.md new file mode 100644 index 00000000..eee29707 --- /dev/null +++ b/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ᵗ`). diff --git a/docs/src/dynamics.md b/docs/src/dynamics.md new file mode 100644 index 00000000..9059fc99 --- /dev/null +++ b/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 +``` diff --git a/docs/src/index.md b/docs/src/index.md index 4bf6c8ba..39fcc78e 100644 --- a/docs/src/index.md +++ b/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 diff --git a/examples/atmos_convection.jl b/examples/atmos_convection.jl new file mode 100644 index 00000000..2bafd187 --- /dev/null +++ b/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) diff --git a/examples/mountain_wave.jl b/examples/mountain_wave.jl new file mode 100644 index 00000000..20394d56 --- /dev/null +++ b/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 diff --git a/examples/test/runtests.jl b/examples/test/runtests.jl deleted file mode 100644 index 59918ade..00000000 --- a/examples/test/runtests.jl +++ /dev/null @@ -1,3 +0,0 @@ -using Test - -@test 1==1 \ No newline at end of file diff --git a/examples/thermal_bubble.jl b/examples/thermal_bubble.jl index 595f50be..4b94c176 100644 --- a/examples/thermal_bubble.jl +++ b/examples/thermal_bubble.jl @@ -35,8 +35,7 @@ buoyancy = MoistAirBuoyancy(; reference_constants) advection = WENO(order=5) # Create the model -#model = NonhydrostaticModel(; grid, advection, buoyancy, tracers = (:θ, :q)) -model = AtmosphereModel(grid; advection) #, buoyancy, tracers = (:θ, :q)) +model = AtmosphereModel(grid; advection) # Thermal bubble parameters x₀ = Lx / 2 # Center of bubble in x @@ -45,7 +44,7 @@ r₀ = 2e3 # Initial radius of bubble (2 km) Δθ = 2 # Potential temperature perturbation (ᵒK) # Background stratification -N² = 1e-6 # Brunt-Väisälä frequency squared (s⁻²) +N² = 1e-8 # Brunt-Väisälä frequency squared (s⁻²) dθdz = N² * θ₀ / 9.81 # Background potential temperature gradient # Initial conditions @@ -67,16 +66,16 @@ conjure_time_step_wizard!(simulation, cfl=0.7) # Progress monitoring function progress(sim) - θ = sim.model.tracers.θ + ρe = sim.model.energy u, w = sim.model.velocities - θ_max = maximum(θ) - θ_min = minimum(θ) + ρe_max = maximum(ρe) + ρe_min = minimum(ρe) u_max = maximum(abs, u) w_max = maximum(abs, w) - msg = @sprintf("Iter: %d, t: %s, Δt: %s, extrema(θ): (%.2f, %.2f) K, max|u|: %.2f m/s, max|w|: %.2f m/s", - iteration(sim), prettytime(sim), prettytime(sim.Δt), θ_min, θ_max, u_max, w_max) + msg = @sprintf("Iter: %d, t: %s, Δt: %s, extrema(ρe): (%.2f, %.2f) J/kg, max|u|: %.2f m/s, max|w|: %.2f m/s", + iteration(sim), prettytime(sim), prettytime(sim.Δt), ρe_min, ρe_max, u_max, w_max) @info msg return nothing @@ -92,12 +91,13 @@ u, v, w = model.velocities ζ = ∂x(w) - ∂z(u) # Temperature perturbation from background state -θ_bg_field = Field{Nothing, Nothing, Center}(grid) -set!(θ_bg_field, z -> θ₀ + dθdz * z) -θ′ = model.tracers.θ - θ_bg_field +ρE = Field{Nothing, Nothing, Center}(grid) +set!(E, Field(Average(model.energy, dims=(1, 2)))) +ρe′ = model.energy - ρE +ρe = model.energy +T = model.temperature -#outputs = merge(model.velocities, model.tracers, (; T, ζ, θ′)) -outputs = merge(model.velocities, model.tracers, (; ζ, θ′)) +outputs = merge(model.velocities, model.tracers, (; ζ, ρe′, ρe, T)) filename = "thermal_bubble_$(Nx)x$(Nz).jld2" writer = JLD2Writer(model, outputs; filename, @@ -117,22 +117,22 @@ if get(ENV, "CI", "false") == "false" @info "Creating visualization..." # Read the output data - θt = FieldTimeSeries(filename, "θ") + ρet = FieldTimeSeries(filename, "ρe") Tt = FieldTimeSeries(filename, "T") ut = FieldTimeSeries(filename, "u") wt = FieldTimeSeries(filename, "w") ζt = FieldTimeSeries(filename, "ζ") - θ′t = FieldTimeSeries(filename, "θ′") + ρe′t = FieldTimeSeries(filename, "ρe′") - times = θt.times - Nt = length(θt) + times = ρet.times + Nt = length(ρet) # Create visualization - expanded to 6 panels fig = Figure(size=(1500, 1000), fontsize=12) # Subplot layout - 3 rows, 2 columns - axθ = Axis(fig[1, 1], xlabel="x (km)", ylabel="z (km)", title="Potential Temperature θ (ᵒK)") - axθ′ = Axis(fig[1, 2], xlabel="x (km)", ylabel="z (km)", title="Temperature Perturbation θ′ (ᵒK)") + axe = Axis(fig[1, 1], xlabel="x (km)", ylabel="z (km)", title="Energy ρe (J / kg)") + axe′ = Axis(fig[1, 2], xlabel="x (km)", ylabel="z (km)", title="Energy Perturbation ρe′ (J / kg)") axT = Axis(fig[2, 1], xlabel="x (km)", ylabel="z (km)", title="Temperature T (ᵒK)") axζ = Axis(fig[2, 2], xlabel="x (km)", ylabel="z (km)", title="Vorticity ζ (s⁻¹)") axu = Axis(fig[3, 1], xlabel="x (km)", ylabel="z (km)", title="Horizontal Velocity u (m/s)") @@ -143,46 +143,46 @@ if get(ENV, "CI", "false") == "false" n = slider.value # Observable fields - θn = @lift interior(θt[$n], :, 1, :) - θ′n = @lift interior(θ′t[$n], :, 1, :) + ρen = @lift interior(ρet[$n], :, 1, :) + ρe′n = @lift interior(ρe′t[$n], :, 1, :) Tn = @lift interior(Tt[$n], :, 1, :) ζn = @lift interior(ζt[$n], :, 1, :) un = @lift interior(ut[$n], :, 1, :) wn = @lift interior(wt[$n], :, 1, :) # Grid coordinates - x = xnodes(θt) ./ 1000 # Convert to km - z = znodes(θt) ./ 1000 # Convert to km + x = xnodes(ρet) ./ 1000 # Convert to km + z = znodes(ρet) ./ 1000 # Convert to km # Title with time title = @lift "Thermal Bubble Evolution - t = $(prettytime(times[$n]))" fig[0, :] = Label(fig, title, fontsize=16, tellwidth=false) # Create heatmaps - θ_range = (minimum(θt), maximum(θt)) - θ′_range = (minimum(θ′t), maximum(θ′t)) + ρe_range = (minimum(ρet), maximum(ρet)) + ρe′_range = (minimum(ρe′t), maximum(ρe′t)) T_range = (minimum(Tt), maximum(Tt)) ζ_range = maximum(abs, ζt) u_range = maximum(abs, ut) w_range = maximum(abs, wt) - hmθ = heatmap!(axθ, x, z, θn, colorrange=θ_range, colormap=:thermal) - hmθ′ = heatmap!(axθ′, x, z, θ′n, colorrange=θ′_range, colormap=:balance) + hme = heatmap!(axe, x, z, ρen, colorrange=ρe_range, colormap=:thermal) + hme′ = heatmap!(axe′, x, z, ρe′n, colorrange=ρe′_range, colormap=:balance) hmT = heatmap!(axT, x, z, Tn, colorrange=T_range, colormap=:thermal) hmζ = heatmap!(axζ, x, z, ζn, colorrange=(-ζ_range, ζ_range), colormap=:balance) hmu = heatmap!(axu, x, z, un, colorrange=(-u_range, u_range), colormap=:balance) hmw = heatmap!(axw, x, z, wn, colorrange=(-w_range, w_range), colormap=:balance) # Add colorbars - Colorbar(fig[1, 3], hmθ, label="θ (ᵒK)", vertical=true) - Colorbar(fig[1, 4], hmθ′, label="θ′ (ᵒK)", vertical=true) + Colorbar(fig[1, 3], hme, label="ρe (J/kg)", vertical=true) + Colorbar(fig[1, 4], hme′, label="ρe′ (J/kg)", vertical=true) Colorbar(fig[2, 3], hmT, label="T (ᵒK)", vertical=true) Colorbar(fig[2, 4], hmζ, label="ζ (s⁻¹)", vertical=true) Colorbar(fig[3, 3], hmu, label="u (m/s)", vertical=true) Colorbar(fig[3, 4], hmw, label="w (m/s)", vertical=true) # Set axis limits - for ax in [axθ, axθ′, axT, axζ, axu, axw] + for ax in [axe, axe′, axT, axζ, axu, axw] xlims!(ax, 0, Lx/1000) ylims!(ax, 0, Lz/1000) end @@ -197,6 +197,7 @@ if get(ENV, "CI", "false") == "false" @info "Drawing frame $nn of $Nt..." n[] = nn end + @info "Saved animation to thermal_bubble.mp4" end diff --git a/src/AtmosphereModels/AtmosphereModels.jl b/src/AtmosphereModels/AtmosphereModels.jl index 68752a8f..62ce1921 100644 --- a/src/AtmosphereModels/AtmosphereModels.jl +++ b/src/AtmosphereModels/AtmosphereModels.jl @@ -5,7 +5,8 @@ export AtmosphereModel, AnelasticFormulation include("atmosphere_model.jl") include("anelastic_formulation.jl") include("saturation_adjustment.jl") -include("update_hydrostatic_pressure.jl") +# include("update_hydrostatic_pressure.jl") +include("dynamics_kernel_functions.jl") include("update_atmosphere_model_state.jl") include("set_atmosphere_model.jl") diff --git a/src/AtmosphereModels/anelastic_formulation.jl b/src/AtmosphereModels/anelastic_formulation.jl index a317920f..2045d97c 100644 --- a/src/AtmosphereModels/anelastic_formulation.jl +++ b/src/AtmosphereModels/anelastic_formulation.jl @@ -7,9 +7,10 @@ using ..Thermodynamics: mixture_heat_capacity, dry_air_gas_constant +using Oceananigans: Oceananigans using Oceananigans.Architectures: architecture using Oceananigans.Grids: inactive_cell, prettysummary -using Oceananigans.Operators: Δzᵃᵃᶜ, Δzᵃᵃᶠ, divᶜᶜᶜ +using Oceananigans.Operators: Δzᵃᵃᶜ, Δzᵃᵃᶠ, divᶜᶜᶜ, Δzᶜᶜᶜ using Oceananigans.Solvers: solve! using KernelAbstractions: @kernel, @index @@ -38,7 +39,7 @@ end Base.show(io::IO, formulation::AnelasticFormulation) = print(io, "AnelasticFormulation") -field_names(::AnelasticFormulation, tracer_names) = (:ρu, :ρv, :ρw, :ρe, :ρq, tracer_names...) +field_names(::AnelasticFormulation, tracer_names) = (:ρu, :ρv, :ρw, :ρe, :ρqᵗ, tracer_names...) struct AnelasticThermodynamicState{FT} potential_temperature :: FT @@ -110,7 +111,7 @@ function collect_prognostic_fields(::AnelasticFormulation, condensates, tracers) - thermodynamic_variables = (ρe=energy, ρq=absolute_humidity) + thermodynamic_variables = (ρe=energy, ρqᵗ=absolute_humidity) return merge(momentum, thermodynamic_variables, condensates, tracers) end @@ -144,7 +145,15 @@ tridiagonal_direction(formulation::AnelasticTridiagonalSolverFormulation) = ZDir function formulation_pressure_solver(anelastic_formulation::AnelasticFormulation, grid) reference_density = anelastic_formulation.reference_density tridiagonal_formulation = AnelasticTridiagonalSolverFormulation(reference_density) - solver = FourierTridiagonalPoissonSolver(grid; tridiagonal_formulation) + + solver = if grid isa Oceananigans.ImmersedBoundaries.ImmersedBoundaryGrid + # With this method, we are using an approximate solver that + # will produce a divergent velocity field near terrain. + FourierTridiagonalPoissonSolver(grid.underlying_grid; tridiagonal_formulation) + else # the solver is exact + FourierTridiagonalPoissonSolver(grid; tridiagonal_formulation) + end + return solver end @@ -160,22 +169,22 @@ end @kernel function _compute_anelastic_main_diagonal!(D, grid, λx, λy, reference_density) i, j = @index(Global, NTuple) Nz = size(grid, 3) - ρʳ = reference_density + ρᵣ = reference_density # Using a homogeneous Neumann (zero Gradient) boundary condition: @inbounds begin - ρ¹ = ρʳ[1, 1, 1] - ρᴺ = ρʳ[1, 1, Nz] - ρ̄² = ℑzᵃᵃᶠ(i, j, 2, grid, ρʳ) - ρ̄ᴺ = ℑzᵃᵃᶠ(i, j, Nz, grid, ρʳ) + ρ¹ = ρᵣ[1, 1, 1] + ρᴺ = ρᵣ[1, 1, Nz] + ρ̄² = ℑzᵃᵃᶠ(i, j, 2, grid, ρᵣ) + ρ̄ᴺ = ℑzᵃᵃᶠ(i, j, Nz, grid, ρᵣ) D[i, j, 1] = - ρ̄² / Δzᵃᵃᶠ(i, j, 2, grid) - ρ¹ * Δzᵃᵃᶜ(i, j, 1, grid) * (λx[i] + λy[j]) D[i, j, Nz] = - ρ̄ᴺ / Δzᵃᵃᶠ(i, j, Nz, grid) - ρᴺ * Δzᵃᵃᶜ(i, j, Nz, grid) * (λx[i] + λy[j]) for k in 2:Nz-1 - ρᵏ = ρʳ[1, 1, k] - ρ̄⁺ = ℑzᵃᵃᶠ(i, j, k+1, grid, ρʳ) - ρ̄ᵏ = ℑzᵃᵃᶠ(i, j, k, grid, ρʳ) + ρᵏ = ρᵣ[1, 1, k] + ρ̄⁺ = ℑzᵃᵃᶠ(i, j, k+1, grid, ρᵣ) + ρ̄ᵏ = ℑzᵃᵃᶠ(i, j, k, grid, ρᵣ) D[i, j, k] = - (ρ̄⁺ / Δzᵃᵃᶠ(i, j, k+1, grid) + ρ̄ᵏ / Δzᵃᵃᶠ(i, j, k, grid)) - ρᵏ * Δzᵃᵃᶜ(i, j, k, grid) * (λx[i] + λy[j]) end @@ -203,7 +212,7 @@ function compute_pressure_correction!(model::AnelasticModel, Δt) foreach(mask_immersed_field!, model.momentum) fill_halo_regions!(model.momentum, model.clock, fields(model)) - ρʳ = model.formulation.reference_density + ρᵣ = model.formulation.reference_density ρŨ = model.momentum solver = model.pressure_solver pₙ = model.nonhydrostatic_pressure @@ -256,15 +265,15 @@ Update the predictor momentum (ρu, ρv, ρw) with the non-hydrostatic pressure u^{n+1} = u^n - δₓp_{NH} / Δx * Δt """ -@kernel function _pressure_correct_momentum!(M, grid, Δt, αʳ_pₙ, ρʳ) +@kernel function _pressure_correct_momentum!(M, grid, Δt, αᵣ_pₙ, ρᵣ) i, j, k = @index(Global, NTuple) - ρᶠ = ℑzᵃᵃᶠ(i, j, k, grid, ρʳ) - ρᶜ = @inbounds ρʳ[i, j, k] + ρᶠ = ℑzᵃᵃᶠ(i, j, k, grid, ρᵣ) + ρᶜ = @inbounds ρᵣ[i, j, k] - @inbounds M.ρu[i, j, k] -= ρᶜ * Δt * ∂xᶠᶜᶜ(i, j, k, grid, αʳ_pₙ) - @inbounds M.ρv[i, j, k] -= ρᶜ * Δt * ∂yᶜᶠᶜ(i, j, k, grid, αʳ_pₙ) - @inbounds M.ρw[i, j, k] -= ρᶠ * Δt * ∂zᶜᶜᶠ(i, j, k, grid, αʳ_pₙ) + @inbounds M.ρu[i, j, k] -= ρᶜ * Δt * ∂xᶠᶜᶜ(i, j, k, grid, αᵣ_pₙ) + @inbounds M.ρv[i, j, k] -= ρᶜ * Δt * ∂yᶜᶠᶜ(i, j, k, grid, αᵣ_pₙ) + @inbounds M.ρw[i, j, k] -= ρᶠ * Δt * ∂zᶜᶜᶠ(i, j, k, grid, αᵣ_pₙ) end function make_pressure_correction!(model::AnelasticModel, Δt) diff --git a/src/AtmosphereModels/atmosphere_model.jl b/src/AtmosphereModels/atmosphere_model.jl index d0912c22..f6ffdfde 100644 --- a/src/AtmosphereModels/atmosphere_model.jl +++ b/src/AtmosphereModels/atmosphere_model.jl @@ -7,14 +7,16 @@ using ..Thermodynamics: mixture_heat_capacity, dry_air_gas_constant -using Oceananigans: Oceananigans, AbstractModel, Center, CenterField, Clock, Face, Field, Flat, - WENO, XFaceField, YFaceField, ZFaceField +using Oceananigans: AbstractModel, Center, CenterField, Clock, Field +using Oceananigans: WENO, XFaceField, YFaceField, ZFaceField using Oceananigans.Advection: adapt_advection_order using Oceananigans.BoundaryConditions: FieldBoundaryConditions, regularize_field_boundary_conditions using Oceananigans.Grids: ZDirection using Oceananigans.Solvers: FourierTridiagonalPoissonSolver using Oceananigans.TimeSteppers: TimeStepper -using Oceananigans.Utils: launch! +using Oceananigans.Utils: launch!, prettytime, prettykeys + +import Oceananigans.Advection: cell_advection_timescale using KernelAbstractions: @kernel, @index @@ -116,7 +118,7 @@ function AtmosphereModel(grid; boundary_conditions = NamedTuple(), forcing = NamedTuple(), advection = WENO(order=5), - microphysics = WarmPhaseSaturationAdjustment(), + microphysics = nothing, # WarmPhaseSaturationAdjustment(), timestepper = :RungeKutta3) arch = grid.architecture @@ -140,11 +142,11 @@ function AtmosphereModel(grid; advection = adapt_advection_order(advection, grid) if absolute_humidity isa DefaultValue - absolute_humidity = CenterField(grid, boundary_conditions=boundary_conditions.ρq) + absolute_humidity = CenterField(grid, boundary_conditions=boundary_conditions.ρqᵗ) end energy = CenterField(grid, boundary_conditions=boundary_conditions.ρe) - specific_humidity = CenterField(grid, boundary_conditions=boundary_conditions.ρq) + specific_humidity = CenterField(grid, boundary_conditions=boundary_conditions.ρqᵗ) temperature = CenterField(grid) prognostic_fields = collect_prognostic_fields(formulation, @@ -192,8 +194,6 @@ function AtmosphereModel(grid; return model end -using Oceananigans.Utils: prettytime, prettykeys - function Base.summary(model::AtmosphereModel) A = nameof(typeof(model.grid.architecture)) G = nameof(typeof(model.grid)) @@ -215,3 +215,5 @@ function Base.show(io::IO, model::AtmosphereModel) "├── coriolis: ", summary(model.coriolis), "\n", "└── microphysics: ", Mic) end + +cell_advection_timescale(model::AtmosphereModel) = cell_advection_timescale(model.grid, model.velocities) \ No newline at end of file diff --git a/src/AtmosphereModels/dynamics_kernel_functions.jl b/src/AtmosphereModels/dynamics_kernel_functions.jl new file mode 100644 index 00000000..8a454732 --- /dev/null +++ b/src/AtmosphereModels/dynamics_kernel_functions.jl @@ -0,0 +1,120 @@ +using Oceananigans.Advection: div_𝐯u, div_𝐯v, div_𝐯w, div_Uc +using Oceananigans.Coriolis: x_f_cross_U, y_f_cross_U, z_f_cross_U +using Oceananigans.Operators: ∂xᶠᶜᶜ, ∂yᶜᶠᶜ, ∂zᶜᶜᶠ, ℑzᵃᵃᶜ, ℑzᵃᵃᶠ + +##### +##### Some key functions +##### + +@inline function buoyancy(i, j, k, grid, formulation, temperature, specific_humidity, thermo) + α = specific_volume(i, j, k, grid, formulation, temperature, specific_humidity, thermo) + αᵣ = reference_specific_volume(i, j, k, grid, formulation, thermo) + g = thermo.gravitational_acceleration + return g * (α - αᵣ) / αᵣ +end + +@inline function ρ_bᶜᶜᶠ(i, j, k, grid, ρ, T, q, formulation, thermo) + ρᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, ρ) + bᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, buoyancy, formulation, T, q, thermo) + return ρᶜᶜᶠ * bᶜᶜᶠ +end + +@inline function ρ_w_bᶜᶜᶠ(i, j, k, grid, w, ρ, T, q, formulation, thermo) + ρ_b = ρ_bᶜᶜᶠ(i, j, k, grid, ρ, T, q, formulation, thermo) + return @inbounds ρ_b * w[i, j, k] +end + +# Note: these are unused currently +hydrostatic_pressure_gradient_x(i, j, k, grid, pₕ′) = ∂xᶠᶜᶜ(i, j, k, grid, pₕ′) +hydrostatic_pressure_gradient_y(i, j, k, grid, pₕ′) = ∂yᶜᶠᶜ(i, j, k, grid, pₕ′) + +@inline function x_momentum_tendency(i, j, k, grid, + advection, + velocities, + momentum, + coriolis, + clock, + model_fields, + forcing, + reference_density, + hydrostatic_pressure_anomaly) + + return ( - div_𝐯u(i, j, k, grid, advection, velocities, momentum.ρu) + - x_f_cross_U(i, j, k, grid, coriolis, momentum) + # - hydrostatic_pressure_gradient_x(i, j, k, grid, hydrostatic_pressure_anomaly) + + forcing(i, j, k, grid, clock, model_fields)) +end + +@inline function y_momentum_tendency(i, j, k, grid, + advection, + velocities, + momentum, + coriolis, + clock, + model_fields, + forcing, + reference_density, + hydrostatic_pressure_anomaly) + + return ( - div_𝐯v(i, j, k, grid, advection, velocities, momentum.ρv) + - y_f_cross_U(i, j, k, grid, coriolis, momentum) + # - hydrostatic_pressure_gradient_y(i, j, k, grid, hydrostatic_pressure_anomaly) + + forcing(i, j, k, grid, clock, model_fields)) +end + +@inline function z_momentum_tendency(i, j, k, grid, + advection, + velocities, + momentum, + coriolis, + clock, + model_fields, + forcing, + reference_density, + formulation, + temperature, + specific_humidity, + thermo) + + return ( - div_𝐯w(i, j, k, grid, advection, velocities, momentum.ρw) + - z_f_cross_U(i, j, k, grid, coriolis, momentum) + + ρ_bᶜᶜᶠ(i, j, k, grid, reference_density, temperature, specific_humidity, formulation, thermo) + + forcing(i, j, k, grid, clock, model_fields)) +end + +@inline function scalar_tendency(i, j, k, grid, + scalar, + forcing, + advection, + velocities, + clock, + model_fields) + + return ( - div_Uc(i, j, k, grid, advection, velocities, scalar) + + forcing(i, j, k, grid, clock, model_fields)) +end + +@inline function moist_static_energy_tendency(i, j, k, grid, + moist_static_energy, + forcing, + advection, + velocities, + clock, + model_fields, + reference_density, + formulation, + temperature, + specific_humidity, + thermo, + condensates, + microphysics) + + # Compute the buoyancy flux term, ρᵣ w b + buoyancy_flux = ℑzᵃᵃᶜ(i, j, k, grid, ρ_w_bᶜᶜᶠ, velocities.w, reference_density, + temperature, specific_humidity, formulation, thermo) + + return ( - div_Uc(i, j, k, grid, advection, velocities, moist_static_energy) + + buoyancy_flux + # + microphysical_energy_tendency(i, j, k, grid, formulation, microphysics, condensates) + + forcing(i, j, k, grid, clock, model_fields)) +end diff --git a/src/AtmosphereModels/saturation_adjustment.jl b/src/AtmosphereModels/saturation_adjustment.jl index dc0a6282..05cfc214 100644 --- a/src/AtmosphereModels/saturation_adjustment.jl +++ b/src/AtmosphereModels/saturation_adjustment.jl @@ -3,9 +3,9 @@ ##### function condensate_specific_humidity(T, state::AnelasticThermodynamicState, thermo) - qᵛ★ = saturation_specific_humidity(T, state.reference_density, thermo, thermo.liquid) + qᵛ⁺ = saturation_specific_humidity(T, state.reference_density, thermo, thermo.liquid) q = state.specific_humidity - return max(0, q - qᵛ★) + return max(0, q - qᵛ⁺) end @inline function compute_temperature(state::AnelasticThermodynamicState{FT}, thermo) where FT diff --git a/src/AtmosphereModels/set_atmosphere_model.jl b/src/AtmosphereModels/set_atmosphere_model.jl index 505dd67d..3d5eef9a 100644 --- a/src/AtmosphereModels/set_atmosphere_model.jl +++ b/src/AtmosphereModels/set_atmosphere_model.jl @@ -23,24 +23,24 @@ function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) θ = model.temperature # use scratch set!(θ, value) - ρʳ = model.formulation.reference_density + ρᵣ = model.formulation.reference_density cᵖᵈ = model.thermodynamics.dry_air.heat_capacity ϕ = model.energy - value = ρʳ * cᵖᵈ * θ + value = ρᵣ * cᵖᵈ * θ elseif name == :q q = model.specific_humidity set!(q, value) - ρʳ = model.formulation.reference_density + ρᵣ = model.formulation.reference_density ϕ = model.absolute_humidity - value = ρʳ * q + value = ρᵣ * q elseif name ∈ (:u, :v, :w) u = model.velocities[name] set!(u, value) - ρʳ = model.formulation.reference_density + ρᵣ = model.formulation.reference_density ϕ = model.momentum[Symbol(:ρ, name)] - value = ρʳ * u + value = ρᵣ * u end set!(ϕ, value) diff --git a/src/AtmosphereModels/update_atmosphere_model_state.jl b/src/AtmosphereModels/update_atmosphere_model_state.jl index 329cc783..8a96711d 100644 --- a/src/AtmosphereModels/update_atmosphere_model_state.jl +++ b/src/AtmosphereModels/update_atmosphere_model_state.jl @@ -6,6 +6,7 @@ using ..Thermodynamics: using Oceananigans.Architectures: architecture using Oceananigans.BoundaryConditions: fill_halo_regions!, compute_x_bcs!, compute_y_bcs!, compute_z_bcs! using Oceananigans.ImmersedBoundaries: mask_immersed_field! +using Oceananigans.Utils: launch! import Oceananigans: fields, prognostic_fields import Oceananigans.TimeSteppers: update_state!, compute_flux_bc_tendencies! @@ -13,7 +14,7 @@ import Oceananigans.TimeSteppers: update_state!, compute_flux_bc_tendencies! const AnelasticModel = AtmosphereModel{<:AnelasticFormulation} function prognostic_fields(model::AnelasticModel) - thermodynamic_fields = (ρe=model.energy, ρq=model.absolute_humidity) + thermodynamic_fields = (ρe=model.energy, ρqᵗ=model.absolute_humidity) return merge(model.momentum, thermodynamic_fields, model.condensates, model.tracers) end @@ -22,7 +23,7 @@ fields(model::AnelasticModel) = prognostic_fields(model) function update_state!(model::AnelasticModel, callbacks=[]; compute_tendencies=true) fill_halo_regions!(prognostic_fields(model), model.clock, fields(model), async=true) compute_auxiliary_variables!(model) - update_hydrostatic_pressure!(model) + # update_hydrostatic_pressure!(model) compute_tendencies && compute_tendencies!(model) return nothing end @@ -45,6 +46,7 @@ function compute_auxiliary_variables!(model) grid, model.thermodynamics, formulation, + model.microphysics, model.energy, model.absolute_humidity) @@ -75,6 +77,7 @@ end grid, thermo, formulation, + microphysics, energy, absolute_humidity) i, j, k = @index(Global, NTuple) @@ -82,24 +85,18 @@ end 𝒰 = thermodynamic_state(i, j, k, grid, formulation, thermo, energy, absolute_humidity) @inbounds specific_humidity[i, j, k] = 𝒰.specific_humidity - # Saturation adjustment - T = compute_temperature(𝒰, thermo) - @inbounds temperature[i, j, k] = T -end + # Possibly perform saturation adjustment + # Note, we will make this much prettier in the future + T = if isnothing(microphysics) + Π = 𝒰.exner_function + θ = 𝒰.potential_temperature + Π * θ + else + compute_temperature(𝒰, microphysics) + end -#= -@inline function specific_volume(state, ref, thermo) - T = temperature(state, ref, thermo) - Rᵐ = mixture_gas_constant(state.q, thermo) - pᵣ = reference_pressure(state.z, ref, thermo) - return Rᵐ * T / pᵣ + @inbounds temperature[i, j, k] = T end -=# - -using Oceananigans.Advection: div_𝐯u, div_𝐯v, div_𝐯w, div_Uc -using Oceananigans.Coriolis: x_f_cross_U, y_f_cross_U, z_f_cross_U -using Oceananigans.Operators: ∂xᶠᶜᶜ, ∂yᶜᶠᶜ, ∂zᶜᶜᶠ -using Oceananigans.Utils: launch! function compute_tendencies!(model::AnelasticModel) grid = model.grid @@ -131,26 +128,31 @@ function compute_tendencies!(model::AnelasticModel) Gρe = model.timestepper.Gⁿ.ρe ρe = model.energy Fρe = model.forcing.ρe - ρe_args = tuple(ρe, Fρe, scalar_args...) - launch!(arch, grid, :xyz, compute_scalar_tendency!, Gρe, grid, ρe_args) + ρe_args = tuple(ρe, Fρe, scalar_args..., ρᵣ, + model.formulation, model.temperature, + model.specific_humidity, model.thermodynamics, model.condensates, model.microphysics) + launch!(arch, grid, :xyz, compute_moist_static_energy_tendency!, Gρe, grid, ρe_args) ρq = model.absolute_humidity - Gρq = model.timestepper.Gⁿ.ρq - Fρq = model.forcing.ρq - ρq_args = tuple(ρq, Fρq, scalar_args...) - launch!(arch, grid, :xyz, compute_scalar_tendency!, Gρq, grid, ρq_args) + Gρqᵗ = model.timestepper.Gⁿ.ρqᵗ + Fρqᵗ = model.forcing.ρqᵗ + ρq_args = tuple(ρq, Fρqᵗ, scalar_args...) + launch!(arch, grid, :xyz, compute_scalar_tendency!, Gρqᵗ, grid, ρq_args) return nothing end -hydrostatic_pressure_gradient_x(i, j, k, grid, pₕ′) = ∂xᶠᶜᶜ(i, j, k, grid, pₕ′) -hydrostatic_pressure_gradient_y(i, j, k, grid, pₕ′) = ∂yᶜᶠᶜ(i, j, k, grid, pₕ′) - +# See dynamics_kernel_functions.jl @kernel function compute_scalar_tendency!(Gc, grid, args) i, j, k = @index(Global, NTuple) @inbounds Gc[i, j, k] = scalar_tendency(i, j, k, grid, args...) end +@kernel function compute_moist_static_energy_tendency!(Gρe, grid, args) + i, j, k = @index(Global, NTuple) + @inbounds Gρe[i, j, k] = moist_static_energy_tendency(i, j, k, grid, args...) +end + @kernel function compute_x_momentum_tendency!(Gρu, grid, args) i, j, k = @index(Global, NTuple) @inbounds Gρu[i, j, k] = x_momentum_tendency(i, j, k, grid, args...) @@ -166,106 +168,11 @@ end @inbounds Gρw[i, j, k] = z_momentum_tendency(i, j, k, grid, args...) end -@inline function x_momentum_tendency(i, j, k, grid, - advection, - velocities, - momentum, - coriolis, - clock, - model_fields, - forcing, - reference_density, - hydrostatic_pressure_anomaly) - - # Note: independent of x - ρᵣ = @inbounds reference_density[i, j, k] - - return ( - div_𝐯u(i, j, k, grid, advection, velocities, momentum.ρu) - - x_f_cross_U(i, j, k, grid, coriolis, momentum) - # - hydrostatic_pressure_gradient_x(i, j, k, grid, hydrostatic_pressure_anomaly) - + forcing(i, j, k, grid, clock, model_fields)) -end - -@inline function y_momentum_tendency(i, j, k, grid, - advection, - velocities, - momentum, - coriolis, - clock, - model_fields, - forcing, - reference_density, - hydrostatic_pressure_anomaly) - - # Note: independent of y - ρᵣ = @inbounds reference_density[i, j, k] - - return ( - div_𝐯v(i, j, k, grid, advection, velocities, momentum.ρv) - - y_f_cross_U(i, j, k, grid, coriolis, momentum) - # - hydrostatic_pressure_gradient_y(i, j, k, grid, hydrostatic_pressure_anomaly) - + forcing(i, j, k, grid, clock, model_fields)) -end - -@inline function z_momentum_tendency(i, j, k, grid, - advection, - velocities, - momentum, - coriolis, - clock, - model_fields, - forcing, - reference_density, - formulation, - temperature, - specific_humidity, - thermo) - - ρᵣᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, reference_density) - bᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, buoyancy, - formulation, temperature, specific_humidity, thermo) - - return ( - div_𝐯w(i, j, k, grid, advection, velocities, momentum.ρw) - - z_f_cross_U(i, j, k, grid, coriolis, momentum) - + ρᵣᶜᶜᶠ * bᶜᶜᶠ - + forcing(i, j, k, grid, clock, model_fields)) -end - -@inline function scalar_tendency(i, j, k, grid, - scalar, - forcing, - advection, - velocities, - clock, - model_fields) - - return ( - div_Uc(i, j, k, grid, advection, velocities, scalar) - + forcing(i, j, k, grid, clock, model_fields)) -end - -#= -@inline function energy_tendency(i, j, k, grid, - formulation, - energy, - forcing, - advection, - velocities, - condensates, - microphysics - clock, - model_fields) - - return ( - div_Uc(i, j, k, grid, advection, velocities, energy) - + microphysical_energy_tendency(i, j, k, grid, formulation, microphysics, condensates) - + forcing(i, j, k, grid, clock, model_fields)) -end -=# - """ Apply boundary conditions by adding flux divergences to the right-hand-side. """ function compute_flux_bc_tendencies!(model::AtmosphereModel) Gⁿ = model.timestepper.Gⁿ arch = model.architecture - clock = model.clock # Compute boundary flux contributions prognostic_model_fields = prognostic_fields(model) @@ -277,4 +184,4 @@ function compute_flux_bc_tendencies!(model::AtmosphereModel) foreach(q -> compute_z_bcs!(Gⁿ[q], prognostic_model_fields[q], args...), field_indices) return nothing -end +end \ No newline at end of file diff --git a/src/AtmosphereModels/update_hydrostatic_pressure.jl b/src/AtmosphereModels/update_hydrostatic_pressure.jl index d9758a53..fa1ad9bc 100644 --- a/src/AtmosphereModels/update_hydrostatic_pressure.jl +++ b/src/AtmosphereModels/update_hydrostatic_pressure.jl @@ -9,13 +9,6 @@ using Oceananigans.Utils: KernelParameters const c = Center() const f = Face() -@inline function buoyancy(i, j, k, grid, formulation, temperature, specific_humidity, thermo) - α = specific_volume(i, j, k, grid, formulation, temperature, specific_humidity, thermo) - αʳ = reference_specific_volume(i, j, k, grid, formulation, thermo) - g = thermo.gravitational_acceleration - return g * (α - αʳ) / αʳ -end - @kernel function _update_hydrostatic_pressure!(pₕ′, grid, formulation, T, q, thermo) i, j = @index(Global, NTuple) diff --git a/src/Breeze.jl b/src/Breeze.jl index 9e9059c4..1eab2927 100644 --- a/src/Breeze.jl +++ b/src/Breeze.jl @@ -38,6 +38,7 @@ using Oceananigans: Oceananigans, @at, AnisotropicMinimumDissipation, Average, minimum_zspacing, nodes, prettytime, regrid!, run!, set!, time_step!, xnodes, xspacings, ynodes, yspacings, znodes, zspacings, ∂x, ∂y, ∂z + using Oceananigans.Grids: znode export diff --git a/src/MoistAirBuoyancies.jl b/src/MoistAirBuoyancies.jl index 67365a8e..eb2c499f 100644 --- a/src/MoistAirBuoyancies.jl +++ b/src/MoistAirBuoyancies.jl @@ -105,14 +105,14 @@ const c = Center() α = specific_volume(𝒰, mb.reference_constants, mb.thermodynamics) # Compute reference specific volume - αʳ = reference_specific_volume(z, mb.reference_constants, mb.thermodynamics) + αᵣ = reference_specific_volume(z, mb.reference_constants, mb.thermodynamics) g = mb.thermodynamics.gravitational_acceleration # Formulation in terms of base density: # ρ₀ = base_density(mb.reference_constants, mb.thermodynamics) - # return ρ₀ * g * (α - αʳ) + # return ρ₀ * g * (α - αᵣ) - return g * (α - αʳ) / αʳ + return g * (α - αᵣ) / αᵣ end @inline ∂z_b(i, j, k, grid, mb::MoistAirBuoyancy, tracers) = @@ -239,7 +239,7 @@ condensate_specific_humidity(T, state::HeightReferenceThermodynamicState, ref, t # Solve # θ = T/Π ( 1 - ℒ qˡ / (cᵖᵐ T)) -# for temperature T with qˡ = max(0, q - qᵛ★). +# for temperature T with qˡ = max(0, q - qᵛ⁺). # root of: f(T) = T - Π θ - ℒ qˡ / cᵖᵐ """ diff --git a/src/Thermodynamics/atmosphere_thermodynamics.jl b/src/Thermodynamics/atmosphere_thermodynamics.jl index 0f286df7..8d9c8f91 100644 --- a/src/Thermodynamics/atmosphere_thermodynamics.jl +++ b/src/Thermodynamics/atmosphere_thermodynamics.jl @@ -356,6 +356,6 @@ condensate_specific_humidity(T, state, ref, thermo) = condensate_specific_humidity(T, state.q, state.z, ref, thermo) function condensate_specific_humidity(T, q, z, ref, thermo) - qᵛ★ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid) - return max(0, q - qᵛ★) + qᵛ⁺ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid) + return max(0, q - qᵛ⁺) end diff --git a/src/Thermodynamics/reference_states.jl b/src/Thermodynamics/reference_states.jl index 9eb14330..5856740b 100644 --- a/src/Thermodynamics/reference_states.jl +++ b/src/Thermodynamics/reference_states.jl @@ -62,11 +62,11 @@ end end function condensate_specific_humidity(T, q, z, ref::ReferenceStateConstants, thermo) - qᵛ★ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid) - return max(0, q - qᵛ★) + qᵛ⁺ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid) + return max(0, q - qᵛ⁺) end function ice_specific_humidity(T, q, z, ref::ReferenceStateConstants, thermo) - qi★ = saturation_specific_humidity(T, z, ref, thermo, thermo.solid) - return max(0, q - qi★) + qi⁺ = saturation_specific_humidity(T, z, ref, thermo, thermo.solid) + return max(0, q - qi⁺) end diff --git a/src/microphysics.jl b/src/microphysics.jl index 0c807df7..f1f1aa0c 100644 --- a/src/microphysics.jl +++ b/src/microphysics.jl @@ -4,7 +4,7 @@ # Solve # θ = T/Π ( 1 - ℒ qˡ / (cᵖᵐ T)) -# for temperature T with qˡ = max(0, q - qᵛ★). +# for temperature T with qˡ = max(0, q - qᵛ⁺). # root of: f(T) = T² - Π θ T - ℒ qˡ / cᵖᵐ @inline function temperature(state::ThermodynamicState{FT}, ref, thermo) where FT state.θ == 0 && return zero(FT)