From 1d7748dd2fc0005a711325d02fbebdd12cf4663f Mon Sep 17 00:00:00 2001 From: tomaschor Date: Wed, 30 Jul 2025 20:16:14 -0700 Subject: [PATCH 01/18] convert vectors to tuples --- src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl | 4 ++-- src/Oceananigans.jl | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl b/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl index be2629657f1..25f65ed3afd 100644 --- a/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl +++ b/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl @@ -149,8 +149,8 @@ function initialize_boundary_mass_fluxes(velocities::NamedTuple) total_area_matching_scheme_boundaries += boundary_fluxes.top_area end - boundary_fluxes = merge(boundary_fluxes, (; left_matching_scheme_boundaries, - right_matching_scheme_boundaries, + boundary_fluxes = merge(boundary_fluxes, (; left_matching_scheme_boundaries = Tuple(left_matching_scheme_boundaries), + right_matching_scheme_boundaries = Tuple(right_matching_scheme_boundaries), total_area_matching_scheme_boundaries)) if length(left_matching_scheme_boundaries) == 0 && length(right_matching_scheme_boundaries) == 0 diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index d600fad17a8..128b38a647d 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -35,6 +35,7 @@ export # Boundary conditions BoundaryCondition, FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, OpenBoundaryCondition, + PerturbationAdvectionOpenBoundaryCondition, FieldBoundaryConditions, # Fields and field manipulation From 7e3ce33b22d2f9b246277757c4d910b5bbd2c661 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Fri, 8 Aug 2025 15:57:50 -0700 Subject: [PATCH 02/18] remove warning --- .../perturbation_advection_open_boundary_matching_scheme.jl | 2 -- src/Oceananigans.jl | 1 - 2 files changed, 3 deletions(-) diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl index 81312cd1efc..e24502a1c49 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl @@ -68,8 +68,6 @@ function PerturbationAdvectionOpenBoundaryCondition(val, FT = defaults.FloatType outflow_timescale = convert(FT, outflow_timescale) classification = Open(PerturbationAdvection(inflow_timescale, outflow_timescale)) - @warn "`PerturbationAdvection` open boundaries matching scheme is experimental and un-tested/validated" - return BoundaryCondition(classification, val; kwargs...) end diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index c3bf24e01df..4efb7f7c37b 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -35,7 +35,6 @@ export # Boundary conditions BoundaryCondition, FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, OpenBoundaryCondition, - PerturbationAdvectionOpenBoundaryCondition, FieldBoundaryConditions, # Fields and field manipulation From 9c031d5dc54fca25ac989a0f4d34088d8584ea4d Mon Sep 17 00:00:00 2001 From: tomaschor Date: Fri, 8 Aug 2025 16:20:27 -0700 Subject: [PATCH 03/18] matching_scheme -> scheme --- src/BoundaryConditions/BoundaryConditions.jl | 4 +- .../boundary_condition_classifications.jl | 4 +- ...lat_extrapolation_open_boundary_scheme.jl} | 6 +- ...rbation_advection_open_boundary_scheme.jl} | 56 +++---------------- .../boundary_mass_fluxes.jl | 44 +++++++-------- .../open_boundaries/oscillating_flow.jl | 8 +-- 6 files changed, 41 insertions(+), 81 deletions(-) rename src/BoundaryConditions/{flat_extrapolation_open_boundary_matching_scheme.jl => flat_extrapolation_open_boundary_scheme.jl} (97%) rename src/BoundaryConditions/{perturbation_advection_open_boundary_matching_scheme.jl => perturbation_advection_open_boundary_scheme.jl} (65%) diff --git a/src/BoundaryConditions/BoundaryConditions.jl b/src/BoundaryConditions/BoundaryConditions.jl index 12a16a5c7ed..7bd24263334 100644 --- a/src/BoundaryConditions/BoundaryConditions.jl +++ b/src/BoundaryConditions/BoundaryConditions.jl @@ -40,6 +40,6 @@ include("compute_flux_bcs.jl") include("update_boundary_conditions.jl") include("polar_boundary_condition.jl") -include("flat_extrapolation_open_boundary_matching_scheme.jl") -include("perturbation_advection_open_boundary_matching_scheme.jl") +include("flat_extrapolation_open_boundary_scheme.jl") +include("perturbation_advection_open_boundary_scheme.jl") end # module diff --git a/src/BoundaryConditions/boundary_condition_classifications.jl b/src/BoundaryConditions/boundary_condition_classifications.jl index 687d765acb9..3aa581f165f 100644 --- a/src/BoundaryConditions/boundary_condition_classifications.jl +++ b/src/BoundaryConditions/boundary_condition_classifications.jl @@ -60,14 +60,14 @@ Open boundary conditions are used to specify the component of a velocity field n and can also be used to describe nested or linked simulation domains. """ struct Open{MS} <: AbstractBoundaryConditionClassification - matching_scheme::MS + scheme::MS end Open() = Open(nothing) (open::Open)() = open -Adapt.adapt_structure(to, open::Open) = Open(adapt(to, open.matching_scheme)) +Adapt.adapt_structure(to, open::Open) = Open(adapt(to, open.scheme)) """ struct MultiRegionCommunication <: AbstractBoundaryConditionClassification diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl similarity index 97% rename from src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl rename to src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl index 6ab41143a59..8365299050b 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl @@ -8,7 +8,7 @@ Zero gradient perpendicular velocity boundary condition. We find the boundary value by Taylor expanding the gradient at the boundary point (`xᵢ`) to second order: ```math -f′(xᵢ) ≈ f′(xᵢ₋₁) + f′′(xᵢ₋₁)(xᵢ₋₁ - xᵢ) + O(Δx²) = f′(xᵢ₋₁) + f′′(xᵢ₋₁)Δx + O(Δx²), +f′(xᵢ) ≈ f′(xᵢ₋₁) + f′′(xᵢ₋₁)(xᵢ₋₁ - xᵢ) + O(Δx²), ``` where ``Δx=xᵢ₋₁ - xᵢ`` (for simplicity, we will also assume the spacing is constant at all ``i`` for now). @@ -47,7 +47,7 @@ end @inline function relax(l, m, grid, ϕ, bc, clock, model_fields) Δt = clock.last_stage_Δt - τ = bc.classification.matching_scheme.relaxation_timescale + τ = bc.classification.scheme.relaxation_timescale Δt̄ = min(1, Δt / τ) ϕₑₓₜ = getbc(bc, l, m, grid, clock, model_fields) @@ -148,3 +148,5 @@ end return nothing end + + diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl similarity index 65% rename from src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl rename to src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl index e24502a1c49..2f658eeef48 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl @@ -7,60 +7,16 @@ using Oceananigans: defaults For cases where we assume that the internal flow is a small perturbation from an external prescribed or coarser flow, we can split the velocity into background and perturbation components. - -We begin with the equation governing the fluid in the interior: - ∂ₜu + u⋅∇u = −∇P + F, -and note that on the boundary the pressure gradient is zero. -We can then assume that the flow composes of mean (U⃗) and pertubation (u⃗′) components, -and considering the x-component of velocity, we can rewrite the equation as - ∂ₜu₁ = -u₁∂₁u - u₂∂₂u₁ - u₃∂₃u₁ + F₁ ≈ - U₁∂₁u₁′ - U₂∂₂u₁′ - U₃∂₃u₁′ + F. - -Simplify by assuming that U⃗ = Ux̂, an then take a numerical step to find u₁. - -When the boundaries are filled the interior is at time tₙ₊₁ so we can take -a backwards euler step (in the case that the mean flow is boundary normal) on a right boundary: - (Uⁿ⁺¹ - Uⁿ) / Δt + (u′ⁿ⁺¹ - u′ⁿ) / Δt = - Uⁿ⁺¹ (u′ⁿ⁺¹ᵢ - u′ⁿ⁺¹ᵢ₋₁) / Δx + Fᵤ. - -This can not be solved for general forcing, but if we assume the dominant forcing is -relaxation to the mean velocity (i.e. u′→0) then Fᵤ = -u′ / τ then we can find u′ⁿ⁺¹: - u′ⁿ⁺¹ = (uⁿ + Ũu′ⁿ⁺¹ᵢ₋₁ - Uⁿ⁺¹) / (1 + Ũ + Δt/τ), - -where Ũ = U Δt / Δx, then uⁿ⁺¹ is: - uⁿ⁺¹ = (uᵢⁿ + Ũuᵢ₋₁ⁿ⁺¹ + Uⁿ⁺¹τ̃) / (1 + τ̃ + Ũ) - -where τ̃ = Δt/τ. - -The same operation can be repeated for left boundaries. - -The relaxation timescale ``τ`` can be set to different values depending on whether -``U`` points in or out of the domain (`inflow_timescale`/`outflow_timescale`). Since the -scheme is only valid when the flow is directed out of the domain the boundary condition -falls back to relaxation to the prescribed value. By default this happens instantly but -if the direction varies this may not be preferable. It is beneficial to relax the outflow -(i.e. non-zero `outflow_timescale`) to reduce the shock when the flow changes direction -to point into the domain. - -The ideal value of the timescales probably depend on the grid spacing and details of the -boundary flow. """ struct PerturbationAdvection{FT} inflow_timescale :: FT - outflow_timescale :: FT + outflow_timescale :: FT end Adapt.adapt_structure(to, pe::PerturbationAdvection) = PerturbationAdvection(adapt(to, pe.inflow_timescale), adapt(to, pe.outflow_timescale)) -""" - PerturbationAdvectionOpenBoundaryCondition(val, FT = defaults.FloatType; - outflow_timescale = Inf, - inflow_timescale = 0, kwargs...) - -Creates a `PerturbationAdvectionOpenBoundaryCondition` with a given exterior value `val`, to which -the flow is forced with an `outflow_timescale` for outflow and `inflow_timescale` for inflow. For -details about this method, refer to the docstring for `PerturbationAdvection`. -""" function PerturbationAdvectionOpenBoundaryCondition(val, FT = defaults.FloatType; outflow_timescale = Inf, inflow_timescale = 0, kwargs...) @@ -85,9 +41,9 @@ const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, iᴬ, jᴬ, kᴬ) U = max(0, min(1, Δt / ΔX * ūⁿ⁺¹)) - pa = bc.classification.matching_scheme + pa = bc.classification.scheme τ = ifelse(ūⁿ⁺¹ >= 0, pa.outflow_timescale, pa.inflow_timescale) - τ̃ = Δt / τ # last stage Δt normalized by the inflow/output timescale + τ̃ = Δt / τ relaxed_uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + U) uᵢⁿ⁺¹ = ifelse(τ == 0, ūⁿ⁺¹, relaxed_uᵢⁿ⁺¹) @@ -109,9 +65,9 @@ end uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, iᴬ, jᴬ, kᴬ) U = min(0, max(-1, Δt / ΔX * ūⁿ⁺¹)) - pa = bc.classification.matching_scheme + pa = bc.classification.scheme τ = ifelse(ūⁿ⁺¹ <= 0, pa.outflow_timescale, pa.inflow_timescale) - τ̃ = Δt / τ # last stage Δt normalized by the inflow/output timescale + τ̃ = Δt / τ relaxed_u₁ⁿ⁺¹ = (uᵢⁿ - U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ - U) u₁ⁿ⁺¹ = ifelse(τ == 0, ūⁿ⁺¹, relaxed_u₁ⁿ⁺¹) @@ -183,3 +139,5 @@ end return nothing end + + diff --git a/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl b/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl index 25f65ed3afd..a0121a9650c 100644 --- a/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl +++ b/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl @@ -4,7 +4,7 @@ using Oceananigans.Fields: Field, interior, XFaceField, YFaceField, ZFaceField using GPUArraysCore: @allowscalar const OBC = BoundaryCondition{<:Open} # OpenBoundaryCondition -const IOBC = BoundaryCondition{<:Open{<:Nothing}} # "Imposed-velocity" OpenBoundaryCondition (with no matching scheme) +const IOBC = BoundaryCondition{<:Open{<:Nothing}} # "Imposed-velocity" OpenBoundaryCondition (with no scheme) const FIOBC = BoundaryCondition{<:Open{<:Nothing}, <:Number} # "Fixed-imposed-velocity" OpenBoundaryCondition const ZIOBC = BoundaryCondition{<:Open{<:Nothing}, <:Nothing} # "Zero-imposed-velocity" OpenBoundaryCondition (no-inflow) @@ -79,7 +79,7 @@ initialize_boundary_mass_flux(velocity, ::Nothing, side) = NamedTuple() initialize_boundary_mass_flux(velocity, bc, side) = NamedTuple() needs_mass_flux_correction(::IOBC) = false -needs_mass_flux_correction(bc::OBC) = bc.classification.matching_scheme !== nothing +needs_mass_flux_correction(bc::OBC) = bc.classification.scheme !== nothing needs_mass_flux_correction(::Nothing) = false needs_mass_flux_correction(bc) = false @@ -97,63 +97,63 @@ function initialize_boundary_mass_fluxes(velocities::NamedTuple) w_bcs = w.boundary_conditions boundary_fluxes = NamedTuple() - right_matching_scheme_boundaries = Symbol[] - left_matching_scheme_boundaries = Symbol[] - total_area_matching_scheme_boundaries = zero(eltype(u)) + right_scheme_boundaries = Symbol[] + left_scheme_boundaries = Symbol[] + total_area_scheme_boundaries = zero(eltype(u)) # Check west boundary (u velocity) west_flux_and_area = initialize_boundary_mass_flux(u, u_bcs.west, Val(:west)) boundary_fluxes = merge(boundary_fluxes, west_flux_and_area) if needs_mass_flux_correction(u_bcs.west) - push!(left_matching_scheme_boundaries, :west) - total_area_matching_scheme_boundaries += boundary_fluxes.west_area + push!(left_scheme_boundaries, :west) + total_area_scheme_boundaries += boundary_fluxes.west_area end # Check east boundary (u velocity) east_flux_and_area = initialize_boundary_mass_flux(u, u_bcs.east, Val(:east)) boundary_fluxes = merge(boundary_fluxes, east_flux_and_area) if needs_mass_flux_correction(u_bcs.east) - push!(right_matching_scheme_boundaries, :east) - total_area_matching_scheme_boundaries += boundary_fluxes.east_area + push!(right_scheme_boundaries, :east) + total_area_scheme_boundaries += boundary_fluxes.east_area end # Check south boundary (v velocity) south_flux_and_area = initialize_boundary_mass_flux(v, v_bcs.south, Val(:south)) boundary_fluxes = merge(boundary_fluxes, south_flux_and_area) if needs_mass_flux_correction(v_bcs.south) - push!(left_matching_scheme_boundaries, :south) - total_area_matching_scheme_boundaries += boundary_fluxes.south_area + push!(left_scheme_boundaries, :south) + total_area_scheme_boundaries += boundary_fluxes.south_area end # Check north boundary (v velocity) north_flux_and_area = initialize_boundary_mass_flux(v, v_bcs.north, Val(:north)) boundary_fluxes = merge(boundary_fluxes, north_flux_and_area) if needs_mass_flux_correction(v_bcs.north) - push!(right_matching_scheme_boundaries, :north) - total_area_matching_scheme_boundaries += boundary_fluxes.north_area + push!(right_scheme_boundaries, :north) + total_area_scheme_boundaries += boundary_fluxes.north_area end # Check bottom boundary (w velocity) bottom_flux_and_area = initialize_boundary_mass_flux(w, w_bcs.bottom, Val(:bottom)) boundary_fluxes = merge(boundary_fluxes, bottom_flux_and_area) if needs_mass_flux_correction(w_bcs.bottom) - push!(left_matching_scheme_boundaries, :bottom) - total_area_matching_scheme_boundaries += boundary_fluxes.bottom_area + push!(left_scheme_boundaries, :bottom) + total_area_scheme_boundaries += boundary_fluxes.bottom_area end # Check top boundary (w velocity) top_flux_and_area = initialize_boundary_mass_flux(w, w_bcs.top, Val(:top)) boundary_fluxes = merge(boundary_fluxes, top_flux_and_area) if needs_mass_flux_correction(w_bcs.top) - push!(right_matching_scheme_boundaries, :top) - total_area_matching_scheme_boundaries += boundary_fluxes.top_area + push!(right_scheme_boundaries, :top) + total_area_scheme_boundaries += boundary_fluxes.top_area end - boundary_fluxes = merge(boundary_fluxes, (; left_matching_scheme_boundaries = Tuple(left_matching_scheme_boundaries), - right_matching_scheme_boundaries = Tuple(right_matching_scheme_boundaries), - total_area_matching_scheme_boundaries)) + boundary_fluxes = merge(boundary_fluxes, (; left_scheme_boundaries = Tuple(left_scheme_boundaries), + right_scheme_boundaries = Tuple(right_scheme_boundaries), + total_area_scheme_boundaries)) - if length(left_matching_scheme_boundaries) == 0 && length(right_matching_scheme_boundaries) == 0 + if length(left_scheme_boundaries) == 0 && length(right_scheme_boundaries) == 0 return nothing else return boundary_fluxes @@ -225,7 +225,7 @@ function enforce_open_boundary_mass_conservation!(model, boundary_mass_fluxes) u, v, w = model.velocities ∮udA = open_boundary_mass_inflow(model) - A = boundary_mass_fluxes.total_area_matching_scheme_boundaries + A = boundary_mass_fluxes.total_area_scheme_boundaries A⁻¹_∮udA = ∮udA / A diff --git a/validation/open_boundaries/oscillating_flow.jl b/validation/open_boundaries/oscillating_flow.jl index 40ecf57a682..027a2cb3e3b 100644 --- a/validation/open_boundaries/oscillating_flow.jl +++ b/validation/open_boundaries/oscillating_flow.jl @@ -5,7 +5,7 @@ # forcings and boundary conditions originally designed for `v` aere then used for `w` without # modification). # -# This case also has a stretched grid to validate the matching scheme on a stretched grid. +# This case also has a stretched grid to validate the scheme on a stretched grid. using Oceananigans, CairoMakie using Oceananigans.BoundaryConditions: FlatExtrapolationOpenBoundaryCondition, PerturbationAdvectionOpenBoundaryCondition @@ -123,7 +123,7 @@ function run_cylinder(grid, boundary_conditions; plot=true, stop_time = 50, simn end inflow_timescale = outflow_timescale = 1/4 -matching_scheme_name(obc) = string(nameof(typeof(obc.classification.matching_scheme))) +scheme_name(obc) = string(nameof(typeof(obc.classification.scheme))) for grid in (xygrid, xzgrid) u_fe = FlatExtrapolationOpenBoundaryCondition(u∞, parameters = (; U, T), relaxation_timescale = 1) @@ -146,10 +146,10 @@ for grid in (xygrid, xzgrid) for obcs in (feobcs, paobcs,) if grid isa Oceananigans.Grids.ZFlatGrid boundary_conditions = (u = obcs.u, v = obcs.v) - simname = "xy_" * matching_scheme_name(boundary_conditions.u.east) + simname = "xy_" * scheme_name(boundary_conditions.u.east) elseif grid isa Oceananigans.Grids.YFlatGrid boundary_conditions = (u = obcs.u, w = obcs.w) - simname = "xz_" * matching_scheme_name(boundary_conditions.u.east) + simname = "xz_" * scheme_name(boundary_conditions.u.east) end @info "Running $simname" run_cylinder(grid, boundary_conditions, simname = simname, stop_time = T) From 99b507b7bae028743dfca6c9d215131a00d21f32 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Fri, 8 Aug 2025 16:55:03 -0700 Subject: [PATCH 04/18] revert some undesirable changes --- ...flat_extrapolation_open_boundary_scheme.jl | 6 +-- ...urbation_advection_open_boundary_scheme.jl | 52 +++++++++++++++++-- 2 files changed, 49 insertions(+), 9 deletions(-) diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl index 8365299050b..dfa1512defe 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl @@ -8,7 +8,7 @@ Zero gradient perpendicular velocity boundary condition. We find the boundary value by Taylor expanding the gradient at the boundary point (`xᵢ`) to second order: ```math -f′(xᵢ) ≈ f′(xᵢ₋₁) + f′′(xᵢ₋₁)(xᵢ₋₁ - xᵢ) + O(Δx²), +f′(xᵢ) ≈ f′(xᵢ₋₁) + f′′(xᵢ₋₁)(xᵢ₋₁ - xᵢ) + O(Δx²) = f′(xᵢ₋₁) + f′′(xᵢ₋₁)Δx + O(Δx²), ``` where ``Δx=xᵢ₋₁ - xᵢ`` (for simplicity, we will also assume the spacing is constant at all ``i`` for now). @@ -147,6 +147,4 @@ end @inbounds ϕ[i, j, k] = relax(i, j, grid, gradient_free_ϕ, bc, clock, model_fields) return nothing -end - - +end \ No newline at end of file diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl index 2f658eeef48..13412436379 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl @@ -7,16 +7,60 @@ using Oceananigans: defaults For cases where we assume that the internal flow is a small perturbation from an external prescribed or coarser flow, we can split the velocity into background and perturbation components. + +We begin with the equation governing the fluid in the interior: + ∂ₜu + u⋅∇u = −∇P + F, +and note that on the boundary the pressure gradient is zero. +We can then assume that the flow composes of mean (U⃗) and pertubation (u⃗′) components, +and considering the x-component of velocity, we can rewrite the equation as + ∂ₜu₁ = -u₁∂₁u - u₂∂₂u₁ - u₃∂₃u₁ + F₁ ≈ - U₁∂₁u₁′ - U₂∂₂u₁′ - U₃∂₃u₁′ + F. + +Simplify by assuming that U⃗ = Ux̂, an then take a numerical step to find u₁. + +When the boundaries are filled the interior is at time tₙ₊₁ so we can take +a backwards euler step (in the case that the mean flow is boundary normal) on a right boundary: + (Uⁿ⁺¹ - Uⁿ) / Δt + (u′ⁿ⁺¹ - u′ⁿ) / Δt = - Uⁿ⁺¹ (u′ⁿ⁺¹ᵢ - u′ⁿ⁺¹ᵢ₋₁) / Δx + Fᵤ. + +This can not be solved for general forcing, but if we assume the dominant forcing is +relaxation to the mean velocity (i.e. u′→0) then Fᵤ = -u′ / τ then we can find u′ⁿ⁺¹: + u′ⁿ⁺¹ = (uⁿ + Ũu′ⁿ⁺¹ᵢ₋₁ - Uⁿ⁺¹) / (1 + Ũ + Δt/τ), + +where Ũ = U Δt / Δx, then uⁿ⁺¹ is: + uⁿ⁺¹ = (uᵢⁿ + Ũuᵢ₋₁ⁿ⁺¹ + Uⁿ⁺¹τ̃) / (1 + τ̃ + Ũ) + +where τ̃ = Δt/τ. + +The same operation can be repeated for left boundaries. + +The relaxation timescale ``τ`` can be set to different values depending on whether +``U`` points in or out of the domain (`inflow_timescale`/`outflow_timescale`). Since the +scheme is only valid when the flow is directed out of the domain the boundary condition +falls back to relaxation to the prescribed value. By default this happens instantly but +if the direction varies this may not be preferable. It is beneficial to relax the outflow +(i.e. non-zero `outflow_timescale`) to reduce the shock when the flow changes direction +to point into the domain. + +The ideal value of the timescales probably depend on the grid spacing and details of the +boundary flow. """ struct PerturbationAdvection{FT} inflow_timescale :: FT - outflow_timescale :: FT + outflow_timescale :: FT end Adapt.adapt_structure(to, pe::PerturbationAdvection) = PerturbationAdvection(adapt(to, pe.inflow_timescale), adapt(to, pe.outflow_timescale)) +""" + PerturbationAdvectionOpenBoundaryCondition(val, FT = defaults.FloatType; + outflow_timescale = Inf, + inflow_timescale = 0, kwargs...) + +Creates a `PerturbationAdvectionOpenBoundaryCondition` with a given exterior value `val`, to which +the flow is forced with an `outflow_timescale` for outflow and `inflow_timescale` for inflow. For +details about this method, refer to the docstring for `PerturbationAdvection`. +""" function PerturbationAdvectionOpenBoundaryCondition(val, FT = defaults.FloatType; outflow_timescale = Inf, inflow_timescale = 0, kwargs...) @@ -43,7 +87,7 @@ const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} pa = bc.classification.scheme τ = ifelse(ūⁿ⁺¹ >= 0, pa.outflow_timescale, pa.inflow_timescale) - τ̃ = Δt / τ + τ̃ = Δt / τ # last stage Δt normalized by the inflow/output timescale relaxed_uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + U) uᵢⁿ⁺¹ = ifelse(τ == 0, ūⁿ⁺¹, relaxed_uᵢⁿ⁺¹) @@ -67,7 +111,7 @@ end pa = bc.classification.scheme τ = ifelse(ūⁿ⁺¹ <= 0, pa.outflow_timescale, pa.inflow_timescale) - τ̃ = Δt / τ + τ̃ = Δt / τ # last stage Δt normalized by the inflow/output timescale relaxed_u₁ⁿ⁺¹ = (uᵢⁿ - U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ - U) u₁ⁿ⁺¹ = ifelse(τ == 0, ūⁿ⁺¹, relaxed_u₁ⁿ⁺¹) @@ -139,5 +183,3 @@ end return nothing end - - From e2c51a566252e882cdc834f1cd3d2bf81d7d4ab5 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Sat, 9 Aug 2025 15:42:26 -0700 Subject: [PATCH 05/18] mograte to new interface --- src/BoundaryConditions/BoundaryConditions.jl | 1 + src/BoundaryConditions/boundary_condition.jl | 2 +- ...flat_extrapolation_open_boundary_scheme.jl | 13 +++---- ...urbation_advection_open_boundary_scheme.jl | 31 ++++++--------- src/Oceananigans.jl | 1 + test/test_boundary_conditions_integration.jl | 38 +++++++++---------- .../test_conjugate_gradient_poisson_solver.jl | 5 +-- 7 files changed, 40 insertions(+), 51 deletions(-) diff --git a/src/BoundaryConditions/BoundaryConditions.jl b/src/BoundaryConditions/BoundaryConditions.jl index 7bd24263334..38c9c3d7fb9 100644 --- a/src/BoundaryConditions/BoundaryConditions.jl +++ b/src/BoundaryConditions/BoundaryConditions.jl @@ -5,6 +5,7 @@ export BoundaryCondition, getbc, setbc!, PeriodicBoundaryCondition, OpenBoundaryCondition, NoFluxBoundaryCondition, MultiRegionCommunicationBoundaryCondition, FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, DistributedCommunicationBoundaryCondition, + FlatExtrapolation, PerturbationAdvection, validate_boundary_condition_topology, validate_boundary_condition_architecture, FieldBoundaryConditions, compute_x_bcs!, compute_y_bcs!, compute_z_bcs!, diff --git a/src/BoundaryConditions/boundary_condition.jl b/src/BoundaryConditions/boundary_condition.jl index 0d5d877d4e5..696c4c6428a 100644 --- a/src/BoundaryConditions/boundary_condition.jl +++ b/src/BoundaryConditions/boundary_condition.jl @@ -96,7 +96,7 @@ MultiRegionCommunicationBoundaryCondition() = BoundaryCondition(MultiRegionCommu FluxBoundaryCondition(val; kwargs...) = BoundaryCondition(Flux(), val; kwargs...) ValueBoundaryCondition(val; kwargs...) = BoundaryCondition(Value(), val; kwargs...) GradientBoundaryCondition(val; kwargs...) = BoundaryCondition(Gradient(), val; kwargs...) - OpenBoundaryCondition(val; kwargs...) = BoundaryCondition(Open(nothing), val; kwargs...) + OpenBoundaryCondition(val; scheme = nothing, kwargs...) = BoundaryCondition(Open(scheme), val; kwargs...) MultiRegionCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(MultiRegionCommunication(), val; kwargs...) ZipperBoundaryCondition(val; kwargs...) = BoundaryCondition(Zipper(), val; kwargs...) DistributedCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(DistributedCommunication(), val; kwargs...) diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl index dfa1512defe..0676e7d261e 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl @@ -1,4 +1,5 @@ using Oceananigans.Operators: Δxᶜᶜᶜ, Δyᶜᶜᶜ, Δzᶜᶜᶜ +using Oceananigans: defaults """ FlatExtrapolation @@ -33,17 +34,13 @@ of 1/2 changes to ``Δx₋₁/(Δx₋₂ + Δx₋₃)`` instead, i.e.: f(xᵢ) ≈ f(xᵢ₋₂) - (f(xᵢ₋₁) - f(xᵢ₋₃))Δxᵢ₋₁/(Δxᵢ₋₂ + Δxᵢ₋₃) + O(Δx²) ```. """ -struct FlatExtrapolation{FT} - relaxation_timescale :: FT +@kwdef struct FlatExtrapolation{FT} + relaxation_timescale :: FT = Inf end -const FEOBC = BoundaryCondition{<:Open{<:FlatExtrapolation}} - -function FlatExtrapolationOpenBoundaryCondition(val = nothing; relaxation_timescale = Inf, kwargs...) - classification = Open(FlatExtrapolation(relaxation_timescale)) +Adapt.adapt_structure(to, fe::FlatExtrapolation) = FlatExtrapolation(adapt(to, fe.relaxation_timescale)) - return BoundaryCondition(classification, val; kwargs...) -end +const FEOBC = BoundaryCondition{<:Open{<:FlatExtrapolation}} @inline function relax(l, m, grid, ϕ, bc, clock, model_fields) Δt = clock.last_stage_Δt diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl index 13412436379..cbd18885bfb 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl @@ -4,6 +4,10 @@ using Oceananigans: defaults """ PerturbationAdvection +Create a `PerturbationAdvection` scheme to be used with an `OpenBoundaryCondition`. +This scheme will nudge the boundary velocity to the OpenBoundaryCondition's exterior value `val`, +using a time-scale `inflow_timescale` for inflow and `outflow_timescale` for outflow. + For cases where we assume that the internal flow is a small perturbation from an external prescribed or coarser flow, we can split the velocity into background and perturbation components. @@ -48,29 +52,18 @@ struct PerturbationAdvection{FT} outflow_timescale :: FT end -Adapt.adapt_structure(to, pe::PerturbationAdvection) = - PerturbationAdvection(adapt(to, pe.inflow_timescale), - adapt(to, pe.outflow_timescale)) - -""" - PerturbationAdvectionOpenBoundaryCondition(val, FT = defaults.FloatType; - outflow_timescale = Inf, - inflow_timescale = 0, kwargs...) - -Creates a `PerturbationAdvectionOpenBoundaryCondition` with a given exterior value `val`, to which -the flow is forced with an `outflow_timescale` for outflow and `inflow_timescale` for inflow. For -details about this method, refer to the docstring for `PerturbationAdvection`. -""" -function PerturbationAdvectionOpenBoundaryCondition(val, FT = defaults.FloatType; - outflow_timescale = Inf, - inflow_timescale = 0, kwargs...) +function PerturbationAdvection(FT = defaults.FloatType; + outflow_timescale = Inf, + inflow_timescale = 0) inflow_timescale = convert(FT, inflow_timescale) outflow_timescale = convert(FT, outflow_timescale) - classification = Open(PerturbationAdvection(inflow_timescale, outflow_timescale)) - - return BoundaryCondition(classification, val; kwargs...) + return PerturbationAdvection(inflow_timescale, outflow_timescale) end +Adapt.adapt_structure(to, pe::PerturbationAdvection) = + PerturbationAdvection(adapt(to, pe.inflow_timescale), + adapt(to, pe.outflow_timescale)) + const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} @inline function step_right_boundary!(bc::PAOBC, l, m, boundary_indices, boundary_adjacent_indices, diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index 4efb7f7c37b..714c3f2f100 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -35,6 +35,7 @@ export # Boundary conditions BoundaryCondition, FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, OpenBoundaryCondition, + FlatExtrapolation, PerturbationAdvection, FieldBoundaryConditions, # Fields and field manipulation diff --git a/test/test_boundary_conditions_integration.jl b/test/test_boundary_conditions_integration.jl index 100162da02f..4ae842947b4 100644 --- a/test/test_boundary_conditions_integration.jl +++ b/test/test_boundary_conditions_integration.jl @@ -1,8 +1,6 @@ include("dependencies_for_runtests.jl") using Oceananigans.BoundaryConditions: ContinuousBoundaryFunction, - FlatExtrapolationOpenBoundaryCondition, - PerturbationAdvectionOpenBoundaryCondition fill_halo_regions! using Oceananigans: prognostic_fields @@ -108,22 +106,22 @@ function fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT) end left_febc(::Val{1}, grid, loc) = FieldBoundaryConditions(grid, loc, east = OpenBoundaryCondition(1), - west = FlatExtrapolationOpenBoundaryCondition()) + west = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) right_febc(::Val{1}, grid, loc) = FieldBoundaryConditions(grid, loc, west = OpenBoundaryCondition(1), - east = FlatExtrapolationOpenBoundaryCondition()) + east = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) left_febc(::Val{2}, grid, loc) = FieldBoundaryConditions(grid, loc, north = OpenBoundaryCondition(1), - south = FlatExtrapolationOpenBoundaryCondition()) + south = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) right_febc(::Val{2}, grid, loc) = FieldBoundaryConditions(grid, loc, south = OpenBoundaryCondition(1), - north = FlatExtrapolationOpenBoundaryCondition()) + north = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) left_febc(::Val{3}, grid, loc) = FieldBoundaryConditions(grid, loc, top = OpenBoundaryCondition(1), - bottom = FlatExtrapolationOpenBoundaryCondition()) + bottom = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) right_febc(::Val{3}, grid, loc) = FieldBoundaryConditions(grid, loc, bottom = OpenBoundaryCondition(1), - top = FlatExtrapolationOpenBoundaryCondition()) + top = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) end_position(::Val{1}, grid) = (grid.Nx+1, 1, 1) end_position(::Val{2}, grid) = (1, grid.Ny+1, 1) @@ -182,7 +180,7 @@ function test_perturbation_advection_open_boundary_conditions(arch, FT) grid = RectilinearGrid(arch, FT; topology, size = (4, ), x = (0, 4), y = (0, 4), z = (0, 4), halo = (1, )) - obc = PerturbationAdvectionOpenBoundaryCondition(-1, inflow_timescale = 10.0) + obc = OpenBoundaryCondition(-1, scheme = PerturbationAdvection(inflow_timescale = 10.0)) boundary_conditions = wall_normal_boundary_condition(Val(orientation), obc) model = NonhydrostaticModel(; grid, boundary_conditions, timestepper = :QuasiAdamsBashforth2) @@ -195,7 +193,7 @@ function test_perturbation_advection_open_boundary_conditions(arch, FT) @test all(view(parent(u), :, :, :) .== -1) @test all(interior(u) .== -1) - obc = PerturbationAdvectionOpenBoundaryCondition((t) -> 0.1*t, inflow_timescale = 0.01, outflow_timescale = 0.5) + obc = OpenBoundaryCondition((t) -> 0.1*t, scheme = PerturbationAdvection(inflow_timescale = 0.01, outflow_timescale = 0.5)) forcing = velocity_forcing(Val(orientation), Forcing((x, t) -> 0.1)) boundary_conditions = wall_normal_boundary_condition(Val(orientation), obc) @@ -427,27 +425,27 @@ test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType), test_flat_extrapolation_open_boundary_conditions(arch, FT) test_perturbation_advection_open_boundary_conditions(arch, FT) - # Only PerturbationAdvectionOpenBoundaryCondition + # Only PerturbationAdvection OpenBoundaryCondition U₀ = 1 inflow_timescale = 1e-1 outflow_timescale = Inf - u_bcs = FieldBoundaryConditions(west = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale), - east = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale)) + u_bcs = FieldBoundaryConditions(west = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale)), + east = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale))) boundary_conditions = (; u = u_bcs) test_open_boundary_condition_mass_conservation(arch, FT, boundary_conditions) - # Only FlatExtrapolationOpenBoundaryCondition - u_bcs = FieldBoundaryConditions(west = FlatExtrapolationOpenBoundaryCondition(), - east = FlatExtrapolationOpenBoundaryCondition()) + # Only FlatExtrapolation OpenBoundaryCondition + u_bcs = FieldBoundaryConditions(west = OpenBoundaryCondition(U₀; scheme = FlatExtrapolation()), + east = OpenBoundaryCondition(U₀; scheme = FlatExtrapolation())) boundary_conditions = (; u = u_bcs) test_open_boundary_condition_mass_conservation(arch, FT, boundary_conditions) # Mixed open boundary conditions - u_bcs = FieldBoundaryConditions(west = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale), - east = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale)) - v_bcs = FieldBoundaryConditions(south = FlatExtrapolationOpenBoundaryCondition(), - north = FlatExtrapolationOpenBoundaryCondition()) + u_bcs = FieldBoundaryConditions(west = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale)), + east = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale))) + v_bcs = FieldBoundaryConditions(south = OpenBoundaryCondition(0; scheme = FlatExtrapolation()), + north = OpenBoundaryCondition(0; scheme = FlatExtrapolation())) w_bcs = FieldBoundaryConditions(bottom = OpenBoundaryCondition(0), top = OpenBoundaryCondition(0)) diff --git a/test/test_conjugate_gradient_poisson_solver.jl b/test/test_conjugate_gradient_poisson_solver.jl index 1442304b583..600ce7cc0dc 100644 --- a/test/test_conjugate_gradient_poisson_solver.jl +++ b/test/test_conjugate_gradient_poisson_solver.jl @@ -3,7 +3,6 @@ include("dependencies_for_poisson_solvers.jl") using Oceananigans.Solvers: fft_poisson_solver, ConjugateGradientPoissonSolver, DiagonallyDominantPreconditioner using Oceananigans.TimeSteppers: compute_pressure_correction! -using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition using Oceananigans.Grids: XYZRegularRG using LinearAlgebra: norm using Random: seed! @@ -183,8 +182,8 @@ function test_CGSolver_with_immersed_boundary_and_open_boundaries(underlying_gri U = 1 inflow_timescale = 1e-4 outflow_timescale = Inf - u_boundaries = FieldBoundaryConditions(west = PerturbationAdvectionOpenBoundaryCondition(U; inflow_timescale, outflow_timescale), - east = PerturbationAdvectionOpenBoundaryCondition(U; inflow_timescale, outflow_timescale)) + u_boundaries = FieldBoundaryConditions(west = OpenBoundaryCondition(U; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale)), + east = OpenBoundaryCondition(U; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale))) model = NonhydrostaticModel(grid = grid, boundary_conditions = (u = u_boundaries,), From f1b229e4ee0d9f0a4d5bc2d84c9e2784982fe8e0 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Sat, 9 Aug 2025 15:44:50 -0700 Subject: [PATCH 06/18] cleanup --- .../flat_extrapolation_open_boundary_scheme.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl index 0676e7d261e..5646245b185 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl @@ -1,5 +1,4 @@ using Oceananigans.Operators: Δxᶜᶜᶜ, Δyᶜᶜᶜ, Δzᶜᶜᶜ -using Oceananigans: defaults """ FlatExtrapolation From 254b55510aae5b91d9d2522ce7cb567596639693 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Mon, 11 Aug 2025 20:16:45 -0700 Subject: [PATCH 07/18] navid's suggestions --- .../perturbation_advection_open_boundary_scheme.jl | 14 ++++++++------ test/test_boundary_conditions_integration.jl | 2 +- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl index cbd18885bfb..7782e296880 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl @@ -1,8 +1,15 @@ using Oceananigans.Operators: Δxᶠᶜᶜ, Δyᶜᶠᶜ, Δzᶜᶜᶠ, Ax_qᶠᶜᶜ, Ay_qᶜᶠᶜ, Az_qᶜᶜᶠ using Oceananigans: defaults +struct PerturbationAdvection{FT} + inflow_timescale :: FT + outflow_timescale :: FT +end + """ - PerturbationAdvection + PerturbationAdvection(FT = defaults.FloatType; + outflow_timescale = Inf, + inflow_timescale = 0) Create a `PerturbationAdvection` scheme to be used with an `OpenBoundaryCondition`. This scheme will nudge the boundary velocity to the OpenBoundaryCondition's exterior value `val`, @@ -47,11 +54,6 @@ to point into the domain. The ideal value of the timescales probably depend on the grid spacing and details of the boundary flow. """ -struct PerturbationAdvection{FT} - inflow_timescale :: FT - outflow_timescale :: FT -end - function PerturbationAdvection(FT = defaults.FloatType; outflow_timescale = Inf, inflow_timescale = 0) diff --git a/test/test_boundary_conditions_integration.jl b/test/test_boundary_conditions_integration.jl index 4ae842947b4..33664ed532e 100644 --- a/test/test_boundary_conditions_integration.jl +++ b/test/test_boundary_conditions_integration.jl @@ -193,7 +193,7 @@ function test_perturbation_advection_open_boundary_conditions(arch, FT) @test all(view(parent(u), :, :, :) .== -1) @test all(interior(u) .== -1) - obc = OpenBoundaryCondition((t) -> 0.1*t, scheme = PerturbationAdvection(inflow_timescale = 0.01, outflow_timescale = 0.5)) + obc = OpenBoundaryCondition(t -> 0.1*t, scheme = PerturbationAdvection(inflow_timescale = 0.01, outflow_timescale = 0.5)) forcing = velocity_forcing(Val(orientation), Forcing((x, t) -> 0.1)) boundary_conditions = wall_normal_boundary_condition(Val(orientation), obc) From c4bb18bb78b7cabd68e23fd97ee9c587f42bfb10 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Mon, 11 Aug 2025 20:20:49 -0700 Subject: [PATCH 08/18] bump minor version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 76733753160..ae89bb3f937 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Oceananigans" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" authors = ["Climate Modeling Alliance and contributors"] -version = "0.97.7" +version = "0.98.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From 6097a24ddd81141b236b658f94faf6a82d8b9269 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Thu, 14 Aug 2025 17:33:46 -0700 Subject: [PATCH 09/18] remove FlatExtrapolation --- src/BoundaryConditions/BoundaryConditions.jl | 3 +- ...flat_extrapolation_open_boundary_scheme.jl | 146 ------------------ .../boundary_mass_fluxes.jl | 2 +- src/Oceananigans.jl | 2 +- test/test_boundary_conditions_integration.jl | 68 -------- validation/open_boundaries/cylinder.jl | 8 +- .../open_boundaries/oscillating_flow.jl | 13 +- 7 files changed, 8 insertions(+), 234 deletions(-) delete mode 100644 src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl diff --git a/src/BoundaryConditions/BoundaryConditions.jl b/src/BoundaryConditions/BoundaryConditions.jl index 38c9c3d7fb9..f93cfa1dfe8 100644 --- a/src/BoundaryConditions/BoundaryConditions.jl +++ b/src/BoundaryConditions/BoundaryConditions.jl @@ -5,7 +5,7 @@ export BoundaryCondition, getbc, setbc!, PeriodicBoundaryCondition, OpenBoundaryCondition, NoFluxBoundaryCondition, MultiRegionCommunicationBoundaryCondition, FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, DistributedCommunicationBoundaryCondition, - FlatExtrapolation, PerturbationAdvection, + PerturbationAdvection, validate_boundary_condition_topology, validate_boundary_condition_architecture, FieldBoundaryConditions, compute_x_bcs!, compute_y_bcs!, compute_z_bcs!, @@ -41,6 +41,5 @@ include("compute_flux_bcs.jl") include("update_boundary_conditions.jl") include("polar_boundary_condition.jl") -include("flat_extrapolation_open_boundary_scheme.jl") include("perturbation_advection_open_boundary_scheme.jl") end # module diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl deleted file mode 100644 index 5646245b185..00000000000 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl +++ /dev/null @@ -1,146 +0,0 @@ -using Oceananigans.Operators: Δxᶜᶜᶜ, Δyᶜᶜᶜ, Δzᶜᶜᶜ - -""" - FlatExtrapolation - -Zero gradient perpendicular velocity boundary condition. - -We find the boundary value by Taylor expanding the gradient at the boundary point (`xᵢ`) -to second order: -```math -f′(xᵢ) ≈ f′(xᵢ₋₁) + f′′(xᵢ₋₁)(xᵢ₋₁ - xᵢ) + O(Δx²) = f′(xᵢ₋₁) + f′′(xᵢ₋₁)Δx + O(Δx²), -``` -where ``Δx=xᵢ₋₁ - xᵢ`` (for simplicity, we will also assume the spacing is constant at -all ``i`` for now). -We can substitute the gradient at some point ``j`` (``f′(xⱼ)``) with the central -difference approximation: -```math -f′(xⱼ) ≈ (f(xⱼ₊₁) - f(xⱼ₋₁)) / 2Δx, -``` -and the second derivative at some point ``j`` (``f′′(xⱼ)``) can be approximated as: -```math -f′′(xⱼ) ≈ (f′(xⱼ₊₁) - f′(xⱼ₋₁)) / 2Δx = ((f(xⱼ₊₂) - f(xⱼ)) - (f(xⱼ) - f(xⱼ₋₂))) / (2Δx)². -``` -When we then substitute for the boundary adjacent point ``f′′(xᵢ₋₁)`` we know that -``f′(xⱼ₊₁)=f′(xᵢ)=0`` so the Taylor expansion becomes: -```math -f(xᵢ) ≈ f(xᵢ₋₂) - (f(xᵢ₋₁) - f(xᵢ₋₃))/2 + O(Δx²). -``` - -When the grid spacing is not constant the above can be repeated resulting in the factor -of 1/2 changes to ``Δx₋₁/(Δx₋₂ + Δx₋₃)`` instead, i.e.: -```math -f(xᵢ) ≈ f(xᵢ₋₂) - (f(xᵢ₋₁) - f(xᵢ₋₃))Δxᵢ₋₁/(Δxᵢ₋₂ + Δxᵢ₋₃) + O(Δx²) -```. -""" -@kwdef struct FlatExtrapolation{FT} - relaxation_timescale :: FT = Inf -end - -Adapt.adapt_structure(to, fe::FlatExtrapolation) = FlatExtrapolation(adapt(to, fe.relaxation_timescale)) - -const FEOBC = BoundaryCondition{<:Open{<:FlatExtrapolation}} - -@inline function relax(l, m, grid, ϕ, bc, clock, model_fields) - Δt = clock.last_stage_Δt - τ = bc.classification.scheme.relaxation_timescale - - Δt̄ = min(1, Δt / τ) - ϕₑₓₜ = getbc(bc, l, m, grid, clock, model_fields) - - Δϕ = (ϕₑₓₜ - ϕ) * Δt̄ - not_relaxing = isnothing(bc.condition) | !isfinite(clock.last_stage_Δt) - Δϕ = ifelse(not_relaxing, zero(ϕ), Δϕ) - - return ϕ + Δϕ -end - -@inline function _fill_west_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - Δx₁ = Δxᶜᶜᶜ(1, j, k, grid) - Δx₂ = Δxᶜᶜᶜ(2, j, k, grid) - Δx₃ = Δxᶜᶜᶜ(3, j, k, grid) - - spacing_factor = Δx₁ / (Δx₂ + Δx₃) - - gradient_free_ϕ = @inbounds ϕ[3, j, k] - (ϕ[2, j, k] - ϕ[4, j, k]) * spacing_factor - - @inbounds ϕ[1, j, k] = relax(j, k, grid, gradient_free_ϕ, bc, clock, model_fields) - - return nothing -end - -@inline function _fill_east_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - i = grid.Nx + 1 - - Δx₁ = Δxᶜᶜᶜ(i-1, j, k, grid) - Δx₂ = Δxᶜᶜᶜ(i-2, j, k, grid) - Δx₃ = Δxᶜᶜᶜ(i-3, j, k, grid) - - spacing_factor = Δx₁ / (Δx₂ + Δx₃) - - gradient_free_ϕ = @inbounds ϕ[i - 2, j, k] - (ϕ[i - 1, j, k] - ϕ[i - 3, j, k]) * spacing_factor - - @inbounds ϕ[i, j, k] = relax(j, k, grid, gradient_free_ϕ, bc, clock, model_fields) - - return nothing -end - -@inline function _fill_south_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - Δy₁ = Δyᶜᶜᶜ(i, 1, k, grid) - Δy₂ = Δyᶜᶜᶜ(i, 2, k, grid) - Δy₃ = Δyᶜᶜᶜ(i, 3, k, grid) - - spacing_factor = Δy₁ / (Δy₂ + Δy₃) - - gradient_free_ϕ = ϕ[i, 3, k] - (ϕ[i, 2, k] - ϕ[i, 4, k]) * spacing_factor - - @inbounds ϕ[i, 1, k] = relax(i, k, grid, gradient_free_ϕ, bc, clock, model_fields) - - return nothing -end - -@inline function _fill_north_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - j = grid.Ny + 1 - - Δy₁ = Δyᶜᶜᶜ(i, j-1, k, grid) - Δy₂ = Δyᶜᶜᶜ(i, j-2, k, grid) - Δy₃ = Δyᶜᶜᶜ(i, j-3, k, grid) - - spacing_factor = Δy₁ / (Δy₂ + Δy₃) - - gradient_free_ϕ = @inbounds ϕ[i, j - 2, k] - (ϕ[i, j - 1, k] - ϕ[i, j - 3, k]) * spacing_factor - - @inbounds ϕ[i, j, k] = relax(i, k, grid, gradient_free_ϕ, bc, clock, model_fields) - - return nothing -end - -@inline function _fill_bottom_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - Δz₁ = Δzᶜᶜᶜ(i, j, 1, grid) - Δz₂ = Δzᶜᶜᶜ(i, j, 2, grid) - Δz₃ = Δzᶜᶜᶜ(i, j, 3, grid) - - spacing_factor = Δz₁ / (Δz₂ + Δz₃) - - gradient_free_ϕ = @inbounds ϕ[i, j, 3] - (ϕ[i, j, 2] - ϕ[i, j, 4]) * spacing_factor - - @inbounds ϕ[i, j, 1] = relax(i, j, grid, gradient_free_ϕ, bc, clock, model_fields) - - return nothing -end - -@inline function _fill_top_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - k = grid.Nz + 1 - - Δz₁ = Δzᶜᶜᶜ(i, j, k-1, grid) - Δz₂ = Δzᶜᶜᶜ(i, j, k-2, grid) - Δz₃ = Δzᶜᶜᶜ(i, j, k-3, grid) - - spacing_factor = Δz₁ / (Δz₂ + Δz₃) - - gradient_free_ϕ = @inbounds ϕ[i, j, k - 2] - (ϕ[i, j, k - 1] - ϕ[i, j, k - 3]) * spacing_factor - - @inbounds ϕ[i, j, k] = relax(i, j, grid, gradient_free_ϕ, bc, clock, model_fields) - - return nothing -end \ No newline at end of file diff --git a/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl b/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl index a0121a9650c..ee2b1476bc6 100644 --- a/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl +++ b/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl @@ -1,4 +1,4 @@ -using Oceananigans.BoundaryConditions: BoundaryCondition, Open, PerturbationAdvection, FlatExtrapolation +using Oceananigans.BoundaryConditions: BoundaryCondition, Open, PerturbationAdvection using Oceananigans.AbstractOperations: Integral, Ax, Ay, Az, GridMetricOperation using Oceananigans.Fields: Field, interior, XFaceField, YFaceField, ZFaceField using GPUArraysCore: @allowscalar diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index 714c3f2f100..1c53a1fa477 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -35,7 +35,7 @@ export # Boundary conditions BoundaryCondition, FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, OpenBoundaryCondition, - FlatExtrapolation, PerturbationAdvection, + PerturbationAdvection, FieldBoundaryConditions, # Fields and field manipulation diff --git a/test/test_boundary_conditions_integration.jl b/test/test_boundary_conditions_integration.jl index 33664ed532e..7567406ffb6 100644 --- a/test/test_boundary_conditions_integration.jl +++ b/test/test_boundary_conditions_integration.jl @@ -105,62 +105,13 @@ function fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT) return isapprox(mean(b) - mean_b₀, flux * model.clock.time / Lz, atol=1e-6) end -left_febc(::Val{1}, grid, loc) = FieldBoundaryConditions(grid, loc, east = OpenBoundaryCondition(1), - west = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) -right_febc(::Val{1}, grid, loc) = FieldBoundaryConditions(grid, loc, west = OpenBoundaryCondition(1), - east = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) - -left_febc(::Val{2}, grid, loc) = FieldBoundaryConditions(grid, loc, north = OpenBoundaryCondition(1), - south = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) - -right_febc(::Val{2}, grid, loc) = FieldBoundaryConditions(grid, loc, south = OpenBoundaryCondition(1), - north = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) - -left_febc(::Val{3}, grid, loc) = FieldBoundaryConditions(grid, loc, top = OpenBoundaryCondition(1), - bottom = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) - -right_febc(::Val{3}, grid, loc) = FieldBoundaryConditions(grid, loc, bottom = OpenBoundaryCondition(1), - top = OpenBoundaryCondition(1, scheme = FlatExtrapolation())) end_position(::Val{1}, grid) = (grid.Nx+1, 1, 1) end_position(::Val{2}, grid) = (1, grid.Ny+1, 1) end_position(::Val{3}, grid) = (1, 1, grid.Nz+1) -function test_flat_extrapolation_open_boundary_conditions(arch, FT) - clock = Clock(; time = zero(FT)) - - for orientation in 1:3 - topology = tuple(map(n -> ifelse(n == orientation, Bounded, Flat), 1:3)...) - - normal_location = tuple(map(n -> ifelse(n == orientation, Face, Center), 1:3)...) - - grid = RectilinearGrid(arch, FT; topology, size = (16, ), x = (0, 1), y = (0, 1), z = (0, 1)) - - bcs1 = left_febc(Val(orientation), grid, normal_location) - bcs2 = right_febc(Val(orientation), grid, normal_location) - - u1 = Field{normal_location...}(grid; boundary_conditions = bcs1) - u2 = Field{normal_location...}(grid; boundary_conditions = bcs2) - - set!(u1, (X, ) -> 2-X) - set!(u2, (X, ) -> 1+X) - - fill_halo_regions!(u1, clock, (); fill_boundary_normal_velocities = false) - fill_halo_regions!(u2, clock, (); fill_boundary_normal_velocities = false) - # we can stop the wall normal halos being filled after the pressure solve - this serves more as a test of the general OBC stuff - @test interior(u1, 1, 1, 1) .== 2 - @test interior(u2, end_position(Val(orientation), grid)...) .== 2 - - fill_halo_regions!(u1, clock, ()) - fill_halo_regions!(u2, clock, ()) - - # now they should be filled - @test interior(u1, 1, 1, 1) .≈ 1.8125 - @test interior(u2, end_position(Val(orientation), grid)...) .≈ 1.8125 - end -end wall_normal_boundary_condition(::Val{1}, obc) = (; u = FieldBoundaryConditions(east = obc, west = obc)) wall_normal_boundary_condition(::Val{2}, obc) = (; v = FieldBoundaryConditions(south = obc, north = obc)) @@ -422,7 +373,6 @@ test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType), for arch in archs, FT in (Float64,) #float_types A = typeof(arch) @info " Testing open boundary conditions [$A, $FT]..." - test_flat_extrapolation_open_boundary_conditions(arch, FT) test_perturbation_advection_open_boundary_conditions(arch, FT) # Only PerturbationAdvection OpenBoundaryCondition @@ -434,24 +384,6 @@ test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType), east = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale))) boundary_conditions = (; u = u_bcs) test_open_boundary_condition_mass_conservation(arch, FT, boundary_conditions) - - # Only FlatExtrapolation OpenBoundaryCondition - u_bcs = FieldBoundaryConditions(west = OpenBoundaryCondition(U₀; scheme = FlatExtrapolation()), - east = OpenBoundaryCondition(U₀; scheme = FlatExtrapolation())) - boundary_conditions = (; u = u_bcs) - test_open_boundary_condition_mass_conservation(arch, FT, boundary_conditions) - - # Mixed open boundary conditions - u_bcs = FieldBoundaryConditions(west = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale)), - east = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale))) - v_bcs = FieldBoundaryConditions(south = OpenBoundaryCondition(0; scheme = FlatExtrapolation()), - north = OpenBoundaryCondition(0; scheme = FlatExtrapolation())) - w_bcs = FieldBoundaryConditions(bottom = OpenBoundaryCondition(0), - top = OpenBoundaryCondition(0)) - - boundary_conditions = (; u = u_bcs, v = v_bcs, w = w_bcs) - test_open_boundary_condition_mass_conservation(arch, FT, boundary_conditions) - end end end diff --git a/validation/open_boundaries/cylinder.jl b/validation/open_boundaries/cylinder.jl index 696ecae2050..bd9b60bb1a7 100644 --- a/validation/open_boundaries/cylinder.jl +++ b/validation/open_boundaries/cylinder.jl @@ -6,8 +6,7 @@ using Oceananigans.Solvers: FFTBasedPoissonSolver using Printf using CUDA -using Oceananigans.BoundaryConditions: FlatExtrapolationOpenBoundaryCondition, - PerturbationAdvectionOpenBoundaryCondition +using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition # there is some problem using ConjugateGradientPoissonSolver with TimeInterval because the timestep can go really small # so while I identify the issue I'm using IterationInterval and a fixed timestep @@ -214,14 +213,11 @@ end u∞ = 1 -feobcs = (west = OpenBoundaryCondition(u∞), east = FlatExtrapolationOpenBoundaryCondition(),) - paobcs = (west = PerturbationAdvectionOpenBoundaryCondition(u∞; inflow_timescale = 0.1, outflow_timescale = 0.1), east = PerturbationAdvectionOpenBoundaryCondition(u∞; inflow_timescale = 1/4, outflow_timescale = Inf), ) -obcs = (; flat_extrapolation=feobcs, - perturbation_advection=paobcs) +obcs = (; perturbation_advection=paobcs,) for (obc_name, obc) in pairs(obcs) @info "Running $(obc_name)" diff --git a/validation/open_boundaries/oscillating_flow.jl b/validation/open_boundaries/oscillating_flow.jl index 027a2cb3e3b..6573b181dbf 100644 --- a/validation/open_boundaries/oscillating_flow.jl +++ b/validation/open_boundaries/oscillating_flow.jl @@ -1,6 +1,6 @@ # This validation script shows open boundaries working in a simple case where a flow past a 2D # cylinder oscillates in two directions. All boundaries have the same -# `FlatExtrapolationOpenBoundaryCondition`s. This is similar to a more realistic case where we know +# open boundary conditions. This is similar to a more realistic case where we know # some arbitary external conditions. First we test an xy flow and then we test an xz flow (the # forcings and boundary conditions originally designed for `v` aere then used for `w` without # modification). @@ -8,7 +8,7 @@ # This case also has a stretched grid to validate the scheme on a stretched grid. using Oceananigans, CairoMakie -using Oceananigans.BoundaryConditions: FlatExtrapolationOpenBoundaryCondition, PerturbationAdvectionOpenBoundaryCondition +using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition @kwdef struct Cylinder{FT} D :: FT = 1.0 @@ -126,14 +126,7 @@ inflow_timescale = outflow_timescale = 1/4 scheme_name(obc) = string(nameof(typeof(obc.classification.scheme))) for grid in (xygrid, xzgrid) - u_fe = FlatExtrapolationOpenBoundaryCondition(u∞, parameters = (; U, T), relaxation_timescale = 1) - v_fe = FlatExtrapolationOpenBoundaryCondition(v∞, parameters = (; U, T), relaxation_timescale = 1) - w_fe = FlatExtrapolationOpenBoundaryCondition(v∞, parameters = (; U, T), relaxation_timescale = 1) - u_boundaries_fe = FieldBoundaryConditions(west = u_fe, east = u_fe) - v_boundaries_fe = FieldBoundaryConditions(south = v_fe, north = v_fe) - w_boundaries_fe = FieldBoundaryConditions(bottom = w_fe, top = w_fe) - feobcs = (u = u_boundaries_fe, v = v_boundaries_fe, w = w_boundaries_fe) u_boundaries_pa = FieldBoundaryConditions(west = PerturbationAdvectionOpenBoundaryCondition(u∞; parameters = (; U, T), inflow_timescale, outflow_timescale), east = PerturbationAdvectionOpenBoundaryCondition(u∞; parameters = (; U, T), inflow_timescale, outflow_timescale)) @@ -143,7 +136,7 @@ for grid in (xygrid, xzgrid) top = PerturbationAdvectionOpenBoundaryCondition(v∞; parameters = (; U, T), inflow_timescale, outflow_timescale)) paobcs = (u = u_boundaries_pa, v = v_boundaries_pa, w = w_boundaries_pa) - for obcs in (feobcs, paobcs,) + for obcs in (paobcs,) if grid isa Oceananigans.Grids.ZFlatGrid boundary_conditions = (u = obcs.u, v = obcs.v) simname = "xy_" * scheme_name(boundary_conditions.u.east) From b0f1da6a23afb356fea94cbb98249b005dfc438e Mon Sep 17 00:00:00 2001 From: tomaschor Date: Thu, 14 Aug 2025 17:39:19 -0700 Subject: [PATCH 10/18] migrate validation to new interface --- .../open_boundaries/2d_flow_over_bathymetry.jl | 6 +++--- validation/open_boundaries/cylinder.jl | 6 +++--- validation/open_boundaries/oscillating_flow.jl | 16 +++++++--------- 3 files changed, 13 insertions(+), 15 deletions(-) diff --git a/validation/open_boundaries/2d_flow_over_bathymetry.jl b/validation/open_boundaries/2d_flow_over_bathymetry.jl index 29b1e9e8b5e..f556526d65c 100644 --- a/validation/open_boundaries/2d_flow_over_bathymetry.jl +++ b/validation/open_boundaries/2d_flow_over_bathymetry.jl @@ -1,6 +1,6 @@ using Oceananigans using Oceananigans.Units -using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition, OpenBoundaryCondition +using Oceananigans.BoundaryConditions: PerturbationAdvection, OpenBoundaryCondition using Oceananigans.Diagnostics: AdvectiveCFL using Oceananigans.Solvers: ConjugateGradientPoissonSolver, FFTBasedPoissonSolver, DiagonallyDominantPreconditioner, fft_poisson_solver @@ -54,8 +54,8 @@ function create_mass_conservation_simulation(; #+++ Set boundary conditions if use_open_boundary_condition u_boundary_conditions = FieldBoundaryConditions( - west = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale), - east = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale) + west = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale)), + east = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale)) ) boundary_conditions = (; u = u_boundary_conditions) else diff --git a/validation/open_boundaries/cylinder.jl b/validation/open_boundaries/cylinder.jl index bd9b60bb1a7..78b306b2bd6 100644 --- a/validation/open_boundaries/cylinder.jl +++ b/validation/open_boundaries/cylinder.jl @@ -6,7 +6,7 @@ using Oceananigans.Solvers: FFTBasedPoissonSolver using Printf using CUDA -using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition +using Oceananigans.BoundaryConditions: PerturbationAdvection # there is some problem using ConjugateGradientPoissonSolver with TimeInterval because the timestep can go really small # so while I identify the issue I'm using IterationInterval and a fixed timestep @@ -213,8 +213,8 @@ end u∞ = 1 -paobcs = (west = PerturbationAdvectionOpenBoundaryCondition(u∞; inflow_timescale = 0.1, outflow_timescale = 0.1), - east = PerturbationAdvectionOpenBoundaryCondition(u∞; inflow_timescale = 1/4, outflow_timescale = Inf), +paobcs = (west = OpenBoundaryCondition(u∞; scheme = PerturbationAdvection(inflow_timescale = 0.1, outflow_timescale = 0.1)), + east = OpenBoundaryCondition(u∞; scheme = PerturbationAdvection(inflow_timescale = 1/4, outflow_timescale = Inf)), ) obcs = (; perturbation_advection=paobcs,) diff --git a/validation/open_boundaries/oscillating_flow.jl b/validation/open_boundaries/oscillating_flow.jl index 6573b181dbf..dd1036c73d2 100644 --- a/validation/open_boundaries/oscillating_flow.jl +++ b/validation/open_boundaries/oscillating_flow.jl @@ -8,7 +8,7 @@ # This case also has a stretched grid to validate the scheme on a stretched grid. using Oceananigans, CairoMakie -using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition +using Oceananigans.BoundaryConditions: PerturbationAdvection @kwdef struct Cylinder{FT} D :: FT = 1.0 @@ -126,14 +126,12 @@ inflow_timescale = outflow_timescale = 1/4 scheme_name(obc) = string(nameof(typeof(obc.classification.scheme))) for grid in (xygrid, xzgrid) - - - u_boundaries_pa = FieldBoundaryConditions(west = PerturbationAdvectionOpenBoundaryCondition(u∞; parameters = (; U, T), inflow_timescale, outflow_timescale), - east = PerturbationAdvectionOpenBoundaryCondition(u∞; parameters = (; U, T), inflow_timescale, outflow_timescale)) - v_boundaries_pa = FieldBoundaryConditions(south = PerturbationAdvectionOpenBoundaryCondition(v∞; parameters = (; U, T), inflow_timescale, outflow_timescale), - north = PerturbationAdvectionOpenBoundaryCondition(v∞; parameters = (; U, T), inflow_timescale, outflow_timescale)) - w_boundaries_pa = FieldBoundaryConditions(bottom = PerturbationAdvectionOpenBoundaryCondition(v∞; parameters = (; U, T), inflow_timescale, outflow_timescale), - top = PerturbationAdvectionOpenBoundaryCondition(v∞; parameters = (; U, T), inflow_timescale, outflow_timescale)) + u_boundaries_pa = FieldBoundaryConditions(west = OpenBoundaryCondition(u∞; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale), parameters = (; U, T)), + east = OpenBoundaryCondition(u∞; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale), parameters = (; U, T))) + v_boundaries_pa = FieldBoundaryConditions(south = OpenBoundaryCondition(v∞; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale), parameters = (; U, T)), + north = OpenBoundaryCondition(v∞; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale), parameters = (; U, T))) + w_boundaries_pa = FieldBoundaryConditions(bottom = OpenBoundaryCondition(v∞; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale), parameters = (; U, T)), + top = OpenBoundaryCondition(v∞; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale), parameters = (; U, T))) paobcs = (u = u_boundaries_pa, v = v_boundaries_pa, w = w_boundaries_pa) for obcs in (paobcs,) From d43b56befffd7ec99b8731e4c02d8e344f785413 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Fri, 15 Aug 2025 09:05:42 -0700 Subject: [PATCH 11/18] slightly better formatting --- validation/open_boundaries/cylinder.jl | 5 ++--- validation/open_boundaries/oscillating_flow.jl | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/validation/open_boundaries/cylinder.jl b/validation/open_boundaries/cylinder.jl index 78b306b2bd6..cf978e4d007 100644 --- a/validation/open_boundaries/cylinder.jl +++ b/validation/open_boundaries/cylinder.jl @@ -214,10 +214,9 @@ end u∞ = 1 paobcs = (west = OpenBoundaryCondition(u∞; scheme = PerturbationAdvection(inflow_timescale = 0.1, outflow_timescale = 0.1)), - east = OpenBoundaryCondition(u∞; scheme = PerturbationAdvection(inflow_timescale = 1/4, outflow_timescale = Inf)), - ) + east = OpenBoundaryCondition(u∞; scheme = PerturbationAdvection(inflow_timescale = 1/4, outflow_timescale = Inf))) -obcs = (; perturbation_advection=paobcs,) +obcs = (; perturbation_advection=paobcs) for (obc_name, obc) in pairs(obcs) @info "Running $(obc_name)" diff --git a/validation/open_boundaries/oscillating_flow.jl b/validation/open_boundaries/oscillating_flow.jl index dd1036c73d2..d9da4f535e7 100644 --- a/validation/open_boundaries/oscillating_flow.jl +++ b/validation/open_boundaries/oscillating_flow.jl @@ -134,7 +134,7 @@ for grid in (xygrid, xzgrid) top = OpenBoundaryCondition(v∞; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale), parameters = (; U, T))) paobcs = (u = u_boundaries_pa, v = v_boundaries_pa, w = w_boundaries_pa) - for obcs in (paobcs,) + for obcs in tuple(paobcs) if grid isa Oceananigans.Grids.ZFlatGrid boundary_conditions = (u = obcs.u, v = obcs.v) simname = "xy_" * scheme_name(boundary_conditions.u.east) From e46c6e23d27296d69d0d6261a3871400fe82915d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1s=20Chor?= Date: Wed, 20 Aug 2025 08:55:13 -0700 Subject: [PATCH 12/18] Update src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl --- src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl b/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl index ee2b1476bc6..9715f803be6 100644 --- a/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl +++ b/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl @@ -153,7 +153,7 @@ function initialize_boundary_mass_fluxes(velocities::NamedTuple) right_scheme_boundaries = Tuple(right_scheme_boundaries), total_area_scheme_boundaries)) - if length(left_scheme_boundaries) == 0 && length(right_scheme_boundaries) == 0 + if length(boundary_fluxes.left_scheme_boundaries) == 0 && length(boundary_fluxes.right_scheme_boundaries) == 0 return nothing else return boundary_fluxes From 9f9b1c48701036881e96eb85866dcee4d9e3527f Mon Sep 17 00:00:00 2001 From: tomaschor Date: Tue, 26 Aug 2025 10:58:05 -0700 Subject: [PATCH 13/18] rename perturbation advection file --- src/BoundaryConditions/BoundaryConditions.jl | 2 +- ...open_boundary_scheme.jl => perturbation_advection_scheme.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/BoundaryConditions/{perturbation_advection_open_boundary_scheme.jl => perturbation_advection_scheme.jl} (100%) diff --git a/src/BoundaryConditions/BoundaryConditions.jl b/src/BoundaryConditions/BoundaryConditions.jl index f93cfa1dfe8..20d9ee2c97e 100644 --- a/src/BoundaryConditions/BoundaryConditions.jl +++ b/src/BoundaryConditions/BoundaryConditions.jl @@ -41,5 +41,5 @@ include("compute_flux_bcs.jl") include("update_boundary_conditions.jl") include("polar_boundary_condition.jl") -include("perturbation_advection_open_boundary_scheme.jl") +include("perturbation_advection_open_scheme.jl") end # module diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl b/src/BoundaryConditions/perturbation_advection_scheme.jl similarity index 100% rename from src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl rename to src/BoundaryConditions/perturbation_advection_scheme.jl From 3b9c451009c495623c57a081680e327861a782e5 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Tue, 26 Aug 2025 11:12:13 -0700 Subject: [PATCH 14/18] oops, typo --- src/BoundaryConditions/BoundaryConditions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BoundaryConditions/BoundaryConditions.jl b/src/BoundaryConditions/BoundaryConditions.jl index 20d9ee2c97e..bbfd2551047 100644 --- a/src/BoundaryConditions/BoundaryConditions.jl +++ b/src/BoundaryConditions/BoundaryConditions.jl @@ -41,5 +41,5 @@ include("compute_flux_bcs.jl") include("update_boundary_conditions.jl") include("polar_boundary_condition.jl") -include("perturbation_advection_open_scheme.jl") +include("perturbation_advection_scheme.jl") end # module From 8ec408430a33b7843677aeb3c17ab7ce207e311e Mon Sep 17 00:00:00 2001 From: tomaschor Date: Tue, 26 Aug 2025 11:45:01 -0700 Subject: [PATCH 15/18] change file name to perturbation_advection.jl --- src/BoundaryConditions/BoundaryConditions.jl | 2 +- ...rturbation_advection_scheme.jl => perturbation_advection.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/BoundaryConditions/{perturbation_advection_scheme.jl => perturbation_advection.jl} (100%) diff --git a/src/BoundaryConditions/BoundaryConditions.jl b/src/BoundaryConditions/BoundaryConditions.jl index bbfd2551047..1d00223a66f 100644 --- a/src/BoundaryConditions/BoundaryConditions.jl +++ b/src/BoundaryConditions/BoundaryConditions.jl @@ -41,5 +41,5 @@ include("compute_flux_bcs.jl") include("update_boundary_conditions.jl") include("polar_boundary_condition.jl") -include("perturbation_advection_scheme.jl") +include("perturbation_advection.jl") end # module diff --git a/src/BoundaryConditions/perturbation_advection_scheme.jl b/src/BoundaryConditions/perturbation_advection.jl similarity index 100% rename from src/BoundaryConditions/perturbation_advection_scheme.jl rename to src/BoundaryConditions/perturbation_advection.jl From 3eb7fa3f15f10c01cf730fcf4dcd4a572203d6e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1s=20Chor?= Date: Wed, 27 Aug 2025 07:18:16 -0700 Subject: [PATCH 16/18] Update src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl Co-authored-by: Simone Silvestri --- src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl b/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl index 9715f803be6..1914bf07a08 100644 --- a/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl +++ b/src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl @@ -79,7 +79,7 @@ initialize_boundary_mass_flux(velocity, ::Nothing, side) = NamedTuple() initialize_boundary_mass_flux(velocity, bc, side) = NamedTuple() needs_mass_flux_correction(::IOBC) = false -needs_mass_flux_correction(bc::OBC) = bc.classification.scheme !== nothing +needs_mass_flux_correction(::OBC) = true needs_mass_flux_correction(::Nothing) = false needs_mass_flux_correction(bc) = false From 60a6f8b2d613dafabf22b7cfafb5a7b0016ed1be Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Mon, 1 Sep 2025 16:28:35 -0800 Subject: [PATCH 17/18] Update cubed_sphere_grid.jl --- src/MultiRegion/cubed_sphere_grid.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MultiRegion/cubed_sphere_grid.jl b/src/MultiRegion/cubed_sphere_grid.jl index c67cf1be7cf..66dbf4ab9d6 100644 --- a/src/MultiRegion/cubed_sphere_grid.jl +++ b/src/MultiRegion/cubed_sphere_grid.jl @@ -147,8 +147,8 @@ julia> using Oceananigans julia> using Oceananigans.MultiRegion: ConformalCubedSphereGrid julia> grid = ConformalCubedSphereGrid(panel_size=(12, 12, 1), z=(-1, 0), radius=1) -ConformalCubedSphereGrid{Float64, Oceananigans.Grids.FullyConnected, Oceananigans.Grids.FullyConnected, Bounded} partitioned on CPU(): -├── grids: 12×12×1 OrthogonalSphericalShellGrid{Float64, Oceananigans.Grids.FullyConnected, Oceananigans.Grids.FullyConnected, Bounded} on CPU with 3×3×3 halo and with precomputed metrics +ConformalCubedSphereGrid{Float64, FullyConnected, FullyConnected, Bounded} partitioned on CPU(): +├── grids: 12×12×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 3×3×3 halo and with precomputed metrics ├── partitioning: CubedSpherePartition with (1 region in each panel) ├── connectivity: CubedSphereConnectivity └── devices: (CPU(), CPU(), CPU(), CPU(), CPU(), CPU()) From 4e0cbef75ca848eef5b63f7218e161d110fd4177 Mon Sep 17 00:00:00 2001 From: tomaschor Date: Mon, 1 Sep 2025 20:21:14 -0700 Subject: [PATCH 18/18] fix doctest --- src/MultiRegion/cubed_sphere_grid.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MultiRegion/cubed_sphere_grid.jl b/src/MultiRegion/cubed_sphere_grid.jl index 66dbf4ab9d6..a77105c811b 100644 --- a/src/MultiRegion/cubed_sphere_grid.jl +++ b/src/MultiRegion/cubed_sphere_grid.jl @@ -148,8 +148,8 @@ julia> using Oceananigans.MultiRegion: ConformalCubedSphereGrid julia> grid = ConformalCubedSphereGrid(panel_size=(12, 12, 1), z=(-1, 0), radius=1) ConformalCubedSphereGrid{Float64, FullyConnected, FullyConnected, Bounded} partitioned on CPU(): -├── grids: 12×12×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 3×3×3 halo and with precomputed metrics -├── partitioning: CubedSpherePartition with (1 region in each panel) +├── region_grids: 12×12×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 3×3×3 halo and with precomputed metrics +├── partition: CubedSpherePartition with (1 region in each panel) ├── connectivity: CubedSphereConnectivity └── devices: (CPU(), CPU(), CPU(), CPU(), CPU(), CPU()) ```