From 7a4b3f04e402be70d45fcb775a4dedef087f3bb0 Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Sun, 5 May 2024 14:49:26 -0600 Subject: [PATCH] (0.91.0) Make hydrostatic pressure anomaly optional in `NonhydrostaticModel` (#3574) * Add kwarg for hydrostatic pressure anomaly, make nothing by default * Modify regression tests to separate hydrostatic pressure * Update docs on pressure decomposition plus nonhydrostatic model docs * Bugfix plus discontinue support for "trying" to set boundary conditions on pressure * Fix regression tests * Bump minor version * Use nonhydrostatic pressure for aux data * Add some validation of the pressure fields * Bugfix * Make sure buoyancy gets added when hydrostatic pressure is omitted * fix typos/phrasing * remove PressureField(s) * Fix AB2 test, which was wrong because update_state computes tendencies now * Relax requirement that w==0 in jld2 output writer test * Try again to get a good comparison * Dont test halo filling for hydrostatic pressure cause who cares --------- Co-authored-by: Navid C. Constantinou --- Project.toml | 2 +- .../pressure_decomposition.md | 43 ++++++--------- src/Fields/Fields.jl | 4 +- src/Fields/field_tuples.jl | 52 ------------------- .../compute_nonhydrostatic_tendencies.jl | 17 ++++-- .../nonhydrostatic_model.jl | 34 ++++++++---- ...onhydrostatic_tendency_kernel_functions.jl | 11 +++- .../update_hydrostatic_pressure.jl | 3 ++ ...n_large_eddy_simulation_regression_test.jl | 4 +- .../rayleigh_benard_regression_test.jl | 24 +++++---- .../thermal_bubble_regression_test.jl | 9 +++- test/test_computed_field.jl | 15 +++--- test/test_distributed_models.jl | 2 +- test/test_jld2_output_writer.jl | 6 ++- test/test_nonhydrostatic_regression.jl | 1 + test/test_time_stepping.jl | 19 ++++--- 16 files changed, 118 insertions(+), 128 deletions(-) diff --git a/Project.toml b/Project.toml index 3fb1256769..552bfcb4b8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Oceananigans" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" authors = ["Climate Modeling Alliance and contributors"] -version = "0.90.14" +version = "0.91.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/docs/src/numerical_implementation/pressure_decomposition.md b/docs/src/numerical_implementation/pressure_decomposition.md index 0a5cf88999..48d6df4008 100644 --- a/docs/src/numerical_implementation/pressure_decomposition.md +++ b/docs/src/numerical_implementation/pressure_decomposition.md @@ -1,62 +1,51 @@ # [Pressure decomposition](@id pressure_decomposition) -In the numerical implementation of the momentum equations, the kinematic pressure ``p`` -is split into "hydrostatic" and "non-hydrostatic" parts via +In the numerical implementation of the momentum equations in the `NonhydrostaticModel`, the kinematic pressure ``p`` +is split into "background" and "dynamic" parts via ```math \begin{equation} \label{eq:pressure} - p(\boldsymbol{x}, t) = p_{\text{total hydrostatic}}(\boldsymbol{x}, t) + p_{\rm{non}}(\boldsymbol{x}, t) \, . + p(\boldsymbol{x}, t) = p_{\text{background}}(\boldsymbol{x}, t) + p'(\boldsymbol{x}, t) \, . \end{equation} ``` -The hydrostatic pressure component in \eqref{eq:pressure} is defined so that the vertical -component of its gradient balances gravity: +The background pressure component in \eqref{eq:pressure} is defined so that the vertical +component of its gradient balances the background density field: ```math \begin{align} - \partial_z p_{\text{total hydrostatic}} & = - g \left ( 1 + \frac{\rho_*}{\rho_0} + \frac{\rho'}{\rho_0} \right ) \, , + \partial_z p_{\text{total hydrostatic}} & = - g \left ( 1 + \frac{\rho_*}{\rho_0} \right ) \, , \end{align} ``` Above, we use the notation introduced in the [Boussinesq approximation](@ref boussinesq_approximation) section. -We can further split the hydrostatic pressure component into +Optionally, we may further decompose the dynamic pressure perturbation ``p'`` into +a "hydrostatic anomaly" and "nonhydrostatic" part: ```math \begin{align} - p_{\text{total hydrostatic}}(\boldsymbol{x}, t) = p_{*}(z) + p_{\rm{hyd}}(\boldsymbol{x}, t) \, , + p'(\boldsymbol{x}, t) = p_{\rm{hyd}}(\boldsymbol(x), t) + p_{\rm{non}}(\boldsymbol{x}, t) \, , \end{align} ``` -i.e., a component that only varies in ``z`` (``p_*``) and a "hydrostatic anomaly" (``p_{\rm{hyd}}``) defined -so that +where ```math \begin{align} - \partial_z p_{*} & = - g \left ( 1 + \frac{\rho_*}{\rho_0} \right ) \, ,\\ - \partial_z p_{\rm{hyd}} & = \underbrace{- g \frac{\rho'}{\rho_0}}_{= b} \, . + \partial_z p_{\rm{hyd}} \equiv \underbrace{- g \frac{\rho'}{\rho_0}}_{= b} \, . \end{align} ``` -Doing so, the gradient of the kinematic pressure becomes: +With this pressure decomposition, the kinematic pressure gradient that appears in the momentum equations +(after we've employed the the [Boussinesq approximation](@ref boussinesq_approximation)) becomes ```math \begin{align} - \boldsymbol{\nabla} p & = \boldsymbol{\nabla} p_{\rm{non}} + \boldsymbol{\nabla}_h p_{\rm{hyd}} + ( \partial_z p_{*} + \partial_z p_{\rm{hyd}} ) \boldsymbol{\hat z}\, , + \boldsymbol{\nabla} p &= - g \frac{\rho}{\rho_0} \hat {\boldsymbol{z}} + \boldsymbol{\nabla} p' + &= - g \frac{\rho}{\rho_0} \hat {\boldsymbol{z}} + \boldsymbol{\nabla} p_{\rm{non}} + \boldsymbol{\nabla}_h p_{\rm{hyd}} \, . \end{align} ``` -where ``\boldsymbol{\nabla}_h \equiv \boldsymbol{\hat x} \partial_x + \boldsymbol{\hat y} \partial_y`` -is the horizontal gradient. +where ``\boldsymbol{\nabla}_h \equiv \boldsymbol{\hat x} \partial_x + \boldsymbol{\hat y} \partial_y``. -Under this pressure decomposition, the kinematic pressure gradient that appears in the momentum equations -(after we've employed the the [Boussinesq approximation](@ref boussinesq_approximation)) combines with -the gravity force to give: - -```math - \begin{align} - \boldsymbol{\nabla} p + g \frac{\rho}{\rho_0} \hat {\boldsymbol{z}} = \boldsymbol{\nabla} p_{\rm{non}} + \boldsymbol{\nabla}_h p_{\rm{hyd}} \, . - \end{align} -``` - -Mathematically, the non-hydrostatic pressure ``p_{\rm{non}}`` enforces the incompressibility constraint. diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl index cae3a72ee9..0c749009b9 100644 --- a/src/Fields/Fields.jl +++ b/src/Fields/Fields.jl @@ -6,7 +6,7 @@ export CenterField, XFaceField, YFaceField, ZFaceField export BackgroundField export interior, data, xnode, ynode, znode, location export set!, compute!, @compute, regrid! -export VelocityFields, TracerFields, tracernames, PressureFields, TendencyFields +export VelocityFields, TracerFields, TendencyFields, tracernames export interpolate using Oceananigans.Architectures @@ -32,7 +32,7 @@ include("broadcasting_abstract_fields.jl") """ field(loc, a, grid) -Build a field from `a` at `loc` and on `grid`. +Build a field from array `a` at `loc` and on `grid`. """ @inline function field(loc, a::AbstractArray, grid) f = Field(loc, grid) diff --git a/src/Fields/field_tuples.jl b/src/Fields/field_tuples.jl index ff426b1c68..debac6aab9 100644 --- a/src/Fields/field_tuples.jl +++ b/src/Fields/field_tuples.jl @@ -203,42 +203,6 @@ TracerFields(::Union{Tuple{}, Nothing}, grid, bcs) = NamedTuple() "Shortcut constructor for empty tracer fields." TracerFields(::NamedTuple{(), Tuple{}}, grid, bcs) = NamedTuple() -##### -##### Pressure fields tuples -##### - -""" - PressureFields(grid, bcs::NamedTuple) - -Return a `NamedTuple` with pressure fields `pHY′` and `pNHS` initialized as -`CenterField`s on `grid`. Boundary conditions `bcs` -may be specified via a named tuple of `FieldBoundaryCondition`s. -""" -function PressureFields(grid, bcs=NamedTuple()) - - default_pressure_boundary_conditions = - (pHY′ = FieldBoundaryConditions(grid, (Center, Center, Center)), - pNHS = FieldBoundaryConditions(grid, (Center, Center, Center))) - - bcs = merge(default_pressure_boundary_conditions, bcs) - - pHY′ = CenterField(grid, boundary_conditions=bcs.pHY′) - pNHS = CenterField(grid, boundary_conditions=bcs.pNHS) - - return (pHY′=pHY′, pNHS=pNHS) -end - -function PressureFields(grid::AbstractGrid{<:Any, <:Any, <:Any, <:Flat}, bcs=NamedTuple()) - default_pressure_boundary_conditions = - (pHY′ = FieldBoundaryConditions(grid, (Center, Center, Center)), - pNHS = FieldBoundaryConditions(grid, (Center, Center, Center))) - - bcs = merge(default_pressure_boundary_conditions, bcs) - pNHS = CenterField(grid, boundary_conditions=bcs.pNHS) - - return (; pHY′=nothing, pNHS=pNHS) -end - """ TendencyFields(grid, tracer_names; u = XFaceField(grid), @@ -268,7 +232,6 @@ end ##### VelocityFields(::Nothing, grid, bcs) = VelocityFields(grid, bcs) -PressureFields(::Nothing, grid, bcs) = PressureFields(grid, bcs) """ VelocityFields(proposed_velocities::NamedTuple{(:u, :v, :w)}, grid, bcs) @@ -302,18 +265,3 @@ function TracerFields(proposed_tracers::NamedTuple, grid, bcs) return NamedTuple{tracer_names}(tracer_fields) end - -""" - PressureFields(proposed_pressures::NamedTuple{(:pHY′, :pNHS)}, grid, bcs) - -Return a `NamedTuple` of pressure fields with, overwriting boundary conditions -in `proposed_tracer_fields` with corresponding fields in the `NamedTuple` `bcs`. -""" -function PressureFields(proposed_pressures::NamedTuple{(:pHY′, :pNHS)}, grid, bcs) - validate_field_tuple_grid("pressures", proposed_pressures, grid) - - pHY′ = CenterField(grid, boundary_conditions=bcs.pHY′, data=proposed_pressures.pHY′.data) - pNHS = CenterField(grid, boundary_conditions=bcs.pNHS, data=proposed_pressures.pNHS.data) - - return (pHY′=pHY′, pNHS=pNHS) -end diff --git a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl index 70bc3b922a..bf840e07b6 100644 --- a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl +++ b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl @@ -85,9 +85,17 @@ function compute_interior_tendency_contributions!(model, kernel_parameters; acti auxiliary_fields, diffusivities) - u_kernel_args = tuple(start_momentum_kernel_args..., u_immersed_bc, end_momentum_kernel_args..., forcings, hydrostatic_pressure, clock) - v_kernel_args = tuple(start_momentum_kernel_args..., v_immersed_bc, end_momentum_kernel_args..., forcings, hydrostatic_pressure, clock) - w_kernel_args = tuple(start_momentum_kernel_args..., w_immersed_bc, end_momentum_kernel_args..., forcings, clock) + u_kernel_args = tuple(start_momentum_kernel_args..., + u_immersed_bc, end_momentum_kernel_args..., + forcings, hydrostatic_pressure, clock) + + v_kernel_args = tuple(start_momentum_kernel_args..., + v_immersed_bc, end_momentum_kernel_args..., + forcings, hydrostatic_pressure, clock) + + w_kernel_args = tuple(start_momentum_kernel_args..., + w_immersed_bc, end_momentum_kernel_args..., + forcings, hydrostatic_pressure, clock) for parameters in kernel_parameters launch!(arch, grid, parameters, compute_Gu!, @@ -104,7 +112,8 @@ function compute_interior_tendency_contributions!(model, kernel_parameters; acti end start_tracer_kernel_args = (advection, closure) - end_tracer_kernel_args = (buoyancy, biogeochemistry, background_fields, velocities, tracers, auxiliary_fields, diffusivities) + end_tracer_kernel_args = (buoyancy, biogeochemistry, background_fields, velocities, + tracers, auxiliary_fields, diffusivities) for tracer_index in 1:length(tracers) @inbounds c_tendency = tendencies[tracer_index + 3] diff --git a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl index 07e0141c1d..b5a22163c3 100644 --- a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl +++ b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl @@ -7,7 +7,7 @@ using Oceananigans.Advection: CenteredSecondOrder using Oceananigans.BuoyancyModels: validate_buoyancy, regularize_buoyancy, SeawaterBuoyancy using Oceananigans.Biogeochemistry: validate_biogeochemistry, AbstractBiogeochemistry, biogeochemical_auxiliary_fields using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions -using Oceananigans.Fields: BackgroundFields, Field, tracernames, VelocityFields, TracerFields, PressureFields +using Oceananigans.Fields: BackgroundFields, Field, tracernames, VelocityFields, TracerFields, CenterField using Oceananigans.Forcings: model_forcing using Oceananigans.Grids: inflate_halo_size, with_halo, architecture using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid @@ -66,7 +66,8 @@ end particles::ParticlesOrNothing = nothing, biogeochemistry::AbstractBGCOrNothing = nothing, velocities = nothing, - pressures = nothing, + nonhydrostatic_pressure = CenterField(grid), + hydrostatic_pressure_anomaly = nothing, diffusivity_fields = nothing, pressure_solver = nothing, immersed_boundary = nothing, @@ -97,7 +98,12 @@ Keyword arguments - `particles`: Lagrangian particles to be advected with the flow. Default: `nothing`. - `biogeochemistry`: Biogeochemical model for `tracers`. - `velocities`: The model velocities. Default: `nothing`. - - `pressures`: Hydrostatic and non-hydrostatic pressure fields. Default: `nothing`. + - `nonhydrostatic_pressure`: The nonhydrostatic pressure field. Default: `CenterField(grid)`. + - `hydrostatic_pressure_anomaly`: An optional field that stores the part of the nonhydrostatic pressure + in hydrostatic balance with the buoyancy field. If `nothing` (default), the anomaly + is not computed. If `CenterField(grid)`, the anomaly is precomputed by + vertically integrating the buoyancy field. In this case, the `nonhydrostatic_pressure` represents + only the part of pressure that deviates from the hydrostatic anomaly. - `diffusivity_fields`: Diffusivity fields. Default: `nothing`. - `pressure_solver`: Pressure solver to be used in the model. If `nothing` (default), the model constructor chooses the default based on the `grid` provide. @@ -119,7 +125,8 @@ function NonhydrostaticModel(; grid, particles::ParticlesOrNothing = nothing, biogeochemistry::AbstractBGCOrNothing = nothing, velocities = nothing, - pressures = nothing, + hydrostatic_pressure_anomaly = nothing, + nonhydrostatic_pressure = CenterField(grid), diffusivity_fields = nothing, pressure_solver = nothing, immersed_boundary = nothing, @@ -129,14 +136,20 @@ function NonhydrostaticModel(; grid, tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example) + # Validate pressure fields + nonhydrostatic_pressure isa Field{Center, Center, Center} || + throw(ArgumentError("nonhydrostatic_pressure must be CenterField(grid).")) + isnothing(hydrostatic_pressure_anomaly) || hydrostatic_pressure_anomaly isa Field{Center, Center, Center} || + throw(ArgumentError("hydrostatic_pressure_anomaly must be `nothing` or `CenterField(grid)`.")) + # We don't support CAKTE for NonhydrostaticModel yet. closure = validate_closure(closure) first_closure = closure isa Tuple ? first(closure) : closure first_closure isa FlavorOfCATKE && - error("CATKEVerticalDiffusivity is not supported for " * - "NonhydrostaticModel --- yet!") + error("CATKEVerticalDiffusivity is not supported for NonhydrostaticModel --- yet!") - tracers, auxiliary_fields = validate_biogeochemistry(tracers, merge(auxiliary_fields, biogeochemical_auxiliary_fields(biogeochemistry)), biogeochemistry, grid, clock) + all_auxiliary_fields = merge(auxiliary_fields, biogeochemical_auxiliary_fields(biogeochemistry)) + tracers, auxiliary_fields = validate_biogeochemistry(tracers, all_auxiliary_fields, biogeochemistry, grid, clock) validate_buoyancy(buoyancy, tracernames(tracers)) buoyancy = regularize_buoyancy(buoyancy) @@ -154,12 +167,12 @@ function NonhydrostaticModel(; grid, # First, we extract boundary conditions that are embedded within any _user-specified_ field tuples: embedded_boundary_conditions = merge(extract_boundary_conditions(velocities), extract_boundary_conditions(tracers), - extract_boundary_conditions(pressures), extract_boundary_conditions(diffusivity_fields)) # Next, we form a list of default boundary conditions: prognostic_field_names = (:u, :v, :w, tracernames(tracers)..., keys(auxiliary_fields)...) - default_boundary_conditions = NamedTuple{prognostic_field_names}(FieldBoundaryConditions() for name in prognostic_field_names) + default_boundary_conditions = NamedTuple{prognostic_field_names}(FieldBoundaryConditions() + for name in prognostic_field_names) # Finally, we merge specified, embedded, and default boundary conditions. Specified boundary conditions # have precedence, followed by embedded, followed by default. @@ -172,7 +185,7 @@ function NonhydrostaticModel(; grid, # Either check grid-correctness, or construct tuples of fields velocities = VelocityFields(velocities, grid, boundary_conditions) tracers = TracerFields(tracers, grid, boundary_conditions) - pressures = PressureFields(pressures, grid, boundary_conditions) + pressures = (pNHS=nonhydrostatic_pressure, pHY′=hydrostatic_pressure_anomaly) diffusivity_fields = DiffusivityFields(diffusivity_fields, grid, tracernames(tracers), boundary_conditions, closure) if isnothing(pressure_solver) @@ -223,3 +236,4 @@ end (u = SumOfArrays{2}(m.velocities.u, m.background_fields.velocities.u), v = SumOfArrays{2}(m.velocities.v, m.background_fields.velocities.v), w = SumOfArrays{2}(m.velocities.w, m.background_fields.velocities.w)) + diff --git a/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl b/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl index e64e110812..cefb5eadc3 100644 --- a/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl +++ b/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl @@ -70,13 +70,13 @@ pressure anomaly. return ( - div_𝐯u(i, j, k, grid, advection, total_velocities, velocities.u) - div_𝐯u(i, j, k, grid, advection, velocities, background_fields.velocities.u) + + x_dot_g_bᶠᶜᶜ(i, j, k, grid, buoyancy, tracers) - x_f_cross_U(i, j, k, grid, coriolis, velocities) - hydrostatic_pressure_gradient_x(i, j, k, grid, hydrostatic_pressure) - ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure, diffusivities, clock, model_fields, buoyancy) - immersed_∂ⱼ_τ₁ⱼ(i, j, k, grid, velocities, u_immersed_bc, closure, diffusivities, clock, model_fields) + x_curl_Uˢ_cross_U(i, j, k, grid, stokes_drift, velocities, clock.time) + ∂t_uˢ(i, j, k, grid, stokes_drift, clock.time) - + x_dot_g_bᶠᶜᶜ(i, j, k, grid, buoyancy, tracers) + forcings.u(i, j, k, grid, clock, model_fields)) end @@ -133,16 +133,21 @@ pressure anomaly. return ( - div_𝐯v(i, j, k, grid, advection, total_velocities, velocities.v) - div_𝐯v(i, j, k, grid, advection, velocities, background_fields.velocities.v) + + y_dot_g_bᶜᶠᶜ(i, j, k, grid, buoyancy, tracers) - y_f_cross_U(i, j, k, grid, coriolis, velocities) - hydrostatic_pressure_gradient_y(i, j, k, grid, hydrostatic_pressure) - ∂ⱼ_τ₂ⱼ(i, j, k, grid, closure, diffusivities, clock, model_fields, buoyancy) - immersed_∂ⱼ_τ₂ⱼ(i, j, k, grid, velocities, v_immersed_bc, closure, diffusivities, clock, model_fields) + y_curl_Uˢ_cross_U(i, j, k, grid, stokes_drift, velocities, clock.time) + ∂t_vˢ(i, j, k, grid, stokes_drift, clock.time) - + y_dot_g_bᶜᶠᶜ(i, j, k, grid, buoyancy, tracers) + forcings.v(i, j, k, grid, clock, model_fields)) end +# Only add buoyancy if the hydrostatic pressure isa Nothing +@inline maybe_z_dot_g_bᶜᶜᶠ(i, j, k, grid, hydrostatic_pressure, buoyancy, tracers) = zero(grid) +@inline maybe_z_dot_g_bᶜᶜᶠ(i, j, k, grid, ::Nothing, buoyancy, tracers) = + z_dot_g_bᶜᶜᶠ(i, j, k, grid, buoyancy, tracers) + """ $(SIGNATURES) @@ -181,6 +186,7 @@ velocity components, tracer fields, and precalculated diffusivities where applic auxiliary_fields, diffusivities, forcings, + hydrostatic_pressure, clock) model_fields = merge(velocities, tracers, auxiliary_fields) @@ -193,6 +199,7 @@ velocity components, tracer fields, and precalculated diffusivities where applic return ( - div_𝐯w(i, j, k, grid, advection, total_velocities, velocities.w) - div_𝐯w(i, j, k, grid, advection, velocities, background_fields.velocities.w) + + maybe_z_dot_g_bᶜᶜᶠ(i, j, k, grid, hydrostatic_pressure, buoyancy, tracers) - z_f_cross_U(i, j, k, grid, coriolis, velocities) - ∂ⱼ_τ₃ⱼ(i, j, k, grid, closure, diffusivities, clock, model_fields, buoyancy) - immersed_∂ⱼ_τ₃ⱼ(i, j, k, grid, velocities, w_immersed_bc, closure, diffusivities, clock, model_fields) diff --git a/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl b/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl index b6b53a86e3..0d02570020 100644 --- a/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl +++ b/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl @@ -34,6 +34,9 @@ update_hydrostatic_pressure!(pHY′, arch, ibg::PCBIBG, buoyancy, tracers; param update_hydrostatic_pressure!(pHY′, arch, grid, buoyancy, tracers; parameters = p_kernel_parameters(grid)) = launch!(arch, grid, parameters, _update_hydrostatic_pressure!, pHY′, grid, buoyancy, tracers) +update_hydrostatic_pressure!(::Nothing, arch, grid, args...; kw...) = nothing +update_hydrostatic_pressure!(::Nothing, arch, ::PCBIBG, args...; kw...) = nothing + # extend p kernel to compute also the boundaries @inline function p_kernel_parameters(grid) Nx, Ny, _ = size(grid) diff --git a/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl b/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl index 7eeed6fa85..395d8ecb60 100644 --- a/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl +++ b/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl @@ -30,11 +30,11 @@ function run_ocean_large_eddy_simulation_regression_test(arch, grid_type, closur equation_of_state = LinearEquationOfState(thermal_expansion=2e-4, haline_contraction=8e-4) # Model instantiation - model = NonhydrostaticModel(; grid, + model = NonhydrostaticModel(; grid, closure, coriolis = FPlane(f=1e-4), buoyancy = SeawaterBuoyancy(; equation_of_state), tracers = (:T, :S), - closure = closure, + hydrostatic_pressure_anomaly = CenterField(grid), boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs)) # The type of the underlying data, not the offset array. diff --git a/test/regression_tests/rayleigh_benard_regression_test.jl b/test/regression_tests/rayleigh_benard_regression_test.jl index 40be053dcf..2f64eb30aa 100644 --- a/test/regression_tests/rayleigh_benard_regression_test.jl +++ b/test/regression_tests/rayleigh_benard_regression_test.jl @@ -34,19 +34,25 @@ function run_rayleigh_benard_regression_test(arch, grid_type) # Force salinity as a passive tracer (βS=0) c★(x, z) = exp(4z) * sin(2π/Lx * x) - Fc(i, j, k, grid, clock, model_fields) = 1/10 * (c★(xnode(i, grid, Center()), znode(k, grid, Center())) - model_fields.c[i, j, k]) + + function Fc(i, j, k, grid, clock, model_fields) + x = xnode(i, grid, Center()) + z = znode(k, grid, Center()) + return 1/10 * (c★(x, z) - model_fields.c[i, j, k]) + end + + cforcing = Forcing(Fc, discrete_form=true) bbcs = FieldBoundaryConditions(top = BoundaryCondition(Value, 0.0), bottom = BoundaryCondition(Value, Δb)) - model = NonhydrostaticModel( - grid = grid, - closure = ScalarDiffusivity(ν=ν, κ=κ), - tracers = (:b, :c), - buoyancy = Buoyancy(model=BuoyancyTracer()), - boundary_conditions = (b=bbcs,), - forcing = (c=Forcing(Fc, discrete_form=true),) - ) + model = NonhydrostaticModel(; grid, + closure = ScalarDiffusivity(ν=ν, κ=κ), + tracers = (:b, :c), + buoyancy = Buoyancy(model=BuoyancyTracer()), + boundary_conditions = (; b=bbcs), + hydrostatic_pressure_anomaly = CenterField(grid), + forcing = (; c=cforcing)) # Lz/Nz will work for both the :regular and :vertically_unstretched grids. Δt = 0.01 * min(model.grid.Δxᶜᵃᵃ, model.grid.Δyᵃᶜᵃ, Lz/Nz)^2 / ν diff --git a/test/regression_tests/thermal_bubble_regression_test.jl b/test/regression_tests/thermal_bubble_regression_test.jl index c75bb6cbda..ad6de22327 100644 --- a/test/regression_tests/thermal_bubble_regression_test.jl +++ b/test/regression_tests/thermal_bubble_regression_test.jl @@ -11,8 +11,13 @@ function run_thermal_bubble_regression_test(arch, grid_type) end closure = ScalarDiffusivity(ν=4e-2, κ=4e-2) - model = NonhydrostaticModel(grid=grid, closure=closure, coriolis=FPlane(f=1e-4), - buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) + + model = NonhydrostaticModel(; grid, closure, + coriolis = FPlane(f=1e-4), + buoyancy = SeawaterBuoyancy(), + hydrostatic_pressure_anomaly = CenterField(grid), + tracers = (:T, :S)) + simulation = Simulation(model, Δt=6, stop_iteration=10) model.tracers.T.data.parent .= 9.85 diff --git a/test/test_computed_field.jl b/test/test_computed_field.jl index 204b4b6540..2f76ceaf23 100644 --- a/test/test_computed_field.jl +++ b/test/test_computed_field.jl @@ -17,7 +17,7 @@ end function compute_unary(unary, model) set!(model; S=π) T, S = model.tracers - @compute uS = Field(unary(S), data=model.pressures.pHY′.data) + @compute uS = Field(unary(S), data=model.pressures.pNHS.data) result = Array(interior(uS)) return all(result .≈ unary(eltype(model.grid)(π))) end @@ -25,7 +25,7 @@ end function compute_plus(model) set!(model; S=π, T=42) T, S = model.tracers - @compute ST = Field(S + T, data=model.pressures.pHY′.data) + @compute ST = Field(S + T, data=model.pressures.pNHS.data) result = Array(interior(ST)) return all(result .≈ eltype(model.grid)(π + 42)) end @@ -42,7 +42,7 @@ end function compute_minus(model) set!(model; S=π, T=42) T, S = model.tracers - @compute ST = Field(S - T, data=model.pressures.pHY′.data) + @compute ST = Field(S - T, data=model.pressures.pNHS.data) result = Array(interior(ST)) return all(result .≈ eltype(model.grid)(π - 42)) end @@ -50,7 +50,7 @@ end function compute_times(model) set!(model; S=π, T=42) T, S = model.tracers - @compute ST = Field(S * T, data=model.pressures.pHY′.data) + @compute ST = Field(S * T, data=model.pressures.pNHS.data) result = Array(interior(ST)) return all(result .≈ eltype(model.grid)(π * 42)) end @@ -62,7 +62,7 @@ function compute_kinetic_energy(model) set!(w, 3) kinetic_energy_operation = @at (Center, Center, Center) (u^2 + v^2 + w^2) / 2 - @compute kinetic_energy = Field(kinetic_energy_operation, data=model.pressures.pHY′.data) + @compute kinetic_energy = Field(kinetic_energy_operation, data=model.pressures.pNHS.data) return all(interior(kinetic_energy, 2:3, 2:3, 2:3) .≈ 7) end @@ -412,7 +412,8 @@ for arch in archs @test begin @inline auxiliary_fields_kernel_function(i, j, k, grid, auxiliary_fields) = 1.0 - op = KernelFunctionOperation{Center, Center, Center}(auxiliary_fields_kernel_function, grid, model.auxiliary_fields) + op = KernelFunctionOperation{Center, Center, Center}(auxiliary_fields_kernel_function, grid, + model.auxiliary_fields) f = Field(op) compute!(f) f isa Field && f.operand === op @@ -471,7 +472,7 @@ for arch in archs set!(model; S=π, T=42) T, S = model.tracers - @compute ST = Field(S + T, data=model.pressures.pHY′.data) + @compute ST = Field(S + T, data=model.pressures.pNHS.data) Nx, Ny, Nz = size(model.grid) Hx, Hy, Hz = halo_size(model.grid) diff --git a/test/test_distributed_models.jl b/test/test_distributed_models.jl index 9ddfe6903c..1648d4b155 100644 --- a/test/test_distributed_models.jl +++ b/test/test_distributed_models.jl @@ -362,7 +362,7 @@ function test_triply_periodic_halo_communication_with_141_ranks(halo, child_arch grid = RectilinearGrid(arch, topology=(Periodic, Periodic, Periodic), size=(8, 8, 8), extent=(1, 2, 3), halo=halo) model = NonhydrostaticModel(grid=grid) - for field in merge(fields(model), model.pressures) + for field in (fields(model)..., model.pressures.pNHS) fill!(field, arch.local_rank) fill_halo_regions!(field) diff --git a/test/test_jld2_output_writer.jl b/test/test_jld2_output_writer.jl index 3fdba376e3..535258586e 100644 --- a/test/test_jld2_output_writer.jl +++ b/test/test_jld2_output_writer.jl @@ -175,8 +175,10 @@ function test_jld2_time_averaging_of_horizontal_averages(model) FT = eltype(model.grid) - @test wu == zero(FT) - @test wT == zero(FT) + # Note: w is not identically 0 because T = 4 introduces a buoyancy term that is + # subsequently cancelled by a large scale pressure field. + @test abs(wu) < eps(FT) + @test abs(wT) < eps(FT) @test uv == FT(2) return nothing diff --git a/test/test_nonhydrostatic_regression.jl b/test/test_nonhydrostatic_regression.jl index f24e7540ac..44be406b4e 100644 --- a/test/test_nonhydrostatic_regression.jl +++ b/test/test_nonhydrostatic_regression.jl @@ -2,6 +2,7 @@ include("dependencies_for_runtests.jl") include("data_dependencies.jl") using Oceananigans.Grids: topology, XRegularLLG, YRegularLLG, ZRegularLLG +using Oceananigans.Fields: CenterField function get_fields_from_checkpoint(filename) file = jldopen(filename) diff --git a/test/test_time_stepping.jl b/test/test_time_stepping.jl index 4c6399e97e..9b07247248 100644 --- a/test/test_time_stepping.jl +++ b/test/test_time_stepping.jl @@ -87,27 +87,32 @@ function time_stepping_works_with_nonlinear_eos(arch, FT, EOS) return true # Test that no errors/crashes happen when time stepping. end +@inline add_ones(args...) = 1 + function run_first_AB2_time_step_tests(arch, FT) - add_ones(args...) = 1.0 # Weird grid size to catch https://github.com/CliMA/Oceananigans.jl/issues/780 grid = RectilinearGrid(arch, FT, size=(13, 17, 19), extent=(1, 2, 3)) - model = NonhydrostaticModel(grid=grid, forcing=(T=add_ones,), - buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) - time_step!(model, 1, euler=true) + model = NonhydrostaticModel(; grid, + forcing = (; T=add_ones), + buoyancy = SeawaterBuoyancy(), + tracers = (:T, :S)) - # Test that GT = 1, T = 1 after 1 time step and that AB2 actually reduced to forward Euler. + # Test that GT = 1 after model construction (note: this computes tendencies) @test all(interior(model.timestepper.Gⁿ.u) .≈ 0) @test all(interior(model.timestepper.Gⁿ.v) .≈ 0) @test all(interior(model.timestepper.Gⁿ.w) .≈ 0) - @test all(interior(model.timestepper.Gⁿ.T) .≈ 1.0) + @test all(interior(model.timestepper.Gⁿ.T) .≈ 1) @test all(interior(model.timestepper.Gⁿ.S) .≈ 0) + # Test that T = 1 after 1 time step and that AB2 actually reduced to forward Euler. + Δt = 1 + time_step!(model, Δt, euler=true) @test all(interior(model.velocities.u) .≈ 0) @test all(interior(model.velocities.v) .≈ 0) @test all(interior(model.velocities.w) .≈ 0) - @test all(interior(model.tracers.T) .≈ 1.0) + @test all(interior(model.tracers.T) .≈ 1) @test all(interior(model.tracers.S) .≈ 0) return nothing