diff --git a/docs/src/momentum_equation.md b/docs/src/momentum_equation.md index 8b041de0..d12b8c2f 100644 --- a/docs/src/momentum_equation.md +++ b/docs/src/momentum_equation.md @@ -66,8 +66,9 @@ KernelFunctionOperation at (Center, Face, Center) The combined RHS — `Tendency(model)` — wraps Oceananigans' internal `*_velocity_tendency` kernel directly. The sum of the individual term diagnostics, taken with the correct -signs from Oceananigans' source, equals `Tendency` to machine precision. For -`NonhydrostaticModel` the u/v formula is +signs from Oceananigans' source, equals `Tendency` to machine precision. + +### `NonhydrostaticModel` u/v budget ```math \text{Tendency} = -\text{Advection} + \text{BuoyancyAcceleration} @@ -77,30 +78,53 @@ signs from Oceananigans' source, equals `Tendency` to machine precision. For +\text{StokesShear} + \text{StokesTendency} + \text{Forcing} ``` -The w-momentum equation is similar but with one subtlety: Oceananigans' `w_velocity_tendency` -has no explicit `-\partial_z p` line. The vertical hydrostatic balance is treated as a -property of the pressure-projection step, and whether the buoyancy term appears in the -tendency depends on whether the model uses hydrostatic/nonhydrostatic pressure splitting: +### `NonhydrostaticModel` w budget + +The w-momentum equation has one subtlety: Oceananigans' `w_velocity_tendency` has no +explicit `-\partial_z p` line. The vertical hydrostatic balance is treated as a property +of the pressure-projection step, and whether the buoyancy term appears in the tendency +depends on whether the model uses hydrostatic/nonhydrostatic pressure splitting: | `hydrostatic_pressure_anomaly` | Buoyancy term in `Tendency` | Budget formula | | --- | --- | --- | | `Field` (default, splitting on) | absorbed by ``pHY'`` | `-ADV - COR - TVISC + SS + ST + FORC` | | `nothing` (splitting off) | included as ``+b`` | `-ADV + BUOY - COR - TVISC + SS + ST + FORC` | +### `HydrostaticFreeSurfaceModel` u/v budget + +HFS uses `hydrostatic_free_surface_{u,v}_velocity_tendency`, which differs from the NH +tendency: there are no Stokes terms, no buoyancy term (absorbed into the hydrostatic +pressure anomaly `pHY′`), and the free-surface contribution to the horizontal pressure +gradient appears as a separate `BarotropicPressureGradient` term: + +```math +\text{Tendency} = -\text{Advection} + -\text{BarotropicPressureGradient} + -\text{CoriolisAcceleration} + -\text{PressureGradient} + -\text{TotalViscousDissipation} + +\text{Forcing} +``` + +For `ImplicitFreeSurface` (the default) and `SplitExplicitFreeSurface`, +`BarotropicPressureGradient` returns zero (the contribution is handled inside the pressure +solve), so the budget closes without it. For `ExplicitFreeSurface` it returns `g ∂_i η` +and must be included for closure. + ## `HydrostaticFreeSurfaceModel` caveats - `WMomentumEquation.Tendency`, `WMomentumEquation.Forcing(model, Val(:w))`, and the U/V/W `StokesShear` and `StokesTendency` throw an `ArgumentError` on `HydrostaticFreeSurfaceModel`: `w` is diagnosed from continuity rather than evolved by a prognostic equation, and HFS has no `stokes_drift` field. -- `Advection(::HydrostaticFreeSurfaceModel)` wraps `div_𝐯u`/`div_𝐯v` (the flux-form - kernel). HFS's default `VectorInvariant` momentum advection uses `U_dot_∇u` instead, so - to use the `Advection` diagnostic with HFS you must build the model with a flux-form - scheme, e.g. `momentum_advection = Centered()`. -- For HFS, `Tendency` wraps `hydrostatic_free_surface_*_velocity_tendency`, which differs - from the NH tendency: it has no Stokes terms, includes a barotropic free-surface - pressure gradient, and a grid-slope contribution. Budget closure for HFS is not yet - supported by this module — it would require additional diagnostics for those terms. +- `Advection(::HydrostaticFreeSurfaceModel)` (and the explicit-args form + `Advection(::HydrostaticFreeSurfaceModel, u, v, w, scheme)`) both wrap + `U_dot_∇u`/`U_dot_∇v` — the vector-invariant kernel that HFS's tendency actually uses. + This works with HFS's default `VectorInvariant` momentum advection as well as any + flux-form scheme. The diagnostic returns a KFO whose kernel function is + `U_dot_∇u`/`U_dot_∇v`, whereas the NH `Advection` diagnostic returns a KFO whose kernel + function is `div_𝐯u`/`div_𝐯v` — both pass the `isa Advection` type check via the + `Advection` Union alias. ## `Forcing` dispatch quirk @@ -125,6 +149,7 @@ Oceanostics.UMomentumEquation.Advection Oceanostics.UMomentumEquation.BuoyancyAcceleration Oceanostics.UMomentumEquation.CoriolisAcceleration Oceanostics.UMomentumEquation.PressureGradient +Oceanostics.UMomentumEquation.BarotropicPressureGradient Oceanostics.UMomentumEquation.ViscousDissipation Oceanostics.UMomentumEquation.ImmersedViscousDissipation Oceanostics.UMomentumEquation.TotalViscousDissipation @@ -141,6 +166,7 @@ Oceanostics.VMomentumEquation.Advection Oceanostics.VMomentumEquation.BuoyancyAcceleration Oceanostics.VMomentumEquation.CoriolisAcceleration Oceanostics.VMomentumEquation.PressureGradient +Oceanostics.VMomentumEquation.BarotropicPressureGradient Oceanostics.VMomentumEquation.ViscousDissipation Oceanostics.VMomentumEquation.ImmersedViscousDissipation Oceanostics.VMomentumEquation.TotalViscousDissipation diff --git a/src/Oceanostics.jl b/src/Oceanostics.jl index 43cb5253..74b74374 100755 --- a/src/Oceanostics.jl +++ b/src/Oceanostics.jl @@ -14,13 +14,13 @@ export TracerAdvection, TracerDiffusion, TracerImmersedDiffusion, TracerTotalDif #--- #+++ UMomentumEquation exports -export UAdvection, UBuoyancyAcceleration, UCoriolisAcceleration, UPressureGradient, +export UAdvection, UBuoyancyAcceleration, UCoriolisAcceleration, UPressureGradient, UBarotropicPressureGradient, UViscousDissipation, UImmersedViscousDissipation, UTotalViscousDissipation, UStokesShear, UStokesTendency, UForcing, UTendency #--- #+++ VMomentumEquation exports -export VAdvection, VBuoyancyAcceleration, VCoriolisAcceleration, VPressureGradient, +export VAdvection, VBuoyancyAcceleration, VCoriolisAcceleration, VPressureGradient, VBarotropicPressureGradient, VViscousDissipation, VImmersedViscousDissipation, VTotalViscousDissipation, VStokesShear, VStokesTendency, VForcing, VTendency #--- diff --git a/src/UMomentumEquation.jl b/src/UMomentumEquation.jl index 708dee33..a36b8b56 100644 --- a/src/UMomentumEquation.jl +++ b/src/UMomentumEquation.jl @@ -5,19 +5,20 @@ using Oceananigans: fields, Face, Center, KernelFunctionOperation, AbstractModel using Oceananigans.Models: HydrostaticFreeSurfaceModel using Oceananigans.Models.NonhydrostaticModels: u_velocity_tendency using Oceananigans.Models.HydrostaticFreeSurfaceModels: hydrostatic_free_surface_u_velocity_tendency -using Oceananigans.Advection: div_𝐯u +using Oceananigans.Advection: div_𝐯u, U_dot_∇u using Oceananigans.BuoyancyFormulations: x_dot_g_bᶠᶜᶜ using Oceananigans.Coriolis: x_f_cross_U using Oceananigans.TurbulenceClosures: ∂ⱼ_τ₁ⱼ, immersed_∂ⱼ_τ₁ⱼ using Oceananigans.StokesDrifts: x_curl_Uˢ_cross_U, ∂t_uˢ using Oceananigans.Operators: ∂xᶠᶜᶜ +using Oceananigans.Models.HydrostaticFreeSurfaceModels: explicit_barotropic_pressure_x_gradient using Oceanostics: validate_location, CustomKFO -export Advection, BuoyancyAcceleration, CoriolisAcceleration, PressureGradient, +export Advection, BuoyancyAcceleration, CoriolisAcceleration, PressureGradient, BarotropicPressureGradient, ViscousDissipation, ImmersedViscousDissipation, TotalViscousDissipation, StokesShear, StokesTendency, Forcing, Tendency, - UAdvection, UBuoyancyAcceleration, UCoriolisAcceleration, UPressureGradient, + UAdvection, UBuoyancyAcceleration, UCoriolisAcceleration, UPressureGradient, UBarotropicPressureGradient, UViscousDissipation, UImmersedViscousDissipation, UTotalViscousDissipation, UStokesShear, UStokesTendency, UForcing, UTendency @@ -30,11 +31,13 @@ export Advection, BuoyancyAcceleration, CoriolisAcceleration, PressureGradient, @inline hydrostatic_pressure_gradient_x(i, j, k, grid, hydrostatic_pressure) = ∂xᶠᶜᶜ(i, j, k, grid, hydrostatic_pressure) @inline hydrostatic_pressure_gradient_x(i, j, k, grid, ::Nothing) = zero(grid) -# Type aliases for major functions -const Advection = CustomKFO{<:typeof(div_𝐯u)} +# Type aliases for major functions. `Advection` is a Union so that the NH flux-form +# (`div_𝐯u`) and the HFS vector-invariant (`U_dot_∇u`) wrappers share the same alias. +const Advection = Union{CustomKFO{<:typeof(div_𝐯u)}, CustomKFO{<:typeof(U_dot_∇u)}} const BuoyancyAcceleration = CustomKFO{<:typeof(x_dot_g_bᶠᶜᶜ)} const CoriolisAcceleration = CustomKFO{<:typeof(x_f_cross_U)} const PressureGradient = CustomKFO{<:typeof(hydrostatic_pressure_gradient_x)} +const BarotropicPressureGradient = CustomKFO{<:typeof(explicit_barotropic_pressure_x_gradient)} const ViscousDissipation = CustomKFO{<:typeof(∂ⱼ_τ₁ⱼ)} const ImmersedViscousDissipation = CustomKFO{<:typeof(immersed_∂ⱼ_τ₁ⱼ)} const TotalViscousDissipation = CustomKFO{<:typeof(total_∂ⱼ_τ₁ⱼ)} @@ -42,12 +45,13 @@ const StokesShear = CustomKFO{<:typeof(x_curl_Uˢ_cross_U)} const StokesTendency = CustomKFO{<:typeof(∂t_uˢ)} const Forcing = KernelFunctionOperation{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any} const Tendency = Union{CustomKFO{<:typeof(u_velocity_tendency)}, - CustomKFO{<:typeof(hydrostatic_free_surface_u_velocity_tendency)}} + CustomKFO{<:typeof(hydrostatic_free_surface_u_velocity_tendency)}} const UAdvection = Advection const UBuoyancyAcceleration = BuoyancyAcceleration const UCoriolisAcceleration = CoriolisAcceleration const UPressureGradient = PressureGradient +const UBarotropicPressureGradient = BarotropicPressureGradient const UViscousDissipation = ViscousDissipation const UImmersedViscousDissipation = ImmersedViscousDissipation const UTotalViscousDissipation = TotalViscousDissipation @@ -60,11 +64,16 @@ const UTendency = Tendency """ $(SIGNATURES) -Calculate the advection of u-momentum as +Calculate the advection of u-momentum. For `NonhydrostaticModel` this wraps the flux-form +divergence kernel [`div_𝐯u`](https://clima.github.io/OceananigansDocumentation/stable/appendix/library/#Oceananigans.Advection.div_𝐯u-NTuple{7,%20Any}) +as ADV = ∂ⱼ (uⱼ u) -using Oceananigans' kernel [`div_𝐯u`.](https://clima.github.io/OceananigansDocumentation/stable/appendix/library/#Oceananigans.Advection.div_𝐯u-NTuple{7,%20Any}) +For `HydrostaticFreeSurfaceModel` this wraps the vector-invariant kernel `U_dot_∇u` (which +is what `hydrostatic_free_surface_u_velocity_tendency` actually uses) as + + ADV = (u · ∇) u ```jldoctest julia> using Oceananigans, Oceanostics @@ -86,8 +95,18 @@ function Advection(model, u, v, w, advection_scheme; location = (Face, Center, C return KernelFunctionOperation{Face, Center, Center}(div_𝐯u, model.grid, advection_scheme, total_velocities, u) end +function Advection(model::HydrostaticFreeSurfaceModel, advection_scheme, velocities; location = (Face, Center, Center)) + validate_location(location, "Advection", (Face, Center, Center)) + return KernelFunctionOperation{Face, Center, Center}(U_dot_∇u, model.grid, advection_scheme, velocities) +end + +# HFS-typed 5-arg overload so that explicit calls `Advection(hfs_model, u, v, w, scheme)` also use +# `U_dot_∇u` rather than dispatching to the generic NH `div_𝐯u` form. +Advection(model::HydrostaticFreeSurfaceModel, u, v, w, advection_scheme; kwargs...) = + Advection(model, advection_scheme, (; u, v, w); kwargs...) + Advection(model; kwargs...) = Advection(model, model.velocities..., model.advection; kwargs...) -Advection(model::HydrostaticFreeSurfaceModel; kwargs...) = Advection(model, model.velocities..., model.advection.momentum; kwargs...) +Advection(model::HydrostaticFreeSurfaceModel; kwargs...) = Advection(model, model.advection.momentum, model.velocities; kwargs...) #--- #+++ Buoyancy acceleration @@ -198,6 +217,41 @@ function PressureGradient(model; kwargs...) end #--- +#+++ Barotropic pressure gradient +""" + $(SIGNATURES) + +Calculate the explicit barotropic pressure gradient in the x-direction. For +`HydrostaticFreeSurfaceModel` this is the contribution of the free surface ``η`` to the +horizontal pressure gradient that's treated explicitly during time-stepping (i.e., not +solved implicitly). For `ExplicitFreeSurface` it equals `g ∂x η`; for `ImplicitFreeSurface` +or `SplitExplicitFreeSurface` the contribution is handled inside the pressure solve so +this kernel returns zero. For `NonhydrostaticModel` (`free_surface == nothing`) it also +returns zero. This term completes the u-momentum budget closure on HFS. + +```jldoctest +julia> using Oceananigans, Oceanostics + +julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1)); + +julia> model = NonhydrostaticModel(grid); + +julia> BARO = UMomentumEquation.BarotropicPressureGradient(model) +KernelFunctionOperation at (Face, Center, Center) +├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── kernel_function: explicit_barotropic_pressure_x_gradient (generic function with 4 methods) +└── arguments: ("Nothing",) +``` +""" +function BarotropicPressureGradient(model, free_surface; location = (Face, Center, Center)) + validate_location(location, "BarotropicPressureGradient", (Face, Center, Center)) + return KernelFunctionOperation{Face, Center, Center}(explicit_barotropic_pressure_x_gradient, model.grid, free_surface) +end + +BarotropicPressureGradient(model; kwargs...) = + BarotropicPressureGradient(model, hasfield(typeof(model), :free_surface) ? model.free_surface : nothing; kwargs...) +#--- + #+++ Viscous dissipation """ $(SIGNATURES) diff --git a/src/VMomentumEquation.jl b/src/VMomentumEquation.jl index 1e4260ba..ec03a599 100644 --- a/src/VMomentumEquation.jl +++ b/src/VMomentumEquation.jl @@ -5,19 +5,20 @@ using Oceananigans: fields, Face, Center, KernelFunctionOperation, AbstractModel using Oceananigans.Models: HydrostaticFreeSurfaceModel using Oceananigans.Models.NonhydrostaticModels: v_velocity_tendency using Oceananigans.Models.HydrostaticFreeSurfaceModels: hydrostatic_free_surface_v_velocity_tendency -using Oceananigans.Advection: div_𝐯v +using Oceananigans.Advection: div_𝐯v, U_dot_∇v using Oceananigans.BuoyancyFormulations: y_dot_g_bᶜᶠᶜ using Oceananigans.Coriolis: y_f_cross_U using Oceananigans.TurbulenceClosures: ∂ⱼ_τ₂ⱼ, immersed_∂ⱼ_τ₂ⱼ using Oceananigans.StokesDrifts: y_curl_Uˢ_cross_U, ∂t_vˢ using Oceananigans.Operators: ∂yᶜᶠᶜ +using Oceananigans.Models.HydrostaticFreeSurfaceModels: explicit_barotropic_pressure_y_gradient using Oceanostics: validate_location, CustomKFO -export Advection, BuoyancyAcceleration, CoriolisAcceleration, PressureGradient, +export Advection, BuoyancyAcceleration, CoriolisAcceleration, PressureGradient, BarotropicPressureGradient, ViscousDissipation, ImmersedViscousDissipation, TotalViscousDissipation, StokesShear, StokesTendency, Forcing, Tendency, - VAdvection, VBuoyancyAcceleration, VCoriolisAcceleration, VPressureGradient, + VAdvection, VBuoyancyAcceleration, VCoriolisAcceleration, VPressureGradient, VBarotropicPressureGradient, VViscousDissipation, VImmersedViscousDissipation, VTotalViscousDissipation, VStokesShear, VStokesTendency, VForcing, VTendency @@ -31,10 +32,11 @@ export Advection, BuoyancyAcceleration, CoriolisAcceleration, PressureGradient, @inline hydrostatic_pressure_gradient_y(i, j, k, grid, ::Nothing) = zero(grid) # Type aliases for major functions -const Advection = CustomKFO{<:typeof(div_𝐯v)} +const Advection = Union{CustomKFO{<:typeof(div_𝐯v)}, CustomKFO{<:typeof(U_dot_∇v)}} const BuoyancyAcceleration = CustomKFO{<:typeof(y_dot_g_bᶜᶠᶜ)} const CoriolisAcceleration = CustomKFO{<:typeof(y_f_cross_U)} const PressureGradient = CustomKFO{<:typeof(hydrostatic_pressure_gradient_y)} +const BarotropicPressureGradient = CustomKFO{<:typeof(explicit_barotropic_pressure_y_gradient)} const ViscousDissipation = CustomKFO{<:typeof(∂ⱼ_τ₂ⱼ)} const ImmersedViscousDissipation = CustomKFO{<:typeof(immersed_∂ⱼ_τ₂ⱼ)} const TotalViscousDissipation = CustomKFO{<:typeof(total_∂ⱼ_τ₂ⱼ)} @@ -42,12 +44,13 @@ const StokesShear = CustomKFO{<:typeof(y_curl_Uˢ_cross_U)} const StokesTendency = CustomKFO{<:typeof(∂t_vˢ)} const Forcing = KernelFunctionOperation{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any} const Tendency = Union{CustomKFO{<:typeof(v_velocity_tendency)}, - CustomKFO{<:typeof(hydrostatic_free_surface_v_velocity_tendency)}} + CustomKFO{<:typeof(hydrostatic_free_surface_v_velocity_tendency)}} const VAdvection = Advection const VBuoyancyAcceleration = BuoyancyAcceleration const VCoriolisAcceleration = CoriolisAcceleration const VPressureGradient = PressureGradient +const VBarotropicPressureGradient = BarotropicPressureGradient const VViscousDissipation = ViscousDissipation const VImmersedViscousDissipation = ImmersedViscousDissipation const VTotalViscousDissipation = TotalViscousDissipation @@ -60,11 +63,16 @@ const VTendency = Tendency """ $(SIGNATURES) -Calculate the advection of v-momentum as +Calculate the advection of v-momentum. For `NonhydrostaticModel` this wraps the flux-form +divergence kernel [`div_𝐯v`](https://clima.github.io/OceananigansDocumentation/stable/appendix/library/#Oceananigans.Advection.div_𝐯v-NTuple{7,%20Any}) +as ADV = ∂ⱼ (uⱼ v) -using Oceananigans' kernel [`div_𝐯v`.](https://clima.github.io/OceananigansDocumentation/stable/appendix/library/#Oceananigans.Advection.div_𝐯v-NTuple{7,%20Any}) +For `HydrostaticFreeSurfaceModel` this wraps the vector-invariant kernel `U_dot_∇v` (which +is what `hydrostatic_free_surface_v_velocity_tendency` actually uses) as + + ADV = (u · ∇) v ```jldoctest julia> using Oceananigans, Oceanostics @@ -86,8 +94,18 @@ function Advection(model, u, v, w, advection_scheme; location = (Center, Face, C return KernelFunctionOperation{Center, Face, Center}(div_𝐯v, model.grid, advection_scheme, total_velocities, v) end +function Advection(model::HydrostaticFreeSurfaceModel, advection_scheme, velocities; location = (Center, Face, Center)) + validate_location(location, "Advection", (Center, Face, Center)) + return KernelFunctionOperation{Center, Face, Center}(U_dot_∇v, model.grid, advection_scheme, velocities) +end + +# HFS-typed 5-arg overload so that explicit calls `Advection(hfs_model, u, v, w, scheme)` also use +# `U_dot_∇v` rather than dispatching to the generic NH `div_𝐯v` form. +Advection(model::HydrostaticFreeSurfaceModel, u, v, w, advection_scheme; kwargs...) = + Advection(model, advection_scheme, (; u, v, w); kwargs...) + Advection(model; kwargs...) = Advection(model, model.velocities..., model.advection; kwargs...) -Advection(model::HydrostaticFreeSurfaceModel; kwargs...) = Advection(model, model.velocities..., model.advection.momentum; kwargs...) +Advection(model::HydrostaticFreeSurfaceModel; kwargs...) = Advection(model, model.advection.momentum, model.velocities; kwargs...) #--- #+++ Buoyancy acceleration @@ -198,6 +216,41 @@ function PressureGradient(model; kwargs...) end #--- +#+++ Barotropic pressure gradient +""" + $(SIGNATURES) + +Calculate the explicit barotropic pressure gradient in the y-direction. For +`HydrostaticFreeSurfaceModel` this is the contribution of the free surface ``η`` to the +horizontal pressure gradient that's treated explicitly during time-stepping (i.e., not +solved implicitly). For `ExplicitFreeSurface` it equals `g ∂y η`; for `ImplicitFreeSurface` +or `SplitExplicitFreeSurface` the contribution is handled inside the pressure solve so +this kernel returns zero. For `NonhydrostaticModel` (`free_surface == nothing`) it also +returns zero. This term completes the v-momentum budget closure on HFS. + +```jldoctest +julia> using Oceananigans, Oceanostics + +julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1)); + +julia> model = NonhydrostaticModel(grid); + +julia> BARO = VMomentumEquation.BarotropicPressureGradient(model) +KernelFunctionOperation at (Center, Face, Center) +├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── kernel_function: explicit_barotropic_pressure_y_gradient (generic function with 4 methods) +└── arguments: ("Nothing",) +``` +""" +function BarotropicPressureGradient(model, free_surface; location = (Center, Face, Center)) + validate_location(location, "BarotropicPressureGradient", (Center, Face, Center)) + return KernelFunctionOperation{Center, Face, Center}(explicit_barotropic_pressure_y_gradient, model.grid, free_surface) +end + +BarotropicPressureGradient(model; kwargs...) = + BarotropicPressureGradient(model, hasfield(typeof(model), :free_surface) ? model.free_surface : nothing; kwargs...) +#--- + #+++ Viscous dissipation """ $(SIGNATURES) diff --git a/test/test_u_momentum_diagnostics.jl b/test/test_u_momentum_diagnostics.jl index 84c8dbde..5563bc48 100644 --- a/test/test_u_momentum_diagnostics.jl +++ b/test/test_u_momentum_diagnostics.jl @@ -109,6 +109,20 @@ function test_u_momentum_terms(model) @test PRES isa UPressureGradient @test PRES_field isa Field + # Test BarotropicPressureGradient + free_surface = hasfield(typeof(model), :free_surface) ? model.free_surface : nothing + BARO = UMomentumEquation.BarotropicPressureGradient(model, free_surface) + BARO_field = Field(BARO) + @test BARO isa UMomentumEquation.BarotropicPressureGradient + @test BARO isa UBarotropicPressureGradient + @test BARO_field isa Field + + BARO = UMomentumEquation.BarotropicPressureGradient(model) + BARO_field = Field(BARO) + @test BARO isa UMomentumEquation.BarotropicPressureGradient + @test BARO isa UBarotropicPressureGradient + @test BARO_field isa Field + # Test ViscousDissipation VISC = UMomentumEquation.ViscousDissipation(model, model.closure, model.closure_fields, model.clock, fields(model), model.buoyancy) VISC_field = Field(VISC) @@ -228,6 +242,9 @@ function test_u_momentum_field_locations(model) PRES = UMomentumEquation.PressureGradient(model) @test location(PRES) == (Face, Center, Center) + BARO = UMomentumEquation.BarotropicPressureGradient(model) + @test location(BARO) == (Face, Center, Center) + VISC = UMomentumEquation.ViscousDissipation(model) @test location(VISC) == (Face, Center, Center) @@ -290,12 +307,48 @@ function test_u_momentum_budget_closure(grid) return nothing end +function test_u_momentum_hfs_budget_closure(grid) + # Build an HFS model with every term active so the budget exercises the full RHS of + # hydrostatic_free_surface_u_velocity_tendency. The default `momentum_advection` + # (`VectorInvariant`) is used — the diagnostic `Advection` dispatches on the model type + # to wrap `U_dot_∇u` with whatever scheme the model carries. + model = HydrostaticFreeSurfaceModel(grid; tracers = :b, + buoyancy = BuoyancyTracer(), + coriolis = FPlane(f = 1e-4), + closure = ScalarDiffusivity(ν = 1e-4, κ = 1e-4), + forcing = (; u = Forcing((x, y, z, t) -> cos(t)))) + set!(model, u = (x, y, z) -> sin(2π*x) * cos(2π*y) * exp(z), + v = (x, y, z) -> cos(2π*x) * sin(2π*y) * exp(z), + b = (x, y, z) -> sin(2π*z)) + update_state!(model) # populates model.pressure.pHY′ from b + + # Reconstruct G_u for HFS: + # G_u = -ADV - BARO - COR - PRES - TVISC + FORCING + # HFS has no Stokes terms and no BuoyancyAcceleration term (buoyancy is absorbed into pHY′). + # BARO is the explicit barotropic free-surface gradient; on ImplicitFreeSurface it returns + # zero (the contribution is handled inside the pressure solve) so the budget closes + # without it, but including it keeps the formula correct for ExplicitFreeSurface too. + ADV = UMomentumEquation.Advection(model) + COR = UMomentumEquation.CoriolisAcceleration(model) + PRES = UMomentumEquation.PressureGradient(model) + BARO = UMomentumEquation.BarotropicPressureGradient(model) + TVISC = UMomentumEquation.TotalViscousDissipation(model) + FORC = UMomentumEquation.Forcing(model) + TEND = UMomentumEquation.Tendency(model) + + budget = Field(-ADV - BARO - COR - PRES - TVISC + FORC) + tend = Field(TEND) + @test interior(budget) ≈ interior(tend) + return nothing +end + function test_u_momentum_location_validation(model) # Test that invalid locations throw errors @test_throws ArgumentError UMomentumEquation.Advection(model; location = (Center, Center, Center)) @test_throws ArgumentError UMomentumEquation.BuoyancyAcceleration(model; location = (Center, Face, Center)) @test_throws ArgumentError UMomentumEquation.CoriolisAcceleration(model; location = (Center, Center, Face)) @test_throws ArgumentError UMomentumEquation.PressureGradient(model; location = (Center, Center, Center)) + @test_throws ArgumentError UMomentumEquation.BarotropicPressureGradient(model; location = (Center, Center, Center)) @test_throws ArgumentError UMomentumEquation.ViscousDissipation(model; location = (Center, Center, Center)) @test_throws ArgumentError UMomentumEquation.ImmersedViscousDissipation(model; location = (Center, Center, Center)) @test_throws ArgumentError UMomentumEquation.TotalViscousDissipation(model; location = (Center, Center, Center)) @@ -316,6 +369,7 @@ end @test UBuoyancyAcceleration !== VBuoyancyAcceleration && UBuoyancyAcceleration !== WBuoyancyAcceleration @test UCoriolisAcceleration !== VCoriolisAcceleration && UCoriolisAcceleration !== WCoriolisAcceleration @test UPressureGradient !== VPressureGradient # W has no PressureGradient + @test UBarotropicPressureGradient !== VBarotropicPressureGradient # W has no BarotropicPressureGradient @test UViscousDissipation !== VViscousDissipation && UViscousDissipation !== WViscousDissipation @test UImmersedViscousDissipation !== VImmersedViscousDissipation && UImmersedViscousDissipation !== WImmersedViscousDissipation @test UTotalViscousDissipation !== VTotalViscousDissipation && UTotalViscousDissipation !== WTotalViscousDissipation @@ -334,12 +388,12 @@ end @info " with $grid_class" for model_type in model_types @info " with $model_type" - # HFS defaults to VectorInvariant momentum advection, which uses U_dot_∇u - # rather than div_𝐯u. Force the flux-form scheme so the Advection diagnostic - # (which wraps div_𝐯u) is well-defined. HFS has no stokes_drift field, so we - # pass the Stokes-drift only to the NH model. + # HFS uses the default `VectorInvariant` momentum advection. The diagnostic + # `Advection` dispatches on the model type to wrap `U_dot_∇u` with whatever + # scheme the model carries, so no override is needed. HFS has no + # `stokes_drift` field, so we pass the Stokes-drift only to the NH model. model = model_type === HydrostaticFreeSurfaceModel ? - model_type(grid; model_kwargs..., momentum_advection=Centered()) : + model_type(grid; model_kwargs...) : model_type(grid; nh_model_kwargs...) @info " Testing u-momentum terms" @@ -355,4 +409,7 @@ end @info " Testing u-momentum budget closure on NonhydrostaticModel" test_u_momentum_budget_closure(underlying_regular_grid) + + @info " Testing u-momentum budget closure on HydrostaticFreeSurfaceModel" + test_u_momentum_hfs_budget_closure(underlying_regular_grid) end diff --git a/test/test_v_momentum_diagnostics.jl b/test/test_v_momentum_diagnostics.jl index c1614183..33939e38 100644 --- a/test/test_v_momentum_diagnostics.jl +++ b/test/test_v_momentum_diagnostics.jl @@ -109,6 +109,20 @@ function test_v_momentum_terms(model) @test PRES isa VPressureGradient @test PRES_field isa Field + # Test BarotropicPressureGradient + free_surface = hasfield(typeof(model), :free_surface) ? model.free_surface : nothing + BARO = VMomentumEquation.BarotropicPressureGradient(model, free_surface) + BARO_field = Field(BARO) + @test BARO isa VMomentumEquation.BarotropicPressureGradient + @test BARO isa VBarotropicPressureGradient + @test BARO_field isa Field + + BARO = VMomentumEquation.BarotropicPressureGradient(model) + BARO_field = Field(BARO) + @test BARO isa VMomentumEquation.BarotropicPressureGradient + @test BARO isa VBarotropicPressureGradient + @test BARO_field isa Field + # Test ViscousDissipation VISC = VMomentumEquation.ViscousDissipation(model, model.closure, model.closure_fields, model.clock, fields(model), model.buoyancy) VISC_field = Field(VISC) @@ -228,6 +242,9 @@ function test_v_momentum_field_locations(model) PRES = VMomentumEquation.PressureGradient(model) @test location(PRES) == (Center, Face, Center) + BARO = VMomentumEquation.BarotropicPressureGradient(model) + @test location(BARO) == (Center, Face, Center) + VISC = VMomentumEquation.ViscousDissipation(model) @test location(VISC) == (Center, Face, Center) @@ -290,12 +307,48 @@ function test_v_momentum_budget_closure(grid) return nothing end +function test_v_momentum_hfs_budget_closure(grid) + # Build an HFS model with every term active so the budget exercises the full RHS of + # hydrostatic_free_surface_v_velocity_tendency. The default `momentum_advection` + # (`VectorInvariant`) is used — the diagnostic `Advection` dispatches on the model type + # to wrap `U_dot_∇v` with whatever scheme the model carries. + model = HydrostaticFreeSurfaceModel(grid; tracers = :b, + buoyancy = BuoyancyTracer(), + coriolis = FPlane(f = 1e-4), + closure = ScalarDiffusivity(ν = 1e-4, κ = 1e-4), + forcing = (; v = Forcing((x, y, z, t) -> cos(t)))) + set!(model, u = (x, y, z) -> sin(2π*x) * cos(2π*y) * exp(z), + v = (x, y, z) -> cos(2π*x) * sin(2π*y) * exp(z), + b = (x, y, z) -> sin(2π*z)) + update_state!(model) # populates model.pressure.pHY′ from b + + # Reconstruct G_v for HFS: + # G_v = -ADV - BARO - COR - PRES - TVISC + FORCING + # HFS has no Stokes terms and no BuoyancyAcceleration term (buoyancy is absorbed into pHY′). + # BARO is the explicit barotropic free-surface gradient; on ImplicitFreeSurface it returns + # zero (the contribution is handled inside the pressure solve) so the budget closes + # without it, but including it keeps the formula correct for ExplicitFreeSurface too. + ADV = VMomentumEquation.Advection(model) + COR = VMomentumEquation.CoriolisAcceleration(model) + PRES = VMomentumEquation.PressureGradient(model) + BARO = VMomentumEquation.BarotropicPressureGradient(model) + TVISC = VMomentumEquation.TotalViscousDissipation(model) + FORC = VMomentumEquation.Forcing(model, Val(:v)) + TEND = VMomentumEquation.Tendency(model) + + budget = Field(-ADV - BARO - COR - PRES - TVISC + FORC) + tend = Field(TEND) + @test interior(budget) ≈ interior(tend) + return nothing +end + function test_v_momentum_location_validation(model) # Test that invalid locations throw errors @test_throws ArgumentError VMomentumEquation.Advection(model; location = (Center, Center, Center)) @test_throws ArgumentError VMomentumEquation.BuoyancyAcceleration(model; location = (Face, Center, Center)) @test_throws ArgumentError VMomentumEquation.CoriolisAcceleration(model; location = (Center, Center, Face)) @test_throws ArgumentError VMomentumEquation.PressureGradient(model; location = (Center, Center, Center)) + @test_throws ArgumentError VMomentumEquation.BarotropicPressureGradient(model; location = (Center, Center, Center)) @test_throws ArgumentError VMomentumEquation.ViscousDissipation(model; location = (Center, Center, Center)) @test_throws ArgumentError VMomentumEquation.ImmersedViscousDissipation(model; location = (Center, Center, Center)) @test_throws ArgumentError VMomentumEquation.TotalViscousDissipation(model; location = (Center, Center, Center)) @@ -313,12 +366,12 @@ end @info " with $grid_class" for model_type in model_types @info " with $model_type" - # HFS defaults to VectorInvariant momentum advection, which uses U_dot_∇u - # rather than div_𝐯v. Force the flux-form scheme so the Advection diagnostic - # (which wraps div_𝐯v) is well-defined. HFS has no stokes_drift field, so we - # pass the Stokes-drift only to the NH model. + # HFS uses the default `VectorInvariant` momentum advection. The diagnostic + # `Advection` dispatches on the model type to wrap `U_dot_∇v` with whatever + # scheme the model carries, so no override is needed. HFS has no + # `stokes_drift` field, so we pass the Stokes-drift only to the NH model. model = model_type === HydrostaticFreeSurfaceModel ? - model_type(grid; model_kwargs..., momentum_advection=Centered()) : + model_type(grid; model_kwargs...) : model_type(grid; nh_model_kwargs...) @info " Testing v-momentum terms" @@ -334,4 +387,7 @@ end @info " Testing v-momentum budget closure on NonhydrostaticModel" test_v_momentum_budget_closure(underlying_regular_grid) + + @info " Testing v-momentum budget closure on HydrostaticFreeSurfaceModel" + test_v_momentum_hfs_budget_closure(underlying_regular_grid) end