From 6e39d3fcc098c69ac207cc21be759cf6bd3ec604 Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Thu, 15 Jul 2021 07:03:48 -0400 Subject: [PATCH] Major refactor of `FieldBoundaryConditions` (#1843) * Reworking FieldBoundaryConditions * Refactors FieldBoundaryConditions * Fixes many tests * More test fixes * Update docs/src/model_setup/boundary_conditions.md Co-authored-by: Navid C. Constantinou * Nuke AuxiliaryFieldBoundaryConditions; just use special FieldBoundaryConditions * No more x, y, z, right, left... * Banish DiffusivityBoundaryConditions * Banish CoordinateBoundaryConditions * Fixes boundary condition test with ContinuousBoundaryFunction * Adds better default bcs for VelocityFields and TracerFields * New syntax! * Iterate better over boundary conditions * Fixes bug in field_boundary_conditions * Fix TracerFields implementation * Were not ready for zero flux boundary conditions... * Fixes bottom-top mix up * Fixes broken ref to Nx in Poisson solver test * Default for ZFaceField was changed so now broadcasting tests need to be changed too * Fix regularize_boundary_condition for CubedSpheres * West vs east * Fix issues with TKEBasedVerticalDiffusivity * Fiddling with AMD boundary conditions * Adds regularize_boundary_conditions for regular bcs on ConformalCubedSphereGrid * Fixes bug in determining tracer diffusivity bcs with AMD * Adds location to FieldBoundaryConditions for broadcasting test * Properly inject CubedSphere exchange bcs during regularization * Adds a method for nested NamedTuple boundary conditions * Fixes VelocityFields constructor to use valid regularized boundary conditions * Updates to new syntax for stratified couette flow * Elide tendency fields when velocities are prescribed * Import regularize_boundary_conditions into field_tuples * Try 2 with poisson solver test refactor * Simplifies default calculation * Fixes ommission of tracer tendencies when using prescribed velocities * Import regularize_field_boundary_conditions into poisson solver tests * Eta is prognostic, not w * Adds missing comma in show methods * Fix doctests * Fix doctests * Fix doctests * Fix docstring in Field constructor * Use FluxBoundaryCondition * Where bc? * Ancient issues with library.md * Comments out library to see what happens * Deletes misleading and impertinent comments * Bump to 0.59.0 * Comments out appendix in docs * Cleans up BoundaryConditions code + docstrings * Docs despration * update BoundaryConditions pages * can't document call syntax; see https://github.com/JuliaDocs/Documenter.jl/issues/228 * back to including everything * improves some docstrings * Updates validation experiments * Fixes typos in TKEBasedVerticalMixing validation Co-authored-by: Navid C. Constantinou --- Project.toml | 2 +- docs/make.jl | 5 +- docs/make_example.jl | 6 +- docs/src/appendix/benchmarks.md | 2 +- docs/src/appendix/fractional_step.md | 9 - docs/src/appendix/library.md | 13 +- docs/src/model_setup/boundary_conditions.md | 87 +++-- docs/src/model_setup/tracers.md | 5 +- docs/src/simulation_tips.md | 23 +- examples/convecting_plankton.jl | 4 +- examples/eady_turbulence.jl | 4 +- examples/langmuir_turbulence.jl | 6 +- examples/ocean_wind_mixing_and_convection.jl | 9 +- src/BoundaryConditions/BoundaryConditions.jl | 13 +- src/BoundaryConditions/apply_flux_bcs.jl | 10 +- src/BoundaryConditions/boundary_condition.jl | 9 +- .../continuous_boundary_function.jl | 73 ++-- .../coordinate_boundary_conditions.jl | 12 - .../discrete_boundary_function.jl | 2 +- .../field_boundary_conditions.jl | 326 +++++++----------- src/BoundaryConditions/fill_halo_regions.jl | 2 +- .../show_boundary_conditions.jl | 41 ++- src/CubedSpheres/CubedSpheres.jl | 18 +- src/CubedSpheres/cubed_sphere_exchange_bcs.jl | 15 +- src/Distributed/halo_communication_bcs.jl | 21 +- src/Fields/computed_field.jl | 22 +- src/Fields/field.jl | 94 ++--- src/Fields/field_tuples.jl | 70 ++-- src/Fields/kernel_computed_field.jl | 4 +- .../HydrostaticFreeSurfaceModels.jl | 3 +- .../hydrostatic_free_surface_field_tuples.jl | 19 + .../hydrostatic_free_surface_model.jl | 27 +- ...ydrostatic_free_surface_tendency_fields.jl | 15 - ...ydrostatic_free_surface_velocity_fields.jl | 18 - .../prescribed_hydrostatic_velocity_fields.jl | 5 + .../vertical_vorticity_field.jl | 10 +- .../incompressible_model.jl | 57 +-- .../ShallowWaterModels/shallow_water_model.jl | 30 +- src/Oceananigans.jl | 4 +- src/OutputWriters/output_writer_utils.jl | 16 +- .../anisotropic_minimum_dissipation.jl | 54 +-- .../leith_enstrophy_diffusivity.jl | 7 +- .../smagorinsky_lilly.jl | 7 +- .../tke_based_vertical_diffusivity.jl | 50 ++- ...n_large_eddy_simulation_regression_test.jl | 7 +- .../rayleigh_benard_regression_test.jl | 4 +- test/test_boundary_conditions.jl | 286 ++++++++------- test/test_boundary_conditions_integration.jl | 66 ++-- test/test_broadcasting.jl | 5 +- test/test_checkpointer.jl | 44 +-- test/test_distributed_poisson_solvers.jl | 14 +- test/test_dynamics.jl | 4 +- test/test_operators.jl | 3 +- test/test_output_readers.jl | 6 +- test/test_poisson_solvers.jl | 27 +- test/test_time_stepping.jl | 2 +- validation/barotropic_gyre/barotropic_gyre.jl | 17 +- validation/convergence_tests/Manifest.toml | 172 +++++---- .../run_forced_fixed_slip.jl | 2 + .../convergence_tests/run_forced_free_slip.jl | 2 + .../convergence_tests/run_taylor_green.jl | 2 + .../src/ForcedFlowFixedSlip.jl | 43 ++- .../lid_driven_cavity/lid_driven_cavity.jl | 12 +- .../stratified_couette_flow.jl | 12 +- .../convective_adjustment_free_convection.jl | 2 +- .../tke_based_free_convection.jl | 9 +- .../tke_based_wind_mixing.jl | 8 +- 67 files changed, 924 insertions(+), 1054 deletions(-) delete mode 100644 src/BoundaryConditions/coordinate_boundary_conditions.jl create mode 100644 src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_field_tuples.jl delete mode 100644 src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_fields.jl delete mode 100644 src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_velocity_fields.jl diff --git a/Project.toml b/Project.toml index b168d9fae4..706dc03232 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Oceananigans" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" -version = "0.58.9" +version = "0.59.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/docs/make.jl b/docs/make.jl index ed74fe1555..540d4a6ab6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -116,7 +116,6 @@ numerical_pages = [ "Large eddy simulation" => "numerical_implementation/large_eddy_simulation.md" ] - appendix_pages = [ "Staggered grid" => "appendix/staggered_grid.md", "Fractional step method" => "appendix/fractional_step.md", @@ -161,7 +160,7 @@ makedocs(bib, doctest = true, strict = true, clean = true, - checkdocs = :none # Should fix our docstring so we can use checkdocs=:exports with strict=true. + checkdocs = :none # Should fix our docstring so we can use checkdocs=:exports with strict=true. ) deploydocs( @@ -173,4 +172,4 @@ deploydocs( @info "Cleaning up temporary .jld2 and .nc files created by doctests..." for file in vcat(glob("docs/*.jld2"), glob("docs/*.nc")) rm(file) -end \ No newline at end of file +end diff --git a/docs/make_example.jl b/docs/make_example.jl index 93b2e0b350..47f942d4ea 100644 --- a/docs/make_example.jl +++ b/docs/make_example.jl @@ -53,9 +53,9 @@ end ##### example_pages = [ - #"Internal wave" => "generated/internal_wave.md", - #"Geostrophic adjustment" => "generated/geostrophic_adjustment.md", - "Shallow water Bickley jet" => "generated/shallow_water_Bickley_jet.md" + #"Internal wave" => "generated/internal_wave.md", + #"Geostrophic adjustment" => "generated/geostrophic_adjustment.md", + "Shallow water Bickley jet" => "generated/shallow_water_Bickley_jet.md" ] pages = [ diff --git a/docs/src/appendix/benchmarks.md b/docs/src/appendix/benchmarks.md index 62aa00c51e..9dee4151a0 100644 --- a/docs/src/appendix/benchmarks.md +++ b/docs/src/appendix/benchmarks.md @@ -10,7 +10,7 @@ format the benchmark results. This is a benchmark of a simple "static ocean" configuration. The time stepping and Poisson solver still takes the same amount of time whether the ocean is static or active, so it is -quite indicative of actual performance. It tests the performance of a bare-bones +indicative of actual performance. It tests the performance of a bare-bones horizontally-periodic model with `topology = (Periodic, Periodic, Bounded)`. ``` diff --git a/docs/src/appendix/fractional_step.md b/docs/src/appendix/fractional_step.md index 8740d10370..f378d1d8a1 100644 --- a/docs/src/appendix/fractional_step.md +++ b/docs/src/appendix/fractional_step.md @@ -35,12 +35,3 @@ methods, which is that "*while the velocity can be reliably computed to second-o in time and space, the pressure is typically only first-order accurate in the ``L_\infty``-norm.*" The numerical boundary conditions must be carefully accounted for to ensure the second-order accuracy promised by the fractional step methods. - -We are currently investigating whether our projection method is indeed second-order accurate -in both velocity and pressure. However, it may not matter too much for simulating high Reynolds -number geophysical fluids as [Brown01](@cite) conclude that "*Quite often, semi-implicit projection -methods are applied to problems in which the viscosity is small. Since the predicted first-order -errors in the pressure are scaled by ``\nu``, it is not clear whether the improved pressure-update -formula is beneficial in such situations. [...] Finally, in some applications of projection methods, -second-order accuracy in the pressure may not be relevant or in some cases even possible due -to the treatment of other terms in the equations.*" diff --git a/docs/src/appendix/library.md b/docs/src/appendix/library.md index 6ac92be0bd..60a943a7c4 100644 --- a/docs/src/appendix/library.md +++ b/docs/src/appendix/library.md @@ -36,15 +36,17 @@ Pages = [ "BoundaryConditions/BoundaryConditions.jl", "BoundaryConditions/boundary_condition_classifications.jl", "BoundaryConditions/boundary_condition.jl", - "BoundaryConditions/coordinate_boundary_conditions.jl", + "BoundaryConditions/discrete_boundary_function.jl", + "BoundaryConditions/continuous_boundary_function.jl", "BoundaryConditions/field_boundary_conditions.jl", - "BoundaryConditions/boundary_function.jl", - "BoundaryConditions/parameterized_boundary_condition.jl", "BoundaryConditions/show_boundary_conditions.jl", "BoundaryConditions/fill_halo_regions.jl", + "BoundaryConditions/fill_halo_regions_value_gradient.jl", + "BoundaryConditions/fill_halo_regions_open.jl", + "BoundaryConditions/fill_halo_regions_periodic.jl", + "BoundaryConditions/fill_halo_regions_flux.jl", + "BoundaryConditions/fill_halo_regions_nothing.jl", "BoundaryConditions/apply_flux_bcs.jl", - "BoundaryConditions/apply_value_gradient_bcs.jl", - "BoundaryConditions/apply_normal_flow_bcs.jl" ] ``` @@ -107,7 +109,6 @@ Pages = [ "Fields/field.jl", "Fields/set!.jl", "Fields/show_fields.jl", - "Fields/field_utils.jl" ] ``` diff --git a/docs/src/model_setup/boundary_conditions.md b/docs/src/model_setup/boundary_conditions.md index 1a8d680eff..4fe7dffa00 100644 --- a/docs/src/model_setup/boundary_conditions.md +++ b/docs/src/model_setup/boundary_conditions.md @@ -45,10 +45,26 @@ Here's a short list of useful tips for defining and prescribing boundary conditi If there is a known `diffusivity`, you can express `FluxBoundaryCondition(flux)` using `GradientBoundaryCondition(-flux / diffusivity)` (aka "Neumann" boundary condition). But when `diffusivity` is not known or is variable (as for large eddy simulation, for example), it's convenient and more straightforward to apply `FluxBoundaryCondition`. +## Default boundary conditions + +By default, periodic boundary conditions are applied on all fields along periodic dimensions. Otherwise tracers +get no-flux boundary conditions and velocities get free-slip and no normal flow boundary conditions. + +## Boundary condition structures + +Oceananigans uses a hierarchical structure to express boundary conditions: + +1. Each boundary has one [`BoundaryCondition`](@ref) +2. Each field has seven [`BoundaryCondition`](@ref) (`west`, `east`, `south`, `north`, `bottom`, `top` and + and an additional experimental condition for `immersed` boundaries) +3. A set of `FieldBoundaryConditions`, up to one for each field, are grouped into a `NamedTuple` and passed + to the model constructor. + ## Specifying boundary conditions for a model Boundary conditions are defined at model construction time by passing a `NamedTuple` of `FieldBoundaryConditions` -specifying non-default boundary conditions for velocities (``u``, ``v``, ``w``) and tracers. +specifying non-default boundary conditions for fields such as velocities and tracers. + Fields for which boundary conditions are not specified are assigned a default boundary conditions. A few illustrations are provided below. See the examples for @@ -59,7 +75,6 @@ further illustrations of boundary condition specification. ```@meta DocTestSetup = quote using Oceananigans - using Oceananigans.BoundaryConditions using Random Random.seed!(1234) @@ -73,7 +88,7 @@ In this section we illustrate usage of the different [`BoundaryCondition`](@ref) ```jldoctest julia> constant_T_bc = ValueBoundaryCondition(20.0) -BoundaryCondition: type=Value, condition=20.0 +BoundaryCondition: classification=Value, condition=20.0 ``` A constant [`Value`](@ref) boundary condition can be used to specify constant tracer (such as temperature), @@ -90,7 +105,7 @@ julia> ρ₀ = 1027; # Reference density [kg/m³] julia> τₓ = 0.08; # Wind stress [N/m²] julia> wind_stress_bc = FluxBoundaryCondition(-τₓ/ρ₀) -BoundaryCondition: type=Flux, condition=-7.789678675754625e-5 +BoundaryCondition: classification=Flux, condition=-7.789678675754625e-5 ``` A constant [`Flux`](@ref) boundary condition can be imposed on tracers and tangential velocity components @@ -115,7 +130,7 @@ Boundary conditions may be specified by functions, julia> @inline surface_flux(x, y, t) = cos(2π * x) * cos(t); julia> top_tracer_bc = FluxBoundaryCondition(surface_flux) -BoundaryCondition: type=Flux, condition=surface_flux(x, y, t) in Main at none:1 +BoundaryCondition: classification=Flux, condition=surface_flux(x, y, t) in Main at none:1 ``` !!! info "Boundary condition functions" @@ -137,8 +152,8 @@ Boundary condition functions may be 'parameterized', ```jldoctest julia> @inline wind_stress(x, y, t, p) = - p.τ * cos(p.k * x) * cos(p.ω * t); # function with parameters -julia> top_u_bc = BoundaryCondition(Flux, wind_stress, parameters=(k=4π, ω=3.0, τ=1e-4)) -BoundaryCondition: type=Flux, condition=wind_stress(x, y, t, p) in Main at none:1 +julia> top_u_bc = FluxBoundaryCondition(wind_stress, parameters=(k=4π, ω=3.0, τ=1e-4)) +BoundaryCondition: classification=Flux, condition=wind_stress(x, y, t, p) in Main at none:1 ``` !!! info "Boundary condition functions with parameters" @@ -157,7 +172,7 @@ julia> @inline linear_drag(x, y, t, u) = - 0.2 * u linear_drag (generic function with 1 method) julia> u_bottom_bc = FluxBoundaryCondition(linear_drag, field_dependencies=:u) -BoundaryCondition: type=Flux, condition=linear_drag(x, y, t, u) in Main at none:1 +BoundaryCondition: classification=Flux, condition=linear_drag(x, y, t, u) in Main at none:1 ``` `field_dependencies` specifies the name of the dependent fields either with a `Symbol` or `Tuple` of `Symbol`s. @@ -171,7 +186,7 @@ julia> @inline quadratic_drag(x, y, t, u, v, drag_coeff) = - drag_coeff * u * sq quadratic_drag (generic function with 1 method) julia> u_bottom_bc = FluxBoundaryCondition(quadratic_drag, field_dependencies=(:u, :v), parameters=1e-3) -BoundaryCondition: type=Flux, condition=quadratic_drag(x, y, t, u, v, drag_coeff) in Main at none:1 +BoundaryCondition: classification=Flux, condition=quadratic_drag(x, y, t, u, v, drag_coeff) in Main at none:1 ``` Put differently, `ξ, η, t` come first in the function signature, followed by field dependencies, @@ -189,7 +204,7 @@ using the `discrete_form`. For example: u_bottom_bc = FluxBoundaryCondition(filtered_drag, discrete_form=true) # output -BoundaryCondition: type=Flux, condition=filtered_drag(i, j, grid, clock, model_fields) in Main at none:1 +BoundaryCondition: classification=Flux, condition=filtered_drag(i, j, grid, clock, model_fields) in Main at none:1 ``` !!! info "The 'discrete form' for boundary condition functions" @@ -212,8 +227,8 @@ julia> Cd = 0.2; # drag coefficient julia> @inline linear_drag(i, j, grid, clock, model_fields, Cd) = @inbounds - Cd * model_fields.u[i, j, 1]; -julia> u_bottom_bc = BoundaryCondition(Flux, linear_drag, discrete_form=true, parameters=Cd) -BoundaryCondition: type=Flux, condition=linear_drag(i, j, grid, clock, model_fields, Cd) in Main at none:1 +julia> u_bottom_bc = FluxBoundaryCondition(linear_drag, discrete_form=true, parameters=Cd) +BoundaryCondition: classification=Flux, condition=linear_drag(i, j, grid, clock, model_fields, Cd) in Main at none:1 ``` !!! info "Inlining and avoiding bounds-checking in boundary condition functions" @@ -231,36 +246,34 @@ julia> Nx = Ny = 16; # Number of grid points. julia> Q = randn(Nx, Ny); # temperature flux julia> white_noise_T_bc = FluxBoundaryCondition(Q) -BoundaryCondition: type=Flux, condition=16×16 Matrix{Float64} +BoundaryCondition: classification=Flux, condition=16×16 Matrix{Float64} ``` When running on the GPU, `Q` must be converted to a `CuArray`. ## Building boundary conditions on a field -To create, for example, a set of horizontally-periodic field boundary conditions, write +To create a set of [`FieldBoundaryConditions`](@ref) for a temperature field, +we write ```jldoctest -julia> topology = (Periodic, Periodic, Bounded); - -julia> grid = RegularRectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=topology); - -julia> T_bcs = TracerBoundaryConditions(grid, top = ValueBoundaryCondition(20), - bottom = GradientBoundaryCondition(0.01)) -Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions -├── x: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}} -├── y: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}} -└── z: CoordinateBoundaryConditions{BoundaryCondition{Gradient, Float64}, BoundaryCondition{Value, Int64}} +julia> T_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(20), + bottom = GradientBoundaryCondition(0.01)) +Oceananigans.FieldBoundaryConditions, with boundary conditions +├── west: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition +├── east: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition +├── south: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition +├── north: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition +├── bottom: BoundaryCondition{Gradient, Float64} +├── top: BoundaryCondition{Value, Int64} +└── immersed: BoundaryCondition{Flux, Nothing} ``` -`T_bcs` is a [`FieldBoundaryConditions`](@ref) object for temperature `T` appropriate -for horizontally periodic grid topologies. -The default `Periodic` boundary conditions in ``x`` and ``y`` are inferred from the `topology` of `grid`. +If the grid is, e.g., horizontally-periodic, then each horizontal `DefaultPrognosticFieldBoundaryCondition` +is converted to `PeriodicBoundaryCondition` inside the model's constructor, before assigning the +boundary conditions to temperature `T`. -For ``u``, ``v``, and ``w``, use the -`UVelocityBoundaryConditions` -`VVelocityBoundaryConditions`, and -`WVelocityBoundaryConditions` constructors, respectively. +In general, boundary condition defaults are inferred from the field location and `topology(grid)`. ## Specifying model boundary conditions @@ -274,11 +287,11 @@ julia> topology = (Periodic, Periodic, Bounded); julia> grid = RegularRectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=topology); -julia> u_bcs = UVelocityBoundaryConditions(grid, top = ValueBoundaryCondition(+0.1), - bottom = ValueBoundaryCondition(-0.1)); +julia> u_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(+0.1), + bottom = ValueBoundaryCondition(-0.1)); -julia> T_bcs = TracerBoundaryConditions(grid, top = ValueBoundaryCondition(20), - bottom = GradientBoundaryCondition(0.01)); +julia> T_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(20), + bottom = GradientBoundaryCondition(0.01)); julia> model = IncompressibleModel(grid=grid, boundary_conditions=(u=u_bcs, T=T_bcs)) IncompressibleModel{CPU, Float64}(time = 0 seconds, iteration = 0) @@ -292,13 +305,13 @@ julia> model.velocities.u Field located at (Face, Center, Center) ├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (16, 16, 16) ├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=16, Ny=16, Nz=16) -└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=Value, top=Value) +└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=Value, top=Value, immersed=ZeroFlux julia> model.tracers.T Field located at (Center, Center, Center) ├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (16, 16, 16) ├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=16, Ny=16, Nz=16) -└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=Gradient, top=Value) +└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=Gradient, top=Value, immersed=ZeroFlux ``` Notice that the specified non-default boundary conditions have been applied at diff --git a/docs/src/model_setup/tracers.md b/docs/src/model_setup/tracers.md index 8abc1a3792..b9174d609e 100644 --- a/docs/src/model_setup/tracers.md +++ b/docs/src/model_setup/tracers.md @@ -28,14 +28,13 @@ julia> model.tracers.T Field located at (Center, Center, Center) ├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (64, 64, 64) ├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=64, Ny=64, Nz=64) -└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux) +└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux julia> model.tracers.S Field located at (Center, Center, Center) ├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (64, 64, 64) ├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=64, Ny=64, Nz=64) -└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux) - +└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux ``` Any number of arbitrary tracers can be appended to this list and passed to a model constructor. For example, to evolve diff --git a/docs/src/simulation_tips.md b/docs/src/simulation_tips.md index 65c1863110..410bae11dc 100644 --- a/docs/src/simulation_tips.md +++ b/docs/src/simulation_tips.md @@ -73,19 +73,26 @@ in GPU simulations by following a few simple principles. ### Variables that need to be used in GPU computations need to be defined as constants Any global variable that needs to be accessed by the GPU needs to be a constant or the simulation -will crash. This includes any variables used in forcing functions and boundary conditions. For -example, if you define a boundary condition like the example below and run your simulation on a GPU -you'll get an error. +will crash. This includes any variables that are referenced as global variables in functions +used for forcing of boundary conditions. For example, ```julia -dTdz = 0.01 # K m⁻¹ -T_bcs = TracerBoundaryConditions(grid, - bottom = GradientBoundaryCondition(dTdz)) +T₀ = 20 # ᵒC +surface_temperature(x, y, t) = T₀ * sin(2π / 86400 * t) +T_bcs = FieldBoundaryConditions(bottom = GradientBoundaryCondition(surface_temperature)) ``` -However, if you define `dTdz` as a constant by replacing the first line with `const dTdz = 0.01`, -then (provided everything else is done properly) your run will be successful. +will throw an error if run on the GPU (and will run more slowly than it should on the CPU). +Replacing the first line above with +```julia +const T₀ = 20 # ᵒC +``` + +fixes the issue by indicating to the compiler that `T₀` will not change. + +Note that the _literal_ `2π / 86400` is not an issue -- it's only the +_variable_ `T₀` that must be declared `const`. ### Complex diagnostics using `ComputedField`s may not work on GPUs diff --git a/examples/convecting_plankton.jl b/examples/convecting_plankton.jl index db99d891ef..8998d8ba18 100644 --- a/examples/convecting_plankton.jl +++ b/examples/convecting_plankton.jl @@ -12,7 +12,7 @@ # The phytoplankton in our model are advected, diffuse, grow, and die according to # # ```math -# ∂_t P + \boldsymbol{u ⋅ ∇} P - κ ∇²P = (μ₀ \exp(z / λ) - m) \, P \, , +# ∂_t P + \boldsymbol{v ⋅ ∇} P - κ ∇²P = (μ₀ \exp(z / λ) - m) \, P \, , # ``` # # where ``\boldsymbol{v}`` is the turbulent velocity field, ``κ`` is an isotropic diffusivity, @@ -92,7 +92,7 @@ buoyancy_gradient_bc = GradientBoundaryCondition(N²) # In summary, the buoyancy boundary conditions impose a destabilizing flux # at the top and a stable buoyancy gradient at the bottom: -buoyancy_bcs = TracerBoundaryConditions(grid, top = buoyancy_flux_bc, bottom = buoyancy_gradient_bc) +buoyancy_bcs = FieldBoundaryConditions(top = buoyancy_flux_bc, bottom = buoyancy_gradient_bc) # ## Phytoplankton dynamics: light-dependent growth and uniform mortality # diff --git a/examples/eady_turbulence.jl b/examples/eady_turbulence.jl index 95605a5bfe..d4eef2d63c 100644 --- a/examples/eady_turbulence.jl +++ b/examples/eady_turbulence.jl @@ -167,8 +167,8 @@ cᴰ = 1e-4 # quadratic drag coefficient drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=cᴰ) drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=cᴰ) -u_bcs = UVelocityBoundaryConditions(grid, bottom = drag_bc_u) -v_bcs = VVelocityBoundaryConditions(grid, bottom = drag_bc_v) +u_bcs = FieldBoundaryConditions(bottom = drag_bc_u) +v_bcs = FieldBoundaryConditions(bottom = drag_bc_v) # ## Turbulence closures # diff --git a/examples/langmuir_turbulence.jl b/examples/langmuir_turbulence.jl index 1a90957cc1..b899f32a06 100644 --- a/examples/langmuir_turbulence.jl +++ b/examples/langmuir_turbulence.jl @@ -95,7 +95,7 @@ uˢ(z) = Uˢ * exp(z / vertical_scale) Qᵘ = -3.72e-5 # m² s⁻², surface kinematic momentum flux -u_boundary_conditions = UVelocityBoundaryConditions(grid, top = FluxBoundaryCondition(Qᵘ)) +u_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) # McWilliams et al. (1997) impose a linear buoyancy gradient `N²` at the bottom # along with a weak, destabilizing flux of buoyancy at the surface to faciliate @@ -104,8 +104,8 @@ u_boundary_conditions = UVelocityBoundaryConditions(grid, top = FluxBoundaryCond Qᵇ = 2.307e-9 # m³ s⁻², surface buoyancy flux N² = 1.936e-5 # s⁻², initial and bottom buoyancy gradient -b_boundary_conditions = TracerBoundaryConditions(grid, top = FluxBoundaryCondition(Qᵇ), - bottom = GradientBoundaryCondition(N²)) +b_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), + bottom = GradientBoundaryCondition(N²)) # !!! info "The flux convention in Oceananigans" # Note that Oceananigans uses "positive upward" conventions for all fluxes. In consequence, diff --git a/examples/ocean_wind_mixing_and_convection.jl b/examples/ocean_wind_mixing_and_convection.jl index b759e7f89a..546a6c230f 100644 --- a/examples/ocean_wind_mixing_and_convection.jl +++ b/examples/ocean_wind_mixing_and_convection.jl @@ -90,9 +90,8 @@ Qᵀ = Qʰ / (ρₒ * cᴾ) # K m s⁻¹, surface temperature flux dTdz = 0.01 # K m⁻¹ -T_bcs = TracerBoundaryConditions(grid, - top = FluxBoundaryCondition(Qᵀ), - bottom = GradientBoundaryCondition(dTdz)) +T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵀ), + bottom = GradientBoundaryCondition(dTdz)) # Note that a positive temperature flux at the surface of the ocean # implies cooling. This is because a positive temperature flux implies @@ -111,7 +110,7 @@ Qᵘ = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻² # The boundary conditions on `u` are thus -u_bcs = UVelocityBoundaryConditions(grid, top = FluxBoundaryCondition(Qᵘ)) +u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) # For salinity, `S`, we impose an evaporative flux of the form @@ -130,7 +129,7 @@ evaporation_bc = FluxBoundaryCondition(Qˢ, field_dependencies=:S, parameters=ev # The full salinity boundary conditions are -S_bcs = TracerBoundaryConditions(grid, top=evaporation_bc) +S_bcs = FieldBoundaryConditions(top=evaporation_bc) # ## Model instantiation # diff --git a/src/BoundaryConditions/BoundaryConditions.jl b/src/BoundaryConditions/BoundaryConditions.jl index 0ec1a9eb51..282243053b 100644 --- a/src/BoundaryConditions/BoundaryConditions.jl +++ b/src/BoundaryConditions/BoundaryConditions.jl @@ -2,21 +2,11 @@ module BoundaryConditions export BCType, Flux, Gradient, Value, Open, - BoundaryCondition, getbc, setbc!, - PeriodicBoundaryCondition, OpenBoundaryCondition, NoFluxBoundaryCondition, FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, - - CoordinateBoundaryConditions, - - FieldBoundaryConditions, UVelocityBoundaryConditions, VVelocityBoundaryConditions, - WVelocityBoundaryConditions, TracerBoundaryConditions, PressureBoundaryConditions, - - DiffusivityBoundaryConditions, - + FieldBoundaryConditions, apply_x_bcs!, apply_y_bcs!, apply_z_bcs!, - fill_halo_regions! using CUDA @@ -31,7 +21,6 @@ include("boundary_condition_classifications.jl") include("boundary_condition.jl") include("discrete_boundary_function.jl") include("continuous_boundary_function.jl") -include("coordinate_boundary_conditions.jl") include("field_boundary_conditions.jl") include("show_boundary_conditions.jl") diff --git a/src/BoundaryConditions/apply_flux_bcs.jl b/src/BoundaryConditions/apply_flux_bcs.jl index c655c1953e..16047f519f 100644 --- a/src/BoundaryConditions/apply_flux_bcs.jl +++ b/src/BoundaryConditions/apply_flux_bcs.jl @@ -12,7 +12,7 @@ Apply flux boundary conditions to a field `c` by adding the associated flux dive the source term `Gc` at the left and right. """ apply_x_bcs!(Gc, c, arch, dep, args...) = launch!(arch, c.grid, :yz, _apply_x_bcs!, Gc, instantiated_location(Gc), c.grid, - c.boundary_conditions.x.left, c.boundary_conditions.x.right, args..., + c.boundary_conditions.west, c.boundary_conditions.east, args..., dependencies=dep) """ @@ -22,7 +22,7 @@ Apply flux boundary conditions to a field `c` by adding the associated flux dive the source term `Gc` at the left and right. """ apply_y_bcs!(Gc, c, arch, dep, args...) = launch!(arch, c.grid, :xz, _apply_y_bcs!, Gc, instantiated_location(Gc), c.grid, - c.boundary_conditions.y.left, c.boundary_conditions.y.right, args..., + c.boundary_conditions.south, c.boundary_conditions.north, args..., dependencies=dep) """ @@ -32,7 +32,7 @@ Apply flux boundary conditions to a field `c` by adding the associated flux dive the source term `Gc` at the top and bottom. """ apply_z_bcs!(Gc, c, arch, dep, args...) = launch!(arch, c.grid, :xy, _apply_z_bcs!, Gc, instantiated_location(Gc), c.grid, - c.boundary_conditions.z.left, c.boundary_conditions.z.right, args..., + c.boundary_conditions.bottom, c.boundary_conditions.top, args..., dependencies=dep) apply_x_bcs!(::Nothing, args...) = nothing @@ -92,7 +92,7 @@ end @inline flip(::Face) = Center() """ - apply_x_west_bc!(Gc, west_flux::BC{<:Flux}, j, k, grid, args...) + apply_x_west_bc!(Gc, loc, west_flux::BC{<:Flux}, j, k, grid, args...) Add the flux divergence associated with a west flux boundary condition on `c`. Note that because @@ -125,7 +125,7 @@ end end """ - apply_x_east_bc!(Gc, top_flux::BC{<:Flux}, j, k, grid, args...) + apply_x_east_bc!(Gc, loc, east_flux::BC{<:Flux}, j, k, grid, args...) Add the part of flux divergence associated with a east boundary condition on `c`. Note that because diff --git a/src/BoundaryConditions/boundary_condition.jl b/src/BoundaryConditions/boundary_condition.jl index b970ecf052..b224f4ad8f 100644 --- a/src/BoundaryConditions/boundary_condition.jl +++ b/src/BoundaryConditions/boundary_condition.jl @@ -11,7 +11,7 @@ struct BoundaryCondition{C<:AbstractBoundaryConditionClassification, T} end """ - BoundaryCondition(BC, condition) + BoundaryCondition(Classification::DataType, condition) Construct a boundary condition of type `BC` with a number or array as a `condition`. @@ -20,9 +20,12 @@ Boundary condition types include `Periodic`, `Flux`, `Value`, `Gradient`, and `O BoundaryCondition(Classification::DataType, condition) = BoundaryCondition(Classification(), condition) """ - BoundaryCondition(BC, condition::Function; parameters=nothing, discrete_form=false) + BoundaryCondition(Classification::DataType, condition::Function; + parameters = nothing, + discrete_form = false, + field_dependencies=()) -Construct a boundary condition of type `BC` with a function boundary `condition`. +Construct a boundary condition of type `Classification` with a function boundary `condition`. By default, the function boudnary `condition` is assumed to have the 'continuous form' `condition(ξ, η, t)`, where `t` is time and `ξ` and `η` vary along the boundary. diff --git a/src/BoundaryConditions/continuous_boundary_function.jl b/src/BoundaryConditions/continuous_boundary_function.jl index 8df086b7b1..961ac4c993 100644 --- a/src/BoundaryConditions/continuous_boundary_function.jl +++ b/src/BoundaryConditions/continuous_boundary_function.jl @@ -4,7 +4,7 @@ using Oceananigans.Utils: tupleit, user_function_arguments import Oceananigans: location """ - ContinuousBoundaryFunction{X, Y, Z, I, F, P, D, N, ℑ} <: Function + struct ContinuousBoundaryFunction{X, Y, Z, I, F, P, D, N, ℑ} <: Function A wrapper for the user-defined boundary condition function `func` at location `X, Y, Z`. `I` denotes the boundary-normal index (`I=1` at western boundaries, @@ -21,49 +21,39 @@ struct ContinuousBoundaryFunction{X, Y, Z, I, F, P, D, N, ℑ} <: Function field_dependencies_interp :: ℑ """ Returns a location-less wrapper for `func`, `parameters`, and `field_dependencies`.""" - function ContinuousBoundaryFunction(func, parameters, field_dependencies) - + function ContinuousBoundaryFunction(func::F, parameters::P, field_dependencies) where {F, P} field_dependencies = tupleit(field_dependencies) - - return new{Nothing, Nothing, Nothing, Nothing, - typeof(func), typeof(parameters), - typeof(field_dependencies), Nothing, Nothing}(func, parameters, field_dependencies, nothing, nothing) + D = typeof(field_dependencies) + return new{Nothing, Nothing, Nothing, Nothing, F, P, D, Nothing, Nothing}(func, parameters, field_dependencies, nothing, nothing) end - function ContinuousBoundaryFunction{X, Y, Z, I}(func, parameters, field_dependencies, - field_dependencies_indices, field_dependencies_interp) where {X, Y, Z, I} - return new{X, Y, Z, I, - typeof(func), - typeof(parameters), - typeof(field_dependencies), - typeof(field_dependencies_indices), - typeof(field_dependencies_interp)}(func, - parameters, - field_dependencies, - field_dependencies_indices, - field_dependencies_interp) + function ContinuousBoundaryFunction{X, Y, Z, I}(func::F, + parameters::P, + field_dependencies::D, + field_dependencies_indices::N, + field_dependencies_interp::ℑ) where {X, Y, Z, I, F, P, D, ℑ, N} + + return new{X, Y, Z, I, F, P, D, N, ℑ}(func, parameters, field_dependencies, field_dependencies_indices, field_dependencies_interp) end end -location(bc::ContinuousBoundaryFunction{X, Y, Z}) where {X, Y, Z} = X, Y, Z +location(::ContinuousBoundaryFunction{X, Y, Z}) where {X, Y, Z} = X, Y, Z ##### ##### "Regularization" for IncompressibleModel setup ##### -regularize_boundary_condition(bc, X, Y, Z, I, model_field_names) = bc # fallback - """ regularize_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoundaryFunction}, - X, Y, Z, I, model_field_names) where C + topo, loc, dim, I, prognostic_field_names) where C -Regularizes `bc.condition` for location `X, Y, Z`, boundary index `I`, and `model_field_names`, +Regularizes `bc.condition` for location `loc`, boundary index `I`, and `prognostic_field_names`, returning `BoundaryCondition(C, regularized_condition)`. The regularization of `bc.condition::ContinuousBoundaryFunction` requries -1. Setting the boundary location to `X, Y, Z`. - The boundary-normal direction is tagged with `Nothing` location. +1. Setting the boundary location to `LX, LY, LZ`. + The location in the boundary-normal direction is `Nothing`. 2. Setting the boundary-normal index `I` for indexing into `field_dependencies`. `I` is either `1` (for left boundaries) or @@ -76,12 +66,16 @@ The regularization of `bc.condition::ContinuousBoundaryFunction` requries of the boundary. """ function regularize_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoundaryFunction}, - LX, LY, LZ, I, model_field_names) where C + topo, loc, dim, I, prognostic_field_names) where C + boundary_func = bc.condition + # Set boundary-normal location to Nothing: + LX, LY, LZ = Tuple(i == dim ? Nothing : loc[i] for i = 1:3) + indices, interps = index_and_interp_dependencies(LX, LY, LZ, boundary_func.field_dependencies, - model_field_names) + prognostic_field_names) regularized_boundary_func = ContinuousBoundaryFunction{LX, LY, LZ, I}(boundary_func.func, boundary_func.parameters, @@ -91,31 +85,42 @@ function regularize_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoun return BoundaryCondition(C, regularized_boundary_func) end +@inline boundary_index(i) = ifelse(i == 1, 1, i + 1) # convert near-boundary Center() index into boundary Face() index + ##### ##### Kernel functions ##### +# Return ContinuousBoundaryFunction on east or west boundaries. @inline function (bc::ContinuousBoundaryFunction{Nothing, LY, LZ, i})(j, k, grid, clock, model_fields) where {LY, LZ, i} args = user_function_arguments(i, j, k, grid, model_fields, bc.parameters, bc) - return bc.func(ynode(nothing, LY(), LZ(), i, j, k, grid), - znode(nothing, LY(), LZ(), i, j, k, grid), + i′ = boundary_index(i) + + return bc.func(ynode(Face(), LY(), LZ(), i′, j, k, grid), + znode(Face(), LY(), LZ(), i′, j, k, grid), clock.time, args...) end +# Return ContinuousBoundaryFunction on south or north boundaries. @inline function (bc::ContinuousBoundaryFunction{LX, Nothing, LZ, j})(i, k, grid, clock, model_fields) where {LX, LZ, j} args = user_function_arguments(i, j, k, grid, model_fields, bc.parameters, bc) - return bc.func(xnode(LX(), nothing, LZ(), i, j, k, grid), - znode(LX(), nothing, LZ(), i, j, k, grid), + j′ = boundary_index(j) + + return bc.func(xnode(LX(), Face(), LZ(), i, j′, k, grid), + znode(LX(), Face(), LZ(), i, j′, k, grid), clock.time, args...) end +# Return ContinuousBoundaryFunction on bottom or top boundaries. @inline function (bc::ContinuousBoundaryFunction{LX, LY, Nothing, k})(i, j, grid, clock, model_fields) where {LX, LY, k} args = user_function_arguments(i, j, k, grid, model_fields, bc.parameters, bc) - return bc.func(xnode(LX(), LY(), nothing, i, j, k, grid), - ynode(LY(), LY(), nothing, i, j, k, grid), + k′ = boundary_index(k) + + return bc.func(xnode(LX(), LY(), Face(), i, j, k′, grid), + ynode(LY(), LY(), Face(), i, j, k′, grid), clock.time, args...) end diff --git a/src/BoundaryConditions/coordinate_boundary_conditions.jl b/src/BoundaryConditions/coordinate_boundary_conditions.jl deleted file mode 100644 index f8f7f4aa9a..0000000000 --- a/src/BoundaryConditions/coordinate_boundary_conditions.jl +++ /dev/null @@ -1,12 +0,0 @@ -""" - CoordinateBoundaryConditions(left, right) - -A set of two `BoundaryCondition`s to be applied along a coordinate x, y, or z. - -The `left` boundary condition is applied on the negative or lower side of the coordinate -while the `right` boundary condition is applied on the positive or higher side. -""" -mutable struct CoordinateBoundaryConditions{L, R} - left :: L - right :: R -end diff --git a/src/BoundaryConditions/discrete_boundary_function.jl b/src/BoundaryConditions/discrete_boundary_function.jl index 0bf87c36b7..287df30113 100644 --- a/src/BoundaryConditions/discrete_boundary_function.jl +++ b/src/BoundaryConditions/discrete_boundary_function.jl @@ -1,5 +1,5 @@ """ - DiscreteBoundaryFunction(func, parameters) + struct DiscreteBoundaryFunction{P, F} <: Function A wrapper for boundary condition functions with optional parameters. When `parameters=nothing`, the boundary condition `func` is called with the signature diff --git a/src/BoundaryConditions/field_boundary_conditions.jl b/src/BoundaryConditions/field_boundary_conditions.jl index d53fe85e4b..06ac8bc9bd 100644 --- a/src/BoundaryConditions/field_boundary_conditions.jl +++ b/src/BoundaryConditions/field_boundary_conditions.jl @@ -1,242 +1,150 @@ using Oceananigans.Operators: assumed_field_location -""" - FieldBoundaryConditions +##### +##### Default boundary conditions +##### + +struct DefaultPrognosticFieldBoundaryCondition end + +default_prognostic_field_boundary_condition(::Grids.Periodic, loc) = PeriodicBoundaryCondition() +default_prognostic_field_boundary_condition(::Bounded, ::Center) = NoFluxBoundaryCondition() +default_prognostic_field_boundary_condition(::Bounded, ::Face) = ImpenetrableBoundaryCondition() + +default_prognostic_field_boundary_condition(::Bounded, ::Nothing) = nothing +default_prognostic_field_boundary_condition(::Flat, ::Nothing) = nothing +default_prognostic_field_boundary_condition(::Grids.Periodic, ::Nothing) = nothing +default_prognostic_field_boundary_condition(::Flat, loc) = nothing + +default_auxiliary_field_boundary_condition(topo, loc) = default_prognostic_field_boundary_condition(topo, loc) +default_auxiliary_field_boundary_condition(::Bounded, ::Face) = nothing + +##### +##### Field boundary conditions +##### + +struct FieldBoundaryConditions{W, E, S, N, B, T, I} + west :: W + east :: E + south :: S + north :: N + bottom :: B + top :: T + immersed :: I +end -An alias for `NamedTuple{(:x, :y, :z)}` that represents a set of three `CoordinateBoundaryCondition`s -applied to a field along x, y, and z. """ -const FieldBoundaryConditions = NamedTuple{(:x, :y, :z)} + FieldBoundaryConditions(; kwargs...) -""" - FieldBoundaryConditions(x, y, z) +Returns a template for boundary conditions on prognostic fields. -Construct a `FieldBoundaryConditions` using a `CoordinateBoundaryCondition` for each of the -`x`, `y`, and `z` coordinates. -""" -FieldBoundaryConditions(x, y, z) = FieldBoundaryConditions((x, y, z)) +Keyword arguments specify boundary conditions on the 7 possible boundaries: -""" - DefaultBoundaryCondition(topo, ::Type{Nothing}) + * `west`, left end point in the `x`-direction where `i=1` + * `east`, right end point in the `x`-direction where `i=grid.Nx` + * `south`, left end point in the `y`-direction where `j=1` + * `north`, right end point in the `y`-direction where `j=grid.Ny` + * `bottom`, right end point in the `z`-direction where `k=1` + * `top`, right end point in the `z`-direction where `k=grid.Nz` + * `immersed`, boundary between solid and fluid for immersed boundaries (experimental support only) -Returns nothing. -""" -DefaultBoundaryCondition(topo, ::Type{Nothing}) = nothing -DefaultBoundaryCondition(::Type{Grids.Periodic}, ::Type{Nothing}) = nothing +If a boundary condition is unspecified, the default for prognostic fields +and the topology in the boundary-normal direction is used: + * `PeriodicBoundaryCondition` for `Periodic` directions + * `NoFluxBoundaryCondition` for `Bounded` directions and `Centered`-located fields + * `ImpenetrableBoundaryCondition` for `Bounded` directions and `Face`-located fields + * `nothing` for `Flat` directions and/or `Nothing`-located fields """ - DefaultBoundaryCondition(::Type{Periodic}, loc) +function FieldBoundaryConditions(; west = DefaultPrognosticFieldBoundaryCondition(), + east = DefaultPrognosticFieldBoundaryCondition(), + south = DefaultPrognosticFieldBoundaryCondition(), + north = DefaultPrognosticFieldBoundaryCondition(), + bottom = DefaultPrognosticFieldBoundaryCondition(), + top = DefaultPrognosticFieldBoundaryCondition(), + immersed = NoFluxBoundaryCondition()) -Returns [`PeriodicBoundaryCondition`](@ref). -""" -DefaultBoundaryCondition(::Type{Grids.Periodic}, loc) = PeriodicBoundaryCondition() + return FieldBoundaryConditions(west, east, south, north, bottom, top, immersed) +end """ - DefaultBoundaryCondition(::Type{Flat}, loc) + FieldBoundaryConditions(grid, location; kwargs...) -Returns `nothing`. -""" -DefaultBoundaryCondition(::Type{Flat}, loc) = nothing -DefaultBoundaryCondition(::Type{Flat}, ::Type{Nothing}) = nothing +Returns boundary conditions for auxiliary fields (fields +whose values are derived from a model's prognostic fields) on `grid` +and at `location`. -""" - DefaultBoundaryCondition(::Type{Bounded}, ::Type{Center}) +Keyword arguments specify boundary conditions on the 6 possible boundaries: -Returns [`NoFluxBoundaryCondition`](@ref). -""" -DefaultBoundaryCondition(::Type{Bounded}, ::Type{Center}) = NoFluxBoundaryCondition() + * `west`, left end point in the `x`-direction where `i=1` + * `east`, right end point in the `x`-direction where `i=grid.Nx` + * `south`, left end point in the `y`-direction where `j=1` + * `north`, right end point in the `y`-direction where `j=grid.Ny` + * `bottom`, right end point in the `z`-direction where `k=1` + * `top`, right end point in the `z`-direction where `k=grid.Nz` -""" - DefaultBoundaryCondition(::Type{Bounded}, ::Type{Face}) +If a boundary condition is unspecified, the default for auxiliary fields +and the topology in the boundary-normal direction is used: -Returns [`ImpenetrableBoundaryCondition`](@ref). + * `PeriodicBoundaryCondition` for `Periodic` directions + * `GradientBoundaryCondition(0)` for `Bounded` directions and `Centered`-located fields + * `nothing` for `Bounded` directions and `Face`-located fields + * `nothing` for `Flat` directions and/or `Nothing`-located fields) """ -DefaultBoundaryCondition(::Type{Bounded}, ::Type{Face}) = ImpenetrableBoundaryCondition() - -function validate_bcs(topology, left_bc, right_bc, default_bc, left_name, right_name, dir) - - if topology isa Periodic && (left_bc != default_bc || right_bc != default_bc) - e = "Cannot specify $left_name or $right_name boundary conditions with " * - "$topology topology in $dir-direction." - throw(ArgumentError(e)) - end +function FieldBoundaryConditions(grid, loc; + west = default_auxiliary_field_boundary_condition(topology(grid, 1)(), loc[1]()), + east = default_auxiliary_field_boundary_condition(topology(grid, 1)(), loc[1]()), + south = default_auxiliary_field_boundary_condition(topology(grid, 2)(), loc[2]()), + north = default_auxiliary_field_boundary_condition(topology(grid, 2)(), loc[2]()), + bottom = default_auxiliary_field_boundary_condition(topology(grid, 3)(), loc[3]()), + top = default_auxiliary_field_boundary_condition(topology(grid, 3)(), loc[3]()), + immersed = NoFluxBoundaryCondition()) - return true + return FieldBoundaryConditions(west, east, south, north, bottom, top, immersed) end -""" - FieldBoundaryConditions(grid, loc; east = DefaultBoundaryCondition(topology(grid, 1), loc[1]), - west = DefaultBoundaryCondition(topology(grid, 1], loc[1]), - south = DefaultBoundaryCondition(topology(grid, 2), loc[2]), - north = DefaultBoundaryCondition(topology(grid, 2), loc[2]), - bottom = DefaultBoundaryCondition(topology(grid, 3), loc[3]), - top = DefaultBoundaryCondition(topology(grid, 3), loc[3])) +##### +##### Boundary condition "regularization" +##### -Construct `FieldBoundaryConditions` for a field with location `loc` (a 3-tuple of `Face` or `Center`) -defined on `grid`. +regularize_boundary_condition(::DefaultPrognosticFieldBoundaryCondition, topo, loc, dim, args...) = + default_prognostic_field_boundary_condition(topo[dim](), loc[dim]()) -Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via -keyword arguments: +regularize_boundary_condition(bc, args...) = bc # fallback - * `west` and `east` for the `-x` and `+x` boundary; - * `south` and `north` for the `-y` and `+y` boundary; - * `bottom` and `top` for the `-z` and `+z` boundary. +""" +Compute default boundary conditions and attach field locations to ContinuousBoundaryFunction +boundary conditions for prognostic model field boundary conditions. -Default boundary conditions depend on `topology(grid)` and `loc`. +Note: don't regularize immersed boundary conditions: we don't support ContinuousBoundaryFunction +for immersed boundary conditions. """ -function FieldBoundaryConditions(grid, loc; east = DefaultBoundaryCondition(topology(grid, 1), loc[1]), - west = DefaultBoundaryCondition(topology(grid, 1), loc[1]), - south = DefaultBoundaryCondition(topology(grid, 2), loc[2]), - north = DefaultBoundaryCondition(topology(grid, 2), loc[2]), - bottom = DefaultBoundaryCondition(topology(grid, 3), loc[3]), - top = DefaultBoundaryCondition(topology(grid, 3), loc[3])) - - TX, TY, TZ = topology(grid) - x_default_bc = DefaultBoundaryCondition(topology(grid, 1), loc[1]) - y_default_bc = DefaultBoundaryCondition(topology(grid, 2), loc[2]) - z_default_bc = DefaultBoundaryCondition(topology(grid, 3), loc[3]) - - validate_bcs(TX, west, east, x_default_bc, :west, :east, :x) - validate_bcs(TY, south, north, y_default_bc, :south, :north, :y) - validate_bcs(TZ, bottom, top, z_default_bc, :bottom, :top, :z) - - x = CoordinateBoundaryConditions(west, east) - y = CoordinateBoundaryConditions(south, north) - z = CoordinateBoundaryConditions(bottom, top) - - return FieldBoundaryConditions(x, y, z) -end +function regularize_field_boundary_conditions(bcs::FieldBoundaryConditions, grid, field_name, prognostic_field_names=nothing) -""" - UVelocityBoundaryConditions(grid; east = DefaultBoundaryCondition(topology(grid, 1), Face), - west = DefaultBoundaryCondition(topology(grid, 1), Face), - south = DefaultBoundaryCondition(topology(grid, 2), Center), - north = DefaultBoundaryCondition(topology(grid, 2), Center), - bottom = DefaultBoundaryCondition(topology(grid, 3), Center), - top = DefaultBoundaryCondition(topology(grid, 3), Center)) - -Construct `FieldBoundaryConditions` for the `u`-velocity field on `grid`. -Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via -keyword arguments: - - * `west` and `east` for the `-x` and `+x` boundary; - * `south` and `north` for the `-y` and `+y` boundary; - * `bottom` and `top` for the `-z` and `+z` boundary. - -Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`. -""" -UVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Face, Center, Center); user_defined_bcs...) + topo = topology(grid) + loc = assumed_field_location(field_name) + + west = regularize_boundary_condition(bcs.west, topo, loc, 1, 1, prognostic_field_names) + east = regularize_boundary_condition(bcs.east, topo, loc, 1, grid.Nx, prognostic_field_names) + south = regularize_boundary_condition(bcs.south, topo, loc, 2, 1, prognostic_field_names) + north = regularize_boundary_condition(bcs.north, topo, loc, 2, grid.Ny, prognostic_field_names) + bottom = regularize_boundary_condition(bcs.bottom, topo, loc, 3, 1, prognostic_field_names) + top = regularize_boundary_condition(bcs.top, topo, loc, 3, grid.Nz, prognostic_field_names) -""" - VVelocityBoundaryConditions(grid; east = DefaultBoundaryCondition(topology(grid, 1), Center), - west = DefaultBoundaryCondition(topology(grid, 1), Center), - south = DefaultBoundaryCondition(topology(grid, 2), Face), - north = DefaultBoundaryCondition(topology(grid, 2), Face), - bottom = DefaultBoundaryCondition(topology(grid, 3), Center), - top = DefaultBoundaryCondition(topology(grid, 3), Center)) - -Construct `FieldBoundaryConditions` for the `v`-velocity field on `grid`. -Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via -keyword arguments: - - * `west` and `east` for the `-x` and `+x` boundary; - * `south` and `north` for the `-y` and `+y` boundary; - * `bottom` and `top` for the `-z` and `+z` boundary. - -Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`. -""" -VVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Center, Face, Center); user_defined_bcs...) + # Eventually we could envision supporting ContinuousForcing-style boundary conditions + # for the immersed boundary condition, which would benefit from regularization. + # But for now we don't regularize the immersed boundary condition. + immersed = bcs.immersed -""" - WVelocityBoundaryConditions(grid; east = DefaultBoundaryCondition(topology(grid, 1), Center), - west = DefaultBoundaryCondition(topology(grid, 1), Center), - south = DefaultBoundaryCondition(topology(grid, 2), Center), - north = DefaultBoundaryCondition(topology(grid, 2), Center), - bottom = DefaultBoundaryCondition(topology(grid, 3), Face), - top = DefaultBoundaryCondition(topology(grid, 3), Face)) - -Construct `FieldBoundaryConditions` for the `w`-velocity field on `grid`. -Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via -keyword arguments: - - * `west` and `east` for the `-x` and `+x` boundary; - * `south` and `north` for the `-y` and `+y` boundary; - * `bottom` and `top` for the `-z` and `+z` boundary. - -Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`. -""" -WVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Center, Center, Face); user_defined_bcs...) - -""" - TracerBoundaryConditions(grid; east = DefaultBoundaryCondition(topology(grid, 1), Center), - west = DefaultBoundaryCondition(topology(grid, 1), Center), - south = DefaultBoundaryCondition(topology(grid, 2), Center), - north = DefaultBoundaryCondition(topology(grid, 2), Center), - bottom = DefaultBoundaryCondition(topology(grid, 3), Center), - top = DefaultBoundaryCondition(topology(grid, 3), Center)) - -Construct `FieldBoundaryConditions` for a tracer field on `grid`. -Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via -keyword arguments: - - * `west` and `east` for the `-x` and `+x` boundary; - * `south` and `north` for the `-y` and `+y` boundary; - * `bottom` and `top` for the `-z` and `+z` boundary. - -Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`. -""" -TracerBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Center, Center, Center); user_defined_bcs...) - -const PressureBoundaryConditions = TracerBoundaryConditions -const DiffusivityBoundaryConditions = TracerBoundaryConditions - -# Here we overload setproperty! and getproperty to permit users to call -# the 'left' and 'right' bcs in the z-direction 'bottom' and 'top' -# and the 'left' and 'right' bcs in the y-direction 'south' and 'north'. -Base.setproperty!(bcs::FieldBoundaryConditions, side::Symbol, bc) = setbc!(bcs, Val(side), bc) - -setbc!(bcs::FieldBoundaryConditions, ::Val{S}, bc) where S = setfield!(bcs, S, bc) -setbc!(bcs::FieldBoundaryConditions, ::Val{:west}, bc) = setfield!(bcs.x, :left, bc) -setbc!(bcs::FieldBoundaryConditions, ::Val{:east}, bc) = setfield!(bcs.x, :right, bc) -setbc!(bcs::FieldBoundaryConditions, ::Val{:south}, bc) = setfield!(bcs.y, :left, bc) -setbc!(bcs::FieldBoundaryConditions, ::Val{:north}, bc) = setfield!(bcs.y, :right, bc) -setbc!(bcs::FieldBoundaryConditions, ::Val{:bottom}, bc) = setfield!(bcs.z, :left, bc) -setbc!(bcs::FieldBoundaryConditions, ::Val{:top}, bc) = setfield!(bcs.z, :right, bc) - -@inline Base.getproperty(bcs::FieldBoundaryConditions, side::Symbol) = getbc(bcs, Val(side)) - -@inline getbc(bcs::FieldBoundaryConditions, ::Val{S}) where S = getfield(bcs, S) -@inline getbc(bcs::FieldBoundaryConditions, ::Val{:west}) = getfield(bcs.x, :left) -@inline getbc(bcs::FieldBoundaryConditions, ::Val{:east}) = getfield(bcs.x, :right) -@inline getbc(bcs::FieldBoundaryConditions, ::Val{:south}) = getfield(bcs.y, :left) -@inline getbc(bcs::FieldBoundaryConditions, ::Val{:north}) = getfield(bcs.y, :right) -@inline getbc(bcs::FieldBoundaryConditions, ::Val{:bottom}) = getfield(bcs.z, :left) -@inline getbc(bcs::FieldBoundaryConditions, ::Val{:top}) = getfield(bcs.z, :right) - -function regularize_field_boundary_conditions(bcs::FieldBoundaryConditions, grid, model_field_names, field_name) - - X, Y, Z = assumed_field_location(field_name) - - east = regularize_boundary_condition(bcs.east, Nothing, Y, Z, grid.Nx, model_field_names) - west = regularize_boundary_condition(bcs.west, Nothing, Y, Z, 1, model_field_names) - north = regularize_boundary_condition(bcs.north, X, Nothing, Z, grid.Ny, model_field_names) - south = regularize_boundary_condition(bcs.south, X, Nothing, Z, 1, model_field_names) - top = regularize_boundary_condition(bcs.top, X, Y, Nothing, grid.Nz, model_field_names) - bottom = regularize_boundary_condition(bcs.bottom, X, Y, Nothing, 1, model_field_names) - - x = CoordinateBoundaryConditions(west, east) - y = CoordinateBoundaryConditions(south, north) - z = CoordinateBoundaryConditions(bottom, top) - - return FieldBoundaryConditions(x, y, z) + return FieldBoundaryConditions(west, east, south, north, bottom, top, immersed) end -function regularize_field_boundary_conditions(boundary_conditions::NamedTuple, grid, model_field_names, field_name=nothing) - boundary_conditions_names = propertynames(boundary_conditions) - boundary_conditions_tuple = Tuple(regularize_field_boundary_conditions(bcs, grid, model_field_names, name) - for (name, bcs) in zip(boundary_conditions_names, boundary_conditions)) - boundary_conditions = NamedTuple{boundary_conditions_names}(boundary_conditions_tuple) - return boundary_conditions -end +regularize_field_boundary_conditions(boundary_conditions::NamedTuple, grid, prognostic_field_names) = + NamedTuple(field_name => regularize_field_boundary_conditions(field_bcs, grid, field_name, prognostic_field_names) + for (field_name, field_bcs) in pairs(boundary_conditions)) + +# For nested NamedTuples of boundary conditions (eg diffusivity boundary conditions) +regularize_field_boundary_conditions(boundary_conditions::NamedTuple, grid, ::Symbol, prognostic_field_names) = + NamedTuple(field_name => regularize_field_boundary_conditions(field_bcs, grid, field_name, prognostic_field_names) + for (field_name, field_bcs) in pairs(boundary_conditions)) -regularize_field_boundary_conditions(::Missing, grid, model_field_names, field_name=nothing) = missing +regularize_field_boundary_conditions(::Missing, grid, field_name, prognostic_field_names=nothing) = missing diff --git a/src/BoundaryConditions/fill_halo_regions.jl b/src/BoundaryConditions/fill_halo_regions.jl index 0e6ed2f0c4..9dbad8017f 100644 --- a/src/BoundaryConditions/fill_halo_regions.jl +++ b/src/BoundaryConditions/fill_halo_regions.jl @@ -7,7 +7,7 @@ using OffsetArrays: OffsetArray fill_halo_regions!(::Nothing, args...) = [] """ - fill_halo_regions!(fields, arch) + fill_halo_regions!(fields::Union{Tuple, NamedTuple}, arch, args...) Fill halo regions for each field in the tuple `fields` according to their boundary conditions, possibly recursing into `fields` if it is a nested tuple-of-tuples. diff --git a/src/BoundaryConditions/show_boundary_conditions.jl b/src/BoundaryConditions/show_boundary_conditions.jl index 479b1bca4c..3d2d4f498a 100644 --- a/src/BoundaryConditions/show_boundary_conditions.jl +++ b/src/BoundaryConditions/show_boundary_conditions.jl @@ -1,13 +1,16 @@ import Base: show import Oceananigans: short_show -bcclassification_str(::FBC) = "Flux" -bcclassification_str(::PBC) = "Periodic" -bcclassification_str(::OBC) = "Open" -bcclassification_str(::VBC) = "Value" -bcclassification_str(::GBC) = "Gradient" -bcclassification_str(::ZFBC) = "ZeroFlux" -bcclassification_str(::Nothing) = "Nothing" +const DFBC = DefaultPrognosticFieldBoundaryCondition + +bc_str(::FBC) = "Flux" +bc_str(::PBC) = "Periodic" +bc_str(::OBC) = "Open" +bc_str(::VBC) = "Value" +bc_str(::GBC) = "Gradient" +bc_str(::ZFBC) = "ZeroFlux" +bc_str(::DFBC) = "Default" +bc_str(::Nothing) = "Nothing" ##### ##### BoundaryCondition @@ -24,21 +27,29 @@ function print_condition(f::Function) end show(io::IO, bc::BoundaryCondition) = - print(io, "BoundaryCondition: type=$(bcclassification_str(bc)), condition=$(print_condition(bc.condition))") + print(io, "BoundaryCondition: classification=$(bc_str(bc)), condition=$(print_condition(bc.condition))") ##### ##### FieldBoundaryConditions ##### short_show(fbcs::FieldBoundaryConditions) = - string("x=(west=$(bcclassification_str(fbcs.x.left)), east=$(bcclassification_str(fbcs.x.right))), ", - "y=(south=$(bcclassification_str(fbcs.y.left)), north=$(bcclassification_str(fbcs.y.right))), ", - "z=(bottom=$(bcclassification_str(fbcs.z.left)), top=$(bcclassification_str(fbcs.z.right)))") + string("west=$(bc_str(fbcs.west)), ", + "east=$(bc_str(fbcs.east)), ", + "south=$(bc_str(fbcs.south)), ", + "north=$(bc_str(fbcs.north)), ", + "bottom=$(bc_str(fbcs.bottom)), ", + "top=$(bc_str(fbcs.top)), ", + "immersed=$(bc_str(fbcs.immersed))") show_field_boundary_conditions(bcs::FieldBoundaryConditions, padding="") = - string("Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions", '\n', - padding, "├── x: ", typeof(bcs.x), '\n', - padding, "├── y: ", typeof(bcs.y), '\n', - padding, "└── z: ", typeof(bcs.z)) + string("Oceananigans.FieldBoundaryConditions, with boundary conditions", '\n', + padding, "├── west: ", typeof(bcs.west), '\n', + padding, "├── east: ", typeof(bcs.east), '\n', + padding, "├── south: ", typeof(bcs.south), '\n', + padding, "├── north: ", typeof(bcs.north), '\n', + padding, "├── bottom: ", typeof(bcs.bottom), '\n', + padding, "├── top: ", typeof(bcs.top), '\n', + padding, "└── immersed: ", typeof(bcs.immersed)) Base.show(io::IO, fieldbcs::FieldBoundaryConditions) = print(io, show_field_boundary_conditions(fieldbcs)) diff --git a/src/CubedSpheres/CubedSpheres.jl b/src/CubedSpheres/CubedSpheres.jl index 9a527a83e8..d9bced66de 100644 --- a/src/CubedSpheres/CubedSpheres.jl +++ b/src/CubedSpheres/CubedSpheres.jl @@ -36,11 +36,23 @@ validate_vertical_velocity_boundary_conditions(w::AbstractCubedSphereField) = import Oceananigans.BoundaryConditions: regularize_field_boundary_conditions -function regularize_field_boundary_conditions(bcs::CubedSphereFaces, grid, model_field_names, field_name) +function regularize_field_boundary_conditions(bcs::CubedSphereFaces, grid, field_name, prognostic_field_names) + + faces = Tuple(regularize_field_boundary_conditions(face_bcs, face_grid, field_name, prognostic_field_names) + for (face_bcs, face_grid) in zip(bcs.faces, grid.faces)) + + return CubedSphereFaces{typeof(faces[1]), typeof(faces)}(faces) +end + +function regularize_field_boundary_conditions(bcs::FieldBoundaryConditions, grid::ConformalCubedSphereGrid, field_name, prognostic_field_names) faces = Tuple( - regularize_field_boundary_conditions(face_bcs, face_grid, model_field_names, field_name) - for (face_bcs, face_grid) in zip(bcs.faces, grid.faces) + inject_cubed_sphere_exchange_boundary_conditions( + regularize_field_boundary_conditions(bcs, face_grid, field_name, prognostic_field_names), + face_index, + grid.face_connectivity + ) + for (face_index, face_grid) in enumerate(grid.faces) ) return CubedSphereFaces{typeof(faces[1]), typeof(faces)}(faces) diff --git a/src/CubedSpheres/cubed_sphere_exchange_bcs.jl b/src/CubedSpheres/cubed_sphere_exchange_bcs.jl index a816b7d0b4..158b64d684 100644 --- a/src/CubedSpheres/cubed_sphere_exchange_bcs.jl +++ b/src/CubedSpheres/cubed_sphere_exchange_bcs.jl @@ -4,13 +4,13 @@ using Oceananigans.BoundaryConditions using Oceananigans.BoundaryConditions: AbstractBoundaryConditionClassification import Base: show -import Oceananigans.BoundaryConditions: bcclassification_str, print_condition +import Oceananigans.BoundaryConditions: bc_str, print_condition struct CubedSphereExchange <: AbstractBoundaryConditionClassification end const CubedSphereExchangeBC = BoundaryCondition{<:CubedSphereExchange} -bcclassification_str(::CubedSphereExchangeBC) ="CubedSphereExchange" +bc_str(::CubedSphereExchangeBC) ="CubedSphereExchange" CubedSphereExchangeBoundaryCondition(val; kwargs...) = BoundaryCondition(CubedSphereExchange, val; kwargs...) @@ -65,10 +65,13 @@ function inject_cubed_sphere_exchange_boundary_conditions(field_bcs, face_number south_exchange_bc = CubedSphereExchangeBoundaryCondition(south_exchange_info) north_exchange_bc = CubedSphereExchangeBoundaryCondition(north_exchange_info) - x_bcs = CoordinateBoundaryConditions(west_exchange_bc, east_exchange_bc) - y_bcs = CoordinateBoundaryConditions(south_exchange_bc, north_exchange_bc) - - return FieldBoundaryConditions(x_bcs, y_bcs, field_bcs.z) + return FieldBoundaryConditions(west_exchange_bc, + east_exchange_bc, + south_exchange_bc, + north_exchange_bc, + field_bcs.bottom, + field_bcs.top, + field_bcs.immersed) end Adapt.adapt_structure(to, ::CubedSphereExchangeInformation) = nothing diff --git a/src/Distributed/halo_communication_bcs.jl b/src/Distributed/halo_communication_bcs.jl index de91d4d8ae..11b9edbf84 100644 --- a/src/Distributed/halo_communication_bcs.jl +++ b/src/Distributed/halo_communication_bcs.jl @@ -1,13 +1,13 @@ using Oceananigans.BoundaryConditions using Oceananigans.BoundaryConditions: AbstractBoundaryConditionClassification -import Oceananigans.BoundaryConditions: bcclassification_str, print_condition +import Oceananigans.BoundaryConditions: bc_str, print_condition struct HaloCommunication <: AbstractBoundaryConditionClassification end const HaloCommunicationBC = BoundaryCondition{<:HaloCommunication} -bcclassification_str(::HaloCommunicationBC) = "HaloCommunication" +bc_str(::HaloCommunicationBC) = "HaloCommunication" HaloCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(HaloCommunication, val; kwargs...) @@ -42,14 +42,13 @@ function inject_halo_communication_boundary_conditions(field_bcs, local_rank, co top_comm_bc = HaloCommunicationBoundaryCondition(top_comm_ranks) bottom_comm_bc = HaloCommunicationBoundaryCondition(bottom_comm_ranks) - x_bcs = CoordinateBoundaryConditions(isnothing(rank_west) ? field_bcs.west : west_comm_bc, - isnothing(rank_east) ? field_bcs.east : east_comm_bc) + west = isnothing(rank_west) ? field_bcs.west : west_comm_bc + east = isnothing(rank_east) ? field_bcs.east : east_comm_bc + south = isnothing(rank_south) ? field_bcs.south : south_comm_bc + north = isnothing(rank_north) ? field_bcs.north : north_comm_bc + bottom = isnothing(rank_bottom) ? field_bcs.bottom : bottom_comm_bc + top = isnothing(rank_top) ? field_bcs.top : top_comm_bc + immersed = field_bcs.immersed - y_bcs = CoordinateBoundaryConditions(isnothing(rank_south) ? field_bcs.south : south_comm_bc, - isnothing(rank_north) ? field_bcs.north : north_comm_bc) - - z_bcs = CoordinateBoundaryConditions(isnothing(rank_bottom) ? field_bcs.bottom : bottom_comm_bc, - isnothing(rank_top) ? field_bcs.top : top_comm_bc) - - return FieldBoundaryConditions(x_bcs, y_bcs, z_bcs) + return FieldBoundaryConditions(west, east, south, north, bottom, top, immersed) end diff --git a/src/Fields/computed_field.jl b/src/Fields/computed_field.jl index aeb1fefdf8..abacf05e27 100644 --- a/src/Fields/computed_field.jl +++ b/src/Fields/computed_field.jl @@ -31,26 +31,6 @@ struct ComputedField{X, Y, Z, S, O, A, D, G, T, C} <: AbstractDataField{X, Y, Z, end end -# We define special default boundary conditions for ComputedField, because currently -# `DefaultBoundaryCondition` uses ImpenetrableBoundaryCondition() in bounded directions -# and for fields on Faces, which is not what we want in general for ComputedFields. -DefaultComputedFieldBoundaryCondition(::Type{Grids.Periodic}, loc) = PeriodicBoundaryCondition() -DefaultComputedFieldBoundaryCondition(::Type{Flat}, loc) = nothing -DefaultComputedFieldBoundaryCondition(::Type{Bounded}, ::Type{Center}) = NoFluxBoundaryCondition() -DefaultComputedFieldBoundaryCondition(::Type{Bounded}, ::Type{Face}) = nothing - -function ComputedFieldBoundaryConditions(grid, loc; - east = DefaultComputedFieldBoundaryCondition(topology(grid, 1), loc[1]), - west = DefaultComputedFieldBoundaryCondition(topology(grid, 1), loc[1]), - south = DefaultComputedFieldBoundaryCondition(topology(grid, 2), loc[2]), - north = DefaultComputedFieldBoundaryCondition(topology(grid, 2), loc[2]), - bottom = DefaultComputedFieldBoundaryCondition(topology(grid, 3), loc[3]), - top = DefaultComputedFieldBoundaryCondition(topology(grid, 3), loc[3])) - - return FieldBoundaryConditions(grid, loc; east=east, west=west, south=south, north=north, bottom=bottom, top=top) -end - - """ ComputedField(operand [, arch=nothing]; data = nothing, recompute_safely = true, boundary_conditions = ComputedFieldBoundaryConditions(operand.grid, location(operand)) @@ -75,7 +55,7 @@ end function ComputedField(LX, LY, LZ, operand, arch, grid; data = nothing, recompute_safely = true, - boundary_conditions = ComputedFieldBoundaryConditions(grid, (LX, LY, LZ))) + boundary_conditions = FieldBoundaryConditions(grid, (LX, LY, LZ))) # Architecturanigans operand_arch = architecture(operand) diff --git a/src/Fields/field.jl b/src/Fields/field.jl index c1cce21ccb..771cd07004 100644 --- a/src/Fields/field.jl +++ b/src/Fields/field.jl @@ -1,24 +1,24 @@ using Adapt -struct Field{X, Y, Z, A, D, G, T, B} <: AbstractDataField{X, Y, Z, A, G, T, 3} +struct Field{LX, LY, LZ, A, D, G, T, B} <: AbstractDataField{LX, LY, LZ, A, G, T, 3} data :: D architecture :: A grid :: G boundary_conditions :: B - function Field{X, Y, Z}(data::D, arch::A, grid::G, bcs::B) where {X, Y, Z, D, A, G, B} + function Field{LX, LY, LZ}(data::D, arch::A, grid::G, bcs::B) where {LX, LY, LZ, D, A, G, B} T = eltype(grid) - return new{X, Y, Z, A, D, G, T, B}(data, arch, grid, bcs) + return new{LX, LY, LZ, A, D, G, T, B}(data, arch, grid, bcs) end end """ - Field(X, Y, Z, [arch = CPU()], grid, - [ bcs = FieldBoundaryConditions(grid, (X, Y, Z)), - data = new_data(eltype(grid), arch, grid, (X, Y, Z))]) + Field(LX, LY, LZ, [arch = CPU()], grid, + [ bcs = FieldBoundaryConditions(grid, (LX, LY, LZ)), + data = new_data(eltype(grid), arch, grid, (LX, LY, LZ))]) Construct a `Field` on `grid` with `data` on architecture `arch` with -boundary conditions `bcs`. Each of `(X, Y, Z)` is either `Center` or `Face` and determines +boundary conditions `bcs`. Each of `(LX, LY, LZ)` is either `Center` or `Face` and determines the field's location in `(x, y, z)`. Example @@ -31,25 +31,25 @@ julia> ω = Field(Face, Face, Center, CPU(), RegularRectilinearGrid(size=(1, 1, Field located at (Face, Face, Center) ├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (1, 1, 1) ├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=1, Ny=1, Nz=1) -└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux) +└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux ``` """ -function Field(X, Y, Z, +function Field(LX, LY, LZ, arch::AbstractArchitecture, grid::AbstractGrid, - bcs = FieldBoundaryConditions(grid, (X, Y, Z)), - data = new_data(eltype(grid), arch, grid, (X, Y, Z))) + bcs = FieldBoundaryConditions(grid, (LX, LY, LZ)), + data = new_data(eltype(grid), arch, grid, (LX, LY, LZ))) - validate_field_data(X, Y, Z, data, grid) + validate_field_data(LX, LY, LZ, data, grid) - return Field{X, Y, Z}(data, arch, grid, bcs) + return Field{LX, LY, LZ}(data, arch, grid, bcs) end # Default CPU architecture -Field(X, Y, Z, grid::AbstractGrid, args...) = Field(X, Y, Z, CPU(), grid, args...) +Field(LX, LY, LZ, grid::AbstractGrid, args...) = Field(LX, LY, LZ, CPU(), grid, args...) # Canonical `similar` for AbstractField (doesn't transfer boundary conditions) -Base.similar(f::AbstractField{X, Y, Z, Arch}) where {X, Y, Z, Arch} = Field(X, Y, Z, Arch(), f.grid) +Base.similar(f::AbstractField{LX, LY, LZ, Arch}) where {LX, LY, LZ, Arch} = Field(LX, LY, LZ, Arch(), f.grid) # Type "destantiation": convert Face() to Face and Center() to Center if needed. destantiate(X) = typeof(X) @@ -68,68 +68,36 @@ Field(L::Tuple, args...) = Field(destantiate.(L)..., args...) ##### """ - CenterField([arch=CPU()], grid - bcs = TracerBoundaryConditions(grid), - data = new_data(eltype(grid), arch, grid, (Center, Center, Center))) + CenterField([arch=CPU()], grid, args...) -Return a `Field{Center, Center, Center}` on architecture `arch` and `grid` containing `data` -with field boundary conditions `bcs`. +Returns `Field{Center, Center, Center}` on `arch`itecture and `grid`. +Additional arguments are passed to the `Field` constructor. """ -function CenterField(arch::AbstractArchitecture, - grid, - bcs = TracerBoundaryConditions(grid), - data = new_data(eltype(grid), arch, grid, (Center, Center, Center))) - - return Field(Center, Center, Center, arch, grid, bcs, data) -end +CenterField(arch::AbstractArchitecture, grid, args...) = Field(Center, Center, Center, arch, grid, args...) """ - XFaceField([arch=CPU()], grid - bcs = UVelocityBoundaryConditions(grid), - data = new_data(eltype(grid), arch, grid, (Face, Center, Center))) + XFaceField([arch=CPU()], grid, args...) -Return a `Field{Face, Center, Center}` on architecture `arch` and `grid` containing `data` -with field boundary conditions `bcs`. +Returns `Field{Face, Center, Center}` on `arch`itecture and `grid`. +Additional arguments are passed to the `Field` constructor. """ -function XFaceField(arch::AbstractArchitecture, - grid, - bcs = UVelocityBoundaryConditions(grid), - data = new_data(eltype(grid), arch, grid, (Face, Center, Center))) - - return Field(Face, Center, Center, arch, grid, bcs, data) -end +XFaceField(arch::AbstractArchitecture, args...) = Field(Face, Center, Center, arch, args...) """ - YFaceField([arch=CPU()], grid, - bcs = VVelocityBoundaryConditions(grid), - data = new_data(eltype(grid), arch, grid, (Center, Face, Center))) + YFaceField([arch=CPU()], grid, args...) -Return a `Field{Center, Face, Center}` on architecture `arch` and `grid` containing `data` -with field boundary conditions `bcs`. +Returns `Field{Center, Face, Center}` on `arch`itecture and `grid`. +Additional arguments are passed to the `Field` constructor. """ -function YFaceField(arch::AbstractArchitecture, - grid, - bcs = VVelocityBoundaryConditions(grid), - data = new_data(eltype(grid), arch, grid, (Center, Face, Center))) - - return Field(Center, Face, Center, arch, grid, bcs, data) -end +YFaceField(arch::AbstractArchitecture, args...) = Field(Center, Face, Center, arch, args...) """ - ZFaceField([arch=CPU()], grid - bcs = WVelocityBoundaryConditions(grid), - data = new_data(eltype(grid), arch, grid, (Center, Center, Face))) + ZFaceField([arch=CPU()], grid, args...) -Return a `Field{Center, Center, Face}` on architecture `arch` and `grid` containing `data` -with field boundary conditions `bcs`. +Returns `Field{Center, Center, Face}` on `arch`itecture and `grid`. +Additional arguments are passed to the `Field` constructor. """ -function ZFaceField(arch::AbstractArchitecture, - grid, - bcs = WVelocityBoundaryConditions(grid), - data = new_data(eltype(grid), arch, grid, (Center, Center, Face))) - - return Field(Center, Center, Face, arch, grid, bcs, data) -end +ZFaceField(arch::AbstractArchitecture, args...) = Field(Center, Center, Face, arch, args...) # Default CPU architectures... CenterField(grid::AbstractGrid, args...) = CenterField(CPU(), grid, args...) diff --git a/src/Fields/field_tuples.jl b/src/Fields/field_tuples.jl index 164c47dc08..2c83ae6a2a 100644 --- a/src/Fields/field_tuples.jl +++ b/src/Fields/field_tuples.jl @@ -1,3 +1,5 @@ +using Oceananigans.BoundaryConditions: FieldBoundaryConditions, regularize_field_boundary_conditions + ##### ##### Velocity fields tuples ##### @@ -9,19 +11,25 @@ Return a NamedTuple with fields `u`, `v`, `w` initialized on the architecture `a and `grid`. Boundary conditions `bcs` may be specified via a named tuple of `FieldBoundaryCondition`s. """ -function VelocityFields(arch, grid, bcs=NamedTuple()) - u_bcs = :u ∈ keys(bcs) ? bcs.u : UVelocityBoundaryConditions(grid) - v_bcs = :v ∈ keys(bcs) ? bcs.v : VVelocityBoundaryConditions(grid) - w_bcs = :w ∈ keys(bcs) ? bcs.w : WVelocityBoundaryConditions(grid) +function VelocityFields(arch, grid, user_bcs = NamedTuple()) + + template = FieldBoundaryConditions() + + default_bcs = ( + u = regularize_field_boundary_conditions(template, grid, :u), + v = regularize_field_boundary_conditions(template, grid, :v), + w = regularize_field_boundary_conditions(template, grid, :w) + ) + + bcs = merge(default_bcs, user_bcs) - u = XFaceField(arch, grid, u_bcs) - v = YFaceField(arch, grid, v_bcs) - w = ZFaceField(arch, grid, w_bcs) + u = XFaceField(arch, grid, bcs.u) + v = YFaceField(arch, grid, bcs.v) + w = ZFaceField(arch, grid, bcs.w) return (u=u, v=v, w=w) end - ##### ##### Tracer fields tuples ##### @@ -33,15 +41,10 @@ Returns a `NamedTuple` with tracer fields specified by `tracer_names` initialize `CenterField`s on the architecture `arch` and `grid`. Boundary conditions `bcs` may be specified via a named tuple of `FieldBoundaryCondition`s. """ -function TracerFields(tracer_names, arch, grid, bcs) - - tracer_fields = - Tuple(c ∈ keys(bcs) ? - CenterField(arch, grid, bcs[c]) : - CenterField(arch, grid, TracerBoundaryConditions(grid)) - for c in tracer_names) - - return NamedTuple{tracer_names}(tracer_fields) +function TracerFields(tracer_names, arch, grid, user_bcs) + default_bcs = NamedTuple(name => FieldBoundaryConditions(grid, (Center, Center, Center)) for name in tracer_names) + bcs = merge(default_bcs, user_bcs) # provided bcs overwrite defaults + return NamedTuple(c => CenterField(arch, grid, bcs[c]) for c in tracer_names) end """ @@ -54,13 +57,8 @@ keyword arguments `kwargs` for each field. This function is used by `OutputWriters.Checkpointer` and `TendencyFields`. ``` """ -function TracerFields(tracer_names, arch, grid; kwargs...) - tracer_fields = - Tuple(c ∈ keys(kwargs) ? kwargs[c] : CenterField(arch, grid, TracerBoundaryConditions(grid)) - for c in tracer_names) - - return NamedTuple{tracer_names}(tracer_fields) -end +TracerFields(tracer_names, arch, grid; kwargs...) = + NamedTuple(c => c ∈ keys(kwargs) ? kwargs[c] : CenterField(arch, grid) for c in tracer_names) # 'Nothing', or empty tracer fields TracerFields(::Union{Tuple{}, Nothing}, arch, grid, bcs) = NamedTuple() @@ -80,11 +78,15 @@ Return a NamedTuple with pressure fields `pHY′` and `pNHS` initialized as be specified via a named tuple of `FieldBoundaryCondition`s. """ function PressureFields(arch, grid, bcs=NamedTuple()) - pHY′_bcs = :pHY′ ∈ keys(bcs) ? bcs[:pHY′] : PressureBoundaryConditions(grid) - pNHS_bcs = :pNHS ∈ keys(bcs) ? bcs[:pNHS] : PressureBoundaryConditions(grid) - pHY′ = CenterField(arch, grid, pHY′_bcs) - pNHS = CenterField(arch, grid, pNHS_bcs) + 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(arch, grid, bcs.pHY′) + pNHS = CenterField(arch, grid, bcs.pNHS) return (pHY′=pHY′, pNHS=pNHS) end @@ -97,9 +99,9 @@ tracer fields), initialized on the architecture `arch` and `grid`. Optional `kwa can be specified to assign data arrays to each tendency field. """ function TendencyFields(arch, grid, tracer_names; - u = XFaceField(arch, grid, UVelocityBoundaryConditions(grid)), - v = YFaceField(arch, grid, VVelocityBoundaryConditions(grid)), - w = ZFaceField(arch, grid, WVelocityBoundaryConditions(grid)), + u = XFaceField(arch, grid), + v = YFaceField(arch, grid), + w = ZFaceField(arch, grid), kwargs...) velocities = (u=u, v=v, w=w) @@ -134,7 +136,7 @@ function VelocityFields(proposed_velocities::NamedTuple{(:u, :v, :w)}, arch, gri end """ - TracerFields(proposed_tracerc::NamedTuple, arch, grid, bcs) + TracerFields(proposed_tracers::NamedTuple, arch, grid, bcs) Returns a `NamedTuple` of tracers, overwriting boundary conditions in `proposed_tracers` with corresponding fields in the `NamedTuple` `bcs`. @@ -144,9 +146,7 @@ function TracerFields(proposed_tracers::NamedTuple, arch, grid, bcs) validate_field_tuple_grid("tracers", proposed_tracers, grid) tracer_names = propertynames(proposed_tracers) - - tracer_fields = Tuple(CenterField(arch, grid, bcs[c], proposed_tracers[c].data) - for c in tracer_names) + tracer_fields = Tuple(CenterField(arch, grid, bcs[c], proposed_tracers[c].data) for c in tracer_names) return NamedTuple{tracer_names}(tracer_fields) end diff --git a/src/Fields/kernel_computed_field.jl b/src/Fields/kernel_computed_field.jl index 1f4aa318ec..c4de51bd67 100644 --- a/src/Fields/kernel_computed_field.jl +++ b/src/Fields/kernel_computed_field.jl @@ -13,7 +13,7 @@ struct KernelComputedField{X, Y, Z, A, S, D, G, T, K, B, F, P} <: AbstractDataFi status :: S function KernelComputedField{X, Y, Z}(kernel::K, arch::A, grid::G; - boundary_conditions::B = ComputedFieldBoundaryConditions(grid, (X, Y, Z)), + boundary_conditions::B = FieldBoundaryConditions(grid, (X, Y, Z)), computed_dependencies = (), parameters::P = nothing, data::D = new_data(arch, grid, (X, Y, Z)), @@ -35,7 +35,7 @@ end """ KernelComputedField(X, Y, Z, kernel, model; - boundary_conditions = ComputedFieldBoundaryConditions(grid, (X, Y, Z)), + boundary_conditions = FieldBoundaryConditions(grid, (X, Y, Z)), computed_dependencies = (), parameters = nothing, data = nothing, diff --git a/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl b/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl index 42a64b4e0c..db7a41e477 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl @@ -33,8 +33,7 @@ include("compute_vertically_integrated_lateral_face_areas.jl") include("compute_vertically_integrated_volume_flux.jl") include("implicit_free_surface_solver.jl") -include("hydrostatic_free_surface_velocity_fields.jl") -include("hydrostatic_free_surface_tendency_fields.jl") +include("hydrostatic_free_surface_field_tuples.jl") include("hydrostatic_free_surface_model.jl") include("show_hydrostatic_free_surface_model.jl") include("set_hydrostatic_free_surface_model.jl") diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_field_tuples.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_field_tuples.jl new file mode 100644 index 0000000000..cf51f6c225 --- /dev/null +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_field_tuples.jl @@ -0,0 +1,19 @@ +using Oceananigans.Grids: Center, Face +using Oceananigans.Fields: XFaceField, YFaceField, ZFaceField, TracerFields + +function HydrostaticFreeSurfaceVelocityFields(::Nothing, arch, grid, clock, bcs=NamedTuple()) + u = XFaceField(arch, grid, bcs.u) + v = YFaceField(arch, grid, bcs.v) + w = ZFaceField(arch, grid) + + return (u=u, v=v, w=w) +end + +function HydrostaticFreeSurfaceTendencyFields(velocities, free_surface, arch, grid, tracer_names) + u = XFaceField(arch, grid) + v = YFaceField(arch, grid) + η = FreeSurfaceDisplacementField(velocities, free_surface, arch, grid) + tracers = TracerFields(tracer_names, arch, grid) + + return merge((u=u, v=v, η=η), tracers) +end diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl index 9a3c053ea6..36e8fa4a54 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl @@ -6,7 +6,7 @@ using Oceananigans: AbstractModel, AbstractOutputWriter, AbstractDiagnostic using Oceananigans.Architectures: AbstractArchitecture, GPU using Oceananigans.Advection: AbstractAdvectionScheme, CenteredSecondOrder using Oceananigans.BuoyancyModels: validate_buoyancy, regularize_buoyancy, SeawaterBuoyancy, g_Earth -using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions, TracerBoundaryConditions +using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions using Oceananigans.Fields: Field, CenterField, tracernames, VelocityFields, TracerFields using Oceananigans.Forcings: model_forcing using Oceananigans.Grids: inflate_halo_size, with_halo, AbstractRectilinearGrid, AbstractCurvilinearGrid, AbstractHorizontallyCurvilinearGrid @@ -24,7 +24,7 @@ validate_tracer_advection(invalid_tracer_advection, grid) = error("$invalid_trac validate_tracer_advection(tracer_advection_tuple::NamedTuple, grid) = CenteredSecondOrder(), tracer_advection_tuple validate_tracer_advection(tracer_advection::AbstractAdvectionScheme, grid) = tracer_advection, NamedTuple() -PressureField(arch, grid) = (pHY′ = CenterField(arch, grid, TracerBoundaryConditions(grid)),) +PressureField(arch, grid) = (; pHY′ = CenterField(arch, grid)) mutable struct HydrostaticFreeSurfaceModel{TS, E, A<:AbstractArchitecture, S, G, T, V, B, R, F, P, U, C, Φ, K, AF} <: AbstractModel{TS} @@ -112,24 +112,33 @@ function HydrostaticFreeSurfaceModel(; grid, momentum_advection = validate_momentum_advection(momentum_advection, grid) tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example) - validate_buoyancy(buoyancy, tracernames(tracers)) + validate_buoyancy(buoyancy, tracernames(tracers)) buoyancy = regularize_buoyancy(buoyancy) - # Recursively "regularize" field-dependent boundary conditions by supplying list of tracer names. - # We also regularize boundary conditions included in velocities, tracers, pressure, and diffusivities. + # Collect boundary conditions for all model prognostic fields and, if specified, some model + # auxiliary fields. Boundary conditions are "regularized" based on the _name_ of the field: + # boundary conditions on u, v are regularized assuming they represent momentum at appropriate + # staggered locations. All other fields are regularized assuming they are tracers. # Note that we do not regularize boundary conditions contained in *tupled* diffusivity fields right now. + # + # 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(pressure), extract_boundary_conditions(diffusivities)) - boundary_conditions = merge(embedded_boundary_conditions, boundary_conditions) + # Next, we form a list of default boundary conditions: + prognostic_field_names = (:u, :v, :η, tracernames(tracers)...) + default_boundary_conditions = NamedTuple{prognostic_field_names}(Tuple(FieldBoundaryConditions() for name in prognostic_field_names)) - model_field_names = (:u, :v, :w, :η, tracernames(tracers)...) - boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, model_field_names) + # Then we merge specified, embedded, and default boundary conditions. Specified boundary conditions + # have precedence, followed by embedded, followed by default. + boundary_conditions = merge(default_boundary_conditions, embedded_boundary_conditions, boundary_conditions) + boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, prognostic_field_names) - # TKEBasedVerticalDiffusivity enforces boundary conditions on TKE: + # Finally, we ensure that closure-specific boundary conditions, such as + # those required by TKEBasedVerticalDiffusivity, are enforced: boundary_conditions = add_closure_specific_boundary_conditions(closure, boundary_conditions, grid, tracernames(tracers), buoyancy) # Ensure `closure` describes all tracers diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_fields.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_fields.jl deleted file mode 100644 index e3c7375ec7..0000000000 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_fields.jl +++ /dev/null @@ -1,15 +0,0 @@ -using Oceananigans.BoundaryConditions: UVelocityBoundaryConditions, VVelocityBoundaryConditions, TracerBoundaryConditions -using Oceananigans.Fields: XFaceField, YFaceField, CenterField, TracerFields - -function HorizontalVelocityFields(velocities, arch, grid) - u = XFaceField(arch, grid, UVelocityBoundaryConditions(grid)) - v = YFaceField(arch, grid, VVelocityBoundaryConditions(grid)) - return u, v -end - -function HydrostaticFreeSurfaceTendencyFields(velocities, free_surface, arch, grid, tracer_names) - u, v = HorizontalVelocityFields(velocities, arch, grid) - η = FreeSurfaceDisplacementField(velocities, free_surface, arch, grid) - tracers = TracerFields(tracer_names, arch, grid) - return merge((u=u, v=v, η=η), tracers) -end diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_velocity_fields.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_velocity_fields.jl deleted file mode 100644 index c924040628..0000000000 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_velocity_fields.jl +++ /dev/null @@ -1,18 +0,0 @@ -using Oceananigans.Grids: Center, Face -using Oceananigans.Fields: XFaceField, YFaceField, ZFaceField -using Oceananigans.BoundaryConditions: UVelocityBoundaryConditions, VVelocityBoundaryConditions, FieldBoundaryConditions - -function HydrostaticFreeSurfaceVelocityFields(::Nothing, arch, grid, clock, bcs=NamedTuple()) - u_bcs = :u ∈ keys(bcs) ? bcs.u : UVelocityBoundaryConditions(grid) - v_bcs = :v ∈ keys(bcs) ? bcs.v : VVelocityBoundaryConditions(grid) - w_bcs = :w ∈ keys(bcs) ? bcs.w : FreeSurfaceWVelocityBoundaryConditions(grid) - - u = XFaceField(arch, grid, u_bcs) - v = YFaceField(arch, grid, v_bcs) - w = ZFaceField(arch, grid, w_bcs) - - return (u=u, v=v, w=w) -end - -FreeSurfaceWVelocityBoundaryConditions(grid; user_defined_bcs...) = - FieldBoundaryConditions(grid, (Center, Center, Face); top=nothing, user_defined_bcs...) diff --git a/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl b/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl index 6f5c342002..2f373d39e6 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl @@ -52,6 +52,11 @@ function HydrostaticFreeSurfaceVelocityFields(velocities::PrescribedVelocityFiel return PrescribedVelocityFields(u, v, w, velocities.parameters) end +function HydrostaticFreeSurfaceTendencyFields(::PrescribedVelocityFields, free_surface, arch, grid, tracer_names) + tracers = TracerFields(tracer_names, arch, grid) + return merge((u = nothing, v = nothing, η = nothing), tracers) +end + @inline fill_halo_regions!(::PrescribedVelocityFields, args...) = nothing @inline fill_halo_regions!(::FunctionField, args...) = nothing diff --git a/src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity_field.jl b/src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity_field.jl index 27a8ae88a3..808fe1dc3c 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity_field.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity_field.jl @@ -1,7 +1,6 @@ using Oceananigans.Grids: Face, Face, Center using Oceananigans.Operators: ζ₃ᶠᶠᵃ -using Oceananigans.Fields: KernelComputedField, ComputedFieldBoundaryConditions - +using Oceananigans.Fields: KernelComputedField using KernelAbstractions: @kernel @kernel function compute_vertical_vorticity!(ζ, grid, u, v) @@ -21,8 +20,5 @@ the vertical vorticity by first integrating the velocity field around the border of the vorticity cell to find the vertical circulation, and then dividing by the area of the vorticity cell to compute vertical vorticity. """ -function VerticalVorticityField(model) - return KernelComputedField(Face, Face, Center, compute_vertical_vorticity!, model, - boundary_conditions = ComputedFieldBoundaryConditions(model.grid, (Face, Face, Center)), - computed_dependencies = (model.velocities.u, model.velocities.v)) -end +VerticalVorticityField(model) = KernelComputedField(Face, Face, Center, compute_vertical_vorticity!, model, + computed_dependencies = (model.velocities.u, model.velocities.v)) diff --git a/src/Models/IncompressibleModels/incompressible_model.jl b/src/Models/IncompressibleModels/incompressible_model.jl index 84512ef656..46fe742c1a 100644 --- a/src/Models/IncompressibleModels/incompressible_model.jl +++ b/src/Models/IncompressibleModels/incompressible_model.jl @@ -22,24 +22,24 @@ const ParticlesOrNothing = Union{Nothing, LagrangianParticles} mutable struct IncompressibleModel{TS, E, A<:AbstractArchitecture, G, T, B, R, SD, U, C, Φ, F, V, S, K, BG, P, I} <: AbstractModel{TS} - architecture :: A # Computer `Architecture` on which `Model` is run - grid :: G # Grid of physical points on which `Model` is solved - clock :: Clock{T} # Tracks iteration number and simulation time of `Model` - advection :: V # Advection scheme for velocities _and_ tracers - buoyancy :: B # Set of parameters for buoyancy model - coriolis :: R # Set of parameters for the background rotation rate of `Model` - stokes_drift :: SD # Set of parameters for surfaces waves via the Craik-Leibovich approximation - forcing :: F # Container for forcing functions defined by the user - closure :: E # Diffusive 'turbulence closure' for all model fields - background_fields :: BG # Background velocity and tracer fields - particles :: P # Particle set for Lagrangian tracking - velocities :: U # Container for velocity fields `u`, `v`, and `w` - tracers :: C # Container for tracer fields - pressures :: Φ # Container for hydrostatic and nonhydrostatic pressure - diffusivities :: K # Container for turbulent diffusivities - timestepper :: TS # Object containing timestepper fields and parameters - pressure_solver :: S # Pressure/Poisson solver - immersed_boundary :: I # Models the physics of immersed boundaries within the grid + architecture :: A # Computer `Architecture` on which `Model` is run + grid :: G # Grid of physical points on which `Model` is solved + clock :: Clock{T} # Tracks iteration number and simulation time of `Model` + advection :: V # Advection scheme for velocities _and_ tracers + buoyancy :: B # Set of parameters for buoyancy model + coriolis :: R # Set of parameters for the background rotation rate of `Model` + stokes_drift :: SD # Set of parameters for surfaces waves via the Craik-Leibovich approximation + forcing :: F # Container for forcing functions defined by the user + closure :: E # Diffusive 'turbulence closure' for all model fields + background_fields :: BG # Background velocity and tracer fields + particles :: P # Particle set for Lagrangian tracking + velocities :: U # Container for velocity fields `u`, `v`, and `w` + tracers :: C # Container for tracer fields + pressures :: Φ # Container for hydrostatic and nonhydrostatic pressure + diffusivities :: K # Container for turbulent diffusivities + timestepper :: TS # Object containing timestepper fields and parameters + pressure_solver :: S # Pressure/Poisson solver + immersed_boundary :: I # Models the physics of immersed boundaries within the grid end """ @@ -109,8 +109,8 @@ function IncompressibleModel(; grid, end tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example) - validate_buoyancy(buoyancy, tracernames(tracers)) + validate_buoyancy(buoyancy, tracernames(tracers)) buoyancy = regularize_buoyancy(buoyancy) # Adjust halos when the advection scheme or turbulence closure requires it. @@ -119,18 +119,26 @@ function IncompressibleModel(; grid, Hx, Hy, Hz = inflate_halo_size(grid.Hx, grid.Hy, grid.Hz, topology(grid), advection, closure) grid = with_halo((Hx, Hy, Hz), grid) - # Recursively "regularize" field-dependent boundary conditions by supplying list of tracer names. - # We also regularize boundary conditions included in velocities, tracers, pressures, and diffusivities. + # Collect boundary conditions for all model prognostic fields and, if specified, some model + # auxiliary fields. Boundary conditions are "regularized" based on the _name_ of the field: + # boundary conditions on u, v, w are regularized assuming they represent momentum at appropriate + # staggered locations. All other fields are regularized assuming they are tracers. # Note that we do not regularize boundary conditions contained in *tupled* diffusivity fields right now. + # + # 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(diffusivities)) - boundary_conditions = merge(embedded_boundary_conditions, boundary_conditions) + # Next, we form a list of default boundary conditions: + prognostic_field_names = (:u, :v, :w, tracernames(tracers)...) + default_boundary_conditions = NamedTuple{prognostic_field_names}(FieldBoundaryConditions() for name in prognostic_field_names) - model_field_names = (:u, :v, :w, tracernames(tracers)...) - boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, model_field_names) + # Finally, we merge specified, embedded, and default boundary conditions. Specified boundary conditions + # have precedence, followed by embedded, followed by default. + boundary_conditions = merge(default_boundary_conditions, embedded_boundary_conditions, boundary_conditions) + boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, prognostic_field_names) # Ensure `closure` describes all tracers closure = with_tracers(tracernames(tracers), closure) @@ -145,6 +153,7 @@ function IncompressibleModel(; grid, pressure_solver = PressureSolver(architecture, grid) end + # Materialize background fields background_fields = BackgroundFields(background_fields, tracernames(tracers), grid, clock) # Instantiate timestepper if not already instantiated diff --git a/src/Models/ShallowWaterModels/shallow_water_model.jl b/src/Models/ShallowWaterModels/shallow_water_model.jl index 80f7ace575..a3610c8e86 100644 --- a/src/Models/ShallowWaterModels/shallow_water_model.jl +++ b/src/Models/ShallowWaterModels/shallow_water_model.jl @@ -3,11 +3,6 @@ using Oceananigans: AbstractModel, AbstractOutputWriter, AbstractDiagnostic using Oceananigans.Architectures: AbstractArchitecture, CPU using Oceananigans.Advection: CenteredSecondOrder using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions - -using Oceananigans.BoundaryConditions: UVelocityBoundaryConditions, - VVelocityBoundaryConditions, - TracerBoundaryConditions - using Oceananigans.Fields: Field, tracernames, TracerFields, XFaceField, YFaceField, CenterField using Oceananigans.Forcings: model_forcing using Oceananigans.Grids: with_halo, topology, inflate_halo_size, halo_size, Flat @@ -17,23 +12,19 @@ using Oceananigans.Utils: tupleit function ShallowWaterTendencyFields(arch, grid, tracer_names) - uh = XFaceField(arch, grid, UVelocityBoundaryConditions(grid)) - vh = YFaceField(arch, grid, VVelocityBoundaryConditions(grid)) - h = CenterField(arch, grid, TracerBoundaryConditions(grid)) + uh = XFaceField(arch, grid) + vh = YFaceField(arch, grid) + h = CenterField(arch, grid) + tracers = TracerFields(tracer_names, arch, grid) return merge((uh=uh, vh=vh, h=h), tracers) end function ShallowWaterSolutionFields(arch, grid, bcs) - - uh_bcs = :uh ∈ keys(bcs) ? bcs.uh : UVelocityBoundaryConditions(grid) - vh_bcs = :vh ∈ keys(bcs) ? bcs.vh : VVelocityBoundaryConditions(grid) - h_bcs = :h ∈ keys(bcs) ? bcs.h : TracerBoundaryConditions(grid) - - uh = XFaceField(arch, grid, uh_bcs) - vh = YFaceField(arch, grid, vh_bcs) - h = CenterField(arch, grid, h_bcs) + uh = XFaceField(arch, grid, bcs.uh) + vh = YFaceField(arch, grid, bcs.vh) + h = CenterField(arch, grid, bcs.h) return (uh=uh, vh=vh, h=h) end @@ -114,8 +105,11 @@ function ShallowWaterModel(; Hx, Hy, Hz = inflate_halo_size(grid.Hx, grid.Hy, 0, topology(grid), advection, closure) grid = with_halo((Hx, Hy, 0), grid) - model_field_names = (:uh, :vh, :h, tracers...) - boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, model_field_names) + prognostic_field_names = (:uh, :vh, :h, tracers...) + default_boundary_conditions = NamedTuple{prognostic_field_names}(Tuple(FieldBoundaryConditions() for name in prognostic_field_names)) + boundary_conditions = merge(default_boundary_conditions, boundary_conditions) + + boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, prognostic_field_names) solution = ShallowWaterSolutionFields(architecture, grid, boundary_conditions) tracers = TracerFields(tracers, architecture, grid, boundary_conditions) diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index e801e7f5fd..a3bdf87a90 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -24,9 +24,7 @@ export # Boundary conditions BoundaryCondition, FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, OpenBoundaryCondition, - CoordinateBoundaryConditions, FieldBoundaryConditions, - UVelocityBoundaryConditions, VVelocityBoundaryConditions, WVelocityBoundaryConditions, - TracerBoundaryConditions, PressureBoundaryConditions, + FieldBoundaryConditions, # Fields and field manipulation Field, CenterField, XFaceField, YFaceField, ZFaceField, diff --git a/src/OutputWriters/output_writer_utils.jl b/src/OutputWriters/output_writer_utils.jl index ea28ab5f0d..35880a979f 100644 --- a/src/OutputWriters/output_writer_utils.jl +++ b/src/OutputWriters/output_writer_utils.jl @@ -1,7 +1,7 @@ using StructArrays: StructArray, replace_storage using Oceananigans.Fields: AbstractField using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid -using Oceananigans.BoundaryConditions: bcclassification_str, CoordinateBoundaryConditions, FieldBoundaryConditions +using Oceananigans.BoundaryConditions: bc_str, FieldBoundaryConditions using Oceananigans.TimeSteppers: QuasiAdamsBashforth2TimeStepper, RungeKutta3TimeStepper using Oceananigans.LagrangianParticleTracking: LagrangianParticles @@ -30,16 +30,16 @@ saveproperty!(file, location, p) = saveproperty!(file, location, p::ImmersedBoundaryGrid) = saveproperty!(file, location, p.grid) # Special saveproperty! so boundary conditions are easily readable outside julia. -function saveproperty!(file, location, cbcs::CoordinateBoundaryConditions) - for endpoint in propertynames(cbcs) - endpoint_bc = getproperty(cbcs, endpoint) - file[location * "/$endpoint/type"] = bcclassification_str(endpoint_bc) +function saveproperty!(file, location, bcs::FieldBoundaryConditions) + for boundary in propertynames(bcs) + bc = getproperty(bcs, endpoint) + file[location * "/$endpoint/type"] = bc_str(bc) - if endpoint_bc.condition isa Function + if bc.condition isa Function @warn "$field.$coord.$endpoint boundary is of type Function and cannot be saved to disk!" - file[location * "/$endpoint/condition"] = missing + file[location * "/$boundary/condition"] = missing else - file[location * "/$endpoint/condition"] = endpoint_bc.condition + file[location * "/$boundary/condition"] = bc.condition end end end diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl b/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl index 5254a7651e..4b3461d1d8 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl @@ -43,18 +43,24 @@ Keyword arguments ================= - `C` : Poincaré constant for both eddy viscosity and eddy diffusivities. `C` is overridden for eddy viscosity or eddy diffusivity if `Cν` or `Cκ` are set, respecitvely. + - `Cν` : Poincaré constant for momentum eddy viscosity. + - `Cκ` : Poincaré constant for tracer eddy diffusivities. If one number or function, the same number or function is applied to all tracers. If a `NamedTuple`, it must possess a field specifying the Poncaré constant for every tracer. + - `Cb` : Buoyancy modification multiplier (`Cb = nothing` turns it off, `Cb = 1` was used by [Abkar16](@cite)). *Note*: that we _do not_ subtract the horizontally-average component before computing this buoyancy modification term. This implementation differs from [Abkar16](@cite)'s proposal and the impact of this approximation has not been tested or validated. + - `ν` : Constant background viscosity for momentum. + - `κ` : Constant background diffusivity for tracer. If a single number, the same background diffusivity is applied to all tracers. If a `NamedTuple`, it must possess a field specifying a background diffusivity for every tracer. + - `time_discretization` : Either `ExplicitTimeDiscretization()` or `VerticallyImplicitTimeDiscretization()`, which integrates the terms involving only z-derivatives in the viscous and diffusive fluxes with an implicit time discretization. @@ -344,51 +350,25 @@ end return cx_ux + cy_uy + cz_uz end -@inline norm_θᵢ²ᶜᶜᶜ(i, j, k, grid, c) = ( - ℑxᶜᵃᵃ(i, j, k, grid, norm_∂x_c², c) - + ℑyᵃᶜᵃ(i, j, k, grid, norm_∂y_c², c) - + ℑzᵃᵃᶜ(i, j, k, grid, norm_∂z_c², c) -) +@inline norm_θᵢ²ᶜᶜᶜ(i, j, k, grid, c) = ℑxᶜᵃᵃ(i, j, k, grid, norm_∂x_c², c) + + ℑyᵃᶜᵃ(i, j, k, grid, norm_∂y_c², c) + + ℑzᵃᵃᶜ(i, j, k, grid, norm_∂z_c², c) ##### ##### DiffusivityFields ##### -function DiffusivityFields(arch, grid, tracer_names, ::AMD; - νₑ = CenterField(arch, grid, DiffusivityBoundaryConditions(grid)), - kwargs...) - - κₑ = TracerDiffusivityFields(arch, grid, tracer_names; kwargs...) - - return (νₑ=νₑ, κₑ=κₑ) -end - -function DiffusivityFields(arch, grid, tracer_names, bcs, ::AMD) - - νₑ_bcs = :νₑ ∈ keys(bcs) ? bcs[:νₑ] : DiffusivityBoundaryConditions(grid) - νₑ = CenterField(arch, grid, νₑ_bcs) - - κₑ = :κₑ ∈ keys(bcs) ? TracerDiffusivityFields(arch, grid, tracer_names, bcs[:κₑ]) : - TracerDiffusivityFields(arch, grid, tracer_names) +function DiffusivityFields(arch, grid, tracer_names, user_bcs, ::AMD) - return (νₑ=νₑ, κₑ=κₑ) -end + default_diffusivity_bcs = FieldBoundaryConditions(grid, (Center, Center, Center)) + default_κₑ_bcs = NamedTuple(c => default_diffusivity_bcs for c in tracer_names) + κₑ_bcs = :κₑ ∈ keys(user_bcs) ? merge(default_κₑ_bcs, user_bcs.κₑ) : default_κₑ_bcs -function TracerDiffusivityFields(arch, grid, tracer_names; kwargs...) + bcs = merge((; νₑ = default_diffusivity_bcs, κₑ = κₑ_bcs), user_bcs) - κ_fields = Tuple(c ∈ keys(kwargs) ? - kwargs[c] : - CenterField(arch, grid, DiffusivityBoundaryConditions(grid)) - for c in tracer_names) + νₑ = CenterField(arch, grid, bcs.νₑ) + κₑ = NamedTuple(c => CenterField(arch, grid, bcs.κₑ[c]) for c in tracer_names) - return NamedTuple{tracer_names}(κ_fields) + return (; νₑ, κₑ) end -function TracerDiffusivityFields(arch, grid, tracer_names, bcs) - - κ_fields = Tuple(c ∈ keys(bcs) ? CenterField(arch, grid, bcs[c]) : - CenterField(arch, grid, DiffusivityBoundaryConditions(grid)) - for c in tracer_names) - - return NamedTuple{tracer_names}(κ_fields) -end diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/leith_enstrophy_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/leith_enstrophy_diffusivity.jl index 95b155e771..6c3b2f4934 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/leith_enstrophy_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/leith_enstrophy_diffusivity.jl @@ -210,7 +210,8 @@ end @inline Δᶠ(i, j, k, grid::RegularRectilinearGrid, ::TwoDimensionalLeith) = sqrt(grid.Δx * grid.Δy) function DiffusivityFields(arch, grid, tracer_names, bcs, ::L2D) - νₑ_bcs = :νₑ ∈ keys(bcs) ? bcs[:νₑ] : DiffusivityBoundaryConditions(grid) - νₑ = CenterField(arch, grid, νₑ_bcs) - return (νₑ = νₑ,) + default_eddy_viscosity_bcs = (; νₑ = FieldBoundaryConditions(grid, (Center, Center, Center))) + bcs = merge(default_eddy_viscosity_bcs, bcs) + νₑ = CenterField(arch, grid, bcs.νₑ) + return (; νₑ) end diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl b/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl index 70f17f91fc..70567ebce5 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl @@ -188,9 +188,10 @@ Base.show(io::IO, closure::SmagorinskyLilly) = ##### function DiffusivityFields(arch, grid, tracer_names, bcs, closure::SmagorinskyLilly) - νₑ_bcs = :νₑ ∈ keys(bcs) ? bcs[:νₑ] : DiffusivityBoundaryConditions(grid) - νₑ = CenterField(arch, grid, νₑ_bcs) + default_eddy_viscosity_bcs = (; νₑ = FieldBoundaryConditions(grid, (Center, Center, Center))) + bcs = merge(default_eddy_viscosity_bcs, bcs) + νₑ = CenterField(arch, grid, bcs.νₑ) # Use AbstractOperations to write eddy diffusivities in terms of # eddy viscosity @@ -206,6 +207,6 @@ function DiffusivityFields(arch, grid, tracer_names, bcs, closure::SmagorinskyLi κₑ = NamedTuple{tracer_names}(Tuple(κₑ_ops)) - return (νₑ=νₑ, κₑ=κₑ) + return (; νₑ, κₑ) end diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/tke_based_vertical_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/tke_based_vertical_diffusivity.jl index 6e5d4d885b..e4e84dd5ba 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/tke_based_vertical_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/tke_based_vertical_diffusivity.jl @@ -1,5 +1,5 @@ using Oceananigans.Architectures: architecture, device_event -using Oceananigans.BoundaryConditions: DefaultBoundaryCondition +using Oceananigans.BoundaryConditions: default_prognostic_field_boundary_condition, FieldBoundaryConditions using Oceananigans.BuoyancyModels: ∂z_b, top_buoyancy_flux using Oceananigans.Operators: ℑzᵃᵃᶜ @@ -243,22 +243,19 @@ end """ Infer tracer boundary conditions from user_bcs and tracer_names. """ function top_tracer_boundary_conditions(grid, tracer_names, user_bcs) - user_bc_names = keys(user_bcs) - default_top_bc = DefaultBoundaryCondition(topology(grid, 3), Center) - - tracer_bcs = Tuple(name ∈ user_bc_names ? user_bcs[name].top : default_top_bc - for name in tracer_names) - - return NamedTuple{tracer_names}(tracer_bcs) + default_tracer_bcs = NamedTuple(c => FieldBoundaryConditions(grid, (Center, Center, Center)) for c in tracer_names) + bcs = merge(default_tracer_bcs, user_bcs) + return NamedTuple(c => bcs[c].top for c in tracer_names) end """ Infer velocity boundary conditions from `user_bcs` and `tracer_names`. """ function top_velocity_boundary_conditions(grid, user_bcs) - user_bc_names = keys(user_bcs) + default_top_bc = default_prognostic_field_boundary_condition(topology(grid, 3)(), Center()) - u_top_bc = :u ∈ user_bc_names ? user_bcs.u.top : DefaultBoundaryCondition(topology(grid, 3), Center) - v_top_bc = :v ∈ user_bc_names ? user_bcs.v.top : DefaultBoundaryCondition(topology(grid, 3), Center) + user_bc_names = keys(user_bcs) + u_top_bc = :u ∈ user_bc_names ? user_bcs.u.top : default_top_bc + v_top_bc = :v ∈ user_bc_names ? user_bcs.v.top : default_top_bc return (u=u_top_bc, v=v_top_bc) end @@ -286,16 +283,15 @@ function add_closure_specific_boundary_conditions(closure::TKEVD, e_bcs = user_bcs[:e] - tke_bcs = TracerBoundaryConditions(grid, - top = top_tke_bc, - bottom = e_bcs.bottom, - north = e_bcs.north, - south = e_bcs.south, - east = e_bcs.east, - west = e_bcs.west) - + tke_bcs = FieldBoundaryConditions(grid, (Center, Center, Center), + top = top_tke_bc, + bottom = e_bcs.bottom, + north = e_bcs.north, + south = e_bcs.south, + east = e_bcs.east, + west = e_bcs.west) else - tke_bcs = TracerBoundaryConditions(grid, top=top_tke_bc) + tke_bcs = FieldBoundaryConditions(grid, (Center, Center, Center), top=top_tke_bc) end new_boundary_conditions = merge(user_bcs, (e = tke_bcs,)) @@ -305,13 +301,15 @@ end function DiffusivityFields(arch, grid, tracer_names, bcs, closure::TKEVD) - Kᵘ_bcs = :Kᵘ ∈ keys(bcs) ? bcs[:Kᵘ] : DiffusivityBoundaryConditions(grid) - Kᶜ_bcs = :Kᶜ ∈ keys(bcs) ? bcs[:Kᶜ] : DiffusivityBoundaryConditions(grid) - Kᵉ_bcs = :Kᵉ ∈ keys(bcs) ? bcs[:Kᵉ] : DiffusivityBoundaryConditions(grid) + default_diffusivity_bcs = (Kᵘ = FieldBoundaryConditions(grid, (Center, Center, Center)), + Kᶜ = FieldBoundaryConditions(grid, (Center, Center, Center)), + Kᵉ = FieldBoundaryConditions(grid, (Center, Center, Center))) + + bcs = merge(default_diffusivity_bcs, bcs) - Kᵘ = CenterField(arch, grid, Kᵘ_bcs) - Kᶜ = CenterField(arch, grid, Kᶜ_bcs) - Kᵉ = CenterField(arch, grid, Kᵉ_bcs) + Kᵘ = CenterField(arch, grid, bcs.Kᵘ) + Kᶜ = CenterField(arch, grid, bcs.Kᶜ) + Kᵉ = CenterField(arch, grid, bcs.Kᵉ) return (; Kᵘ, Kᶜ, Kᵉ) end 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 99d3268687..71fc40d3de 100644 --- a/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl +++ b/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl @@ -23,10 +23,9 @@ function run_ocean_large_eddy_simulation_regression_test(arch, grid_type, closur end # Boundary conditions - u_bcs = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ)) - T_bcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵀ), - bottom = BoundaryCondition(Gradient, ∂T∂z)) - S_bcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, 5e-8)) + u_bcs = FieldBoundaryConditions(top = BoundaryCondition(Flux, Qᵘ)) + T_bcs = FieldBoundaryConditions(top = BoundaryCondition(Flux, Qᵀ), bottom = BoundaryCondition(Gradient, ∂T∂z)) + S_bcs = FieldBoundaryConditions(top = BoundaryCondition(Flux, 5e-8)) # Model instantiation model = IncompressibleModel( diff --git a/test/regression_tests/rayleigh_benard_regression_test.jl b/test/regression_tests/rayleigh_benard_regression_test.jl index c4a6d16bf7..8028169dd0 100644 --- a/test/regression_tests/rayleigh_benard_regression_test.jl +++ b/test/regression_tests/rayleigh_benard_regression_test.jl @@ -36,8 +36,8 @@ function run_rayleigh_benard_regression_test(arch, grid_type) c★(x, z) = exp(4z) * sin(2π/Lx * x) Fc(i, j, k, grid, clock, model_fields) = 1/10 * (c★(xnode(Center(), i, grid), znode(Center(), k, grid)) - model_fields.c[i, j, k]) - bbcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Value, 0.0), - bottom = BoundaryCondition(Value, Δb)) + bbcs = FieldBoundaryConditions(top = BoundaryCondition(Value, 0.0), + bottom = BoundaryCondition(Value, Δb)) model = IncompressibleModel( architecture = arch, diff --git a/test/test_boundary_conditions.jl b/test/test_boundary_conditions.jl index 84f9d3ee78..c733ec04d1 100644 --- a/test/test_boundary_conditions.jl +++ b/test/test_boundary_conditions.jl @@ -1,5 +1,4 @@ -using Oceananigans.BoundaryConditions: PBC, ZFBC, OBC, ContinuousBoundaryFunction, DiscreteBoundaryFunction - +using Oceananigans.BoundaryConditions: PBC, ZFBC, OBC, ContinuousBoundaryFunction, DiscreteBoundaryFunction, regularize_field_boundary_conditions using Oceananigans.Fields: Face, Center simple_bc(ξ, η, t) = exp(ξ) * cos(η) * sin(t) @@ -48,191 +47,202 @@ end ppp_topology = (Periodic, Periodic, Periodic) ppp_grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=ppp_topology) - u_bcs = UVelocityBoundaryConditions(ppp_grid) - v_bcs = VVelocityBoundaryConditions(ppp_grid) - w_bcs = WVelocityBoundaryConditions(ppp_grid) - T_bcs = TracerBoundaryConditions(ppp_grid) + default_bcs = FieldBoundaryConditions() + + u_bcs = regularize_field_boundary_conditions(default_bcs, ppp_grid, :u) + v_bcs = regularize_field_boundary_conditions(default_bcs, ppp_grid, :v) + w_bcs = regularize_field_boundary_conditions(default_bcs, ppp_grid, :w) + T_bcs = regularize_field_boundary_conditions(default_bcs, ppp_grid, :T) @test u_bcs isa FieldBoundaryConditions - @test u_bcs.x.left isa PBC - @test u_bcs.x.right isa PBC - @test u_bcs.y.left isa PBC - @test u_bcs.y.right isa PBC - @test u_bcs.z.left isa PBC - @test u_bcs.z.right isa PBC + @test u_bcs.west isa PBC + @test u_bcs.east isa PBC + @test u_bcs.south isa PBC + @test u_bcs.north isa PBC + @test u_bcs.bottom isa PBC + @test u_bcs.top isa PBC @test v_bcs isa FieldBoundaryConditions - @test v_bcs.x.left isa PBC - @test v_bcs.x.right isa PBC - @test v_bcs.y.left isa PBC - @test v_bcs.y.right isa PBC - @test v_bcs.z.left isa PBC - @test v_bcs.z.right isa PBC + @test v_bcs.west isa PBC + @test v_bcs.east isa PBC + @test v_bcs.south isa PBC + @test v_bcs.north isa PBC + @test v_bcs.bottom isa PBC + @test v_bcs.top isa PBC @test w_bcs isa FieldBoundaryConditions - @test w_bcs.x.left isa PBC - @test w_bcs.x.right isa PBC - @test w_bcs.y.left isa PBC - @test w_bcs.y.right isa PBC - @test w_bcs.z.left isa PBC - @test w_bcs.z.right isa PBC + @test w_bcs.west isa PBC + @test w_bcs.east isa PBC + @test w_bcs.south isa PBC + @test w_bcs.north isa PBC + @test w_bcs.bottom isa PBC + @test w_bcs.top isa PBC @test T_bcs isa FieldBoundaryConditions - @test T_bcs.x.left isa PBC - @test T_bcs.x.right isa PBC - @test T_bcs.y.left isa PBC - @test T_bcs.y.right isa PBC - @test T_bcs.z.left isa PBC - @test T_bcs.z.right isa PBC + @test T_bcs.west isa PBC + @test T_bcs.east isa PBC + @test T_bcs.south isa PBC + @test T_bcs.north isa PBC + @test T_bcs.bottom isa PBC + @test T_bcs.top isa PBC # Doubly periodic. Engineers call this a "Channel geometry". ppb_topology = (Periodic, Periodic, Bounded) ppb_grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=ppb_topology) - u_bcs = UVelocityBoundaryConditions(ppb_grid) - v_bcs = VVelocityBoundaryConditions(ppb_grid) - w_bcs = WVelocityBoundaryConditions(ppb_grid) - T_bcs = TracerBoundaryConditions(ppb_grid) + u_bcs = regularize_field_boundary_conditions(default_bcs, ppb_grid, :u) + v_bcs = regularize_field_boundary_conditions(default_bcs, ppb_grid, :v) + w_bcs = regularize_field_boundary_conditions(default_bcs, ppb_grid, :w) + T_bcs = regularize_field_boundary_conditions(default_bcs, ppb_grid, :T) @test u_bcs isa FieldBoundaryConditions - @test u_bcs.x.left isa PBC - @test u_bcs.x.right isa PBC - @test u_bcs.y.left isa PBC - @test u_bcs.y.right isa PBC - @test u_bcs.z.left isa ZFBC - @test u_bcs.z.right isa ZFBC + @test u_bcs.west isa PBC + @test u_bcs.east isa PBC + @test u_bcs.south isa PBC + @test u_bcs.north isa PBC + @test u_bcs.bottom isa ZFBC + @test u_bcs.top isa ZFBC @test v_bcs isa FieldBoundaryConditions - @test v_bcs.x.left isa PBC - @test v_bcs.x.right isa PBC - @test v_bcs.y.left isa PBC - @test v_bcs.y.right isa PBC - @test v_bcs.z.left isa ZFBC - @test v_bcs.z.right isa ZFBC + @test v_bcs.west isa PBC + @test v_bcs.east isa PBC + @test v_bcs.south isa PBC + @test v_bcs.north isa PBC + @test v_bcs.bottom isa ZFBC + @test v_bcs.top isa ZFBC @test w_bcs isa FieldBoundaryConditions - @test w_bcs.x.left isa PBC - @test w_bcs.x.right isa PBC - @test w_bcs.y.left isa PBC - @test w_bcs.y.right isa PBC - @test w_bcs.z.left isa OBC - @test w_bcs.z.right isa OBC + @test w_bcs.west isa PBC + @test w_bcs.east isa PBC + @test w_bcs.south isa PBC + @test w_bcs.north isa PBC + @test w_bcs.bottom isa OBC + @test w_bcs.top isa OBC @test T_bcs isa FieldBoundaryConditions - @test T_bcs.x.left isa PBC - @test T_bcs.x.right isa PBC - @test T_bcs.y.left isa PBC - @test T_bcs.y.right isa PBC - @test T_bcs.z.left isa ZFBC - @test T_bcs.z.right isa ZFBC + @test T_bcs.west isa PBC + @test T_bcs.east isa PBC + @test T_bcs.south isa PBC + @test T_bcs.north isa PBC + @test T_bcs.bottom isa ZFBC + @test T_bcs.top isa ZFBC # Singly periodic. Oceanographers call this a "Channel", engineers call it a "Pipe" pbb_topology = (Periodic, Bounded, Bounded) pbb_grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=pbb_topology) - u_bcs = UVelocityBoundaryConditions(pbb_grid) - v_bcs = VVelocityBoundaryConditions(pbb_grid) - w_bcs = WVelocityBoundaryConditions(pbb_grid) - T_bcs = TracerBoundaryConditions(pbb_grid) + u_bcs = regularize_field_boundary_conditions(default_bcs, pbb_grid, :u) + v_bcs = regularize_field_boundary_conditions(default_bcs, pbb_grid, :v) + w_bcs = regularize_field_boundary_conditions(default_bcs, pbb_grid, :w) + T_bcs = regularize_field_boundary_conditions(default_bcs, pbb_grid, :T) @test u_bcs isa FieldBoundaryConditions - @test u_bcs.x.left isa PBC - @test u_bcs.x.right isa PBC - @test u_bcs.y.left isa ZFBC - @test u_bcs.y.right isa ZFBC - @test u_bcs.z.left isa ZFBC - @test u_bcs.z.right isa ZFBC + @test u_bcs.west isa PBC + @test u_bcs.east isa PBC + @test u_bcs.south isa ZFBC + @test u_bcs.north isa ZFBC + @test u_bcs.bottom isa ZFBC + @test u_bcs.top isa ZFBC @test v_bcs isa FieldBoundaryConditions - @test v_bcs.x.left isa PBC - @test v_bcs.x.right isa PBC - @test v_bcs.y.left isa OBC - @test v_bcs.y.right isa OBC - @test v_bcs.z.left isa ZFBC - @test v_bcs.z.right isa ZFBC + @test v_bcs.west isa PBC + @test v_bcs.east isa PBC + @test v_bcs.south isa OBC + @test v_bcs.north isa OBC + @test v_bcs.bottom isa ZFBC + @test v_bcs.top isa ZFBC @test w_bcs isa FieldBoundaryConditions - @test w_bcs.x.left isa PBC - @test w_bcs.x.right isa PBC - @test w_bcs.y.left isa ZFBC - @test w_bcs.y.right isa ZFBC - @test w_bcs.z.left isa OBC - @test w_bcs.z.right isa OBC + @test w_bcs.west isa PBC + @test w_bcs.east isa PBC + @test w_bcs.south isa ZFBC + @test w_bcs.north isa ZFBC + @test w_bcs.bottom isa OBC + @test w_bcs.top isa OBC @test T_bcs isa FieldBoundaryConditions - @test T_bcs.x.left isa PBC - @test T_bcs.x.right isa PBC - @test T_bcs.y.left isa ZFBC - @test T_bcs.y.right isa ZFBC - @test T_bcs.z.left isa ZFBC - @test T_bcs.z.right isa ZFBC + @test T_bcs.west isa PBC + @test T_bcs.east isa PBC + @test T_bcs.south isa ZFBC + @test T_bcs.north isa ZFBC + @test T_bcs.bottom isa ZFBC + @test T_bcs.top isa ZFBC # Triply bounded. Oceanographers call this a "Basin", engineers call it a "Box" bbb_topology = (Bounded, Bounded, Bounded) bbb_grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=bbb_topology) - u_bcs = UVelocityBoundaryConditions(bbb_grid) - v_bcs = VVelocityBoundaryConditions(bbb_grid) - w_bcs = WVelocityBoundaryConditions(bbb_grid) - T_bcs = TracerBoundaryConditions(bbb_grid) + u_bcs = regularize_field_boundary_conditions(default_bcs, bbb_grid, :u) + v_bcs = regularize_field_boundary_conditions(default_bcs, bbb_grid, :v) + w_bcs = regularize_field_boundary_conditions(default_bcs, bbb_grid, :w) + T_bcs = regularize_field_boundary_conditions(default_bcs, bbb_grid, :T) @test u_bcs isa FieldBoundaryConditions - @test u_bcs.x.left isa OBC - @test u_bcs.x.right isa OBC - @test u_bcs.y.left isa ZFBC - @test u_bcs.y.right isa ZFBC - @test u_bcs.z.left isa ZFBC - @test u_bcs.z.right isa ZFBC + @test u_bcs.west isa OBC + @test u_bcs.east isa OBC + @test u_bcs.south isa ZFBC + @test u_bcs.north isa ZFBC + @test u_bcs.bottom isa ZFBC + @test u_bcs.top isa ZFBC @test v_bcs isa FieldBoundaryConditions - @test v_bcs.x.left isa ZFBC - @test v_bcs.x.right isa ZFBC - @test v_bcs.y.left isa OBC - @test v_bcs.y.right isa OBC - @test v_bcs.z.left isa ZFBC - @test v_bcs.z.right isa ZFBC + @test v_bcs.west isa ZFBC + @test v_bcs.east isa ZFBC + @test v_bcs.south isa OBC + @test v_bcs.north isa OBC + @test v_bcs.bottom isa ZFBC + @test v_bcs.top isa ZFBC @test w_bcs isa FieldBoundaryConditions - @test w_bcs.x.left isa ZFBC - @test w_bcs.x.right isa ZFBC - @test w_bcs.y.left isa ZFBC - @test w_bcs.y.right isa ZFBC - @test w_bcs.z.left isa OBC - @test w_bcs.z.right isa OBC + @test w_bcs.west isa ZFBC + @test w_bcs.east isa ZFBC + @test w_bcs.south isa ZFBC + @test w_bcs.north isa ZFBC + @test w_bcs.bottom isa OBC + @test w_bcs.top isa OBC @test T_bcs isa FieldBoundaryConditions - @test T_bcs.x.left isa ZFBC - @test T_bcs.x.right isa ZFBC - @test T_bcs.y.left isa ZFBC - @test T_bcs.y.right isa ZFBC - @test T_bcs.z.left isa ZFBC - @test T_bcs.z.right isa ZFBC + @test T_bcs.west isa ZFBC + @test T_bcs.east isa ZFBC + @test T_bcs.south isa ZFBC + @test T_bcs.north isa ZFBC + @test T_bcs.bottom isa ZFBC + @test T_bcs.top isa ZFBC grid = bbb_grid - T_bcs = TracerBoundaryConditions(grid; - east = BoundaryCondition(Value, simple_bc), - west = BoundaryCondition(Value, simple_bc), - bottom = BoundaryCondition(Value, simple_bc), - top = BoundaryCondition(Value, simple_bc), - north = BoundaryCondition(Value, simple_bc), - south = BoundaryCondition(Value, simple_bc) - ) - - # Note that boundary condition setting is valid only in cases where - # boundary conditionsare same type - + T_bcs = FieldBoundaryConditions(grid, (Center, Center, Center), + east = ValueBoundaryCondition(simple_bc), + west = ValueBoundaryCondition(simple_bc), + bottom = ValueBoundaryCondition(simple_bc), + top = ValueBoundaryCondition(simple_bc), + north = ValueBoundaryCondition(simple_bc), + south = ValueBoundaryCondition(simple_bc)) + + @test T_bcs.east.condition isa ContinuousBoundaryFunction + @test T_bcs.west.condition isa ContinuousBoundaryFunction + @test T_bcs.north.condition isa ContinuousBoundaryFunction + @test T_bcs.south.condition isa ContinuousBoundaryFunction + @test T_bcs.top.condition isa ContinuousBoundaryFunction + @test T_bcs.bottom.condition isa ContinuousBoundaryFunction + + @test T_bcs.east.condition.func === simple_bc + @test T_bcs.west.condition.func === simple_bc + @test T_bcs.north.condition.func === simple_bc + @test T_bcs.south.condition.func === simple_bc + @test T_bcs.top.condition.func === simple_bc + @test T_bcs.bottom.condition.func === simple_bc + one_bc = BoundaryCondition(Value, 1.0) - two_bc = BoundaryCondition(Value, 2.0) - T_bcs = TracerBoundaryConditions(grid; - east = one_bc, + T_bcs = FieldBoundaryConditions( east = one_bc, west = one_bc, bottom = one_bc, top = one_bc, north = one_bc, - south = one_bc, - ) + south = one_bc) + + T_bcs = regularize_field_boundary_conditions(T_bcs, grid, :T) @test T_bcs.east === one_bc @test T_bcs.west === one_bc @@ -240,19 +250,5 @@ end @test T_bcs.south === one_bc @test T_bcs.top === one_bc @test T_bcs.bottom === one_bc - - T_bcs.east = two_bc - T_bcs.west = two_bc - T_bcs.north = two_bc - T_bcs.south = two_bc - T_bcs.top = two_bc - T_bcs.bottom = two_bc - - @test T_bcs.east === two_bc - @test T_bcs.west === two_bc - @test T_bcs.north === two_bc - @test T_bcs.south === two_bc - @test T_bcs.top === two_bc - @test T_bcs.bottom === two_bc end end diff --git a/test/test_boundary_conditions_integration.jl b/test/test_boundary_conditions_integration.jl index 1bb653a1c1..55c3915ea5 100644 --- a/test/test_boundary_conditions_integration.jl +++ b/test/test_boundary_conditions_integration.jl @@ -3,12 +3,10 @@ using Oceananigans.BoundaryConditions: ContinuousBoundaryFunction function test_boundary_condition(arch, FT, topo, side, field_name, boundary_condition) grid = RegularRectilinearGrid(FT, size=(1, 1, 1), extent=(1, π, 42), topology=topo) - boundary_condition_kwarg = Dict(side => boundary_condition) - field_boundary_conditions = TracerBoundaryConditions(grid; boundary_condition_kwarg...) - bcs = NamedTuple{(field_name,)}((field_boundary_conditions,)) - - model = IncompressibleModel(grid=grid, architecture=arch, - boundary_conditions=bcs) + boundary_condition_kwarg = (; side => boundary_condition) + field_boundary_conditions = FieldBoundaryConditions(; boundary_condition_kwarg...) + bcs = (; field_name => field_boundary_conditions) + model = IncompressibleModel(grid=grid, architecture=arch, boundary_conditions=bcs) success = try time_step!(model, 1e-16, euler=true) @@ -39,13 +37,9 @@ function test_incompressible_flux_budget(arch, name, side, topo) direction = side ∈ (:west, :south, :bottom) ? 1 : -1 bc_kwarg = Dict(side => BoundaryCondition(Flux, flux * direction)) - FieldBoundaryConditions = name === :u ? UVelocityBoundaryConditions : - name === :v ? VVelocityBoundaryConditions : - name === :w ? WVelocityBoundaryConditions : TracerBoundaryConditions - - field_bcs = FieldBoundaryConditions(grid; bc_kwarg...) + field_bcs = FieldBoundaryConditions(; bc_kwarg...) - model_bcs = NamedTuple{(name,)}((field_bcs,)) + model_bcs = NamedTuple{tuple(name)}(tuple(field_bcs)) model = IncompressibleModel(grid=grid, buoyancy=nothing, boundary_conditions=model_bcs, closure=nothing, architecture=arch, tracers=:c) @@ -77,8 +71,8 @@ function fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT) grid = RegularRectilinearGrid(FT, size=(16, 16, 16), extent=(1, 1, Lz)) - buoyancy_bcs = TracerBoundaryConditions(grid, bottom=BoundaryCondition(Gradient, bz)) - κₑ_bcs = DiffusivityBoundaryConditions(grid, bottom=BoundaryCondition(Value, κ₀)) + buoyancy_bcs = FieldBoundaryConditions(bottom=GradientBoundaryCondition(bz)) + κₑ_bcs = FieldBoundaryConditions(grid, (Center, Center, Center), bottom=ValueBoundaryCondition(κ₀)) model_bcs = (b=buoyancy_bcs, κₑ=(b=κₑ_bcs,)) model = IncompressibleModel( @@ -146,30 +140,26 @@ test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType), # yet. grid = RegularRectilinearGrid(FT, size=(1, 1, 1), extent=(1, π, 42), topology=(Periodic, Bounded, Bounded)) - u_boundary_conditions = UVelocityBoundaryConditions(grid; - bottom = simple_function_bc(Value), - top = simple_function_bc(Value), - north = simple_function_bc(Value), - south = simple_function_bc(Value)) - - v_boundary_conditions = VVelocityBoundaryConditions(grid; - bottom = simple_function_bc(Value), - top = simple_function_bc(Value), - north = simple_function_bc(Open), - south = simple_function_bc(Open)) - - - w_boundary_conditions = VVelocityBoundaryConditions(grid; - bottom = simple_function_bc(Open), - top = simple_function_bc(Open), - north = simple_function_bc(Value), - south = simple_function_bc(Value)) - - T_boundary_conditions = TracerBoundaryConditions(grid; - bottom = simple_function_bc(Value), - top = simple_function_bc(Value), - north = simple_function_bc(Value), - south = simple_function_bc(Value)) + u_boundary_conditions = FieldBoundaryConditions(bottom = simple_function_bc(Value), + top = simple_function_bc(Value), + north = simple_function_bc(Value), + south = simple_function_bc(Value)) + + v_boundary_conditions = FieldBoundaryConditions(bottom = simple_function_bc(Value), + top = simple_function_bc(Value), + north = simple_function_bc(Open), + south = simple_function_bc(Open)) + + + w_boundary_conditions = FieldBoundaryConditions(bottom = simple_function_bc(Open), + top = simple_function_bc(Open), + north = simple_function_bc(Value), + south = simple_function_bc(Value)) + + T_boundary_conditions = FieldBoundaryConditions(bottom = simple_function_bc(Value), + top = simple_function_bc(Value), + north = simple_function_bc(Value), + south = simple_function_bc(Value)) boundary_conditions = (u=u_boundary_conditions, v=v_boundary_conditions, diff --git a/test/test_broadcasting.jl b/test/test_broadcasting.jl index a3f9dd923f..f562fa8775 100644 --- a/test/test_broadcasting.jl +++ b/test/test_broadcasting.jl @@ -27,7 +27,10 @@ three_point_grid = RegularRectilinearGrid(size=(1, 1, 3), extent=(1, 1, 1)) a2 = CenterField(arch, three_point_grid) - b2 = ZFaceField(arch, three_point_grid) + + b2_bcs = FieldBoundaryConditions(grid, (Center, Center, Face), top=OpenBoundaryCondition(0), bottom=OpenBoundaryCondition(0)) + b2 = ZFaceField(arch, three_point_grid, b2_bcs) + b2 .= 1 fill_halo_regions!(b2, arch) # sets b2[1, 1, 1] = b[1, 1, 4] = 0 diff --git a/test/test_checkpointer.jl b/test/test_checkpointer.jl index ba556bdb38..e114881219 100644 --- a/test/test_checkpointer.jl +++ b/test/test_checkpointer.jl @@ -62,19 +62,13 @@ function test_thermal_bubble_checkpointer_output(arch) # Checkpoint should be saved as "checkpoint_iteration5.jld" after the 5th iteration. run!(checkpointed_simulation) # for 5 iterations - # model_kwargs = Dict{Symbol, Any}(:boundary_conditions => SolutionBoundaryConditions(grid)) restored_model = restore_from_checkpoint("checkpoint_iteration5.jld2") - #restored_simulation = Simulation(restored_model, Δt=Δt, stop_iteration=9) - #run!(restored_simulation) - for n in 1:4 update_state!(restored_model) time_step!(restored_model, Δt, euler=false) # time-step for 4 iterations end - # test_model_equality(restored_model, true_model) - ##### ##### Test `set!(model, checkpoint_file)` ##### @@ -124,8 +118,8 @@ function test_checkpoint_output_with_function_bcs(arch) @inline some_flux(x, y, t) = 2x + exp(y) top_u_bc = top_T_bc = FluxBoundaryCondition(some_flux) - u_bcs = UVelocityBoundaryConditions(grid, top=top_u_bc) - T_bcs = TracerBoundaryConditions(grid, top=top_T_bc) + u_bcs = FieldBoundaryConditions(top=top_u_bc) + T_bcs = FieldBoundaryConditions(top=top_T_bc) model = IncompressibleModel(architecture=arch, grid=grid, boundary_conditions=(u=u_bcs, T=T_bcs)) set!(model, u=π/2, v=ℯ, T=Base.MathConstants.γ, S=Base.MathConstants.φ) @@ -170,23 +164,23 @@ function test_checkpoint_output_with_function_bcs(arch) u, v, w = properly_restored_model.velocities T, S = properly_restored_model.tracers - @test u.boundary_conditions.x.left isa PBC - @test u.boundary_conditions.x.right isa PBC - @test u.boundary_conditions.y.left isa PBC - @test u.boundary_conditions.y.right isa PBC - @test u.boundary_conditions.z.left isa ZFBC - @test u.boundary_conditions.z.right isa FBC - @test u.boundary_conditions.z.right.condition isa ContinuousBoundaryFunction - @test u.boundary_conditions.z.right.condition.func(1, 2, 3) == some_flux(1, 2, 3) - - @test T.boundary_conditions.x.left isa PBC - @test T.boundary_conditions.x.right isa PBC - @test T.boundary_conditions.y.left isa PBC - @test T.boundary_conditions.y.right isa PBC - @test T.boundary_conditions.z.left isa ZFBC - @test T.boundary_conditions.z.right isa FBC - @test T.boundary_conditions.z.right.condition isa ContinuousBoundaryFunction - @test T.boundary_conditions.z.right.condition.func(1, 2, 3) == some_flux(1, 2, 3) + @test u.boundary_conditions.west isa PBC + @test u.boundary_conditions.east isa PBC + @test u.boundary_conditions.south isa PBC + @test u.boundary_conditions.north isa PBC + @test u.boundary_conditions.bottom isa ZFBC + @test u.boundary_conditions.top isa FBC + @test u.boundary_conditions.top.condition isa ContinuousBoundaryFunction + @test u.boundary_conditions.top.condition.func(1, 2, 3) == some_flux(1, 2, 3) + + @test T.boundary_conditions.west isa PBC + @test T.boundary_conditions.east isa PBC + @test T.boundary_conditions.south isa PBC + @test T.boundary_conditions.north isa PBC + @test T.boundary_conditions.bottom isa ZFBC + @test T.boundary_conditions.top isa FBC + @test T.boundary_conditions.top.condition isa ContinuousBoundaryFunction + @test T.boundary_conditions.top.condition.func(1, 2, 3) == some_flux(1, 2, 3) # Test that the restored model can be time stepped time_step!(properly_restored_model, 1) diff --git a/test/test_distributed_poisson_solvers.jl b/test/test_distributed_poisson_solvers.jl index 9b0dcea58f..6fc7acc3a1 100644 --- a/test/test_distributed_poisson_solvers.jl +++ b/test/test_distributed_poisson_solvers.jl @@ -1,15 +1,15 @@ function random_divergent_source_term(arch, grid) # Generate right hand side from a random (divergent) velocity field. - Ru = CenterField(arch, grid, UVelocityBoundaryConditions(grid)) - Rv = CenterField(arch, grid, VVelocityBoundaryConditions(grid)) - Rw = CenterField(arch, grid, WVelocityBoundaryConditions(grid)) + Ru = XFaceField(arch, grid) + Rv = YFaceField(arch, grid) + Rw = ZFaceField(arch, grid) U = (u=Ru, v=Rv, w=Rw) Nx, Ny, Nz = size(grid) - set!(Ru, rand(Nx, Ny, Nz)) - set!(Rv, rand(Nx, Ny, Nz)) - set!(Rw, rand(Nx, Ny, Nz)) + set!(Ru, (x, y, z) -> rand()) + set!(Rv, (x, y, z) -> rand()) + set!(Rw, (x, y, z) -> rand()) fill_halo_regions!(Ru, arch) fill_halo_regions!(Rv, arch) @@ -39,7 +39,7 @@ function divergence_free_poisson_solution_triply_periodic(grid_points, ranks) solve_poisson_equation!(solver) - p_bcs = PressureBoundaryConditions(local_grid) + p_bcs = FieldBoundaryConditions(local_grid, (Center, Center, Center)) p_bcs = inject_halo_communication_boundary_conditions(p_bcs, arch.local_rank, arch.connectivity) ϕ = CenterField(child_architecture(arch), local_grid, p_bcs) # "pressure" diff --git a/test/test_dynamics.jl b/test/test_dynamics.jl index 7a5d1dc6d1..c769835880 100644 --- a/test/test_dynamics.jl +++ b/test/test_dynamics.jl @@ -256,7 +256,7 @@ function stratified_fluid_remains_at_rest_with_tilted_gravity_buoyancy_tracer(ar y_bc = GradientBoundaryCondition(N² * g̃[2]) z_bc = GradientBoundaryCondition(N² * g̃[3]) - b_bcs = TracerBoundaryConditions(grid, bottom=z_bc, top=z_bc, south=y_bc, north=y_bc) + b_bcs = FieldBoundaryConditions(bottom=z_bc, top=z_bc, south=y_bc, north=y_bc) model = IncompressibleModel( architecture = arch, @@ -312,7 +312,7 @@ function stratified_fluid_remains_at_rest_with_tilted_gravity_temperature_tracer y_bc = GradientBoundaryCondition(∂T∂z * g̃[2]) z_bc = GradientBoundaryCondition(∂T∂z * g̃[3]) - T_bcs = TracerBoundaryConditions(grid, bottom=z_bc, top=z_bc, south=y_bc, north=y_bc) + T_bcs = FieldBoundaryConditions(bottom=z_bc, top=z_bc, south=y_bc, north=y_bc) model = IncompressibleModel( architecture = arch, diff --git a/test/test_operators.jl b/test/test_operators.jl index 7d9295878b..7b377bbd57 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -132,13 +132,14 @@ end arch = CPU() grid = RegularRectilinearGrid(size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz)) + bcs = FieldBoundaryConditions(grid, (Center, Center, Center)) Hx, Hy, Hz = grid.Hx, grid.Hy, grid.Hz Tx, Ty, Tz = Nx+2Hx, Ny+2Hy, Nz+2Hz A3 = OffsetArray(zeros(Tx, Ty, Tz), 1-Hx:Nx+Hx, 1-Hy:Ny+Hy, 1-Hz:Nz+Hz) @. @views A3[1:Nx, 1:Ny, 1:Nz] = rand() - fill_halo_regions!(A3, TracerBoundaryConditions(grid), arch, grid) + fill_halo_regions!(A3, bcs, arch, grid) # A yz-slice with Nx==1. A2yz = OffsetArray(zeros(1+2Hx, Ty, Tz), 1-Hx:1+Hx, 1-Hy:Ny+Hy, 1-Hz:Nz+Hz) diff --git a/test/test_output_readers.jl b/test/test_output_readers.jl index e7f5b2cb7a..4347f40498 100644 --- a/test/test_output_readers.jl +++ b/test/test_output_readers.jl @@ -10,12 +10,12 @@ using Oceananigans.Fields: location function generate_some_interesting_simulation_data(Nx, Ny, Nz; architecture=CPU()) grid = RegularRectilinearGrid(size=(Nx, Ny, Nz), extent=(64, 64, 32)) - T_bcs = TracerBoundaryConditions(grid, top = FluxBoundaryCondition(5e-5), bottom = GradientBoundaryCondition(0.01)) - u_bcs = UVelocityBoundaryConditions(grid, top = FluxBoundaryCondition(-3e-4)) + T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(5e-5), bottom = GradientBoundaryCondition(0.01)) + u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(-3e-4)) @inline Qˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S evaporation_bc = FluxBoundaryCondition(Qˢ, field_dependencies=:S, parameters=3e-7) - S_bcs = TracerBoundaryConditions(grid, top=evaporation_bc) + S_bcs = FieldBoundaryConditions(top=evaporation_bc) model = IncompressibleModel( architecture = architecture, diff --git a/test/test_poisson_solvers.jl b/test/test_poisson_solvers.jl index 9bc16ad43d..a847851166 100644 --- a/test/test_poisson_solvers.jl +++ b/test/test_poisson_solvers.jl @@ -1,6 +1,7 @@ using Oceananigans.Solvers: solve_for_pressure!, solve_poisson_equation!, set_source_term! using Oceananigans.Solvers: poisson_eigenvalues using Oceananigans.Models.HydrostaticFreeSurfaceModels: _compute_w_from_continuity! +using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions function poisson_solver_instantiates(arch, grid, planner_flag) solver = FFTBasedPoissonSolver(arch, grid, planner_flag) @@ -8,9 +9,14 @@ function poisson_solver_instantiates(arch, grid, planner_flag) end function random_divergent_source_term(arch, grid) - Ru = CenterField(arch, grid, UVelocityBoundaryConditions(grid)) - Rv = CenterField(arch, grid, VVelocityBoundaryConditions(grid)) - Rw = CenterField(arch, grid, WVelocityBoundaryConditions(grid)) + default_bcs = FieldBoundaryConditions() + u_bcs = regularize_field_boundary_conditions(default_bcs, grid, :u) + v_bcs = regularize_field_boundary_conditions(default_bcs, grid, :v) + w_bcs = regularize_field_boundary_conditions(default_bcs, grid, :w) + + Ru = CenterField(arch, grid, u_bcs) + Rv = CenterField(arch, grid, v_bcs) + Rw = CenterField(arch, grid, w_bcs) U = (u=Ru, v=Rv, w=Rw) Nx, Ny, Nz = size(grid) @@ -33,10 +39,15 @@ function random_divergent_source_term(arch, grid) end function random_divergence_free_source_term(arch, grid) + default_bcs = FieldBoundaryConditions() + u_bcs = regularize_field_boundary_conditions(default_bcs, grid, :u) + v_bcs = regularize_field_boundary_conditions(default_bcs, grid, :v) + w_bcs = regularize_field_boundary_conditions(default_bcs, grid, :w) + # Random right hand side - Ru = CenterField(arch, grid, UVelocityBoundaryConditions(grid)) - Rv = CenterField(arch, grid, VVelocityBoundaryConditions(grid)) - Rw = CenterField(arch, grid, WVelocityBoundaryConditions(grid)) + Ru = CenterField(arch, grid, u_bcs) + Rv = CenterField(arch, grid, v_bcs) + Rw = CenterField(arch, grid, w_bcs) U = (u=Ru, v=Rv, w=Rw) Nx, Ny, Nz = size(grid) @@ -75,7 +86,7 @@ function divergence_free_poisson_solution(arch, grid, planner_flag=FFTW.MEASURE) solver = FFTBasedPoissonSolver(arch, grid, planner_flag) R, U = random_divergent_source_term(arch, grid) - p_bcs = PressureBoundaryConditions(grid) + p_bcs = FieldBoundaryConditions(grid, (Center, Center, Center)) ϕ = CenterField(arch, grid, p_bcs) # "kinematic pressure" ∇²ϕ = CenterField(arch, grid, p_bcs) @@ -149,7 +160,7 @@ function vertically_stretched_poisson_solver_correct_answer(FT, arch, topo, Nx, vs_grid = VerticallyStretchedRectilinearGrid(FT; architecture=arch, topology=topo, size=sz, z_faces=zF, xy_intervals...) solver = FourierTridiagonalPoissonSolver(arch, vs_grid) - p_bcs = PressureBoundaryConditions(vs_grid) + p_bcs = FieldBoundaryConditions(vs_grid, (Center, Center, Center)) ϕ = CenterField(arch, vs_grid, p_bcs) # "kinematic pressure" ∇²ϕ = CenterField(arch, vs_grid, p_bcs) diff --git a/test/test_time_stepping.jl b/test/test_time_stepping.jl index b838232ae5..0f0fcced9f 100644 --- a/test/test_time_stepping.jl +++ b/test/test_time_stepping.jl @@ -101,7 +101,7 @@ function incompressible_in_time(arch, grid, Nt, timestepper) grid = model.grid u, v, w = model.velocities - div_U = CenterField(arch, grid, TracerBoundaryConditions(grid)) + div_U = CenterField(arch, grid) # Just add a temperature perturbation so we get some velocity field. CUDA.@allowscalar interior(model.tracers.T)[8:24, 8:24, 8:24] .+= 0.01 diff --git a/validation/barotropic_gyre/barotropic_gyre.jl b/validation/barotropic_gyre/barotropic_gyre.jl index d73d7c3c64..93c848c6d6 100644 --- a/validation/barotropic_gyre/barotropic_gyre.jl +++ b/validation/barotropic_gyre/barotropic_gyre.jl @@ -37,15 +37,6 @@ underlying_grid = RegularLatitudeLongitudeGrid(size = (Nx, Ny, 1), grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBoundary(raster_depth, mask_type=RasterDepthMask())) -solid(x, y, z, i, j, k) = ( - if i > 30 && i < 35; - if j > 42 && j < 48; - return true; - end; - end; - return false; - ) - free_surface = ImplicitFreeSurface(gravitational_acceleration=0.1) # free_surface = ExplicitFreeSurface(gravitational_acceleration=0.1) @@ -73,12 +64,10 @@ v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag, discrete_form = true, parameters = μ) -u_bcs = UVelocityBoundaryConditions(grid, - top = surface_wind_stress_bc, - bottom = u_bottom_drag_bc) +u_bcs = FieldBoundaryConditions(top = surface_wind_stress_bc, + bottom = u_bottom_drag_bc) -v_bcs = VVelocityBoundaryConditions(grid, - bottom = v_bottom_drag_bc) +v_bcs = FieldBoundaryConditions(bottom = v_bottom_drag_bc) @show const νh₀ = 5e3 * (60 / grid.Nx)^2 diff --git a/validation/convergence_tests/Manifest.toml b/validation/convergence_tests/Manifest.toml index 9fd6cef879..d3b5b9dbe0 100644 --- a/validation/convergence_tests/Manifest.toml +++ b/validation/convergence_tests/Manifest.toml @@ -8,9 +8,9 @@ version = "1.0.1" [[Adapt]] deps = ["LinearAlgebra"] -git-tree-sha1 = "ffcfa2d345aaee0ef3d8346a073d5dd03c983ebe" +git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.2.0" +version = "3.3.1" [[ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" @@ -33,40 +33,34 @@ uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" version = "0.4.1" [[CUDA]] -deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CompilerSupportLibraries_jll", "DataStructures", "ExprTools", "GPUArrays", "GPUCompiler", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "MacroTools", "Memoize", "NNlib", "Printf", "Random", "Reexport", "Requires", "SparseArrays", "Statistics", "TimerOutputs"] -git-tree-sha1 = "870e029382294443a6578190e992bf4cbfd34e22" +deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CompilerSupportLibraries_jll", "DataStructures", "ExprTools", "GPUArrays", "GPUCompiler", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "SpecialFunctions", "TimerOutputs"] +git-tree-sha1 = "8ef71bf6d6602cf227196b43650924bf9ef7babc" uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" -version = "2.6.2" +version = "3.3.3" [[ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "0893f8d90331a0f5223c7ef2a8868464394a886c" +git-tree-sha1 = "dcc25ff085cf548bc8befad5ce048391a7c07d40" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "0.9.33" - -[[CodecZlib]] -deps = ["TranscodingStreams", "Zlib_jll"] -git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da" -uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.7.0" +version = "0.10.11" [[ColorTypes]] deps = ["FixedPointNumbers", "Random"] -git-tree-sha1 = "32a2b8af383f11cbb65803883837a149d10dfe8a" +git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597" uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" -version = "0.10.12" +version = "0.11.0" [[Colors]] -deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Reexport"] -git-tree-sha1 = "ac5f2213e56ed8a34a3dd2f681f4df1166b34929" +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" -version = "0.12.6" +version = "0.12.8" [[Compat]] deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "919c7f3151e79ff196add81d7f4e45d91bbf420b" +git-tree-sha1 = "dc7dedc2c2aa9faf59a55c622760a25cbefbe941" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.25.0" +version = "3.31.0" [[CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] @@ -74,9 +68,9 @@ uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" [[Conda]] deps = ["JSON", "VersionParsing"] -git-tree-sha1 = "6231e40619c15148bcb80aa19d731e629877d762" +git-tree-sha1 = "299304989a5e6473d985212c28928899c74e9421" uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d" -version = "1.5.1" +version = "1.5.2" [[DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] @@ -96,14 +90,26 @@ uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" +[[DocStringExtensions]] +deps = ["LibGit2"] +git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.8.5" + [[Downloads]] deps = ["ArgTools", "LibCURL", "NetworkOptions"] uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" [[ExprTools]] -git-tree-sha1 = "10407a39b87f29d47ebaca8edbc75d7c302ff93e" +git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92" uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" -version = "0.1.3" +version = "0.1.6" + +[[FileIO]] +deps = ["Pkg", "Requires", "UUIDs"] +git-tree-sha1 = "256d8e6188f3f1ebfa1a5d17e072a0efafa8c5bf" +uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +version = "1.10.1" [[FixedPointNumbers]] deps = ["Statistics"] @@ -112,16 +118,16 @@ uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" version = "0.8.4" [[GPUArrays]] -deps = ["AbstractFFTs", "Adapt", "LinearAlgebra", "Printf", "Random", "Serialization"] -git-tree-sha1 = "f99a25fe0313121f2f9627002734c7d63b4dd3bd" +deps = ["AbstractFFTs", "Adapt", "LinearAlgebra", "Printf", "Random", "Serialization", "Statistics"] +git-tree-sha1 = "ececbf05f8904c92814bdbd0aafd5540b0bf2e9a" uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" -version = "6.2.0" +version = "7.0.1" [[GPUCompiler]] -deps = ["DataStructures", "ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "Serialization", "TimerOutputs", "UUIDs"] -git-tree-sha1 = "ef2839b063e158672583b9c09d2cf4876a8d3d55" +deps = ["DataStructures", "ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "TimerOutputs", "UUIDs"] +git-tree-sha1 = "e8a09182a4440489e2e3dedff5ad3f6bbe555396" uuid = "61eb1bfa-7361-4325-ad38-22787b887f55" -version = "0.10.0" +version = "0.12.5" [[Glob]] git-tree-sha1 = "4df9f7e06108728ebf00a0a11edee4b29a482bb2" @@ -133,10 +139,16 @@ deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" [[JLD2]] -deps = ["CodecZlib", "DataStructures", "MacroTools", "Mmap", "Pkg", "Printf", "Requires", "UUIDs"] -git-tree-sha1 = "9f2f2f24e60305feb6ae293a617ddf06f429efc3" +deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"] +git-tree-sha1 = "4813826871754cf52607e76ad37acb36ccf52719" uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -version = "0.4.3" +version = "0.4.11" + +[[JLLWrappers]] +deps = ["Preferences"] +git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.3.0" [[JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] @@ -145,10 +157,16 @@ uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" version = "0.21.1" [[LLVM]] -deps = ["CEnum", "Libdl", "Printf", "Unicode"] -git-tree-sha1 = "b616937c31337576360cb9fb872ec7633af7b194" +deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"] +git-tree-sha1 = "1b7ba36ea7aa6fa2278118951bad114fbb8359f2" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "3.6.0" +version = "4.1.0" + +[[LLVMExtra_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "b36c0677a0549c7d1dc8719899a4133abbfacf7d" +uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" +version = "0.0.6+0" [[LaTeXStrings]] git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104" @@ -182,6 +200,12 @@ uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" deps = ["Libdl"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +[[LogExpFunctions]] +deps = ["DocStringExtensions", "LinearAlgebra"] +git-tree-sha1 = "7bd5f6565d80b6bf753738d2bc40a5dfea072070" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.2.5" + [[Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" @@ -199,37 +223,31 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -[[Memoize]] -deps = ["MacroTools"] -git-tree-sha1 = "2b1dfcba103de714d31c033b5dacc2e4a12c7caa" -uuid = "c03570c3-d221-55d1-a50c-7939bbd78826" -version = "0.4.4" - [[Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" [[MozillaCACerts_jll]] uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -[[NNlib]] -deps = ["ChainRulesCore", "Compat", "LinearAlgebra", "Pkg", "Requires", "Statistics"] -git-tree-sha1 = "ab1d43fead2ecb9aa5ae460d3d547c2cf8d89461" -uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" -version = "0.7.17" - [[NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" [[OffsetArrays]] deps = ["Adapt"] -git-tree-sha1 = "b3dfef5f2be7d7eb0e782ba9146a5271ee426e90" +git-tree-sha1 = "2bf78c5fd7fa56d2bbf1efbadd45c1b8789e6f57" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.6.2" +version = "1.10.2" + +[[OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.5+0" [[OrderedCollections]] -git-tree-sha1 = "4fa2ba51070ec13fcc7517db714445b4ab986bdf" +git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.4.0" +version = "1.4.1" [[Parsers]] deps = ["Dates"] @@ -238,18 +256,24 @@ uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" version = "1.1.0" [[Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs"] +deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +[[Preferences]] +deps = ["TOML"] +git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.2.2" + [[Printf]] deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" [[PyCall]] deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"] -git-tree-sha1 = "dd1a970b543bd02efce2984582e996af28cab27f" +git-tree-sha1 = "169bb8ea6b1b143c5cf57df6d34d022a7b60c6db" uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" -version = "1.92.2" +version = "1.92.3" [[PyPlot]] deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"] @@ -265,10 +289,22 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" deps = ["Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +[[Random123]] +deps = ["Libdl", "Random", "RandomNumbers"] +git-tree-sha1 = "0e8b146557ad1c6deb1367655e052276690e71a3" +uuid = "74087812-796a-5b5d-8853-05524746bad3" +version = "1.4.2" + +[[RandomNumbers]] +deps = ["Random", "Requires"] +git-tree-sha1 = "441e6fc35597524ada7f85e13df1f4e10137d16f" +uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143" +version = "1.4.0" + [[Reexport]] -git-tree-sha1 = "57d8440b0c7d98fc4f889e478e80f268d534c9d5" +git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220" uuid = "189a3867-3050-52da-a836-e630ba90ab69" -version = "1.0.0" +version = "1.1.0" [[Requires]] deps = ["UUIDs"] @@ -279,12 +315,6 @@ version = "1.1.3" [[SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" -[[Scratch]] -deps = ["Dates"] -git-tree-sha1 = "ad4b278adb62d185bbcb6864dc24959ab0627bf6" -uuid = "6c6a2e73-6563-6170-7368-637461726353" -version = "1.0.3" - [[Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -299,6 +329,12 @@ uuid = "6462fe0b-24de-5631-8697-dd941f90decc" deps = ["LinearAlgebra", "Random"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +[[SpecialFunctions]] +deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"] +git-tree-sha1 = "a50550fa3164a8c46747e62063b4d774ac1bcf49" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "1.5.1" + [[Statistics]] deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" @@ -316,10 +352,10 @@ deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [[TimerOutputs]] -deps = ["Printf"] -git-tree-sha1 = "32cdbe6cd2d214c25a0b88f985c9e0092877c236" +deps = ["ExprTools", "Printf"] +git-tree-sha1 = "209a8326c4f955e2442c07b56029e88bb48299c7" uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" -version = "0.5.8" +version = "0.5.12" [[TranscodingStreams]] deps = ["Random", "Test"] @@ -346,3 +382,7 @@ uuid = "83775a58-1f1d-513f-b197-d71354ab007a" [[nghttp2_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" + +[[p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" diff --git a/validation/convergence_tests/run_forced_fixed_slip.jl b/validation/convergence_tests/run_forced_fixed_slip.jl index a1f8137a02..c9fbf599ee 100644 --- a/validation/convergence_tests/run_forced_fixed_slip.jl +++ b/validation/convergence_tests/run_forced_fixed_slip.jl @@ -1,3 +1,5 @@ +pushfirst!(LOAD_PATH, joinpath("..", "..")) + using CUDA using Oceananigans diff --git a/validation/convergence_tests/run_forced_free_slip.jl b/validation/convergence_tests/run_forced_free_slip.jl index c0faf52657..f2ef1f8cb7 100644 --- a/validation/convergence_tests/run_forced_free_slip.jl +++ b/validation/convergence_tests/run_forced_free_slip.jl @@ -1,3 +1,5 @@ +pushfirst!(LOAD_PATH, joinpath("..", "..")) + using CUDA using Oceananigans diff --git a/validation/convergence_tests/run_taylor_green.jl b/validation/convergence_tests/run_taylor_green.jl index 2baa363fcf..74dd4c25fb 100644 --- a/validation/convergence_tests/run_taylor_green.jl +++ b/validation/convergence_tests/run_taylor_green.jl @@ -1,3 +1,5 @@ +pushfirst!(LOAD_PATH, joinpath("..", "..")) + using CUDA using Oceananigans diff --git a/validation/convergence_tests/src/ForcedFlowFixedSlip.jl b/validation/convergence_tests/src/ForcedFlowFixedSlip.jl index 9d2ff4334e..592985961c 100644 --- a/validation/convergence_tests/src/ForcedFlowFixedSlip.jl +++ b/validation/convergence_tests/src/ForcedFlowFixedSlip.jl @@ -3,9 +3,6 @@ module ForcedFlowFixedSlip using Printf using Oceananigans -using Oceananigans.BoundaryConditions -using Oceananigans.OutputWriters -using Oceananigans.Fields # Functions that define the forced flow problem @@ -48,17 +45,17 @@ function setup_xy_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir= topology=(Periodic, Bounded, Bounded)) # "Fixed slip" boundary conditions (eg, no-slip on south wall, finite slip on north wall)." - u_bcs = UVelocityBoundaryConditions(grid, north = ValueBoundaryCondition((x, y, t) -> f(x, t)), - south = ValueBoundaryCondition(0)) - - model = IncompressibleModel( architecture = CPU(), - grid = grid, - coriolis = nothing, - buoyancy = nothing, - tracers = nothing, - closure = IsotropicDiffusivity(ν=1), + u_bcs = FieldBoundaryConditions(north = ValueBoundaryCondition((x, y, t) -> f(x, t)), + south = ValueBoundaryCondition(0)) + + model = IncompressibleModel(architecture = CPU(), + grid = grid, + coriolis = nothing, + buoyancy = nothing, + tracers = nothing, + closure = IsotropicDiffusivity(ν=1), boundary_conditions = (u=u_bcs,), - forcing = (u=u_forcing, v=v_forcing)) + forcing = (u=u_forcing, v=v_forcing)) set!(model, u = (x, y, z) -> u(x, y, 0), v = (x, y, z) -> v(x, y, 0)) @@ -94,17 +91,17 @@ function setup_xz_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir= topology=(Periodic, Bounded, Bounded)) # "Fixed slip" boundary conditions (eg, no-slip on bottom and finite slip on top)." - u_bcs = UVelocityBoundaryConditions(grid, top = ValueBoundaryCondition((x, z, t) -> f(x, t)), - bottom = ValueBoundaryCondition(0)) - - model = IncompressibleModel( architecture = CPU(), - grid = grid, - coriolis = nothing, - buoyancy = nothing, - tracers = nothing, - closure = IsotropicDiffusivity(ν=1), + u_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition((x, z, t) -> f(x, t)), + bottom = ValueBoundaryCondition(0)) + + model = IncompressibleModel(architecture = CPU(), + grid = grid, + coriolis = nothing, + buoyancy = nothing, + tracers = nothing, + closure = IsotropicDiffusivity(ν=1), boundary_conditions = (u=u_bcs,), - forcing = (u=u_forcing, w=w_forcing)) + forcing = (u=u_forcing, w=w_forcing)) set!(model, u = (x, y, z) -> u(x, z, 0), w = (x, y, z) -> v(x, z, 0)) diff --git a/validation/lid_driven_cavity/lid_driven_cavity.jl b/validation/lid_driven_cavity/lid_driven_cavity.jl index f755765d23..b3a636f132 100644 --- a/validation/lid_driven_cavity/lid_driven_cavity.jl +++ b/validation/lid_driven_cavity/lid_driven_cavity.jl @@ -9,15 +9,11 @@ function simulate_lid_driven_cavity(; Re, N, end_time) domain = (y=(0, 1), z=(0, 1)) grid = RegularRectilinearGrid(topology=topology, size=(N, N); domain...) - v_bcs = VVelocityBoundaryConditions(grid, - top = ValueBoundaryCondition(1), - bottom = ValueBoundaryCondition(0) - ) + v_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(1), + bottom = ValueBoundaryCondition(0)) - w_bcs = WVelocityBoundaryConditions(grid, - north = ValueBoundaryCondition(0), - south = ValueBoundaryCondition(0) - ) + w_bcs = FieldBoundaryConditions(north = ValueBoundaryCondition(0), + south = ValueBoundaryCondition(0)) model = IncompressibleModel( grid = grid, diff --git a/validation/stratified_couette_flow/stratified_couette_flow.jl b/validation/stratified_couette_flow/stratified_couette_flow.jl index 9faa9b7807..295bf0fbd0 100644 --- a/validation/stratified_couette_flow/stratified_couette_flow.jl +++ b/validation/stratified_couette_flow/stratified_couette_flow.jl @@ -100,14 +100,14 @@ function simulate_stratified_couette_flow(; Nxy, Nz, arch=GPU(), h=1, U_wall=1, grid = RegularRectilinearGrid(size = (Nxy, Nxy, Nz), extent = (4π*h, 2π*h, 2h)) - Tbcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Value, Θ_wall), - bottom = BoundaryCondition(Value, -Θ_wall)) + Tbcs = FieldBoundaryConditions(top = BoundaryCondition(Value, Θ_wall), + bottom = BoundaryCondition(Value, -Θ_wall)) - ubcs = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Value, U_wall), - bottom = BoundaryCondition(Value, -U_wall)) + ubcs = FieldBoundaryConditions(top = BoundaryCondition(Value, U_wall), + bottom = BoundaryCondition(Value, -U_wall)) - vbcs = VVelocityBoundaryConditions(grid, top = BoundaryCondition(Value, 0), - bottom = BoundaryCondition(Value, 0)) + vbcs = FieldBoundaryConditions(top = BoundaryCondition(Value, 0), + bottom = BoundaryCondition(Value, 0)) ##### ##### Non-dimensional model setup diff --git a/validation/vertical_mixing_closures/convective_adjustment_free_convection.jl b/validation/vertical_mixing_closures/convective_adjustment_free_convection.jl index 6718536189..ff8b061003 100644 --- a/validation/vertical_mixing_closures/convective_adjustment_free_convection.jl +++ b/validation/vertical_mixing_closures/convective_adjustment_free_convection.jl @@ -9,7 +9,7 @@ grid = RegularRectilinearGrid(size=32, z=(-64, 0), topology=(Flat, Flat, Bounded closure = ConvectiveAdjustmentVerticalDiffusivity(convective_κz=1, background_κz=1e-5) Qᵇ = 1e-8 -b_bcs = TracerBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵇ)) +b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ)) model = HydrostaticFreeSurfaceModel(grid = grid, tracers = :b, diff --git a/validation/vertical_mixing_closures/tke_based_free_convection.jl b/validation/vertical_mixing_closures/tke_based_free_convection.jl index fee3f23599..e0b576f833 100644 --- a/validation/vertical_mixing_closures/tke_based_free_convection.jl +++ b/validation/vertical_mixing_closures/tke_based_free_convection.jl @@ -18,15 +18,14 @@ w★ = Qᵇ * grid.Δz Qᵉ = - closure.dissipation_parameter * (closure.surface_model.CᵂwΔ * w★ + closure.surface_model.Cᵂu★ * u★) -u_bcs = UVelocityBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵘ)) -v_bcs = VVelocityBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵛ)) -b_bcs = TracerBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵇ)) -tke_bcs = TracerBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵉ)) +u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) +v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ)) +b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ)) model = HydrostaticFreeSurfaceModel(grid = grid, tracers = (:b, :e), buoyancy = BuoyancyTracer(), - boundary_conditions = (b=b_bcs, e=tke_bcs), + boundary_conditions = (; b=b_bcs), closure = closure) N² = 1e-5 diff --git a/validation/vertical_mixing_closures/tke_based_wind_mixing.jl b/validation/vertical_mixing_closures/tke_based_wind_mixing.jl index 55d9152bfb..94e478ab8d 100644 --- a/validation/vertical_mixing_closures/tke_based_wind_mixing.jl +++ b/validation/vertical_mixing_closures/tke_based_wind_mixing.jl @@ -17,16 +17,14 @@ w★³ = Qᵇ * grid.Δz Qᵉ = - closure.dissipation_parameter * (closure.surface_model.CᵂwΔ * w★³ + closure.surface_model.Cᵂu★ * u★^3) -u_bcs = UVelocityBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵘ)) -v_bcs = VVelocityBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵛ)) -b_bcs = TracerBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵇ)) -tke_bcs = TracerBoundaryConditions(grid; top = FluxBoundaryCondition(Qᵉ)) +u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) +v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ)) +b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ)) model = HydrostaticFreeSurfaceModel(grid = grid, tracers = (:b, :e), buoyancy = BuoyancyTracer(), coriolis = FPlane(f=1e-4), - #boundary_conditions = (b=b_bcs, e=tke_bcs, u=u_bcs, v=v_bcs), boundary_conditions = (b=b_bcs, u=u_bcs, v=v_bcs), closure = closure)