From d992d81db35d42ce18af87f52d4596b4df3cefd3 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 18 Nov 2024 17:32:08 +0000 Subject: [PATCH 01/24] fix FEOBC --- ...at_extrapolation_open_boundary_matching_scheme.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl index b91ac36249..9650ddb100 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl @@ -61,7 +61,7 @@ end const c = Center() -@inline function _fill_west_open_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) +@inline function _fill_west_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) Δx₁ = xspacing(1, j, k, grid, c, c, c) Δx₂ = xspacing(2, j, k, grid, c, c, c) Δx₃ = xspacing(3, j, k, grid, c, c, c) @@ -75,7 +75,7 @@ const c = Center() return nothing end -@inline function _fill_east_open_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) +@inline function _fill_east_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) i = grid.Nx + 1 Δx₁ = xspacing(i-1, j, k, grid, c, c, c) @@ -91,7 +91,7 @@ end return nothing end -@inline function _fill_south_open_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) +@inline function _fill_south_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) Δy₁ = yspacing(i, 1, k, grid, c, c, c) Δy₂ = yspacing(i, 2, k, grid, c, c, c) Δy₃ = yspacing(i, 3, k, grid, c, c, c) @@ -105,7 +105,7 @@ end return nothing end -@inline function _fill_north_open_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) +@inline function _fill_north_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) j = grid.Ny + 1 Δy₁ = yspacing(i, j-1, k, grid, c, c, c) @@ -121,7 +121,7 @@ end return nothing end -@inline function _fill_bottom_open_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields) +@inline function _fill_bottom_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields) Δz₁ = zspacing(i, j, 1, grid, c, c, c) Δz₂ = zspacing(i, j, 2, grid, c, c, c) Δz₃ = zspacing(i, j, 3, grid, c, c, c) @@ -135,7 +135,7 @@ end return nothing end -@inline function _fill_top_open_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields) +@inline function _fill_top_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields) k = grid.Nz + 1 Δz₁ = zspacing(i, j, k-1, grid, c, c, c) From 2a374cb23f71cbb02729f91fa4749e4b9b6c6d80 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 3 Dec 2024 23:18:21 -0500 Subject: [PATCH 02/24] x-boundary condition implemented --- src/BoundaryConditions/BoundaryConditions.jl | 1 + ...advection_open_boundary_matching_scheme.jl | 87 ++++++++++++ validation/open_boundaries/cylinder.jl | 132 +++--------------- 3 files changed, 105 insertions(+), 115 deletions(-) create mode 100644 src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl diff --git a/src/BoundaryConditions/BoundaryConditions.jl b/src/BoundaryConditions/BoundaryConditions.jl index 2011bd3df4..2751bb621b 100644 --- a/src/BoundaryConditions/BoundaryConditions.jl +++ b/src/BoundaryConditions/BoundaryConditions.jl @@ -39,4 +39,5 @@ include("apply_flux_bcs.jl") include("update_boundary_conditions.jl") include("flat_extrapolation_open_boundary_matching_scheme.jl") +include("perturbation_advection_open_boundary_matching_scheme.jl") end # module diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl new file mode 100644 index 0000000000..3d62658e53 --- /dev/null +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl @@ -0,0 +1,87 @@ +using Oceananigans.Grids: xspacing +# Immersed boundaries are defined later but we probably need todo this? +#using Oceananigans.ImmersedBoundaries: active_cell + +""" + PerturbationAdvection + +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: +... +see latex document for now + +TODO: check what the coriolis is doing, and check what happens if U is the mean velocity +""" +struct PerturbationAdvection{FT} + inflow_timescale :: FT + outflow_timescale :: FT +end + +Adapt.adapt_structure(to, pe::PerturbationAdvection) = + PerturbationAdvection(adapt(to, pe.outflow_timescale), + adapt(to, pe.inflow_timescale)) + +function PerturbationAdvectionOpenBoundaryCondition(val, FT = Float64; + outflow_timescale = Inf, + inflow_timescale = 300.0, kwargs...) + + 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 + +const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} + +@inline function _fill_east_halo!(j, k, grid, u, bc::PAOBC, loc::Tuple{Face, Any, Any}, clock, model_fields) + i = grid.Nx + 1 + + Δt = clock.last_stage_Δt + + Δt = ifelse(isinf(Δt), 0, Δt) + + Δx = xspacing(i, j, k, grid, loc...) + + ūⁿ⁺¹ = getbc(bc, j, k, grid, clock, model_fields) + + uᵢⁿ = @inbounds u[i, j, k] + u′ᵢ₋₁ⁿ⁺¹ = @inbounds u[i - 1, j, k] - ūⁿ⁺¹ + + U = max(0, min(1, Δt / Δx * ūⁿ⁺¹)) + + τ = ifelse(ūⁿ⁺¹ >= 0, + bc.classification.matching_scheme.outflow_timescale, + bc.classification.matching_scheme.inflow_timescale) + + uᵢⁿ⁺¹ = (uᵢⁿ + U * u′ᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * (Δt / τ + U)) / (1 + Δt / τ + U) + + @inbounds u[i, j, k] = uᵢⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) +end + +@inline function _fill_west_halo!(j, k, grid, u, bc::PAOBC, loc::Tuple{Face, Any, Any}, clock, model_fields) + Δt = clock.last_stage_Δt + + Δt = ifelse(isinf(Δt), 0, Δt) + + Δx = xspacing(1, j, k, grid, loc...) + + ūⁿ⁺¹ = getbc(bc, j, k, grid, clock, model_fields) + + uᵢⁿ = @inbounds u[2, j, k] + u′ᵢ₋₁ⁿ⁺¹ = @inbounds u[0, j, k] - ūⁿ⁺¹ + + U = min(0, max(-1, Δt / Δx * ūⁿ⁺¹)) + + τ = ifelse(ūⁿ⁺¹ <= 0, + bc.classification.matching_scheme.outflow_timescale, + bc.classification.matching_scheme.inflow_timescale) + + τ̃ = min(1, Δt / τ) + + u₁ⁿ⁺¹ = (uᵢⁿ - U * u′ᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * (τ̃ - U)) / (1 + τ̃ - U) + + @inbounds u[1, j, k] = u₁ⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) + @inbounds u[0, j, k] = u₁ⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) +end diff --git a/validation/open_boundaries/cylinder.jl b/validation/open_boundaries/cylinder.jl index 993a41c4c1..0fb08d756b 100644 --- a/validation/open_boundaries/cylinder.jl +++ b/validation/open_boundaries/cylinder.jl @@ -1,127 +1,29 @@ -# This validation script shows open boundaries working in a simple case where the -# flow remains largely unidirectional and so at one end we have no matching scheme -# but just prescribe the inflow. At the other end we then make no assumptions about -# the flow and use a very simple open boundary condition to permit information to -# exit the domain. If, for example, the flow at the prescribed boundary was reversed -# then the model would likely fail. +using Oceananigans +using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition -using Oceananigans, CairoMakie -using Oceananigans.BoundaryConditions: FlatExtrapolationOpenBoundaryCondition +U(y, z, t) = 1#(1 + tanh(3t))/2 +u_east = PerturbationAdvectionOpenBoundaryCondition(U, inflow_timescale = Inf, outflow_timescale = Inf) +u_west = PerturbationAdvectionOpenBoundaryCondition(U, inflow_timescale = 1.0, outflow_timescale = Inf)#OpenBoundaryCondition(U) -@kwdef struct Cylinder{FT} - D :: FT = 1.0 - x₀ :: FT = 0.0 - y₀ :: FT = 0.0 -end +u_bcs = FieldBoundaryConditions(east = u_east, west = u_west) +v_bcs = FieldBoundaryConditions(east = GradientBoundaryCondition(0), west = GradientBoundaryCondition(0)) +w_bcs = FieldBoundaryConditions(east = GradientBoundaryCondition(0), west = GradientBoundaryCondition(0)) -@inline (cylinder::Cylinder)(x, y) = ifelse((x - cylinder.x₀)^2 + (y - cylinder.y₀)^2 < (cylinder.D/2)^2, 1, 0) +grid = RectilinearGrid(topology = (Bounded, Periodic, Bounded), size = (32, 16, 1), x = (-1, 1), y = (-0.5, 0.5), z = (-1, 0)) -architecture = GPU() +circle(x, y) = ifelse((x^2 + y^2) < 0.1^2, 0, -1) -# model parameters -Re = 200 -U = 1 -D = 1. -resolution = D / 40 - -# add extra downstream distance to see if the solution near the cylinder changes -extra_downstream = 0 - -cylinder = Cylinder(; D) - -x = (-5, 5 + extra_downstream) .* D -y = (-5, 5) .* D - -Ny = Int(10 / resolution) -Nx = Ny + Int(extra_downstream / resolution) - -ν = U * D / Re - -closure = ScalarDiffusivity(;ν, κ = ν) - -grid = RectilinearGrid(architecture; topology = (Bounded, Periodic, Flat), size = (Nx, Ny), x, y) - -@inline u(y, t, U) = U * (1 + 0.01 * randn()) - -u_boundaries = FieldBoundaryConditions(east = FlatExtrapolationOpenBoundaryCondition(), - west = OpenBoundaryCondition(u, parameters = U)) - -v_boundaries = FieldBoundaryConditions(east = GradientBoundaryCondition(0), - west = GradientBoundaryCondition(0)) - -Δt = .3 * resolution / U - -u_forcing = Relaxation(; rate = 1 / (2 * Δt), mask = cylinder) -v_forcing = Relaxation(; rate = 1 / (2 * Δt), mask = cylinder) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(circle)) model = NonhydrostaticModel(; grid, - closure, - forcing = (u = u_forcing, v = v_forcing), - boundary_conditions = (u = u_boundaries, v = v_boundaries)) - -@info "Constructed model" - -# initial noise to induce turbulance faster -set!(model, u = U, v = (x, y) -> randn() * U * 0.01) - -@info "Set initial conditions" - -simulation = Simulation(model; Δt = Δt, stop_time = 300) - -wizard = TimeStepWizard(cfl = 0.3) - -simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(100)) - -progress(sim) = @info "$(time(sim)) with Δt = $(prettytime(sim.Δt)) in $(prettytime(sim.run_wall_time))" - -simulation.callbacks[:progress] = Callback(progress, IterationInterval(1000)) - -simulation.output_writers[:velocity] = JLD2OutputWriter(model, model.velocities, - overwrite_existing = true, - filename = "cylinder_$(extra_downstream)_Re_$Re.jld2", - schedule = TimeInterval(1), - with_halos = true) - -run!(simulation) - -# load the results - -u_ts = FieldTimeSeries("cylinder_$(extra_downstream)_Re_$Re.jld2", "u") -v_ts = FieldTimeSeries("cylinder_$(extra_downstream)_Re_$Re.jld2", "v") - -u′, v′, w′ = Oceananigans.Fields.VelocityFields(u_ts.grid) - -ζ = Field((@at (Center, Center, Center) ∂x(v′)) - (@at (Center, Center, Center) ∂y(u′))) - -# there is probably a more memory efficient way todo this - -ζ_ts = zeros(size(grid, 1), size(grid, 2), length(u_ts.times)) # u_ts.grid so its always on cpu - -for n in 1:length(u_ts.times) - set!(u′, u_ts[n]) - set!(v′, v_ts[n]) - compute!(ζ) - ζ_ts[:, :, n] = interior(ζ, :, :, 1) -end - -@info "Loaded results" - -# plot the results - -fig = Figure(size = (600, 600)) - -ax = Axis(fig[1, 1], aspect = DataAspect()) - -xc, yc, zc = nodes(ζ) + boundary_conditions = (u = u_bcs, v = v_bcs, w = w_bcs), + advection = WENO(order = 5)) -n = Observable(1) +Δt = 0.3 * minimum_xspacing(grid) / 2 -ζ_plt = @lift ζ_ts[:, :, $n] +simulation = Simulation(model; Δt, stop_time = 20) -contour!(ax, xc, yc, ζ_plt, levels = [-2, 2], colorrange = (-2, 2), colormap = :roma) +simulation.output_writers[:fields] = JLD2OutputWriter(model, model.velocities; filename = "cylinder_$(grid.underlying_grid.Nz).jld2", schedule = TimeInterval(0.1), overwrite_existing = true) -record(fig, "ζ_Re_$Re.mp4", 1:length(u_ts.times), framerate = 5) do i; - n[] = i - i % 10 == 0 && @info "$(n.val) of $(length(u_ts.times))" -end \ No newline at end of file +run!(simulation) \ No newline at end of file From babe922fac47e20290b48abc0c210a726039b35f Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Wed, 4 Dec 2024 19:11:25 -0500 Subject: [PATCH 03/24] changed relaxation a bit --- ...perturbation_advection_open_boundary_matching_scheme.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl index 3d62658e53..5ad912c1e4 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl @@ -55,7 +55,10 @@ const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} bc.classification.matching_scheme.outflow_timescale, bc.classification.matching_scheme.inflow_timescale) - uᵢⁿ⁺¹ = (uᵢⁿ + U * u′ᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * (Δt / τ + U)) / (1 + Δt / τ + U) + + τ̃ = Δt / τ + + uᵢⁿ⁺¹ = (uᵢⁿ + U * u′ᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * (τ̃ + U)) / (1 + τ̃ + U) @inbounds u[i, j, k] = uᵢⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) end @@ -78,7 +81,7 @@ end bc.classification.matching_scheme.outflow_timescale, bc.classification.matching_scheme.inflow_timescale) - τ̃ = min(1, Δt / τ) + τ̃ = Δt / τ u₁ⁿ⁺¹ = (uᵢⁿ - U * u′ᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * (τ̃ - U)) / (1 + τ̃ - U) From 32cd354fdf6a1e8bcc27e774f7ef323f18c23703 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 6 Dec 2024 14:02:25 -0500 Subject: [PATCH 04/24] cleaning up + validation case --- ...advection_open_boundary_matching_scheme.jl | 14 +- validation/open_boundaries/cylinder.jl | 217 ++++++++++++++++-- 2 files changed, 208 insertions(+), 23 deletions(-) diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl index 5ad912c1e4..3cde0615af 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl @@ -1,5 +1,5 @@ using Oceananigans.Grids: xspacing -# Immersed boundaries are defined later but we probably need todo this? +# Immersed boundaries are defined later but we probably need todo this for when a boundary intersects the bathymetry #using Oceananigans.ImmersedBoundaries: active_cell """ @@ -46,8 +46,8 @@ const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} ūⁿ⁺¹ = getbc(bc, j, k, grid, clock, model_fields) - uᵢⁿ = @inbounds u[i, j, k] - u′ᵢ₋₁ⁿ⁺¹ = @inbounds u[i - 1, j, k] - ūⁿ⁺¹ + uᵢⁿ = @inbounds u[i, j, k] + uᵢ₋₁ⁿ⁺¹ = @inbounds u[i - 1, j, k] U = max(0, min(1, Δt / Δx * ūⁿ⁺¹)) @@ -58,7 +58,7 @@ const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} τ̃ = Δt / τ - uᵢⁿ⁺¹ = (uᵢⁿ + U * u′ᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * (τ̃ + U)) / (1 + τ̃ + U) + uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + U) @inbounds u[i, j, k] = uᵢⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) end @@ -72,8 +72,8 @@ end ūⁿ⁺¹ = getbc(bc, j, k, grid, clock, model_fields) - uᵢⁿ = @inbounds u[2, j, k] - u′ᵢ₋₁ⁿ⁺¹ = @inbounds u[0, j, k] - ūⁿ⁺¹ + uᵢⁿ = @inbounds u[2, j, k] + uᵢ₋₁ⁿ⁺¹ = @inbounds u[0, j, k] U = min(0, max(-1, Δt / Δx * ūⁿ⁺¹)) @@ -83,7 +83,7 @@ end τ̃ = Δt / τ - u₁ⁿ⁺¹ = (uᵢⁿ - U * u′ᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * (τ̃ - U)) / (1 + τ̃ - U) + u₁ⁿ⁺¹ = (uᵢⁿ - U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ - U) @inbounds u[1, j, k] = u₁ⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) @inbounds u[0, j, k] = u₁ⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) diff --git a/validation/open_boundaries/cylinder.jl b/validation/open_boundaries/cylinder.jl index 0fb08d756b..a058320577 100644 --- a/validation/open_boundaries/cylinder.jl +++ b/validation/open_boundaries/cylinder.jl @@ -1,29 +1,214 @@ using Oceananigans +using Oceananigans.Models.NonhydrostaticModels: ConjugateGradientPoissonSolver +using Oceananigans.Solvers: DiagonallyDominantPreconditioner +using Oceananigans.Operators: ℑxyᶠᶜᵃ, ℑxyᶜᶠᵃ +using Oceananigans.Solvers: FFTBasedPoissonSolver +using Printf +using CUDA + using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition -U(y, z, t) = 1#(1 + tanh(3t))/2 +u∞ = 1 +r = 1/2 +arch = GPU() +stop_time = 200 + +cylinder(x, y) = (x^2 + y^2) ≤ r^2 + +""" + drag(modell; bounding_box = (-1, 3, -2, 2), ν = 1e-3) + +Returns the drag within the `bounding_box` computed by: + +```math +\\frac{\\partial \\vec{u}}{\\partial t} + (\\vec{u}\\cdot\\nabla)\\vec{u}=-\\nabla P + \\nabla\\cdot\\vec{\\tau} + \\vec{F},\\newline +\\vec{F}_T=\\int_\\Omega\\vec{F}dV = \\int_\\Omega\\left(\\frac{\\partial \\vec{u}}{\\partial t} + (\\vec{u}\\cdot\\nabla)\\vec{u}+\\nabla P - \\nabla\\cdot\\vec{\\tau}\\right)dV,\\newline +\\vec{F}_T=\\int_\\Omega\\left(\\frac{\\partial \\vec{u}}{\\partial t}\\right)dV + \\oint_{\\partial\\Omega}\\left(\\vec{u}(\\vec{u}\\cdot\\hat{n}) + P\\hat{n} - \\vec{\\tau}\\cdot\\hat{n}\\right)dS,\\newline +\\F_u=\\int_\\Omega\\left(\\frac{\\partial u}{\\partial t}\\right)dV + \\oint_{\\partial\\Omega}\\left(u(\\vec{u}\\cdot\\hat{n}) - \\tau_{xx}\\right)dS + \\int_{\\partial\\Omega}P\\hat{x}\\cdot d\\vec{S},\\newline +F_u=\\int_\\Omega\\left(\\frac{\\partial u}{\\partial t}\\right)dV - \\int_{\\partial\\Omega_1} \\left(u^2 - 2\\nu\\frac{\\partial u}{\\partial x} + P\\right)dS + \\int_{\\partial\\Omega_2}\\left(u^2 - 2\\nu\\frac{\\partial u}{\\partial x}+P\\right)dS - \\int_{\\partial\\Omega_2} uvdS + \\int_{\\partial\\Omega_4} uvdS, +``` +where the bounding box is ``\\Omega`` which is formed from the boundaries ``\\partial\\Omega_{1}``, ``\\partial\\Omega_{2}``, ``\\partial\\Omega_{3}``, and ``\\partial\\Omega_{4}`` +which have outward directed normals ``-\\hat{x}``, ``\\hat{x}``, ``-\\hat{y}``, and ``\\hat{y}`` +""" +function drag(model; + bounding_box = (-1, 3, -2, 2), + ν = 1e-3) + + u, v, _ = model.velocities + + uᶜ = Field(@at (Center, Center, Center) u) + vᶜ = Field(@at (Center, Center, Center) v) + + xc, yc, _ = nodes(uᶜ) + + i₁ = findfirst(xc .> bounding_box[1]) + i₂ = findlast(xc .< bounding_box[2]) + + j₁ = findfirst(yc .> bounding_box[3]) + j₂ = findlast(yc .< bounding_box[4]) + + uₗ² = Field(uᶜ^2, indices = (i₁, j₁:j₂, 1)) + uᵣ² = Field(uᶜ^2, indices = (i₂, j₁:j₂, 1)) + + uvₗ = Field(uᶜ*vᶜ, indices = (i₁:i₂, j₁, 1)) + uvᵣ = Field(uᶜ*vᶜ, indices = (i₁:i₂, j₂, 1)) + + ∂₁uₗ = Field(∂x(uᶜ), indices = (i₁, j₁:j₂, 1)) + ∂₁uᵣ = Field(∂x(uᶜ), indices = (i₂, j₁:j₂, 1)) + + ∂ₜuᶜ = Field(@at (Center, Center, Center) model.timestepper.Gⁿ.u) + + ∂ₜu = Field(∂ₜuᶜ, indices = (i₁:i₂, j₁:j₂, 1)) + + p = model.pressures.pNHS + + ∫∂ₓp = Field(∂x(p), indices = (i₁:i₂, j₁:j₂, 1)) + + a_local = Field(Integral(∂ₜu)) + + a_flux = Field(Integral(uᵣ²)) - Field(Integral(uₗ²)) + Field(Integral(uvᵣ)) - Field(Integral(uvₗ)) + + a_viscous_stress = 2ν * (Field(Integral(∂₁uᵣ)) - Field(Integral(∂₁uₗ))) + + a_pressure = Field(Integral(∫∂ₓp)) + + return a_local + a_flux + a_pressure - a_viscous_stress +end + +Re = 1000 +#= +if Re <= 100 + Ny = 512 + Nx = Ny +elseif Re <= 1000 + Ny = 2^11 + Nx = Ny +elseif Re == 10^4 + Ny = 2^12 + Nx = Ny +elseif Re == 10^5 + Ny = 2^13 + Nx = Ny +elseif Re == 10^6 + Ny = 3/2 * 2^13 |> Int + Nx = Ny +end +=# + +Ny = 2048 +Nx = Ny + +prefix = "flow_around_cylinder_Re$(Re)_Ny$(Ny)" + +ϵ = 0 # break up-down symmetry +x = (-6, 12) # 18 +y = (-6 + ϵ, 6 + ϵ) # 12 +# TODO: temporatily use an iteration interval thingy with a fixed timestep!!! +kw = (; size=(Nx, Ny), x, y, halo=(6, 6), topology=(Bounded, Bounded, Flat)) +grid = RectilinearGrid(arch; kw...) +reduced_precision_grid = RectilinearGrid(arch, Float32; kw...) + +grid = ImmersedBoundaryGrid(grid, GridFittedBoundary(cylinder)) + +advection = Centered(order=2) +closure = ScalarDiffusivity(ν=1/Re) + +no_slip = ValueBoundaryCondition(0) +u_bcs = FieldBoundaryConditions(immersed=no_slip, + east = PerturbationAdvectionOpenBoundaryCondition(u∞; inflow_timescale = 1/4, outflow_timescale = Inf), + west = PerturbationAdvectionOpenBoundaryCondition(u∞; inflow_timescale = 0.1, outflow_timescale = Inf)) +v_bcs = FieldBoundaryConditions(immersed=no_slip, + east = GradientBoundaryCondition(0), + west = ValueBoundaryCondition(0)) +boundary_conditions = (u=u_bcs, v=v_bcs) + +ddp = DiagonallyDominantPreconditioner() +preconditioner = FFTBasedPoissonSolver(reduced_precision_grid) +reltol = abstol = 1e-7 +pressure_solver = ConjugateGradientPoissonSolver(grid, maxiter=10; + reltol, abstol, preconditioner) + +model = NonhydrostaticModel(; grid, pressure_solver, closure, + advection, boundary_conditions) + +@show model + +uᵢ(x, y) = 1e-2 * randn() +vᵢ(x, y) = 1e-2 * randn() +set!(model, u=uᵢ, v=vᵢ) + +Δx = minimum_xspacing(grid) +Δt = max_Δt = 0.002#0.2 * Δx^2 * Re + +simulation = Simulation(model; Δt, stop_time) +#conjure_time_step_wizard!(simulation, cfl=1.0, IterationInterval(3); max_Δt) + +u, v, w = model.velocities +#d = ∂x(u) + ∂y(v) + +# Drag computation +drag_force = drag(model; ν=1/Re) +compute!(drag_force) + +wall_time = Ref(time_ns()) + +function progress(sim) + if pressure_solver isa ConjugateGradientPoissonSolver + pressure_iters = iteration(pressure_solver) + else + pressure_iters = 0 + end + + #compute!(drag_force) + D = CUDA.@allowscalar drag_force[1, 1, 1] + cᴰ = D / (u∞ * r) + vmax = maximum(model.velocities.v) + dmax = 0#maximum(abs, d) + + msg = @sprintf("Iter: %d, time: %.2f, Δt: %.4f, Poisson iters: %d", + iteration(sim), time(sim), sim.Δt, pressure_iters) + + elapsed = 1e-9 * (time_ns() - wall_time[]) + + msg *= @sprintf(", max d: %.2e, max v: %.2e, Cd: %0.2f, wall time: %s", + dmax, vmax, cᴰ, prettytime(elapsed)) + + @info msg + wall_time[] = time_ns() + + return nothing +end -u_east = PerturbationAdvectionOpenBoundaryCondition(U, inflow_timescale = Inf, outflow_timescale = Inf) -u_west = PerturbationAdvectionOpenBoundaryCondition(U, inflow_timescale = 1.0, outflow_timescale = Inf)#OpenBoundaryCondition(U) +add_callback!(simulation, progress, IterationInterval(100)) -u_bcs = FieldBoundaryConditions(east = u_east, west = u_west) -v_bcs = FieldBoundaryConditions(east = GradientBoundaryCondition(0), west = GradientBoundaryCondition(0)) -w_bcs = FieldBoundaryConditions(east = GradientBoundaryCondition(0), west = GradientBoundaryCondition(0)) +ζ = ∂x(v) - ∂y(u) -grid = RectilinearGrid(topology = (Bounded, Periodic, Bounded), size = (32, 16, 1), x = (-1, 1), y = (-0.5, 0.5), z = (-1, 0)) +p = model.pressures.pNHS -circle(x, y) = ifelse((x^2 + y^2) < 0.1^2, 0, -1) +outputs = (; u, v, p, ζ) -grid = ImmersedBoundaryGrid(grid, GridFittedBottom(circle)) +simulation.output_writers[:jld2] = JLD2OutputWriter(model, outputs, + schedule = IterationInterval(Int(2/Δt)),#TimeInterval(0.1), + filename = prefix * "_fields.jld2", + overwrite_existing = true, + with_halos = true) -model = NonhydrostaticModel(; grid, - boundary_conditions = (u = u_bcs, v = v_bcs, w = w_bcs), - advection = WENO(order = 5)) +simulation.output_writers[:drag] = JLD2OutputWriter(model, (; drag_force), + schedule = IterationInterval(Int(0.1/Δt)),#TimeInterval(0.1), + filename = prefix * "_drag.jld2", + overwrite_existing = true, + with_halos = true, + indices = (1, 1, 1)) -Δt = 0.3 * minimum_xspacing(grid) / 2 +run!(simulation) -simulation = Simulation(model; Δt, stop_time = 20) -simulation.output_writers[:fields] = JLD2OutputWriter(model, model.velocities; filename = "cylinder_$(grid.underlying_grid.Nz).jld2", schedule = TimeInterval(0.1), overwrite_existing = true) +# Re = 100 +# with 12m ~0.33 Hz and Cd ~ 1.403 +# with 6m ~0.33 Hz and Cd ~ (higher, don't seem to have computed it!) +# with 18m ~0.32 Hz and Cd ~ 1.28 +# with 30m ~0.34 Hz and Cd ~ 1.37 +# the Strouhal number should be 0.16 to 0.18 which I think means we're pretty close as 1 drag osccilation cycle is 2 sheddings so we have 0.165 ish -run!(simulation) \ No newline at end of file +# Re = 1000 +# with 12m ~ HZ and Cd ~ \ No newline at end of file From 5a0285f20aeb4ad34a211889bf49466621dcf2c0 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 6 Dec 2024 14:18:01 -0500 Subject: [PATCH 05/24] changed to pointwise --- ...apolation_open_boundary_matching_scheme.jl | 38 +++++++++---------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl index 9650ddb100..382d5f685a 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl @@ -1,4 +1,4 @@ -using Oceananigans.Grids: xspacing, yspacing, zspacing +using Oceananigans.Grids: Δxᶜᶜᶜ, Δyᶜᶜᶜ, Δzᶜᶜᶜ """ FlatExtrapolation @@ -62,9 +62,9 @@ end const c = Center() @inline function _fill_west_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - Δx₁ = xspacing(1, j, k, grid, c, c, c) - Δx₂ = xspacing(2, j, k, grid, c, c, c) - Δx₃ = xspacing(3, j, k, grid, c, c, c) + Δx₁ = Δxᶜᶜᶜ(1, j, k, grid) + Δx₂ = Δxᶜᶜᶜ(2, j, k, grid) + Δx₃ = Δxᶜᶜᶜ(3, j, k, grid) spacing_factor = Δx₁ / (Δx₂ + Δx₃) @@ -78,9 +78,9 @@ end @inline function _fill_east_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) i = grid.Nx + 1 - Δx₁ = xspacing(i-1, j, k, grid, c, c, c) - Δx₂ = xspacing(i-2, j, k, grid, c, c, c) - Δx₃ = xspacing(i-3, j, k, grid, c, c, c) + Δ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₃) @@ -92,9 +92,9 @@ end end @inline function _fill_south_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - Δy₁ = yspacing(i, 1, k, grid, c, c, c) - Δy₂ = yspacing(i, 2, k, grid, c, c, c) - Δy₃ = yspacing(i, 3, k, grid, c, c, c) + Δy₁ = Δyᶜᶜᶜ(i, 1, k, grid) + Δy₂ = Δyᶜᶜᶜ(i, 2, k, grid) + Δy₃ = Δyᶜᶜᶜ(i, 3, k, grid) spacing_factor = Δy₁ / (Δy₂ + Δy₃) @@ -108,9 +108,9 @@ end @inline function _fill_north_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) j = grid.Ny + 1 - Δy₁ = yspacing(i, j-1, k, grid, c, c, c) - Δy₂ = yspacing(i, j-2, k, grid, c, c, c) - Δy₃ = yspacing(i, j-3, k, grid, c, c, c) + Δ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₃) @@ -122,9 +122,9 @@ end end @inline function _fill_bottom_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - Δz₁ = zspacing(i, j, 1, grid, c, c, c) - Δz₂ = zspacing(i, j, 2, grid, c, c, c) - Δz₃ = zspacing(i, j, 3, grid, c, c, c) + Δz₁ = Δzᶜᶜᶜ(i, j, 1, grid) + Δz₂ = Δzᶜᶜᶜ(i, j, 2, grid) + Δz₃ = Δzᶜᶜᶜ(i, j, 3, grid) spacing_factor = Δz₁ / (Δz₂ + Δz₃) @@ -138,9 +138,9 @@ end @inline function _fill_top_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields) k = grid.Nz + 1 - Δz₁ = zspacing(i, j, k-1, grid, c, c, c) - Δz₂ = zspacing(i, j, k-2, grid, c, c, c) - Δz₃ = zspacing(i, j, k-3, grid, c, c, c) + Δ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₃) From fef6185cd3fe7746968933e4b4074f61c025bdf2 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 6 Dec 2024 14:25:21 -0500 Subject: [PATCH 06/24] oops --- .../flat_extrapolation_open_boundary_matching_scheme.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl index 382d5f685a..11227437ae 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl @@ -1,4 +1,4 @@ -using Oceananigans.Grids: Δxᶜᶜᶜ, Δyᶜᶜᶜ, Δzᶜᶜᶜ +using Oceananigans.Operators: Δxᶜᶜᶜ, Δyᶜᶜᶜ, Δzᶜᶜᶜ """ FlatExtrapolation From fa5a471dc0930758ad1d26572bdabb7fa9bc8eed Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 6 Dec 2024 15:05:58 -0500 Subject: [PATCH 07/24] fixed open boundary filling + tests --- src/BoundaryConditions/fill_halo_regions.jl | 4 +- .../fill_halo_regions_open.jl | 4 +- ...apolation_open_boundary_matching_scheme.jl | 2 - test/test_boundary_conditions_integration.jl | 67 ++++++++++++++++++- 4 files changed, 72 insertions(+), 5 deletions(-) diff --git a/src/BoundaryConditions/fill_halo_regions.jl b/src/BoundaryConditions/fill_halo_regions.jl index 7ac3e79ea5..5eb29f687e 100644 --- a/src/BoundaryConditions/fill_halo_regions.jl +++ b/src/BoundaryConditions/fill_halo_regions.jl @@ -20,12 +20,14 @@ conditions, possibly recursing into `fields` if it is a nested tuple-of-tuples. # Some fields have `nothing` boundary conditions, such as `FunctionField` and `ZeroField`. fill_halo_regions!(c::OffsetArray, ::Nothing, args...; kwargs...) = nothing +retrieve_bc(bc) = bc + # Returns the boundary conditions a specific side for `FieldBoundaryConditions` inputs and # a tuple of boundary conditions for `NTuple{N, <:FieldBoundaryConditions}` inputs for dir in (:west, :east, :south, :north, :bottom, :top) extract_side_bc = Symbol(:extract_, dir, :_bc) @eval begin - @inline $extract_side_bc(bc) = bc.$dir + @inline $extract_side_bc(bc) = retrieve_bc(bc.$dir) @inline $extract_side_bc(bc::Tuple) = map($extract_side_bc, bc) end end diff --git a/src/BoundaryConditions/fill_halo_regions_open.jl b/src/BoundaryConditions/fill_halo_regions_open.jl index 039de899dd..e19bd3245a 100644 --- a/src/BoundaryConditions/fill_halo_regions_open.jl +++ b/src/BoundaryConditions/fill_halo_regions_open.jl @@ -48,6 +48,8 @@ end @inline retrieve_open_bc(bc::OBC) = bc @inline retrieve_open_bc(bc) = nothing +@inline retrieve_bc(bc::OBC) = nothing + # for regular halo fills, return nothing if the BC is not an OBC @inline left_open_boundary_condition(boundary_condition, loc) = nothing @inline left_open_boundary_condition(boundary_conditions, ::Tuple{Face, Center, Center}) = retrieve_open_bc(boundary_conditions.west) @@ -59,7 +61,7 @@ end @inline right_open_boundary_condition(boundary_conditions, ::Tuple{Center, Face, Center}) = retrieve_open_bc(boundary_conditions.north) @inline right_open_boundary_condition(boundary_conditions, ::Tuple{Center, Center, Face}) = retrieve_open_bc(boundary_conditions.top) -# Opern boundary fill +# Open boundary fill @inline _fill_west_halo!(j, k, grid, c, bc::OBC, loc, args...) = @inbounds c[1, j, k] = getbc(bc, j, k, grid, args...) @inline _fill_east_halo!(j, k, grid, c, bc::OBC, loc, args...) = @inbounds c[grid.Nx + 1, j, k] = getbc(bc, j, k, grid, args...) @inline _fill_south_halo!(i, k, grid, c, bc::OBC, loc, args...) = @inbounds c[i, 1, k] = getbc(bc, i, k, grid, args...) diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl index 11227437ae..1325b26173 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl @@ -59,8 +59,6 @@ end return ϕ + Δϕ end -const c = Center() - @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) diff --git a/test/test_boundary_conditions_integration.jl b/test/test_boundary_conditions_integration.jl index 16efa34729..43952400ce 100644 --- a/test/test_boundary_conditions_integration.jl +++ b/test/test_boundary_conditions_integration.jl @@ -1,6 +1,6 @@ include("dependencies_for_runtests.jl") -using Oceananigans.BoundaryConditions: ContinuousBoundaryFunction +using Oceananigans.BoundaryConditions: ContinuousBoundaryFunction, FlatExtrapolationOpenBoundaryCondition, fill_halo_regions! using Oceananigans: prognostic_fields function test_boundary_condition(arch, FT, topo, side, field_name, boundary_condition) @@ -103,6 +103,63 @@ 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 = FlatExtrapolationOpenBoundaryCondition()) + +right_febc(::Val{1}, grid, loc) = FieldBoundaryConditions(grid, loc, west = OpenBoundaryCondition(1), + east = FlatExtrapolationOpenBoundaryCondition()) + +left_febc(::Val{2}, grid, loc) = FieldBoundaryConditions(grid, loc, north = OpenBoundaryCondition(1), + south = FlatExtrapolationOpenBoundaryCondition()) + +right_febc(::Val{2}, grid, loc) = FieldBoundaryConditions(grid, loc, south = OpenBoundaryCondition(1), + north = FlatExtrapolationOpenBoundaryCondition()) + +left_febc(::Val{3}, grid, loc) = FieldBoundaryConditions(grid, loc, top = OpenBoundaryCondition(1), + bottom = FlatExtrapolationOpenBoundaryCondition()) + +right_febc(::Val{3}, grid, loc) = FieldBoundaryConditions(grid, loc, bottom = OpenBoundaryCondition(1), + top = FlatExtrapolationOpenBoundaryCondition()) + +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 flat_extrapolation_open_boundary_conditions_are_correct(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 + @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 + test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType), float_bc(C, FT, ArrayType), irrational_bc(C, FT, ArrayType), @@ -273,4 +330,12 @@ test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType), @test fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT) end end + + @testset "Open boundary conditions" begin + 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_are_correct(arch, FT) + end + end end From b55395c168e203567a92472c99c10f62462c726b Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 6 Dec 2024 15:33:23 -0500 Subject: [PATCH 08/24] generalised and fixed bug in left boundary --- ...advection_open_boundary_matching_scheme.jl | 75 ++++++++++++------- 1 file changed, 49 insertions(+), 26 deletions(-) diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl index 3cde0615af..e1efa4cf17 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl @@ -1,6 +1,4 @@ using Oceananigans.Grids: xspacing -# Immersed boundaries are defined later but we probably need todo this for when a boundary intersects the bathymetry -#using Oceananigans.ImmersedBoundaries: active_cell """ PerturbationAdvection @@ -35,56 +33,81 @@ end const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} -@inline function _fill_east_halo!(j, k, grid, u, bc::PAOBC, loc::Tuple{Face, Any, Any}, clock, model_fields) - i = grid.Nx + 1 - +@inline function step_right_boundary_node!(bc, l, m, boundary_indices, boundary_adjacent_indices, + grid, u, clock, model_fields, ΔX) Δt = clock.last_stage_Δt Δt = ifelse(isinf(Δt), 0, Δt) - Δx = xspacing(i, j, k, grid, loc...) + ūⁿ⁺¹ = getbc(bc, l, m, grid, clock, model_fields) - ūⁿ⁺¹ = getbc(bc, j, k, grid, clock, model_fields) + uᵢⁿ = @inbounds getindex(u, boundary_indices...) + uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, boundary_adjacent_indices...) - uᵢⁿ = @inbounds u[i, j, k] - uᵢ₋₁ⁿ⁺¹ = @inbounds u[i - 1, j, k] + U = max(0, min(1, Δt / ΔX * ūⁿ⁺¹)) - U = max(0, min(1, Δt / Δx * ūⁿ⁺¹)) - - τ = ifelse(ūⁿ⁺¹ >= 0, - bc.classification.matching_scheme.outflow_timescale, - bc.classification.matching_scheme.inflow_timescale) + pa = bc.classification.matching_scheme + τ = ifelse(ūⁿ⁺¹ >= 0, pa.outflow_timescale, pa.inflow_timescale) τ̃ = Δt / τ uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + U) - @inbounds u[i, j, k] = uᵢⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) + @inbounds setindex!(u, uᵢⁿ⁺¹, boundary_indices...) + + return nothing end -@inline function _fill_west_halo!(j, k, grid, u, bc::PAOBC, loc::Tuple{Face, Any, Any}, clock, model_fields) + +@inline function step_left_boundary_node!(bc, l, m, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, + grid, u, clock, model_fields, ΔX) Δt = clock.last_stage_Δt Δt = ifelse(isinf(Δt), 0, Δt) - Δx = xspacing(1, j, k, grid, loc...) + ūⁿ⁺¹ = getbc(bc, l, m, grid, clock, model_fields) - ūⁿ⁺¹ = getbc(bc, j, k, grid, clock, model_fields) + uᵢⁿ = @inbounds getindex(u, boundary_secret_storage_indices...) + uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, boundary_adjacent_indices...) - uᵢⁿ = @inbounds u[2, j, k] - uᵢ₋₁ⁿ⁺¹ = @inbounds u[0, j, k] + U = min(0, max(-1, Δt / ΔX * ūⁿ⁺¹)) - U = min(0, max(-1, Δt / Δx * ūⁿ⁺¹)) + pa = bc.classification.matching_scheme - τ = ifelse(ūⁿ⁺¹ <= 0, - bc.classification.matching_scheme.outflow_timescale, - bc.classification.matching_scheme.inflow_timescale) + τ = ifelse(ūⁿ⁺¹ <= 0, pa.outflow_timescale, pa.inflow_timescale) τ̃ = Δt / τ u₁ⁿ⁺¹ = (uᵢⁿ - U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ - U) - @inbounds u[1, j, k] = u₁ⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) - @inbounds u[0, j, k] = u₁ⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) + @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_indices...) + @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_secret_storage_indices...) + + return nothing +end + +@inline function _fill_east_halo!(j, k, grid, u, bc::PAOBC, loc::Tuple{Face, Any, Any}, clock, model_fields) + i = grid.Nx + 1 + + boundary_indices = (i, j, k) + boundary_adjacent_indices = (i-1, j, k) + + Δx = xspacing(i, j, k, grid, loc...) + + step_right_boundary_node!(bc, j, k, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, Δx) + + return nothing +end + +@inline function _fill_west_halo!(j, k, grid, u, bc::PAOBC, loc::Tuple{Face, Any, Any}, clock, model_fields) + boundary_indices = (1, j, k) + boundary_adjacent_indices = (2, j, k) + boundary_secret_storage_indices = (0, j, k) + + Δx = xspacing(1, j, k, grid, loc...) + + step_left_boundary_node!(bc, j, k, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, grid, u, clock, model_fields, Δx) + + return nothing end From 1a9fbecfe8e0fd927cf8f9e82d55676c0cdec028 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 6 Dec 2024 15:35:37 -0500 Subject: [PATCH 09/24] bump patch --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index cec2f85661..8451e115c0 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.94.3" +version = "0.94.4" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From e66f3d14aacc6c593d017df5be134571fce136c6 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 6 Dec 2024 15:34:33 -0500 Subject: [PATCH 10/24] Formating --- ...bation_advection_open_boundary_matching_scheme.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl index e1efa4cf17..a89796ca94 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl @@ -33,8 +33,8 @@ end const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} -@inline function step_right_boundary_node!(bc, l, m, boundary_indices, boundary_adjacent_indices, - grid, u, clock, model_fields, ΔX) +@inline function step_right_boundary!(bc, l, m, boundary_indices, boundary_adjacent_indices, + grid, u, clock, model_fields, ΔX) Δt = clock.last_stage_Δt Δt = ifelse(isinf(Δt), 0, Δt) @@ -60,8 +60,8 @@ const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} end -@inline function step_left_boundary_node!(bc, l, m, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, - grid, u, clock, model_fields, ΔX) +@inline function step_left_boundary!(bc, l, m, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, + grid, u, clock, model_fields, ΔX) Δt = clock.last_stage_Δt Δt = ifelse(isinf(Δt), 0, Δt) @@ -95,7 +95,7 @@ end Δx = xspacing(i, j, k, grid, loc...) - step_right_boundary_node!(bc, j, k, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, Δx) + step_right_boundary!(bc, j, k, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, Δx) return nothing end @@ -107,7 +107,7 @@ end Δx = xspacing(1, j, k, grid, loc...) - step_left_boundary_node!(bc, j, k, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, grid, u, clock, model_fields, Δx) + step_left_boundary!(bc, j, k, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, grid, u, clock, model_fields, Δx) return nothing end From 51c9938ec37033891aa4cdc81b4d80ef4492cfa7 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 10 Dec 2024 15:34:58 -0500 Subject: [PATCH 11/24] test problem --- test/test_boundary_conditions_integration.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_boundary_conditions_integration.jl b/test/test_boundary_conditions_integration.jl index 43952400ce..de169e0a2d 100644 --- a/test/test_boundary_conditions_integration.jl +++ b/test/test_boundary_conditions_integration.jl @@ -125,7 +125,7 @@ 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 flat_extrapolation_open_boundary_conditions_are_correct(arch, FT) +function test_flat_extrapolation_open_boundary_conditions(arch, FT) clock = Clock(; time = zero(FT)) for orientation in 1:3 @@ -335,7 +335,7 @@ 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_are_correct(arch, FT) + test_flat_extrapolation_open_boundary_conditions(arch, FT) end end end From 6ff1bc7375606f03d08c17c83f5684e8a9ec0a46 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 10 Dec 2024 17:06:54 -0500 Subject: [PATCH 12/24] bump patch --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8451e115c0..ec1a5f4bda 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.94.4" +version = "0.94.5" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From 62e27b2a4ede9aa9e1a023d65440a699b37eb3a0 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 12 Dec 2024 19:28:22 +0000 Subject: [PATCH 13/24] clean up --- ...rbation_advection_open_boundary_matching_scheme.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl index a89796ca94..0f12b30c50 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl @@ -1,4 +1,4 @@ -using Oceananigans.Grids: xspacing +using Oceananigans.Operators: Δxᶠᶜᶜ, Δyᶜᶠᶜ, Δzᶜᶜᶠ, Ax_qᶠᶜᶜ, Ay_qᶜᶠᶜ, Az_qᶜᶜᶠ """ PerturbationAdvection @@ -59,7 +59,6 @@ const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} return nothing end - @inline function step_left_boundary!(bc, l, m, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, grid, u, clock, model_fields, ΔX) Δt = clock.last_stage_Δt @@ -87,25 +86,25 @@ end return nothing end -@inline function _fill_east_halo!(j, k, grid, u, bc::PAOBC, loc::Tuple{Face, Any, Any}, clock, model_fields) +@inline function _fill_east_halo!(j, k, grid, u, bc::PAOBC, ::Tuple{Face, Any, Any}, clock, model_fields) i = grid.Nx + 1 boundary_indices = (i, j, k) boundary_adjacent_indices = (i-1, j, k) - Δx = xspacing(i, j, k, grid, loc...) + Δx = Δxᶠᶜᶜ(i, j, k, grid) step_right_boundary!(bc, j, k, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, Δx) return nothing end -@inline function _fill_west_halo!(j, k, grid, u, bc::PAOBC, loc::Tuple{Face, Any, Any}, clock, model_fields) +@inline function _fill_west_halo!(j, k, grid, u, bc::PAOBC, ::Tuple{Face, Any, Any}, clock, model_fields) boundary_indices = (1, j, k) boundary_adjacent_indices = (2, j, k) boundary_secret_storage_indices = (0, j, k) - Δx = xspacing(1, j, k, grid, loc...) + Δx = Δxᶠᶜᶜ(1, j, k, grid) step_left_boundary!(bc, j, k, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, grid, u, clock, model_fields, Δx) From a23adae025498a92b1fbfd46e1279d1beae4a1f0 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 12 Dec 2024 19:28:36 +0000 Subject: [PATCH 14/24] boundary normal velocity --- src/Oceananigans.jl | 2 ++ src/boundary_normal_mean_velocity.jl | 47 ++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+) create mode 100644 src/boundary_normal_mean_velocity.jl diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index d49accb100..c8561680b6 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -221,6 +221,8 @@ include("Simulations/Simulations.jl") # Abstractions for distributed and multi-region models include("MultiRegion/MultiRegion.jl") +include("boundary_normal_mean_velocity.jl") + ##### ##### Needed so we can export names from sub-modules at the top-level ##### diff --git a/src/boundary_normal_mean_velocity.jl b/src/boundary_normal_mean_velocity.jl new file mode 100644 index 0000000000..2a98162cbf --- /dev/null +++ b/src/boundary_normal_mean_velocity.jl @@ -0,0 +1,47 @@ +using Adapt +using Oceananigans.AbstractOperations: GridMetricOperation, Ax, Ay, Az +using Oceananigans.BoundaryConditions: BoundaryCondition, Open, PerturbationAdvection + +import Adapt: adapt_structure +import Oceananigans.BoundaryConditions: update_boundary_condition! + +struct BoundaryNormalMeanVelocity{BV} + boundary_velocity :: BV + + BoundaryNormalMeanVelocity(; value::BV = Ref(0.0)) where BV = new{BV}(value) +end + +@inline (value::BoundaryNormalMeanVelocity)(args...) = value.boundary_velocity[] + +Adapt.adapt_structure(to, mo::BoundaryNormalMeanVelocity) = + BoundaryNormalMeanVelocity(adapt(to, mo.mean_outflow_velocity[])) + +const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryNormalMeanVelocity} + +@inline boundary_normal_area(::Union{Val{:west}, Val{:east}}, grid) = GridMetricOperation((Face, Center, Center), Ax, grid) +@inline boundary_normal_area(::Union{Val{:south}, Val{:north}}, grid) = GridMetricOperation((Center, Face, Center), Ay, grid) +@inline boundary_normal_area(::Union{Val{:bottom}, Val{:top}}, grid) = GridMetricOperation((Center, Center, Face), Az, grid) + +@inline boundary_adjacent_index(::Val{:east}, grid) = (size(grid, 1), 1, 1), (2, 3) +@inline boundary_adjacent_index(::Val{:west}, grid) = (2, 1, 1), (2, 3) + +@inline boundary_adjacent_index(::Val{:north}, grid) = (1, size(grid, 2), 1), (1, 3) +@inline boundary_adjacent_index(::Val{:south}, grid) = (1, 2, 1), (1, 3) + +@inline boundary_adjacent_index(::Val{:top}, grid) = (1, 1, size(grid, 3)), (2, 3) +@inline boundary_adjacent_index(::Val{:bottom}, grid) = (1, 1, 2), (2, 3) + +function update_boundary_condition!(bc::MOPABC, side, u, model) + grid = model.grid + val_side = Val(side) + + An = boundary_normal_area(val_side, grid) + + (i, j, k), dims = boundary_adjacent_index(val_side, grid) + + total_area = sum(An; dims)[i, j, k] + + Ū = sum(u * An; dims) / total_area + + bc.condition.boundary_velocity[] = Ū[i, j, k] +end \ No newline at end of file From 5246c9aa953d44ce05fd8a209ee891927ddcbf44 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 12 Dec 2024 19:29:32 +0000 Subject: [PATCH 15/24] oops --- src/boundary_normal_mean_velocity.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/boundary_normal_mean_velocity.jl b/src/boundary_normal_mean_velocity.jl index 2a98162cbf..cf85ecb462 100644 --- a/src/boundary_normal_mean_velocity.jl +++ b/src/boundary_normal_mean_velocity.jl @@ -31,9 +31,8 @@ const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryNorm @inline boundary_adjacent_index(::Val{:top}, grid) = (1, 1, size(grid, 3)), (2, 3) @inline boundary_adjacent_index(::Val{:bottom}, grid) = (1, 1, 2), (2, 3) -function update_boundary_condition!(bc::MOPABC, side, u, model) +function update_boundary_condition!(bc::MOPABC, val_side, u, model) grid = model.grid - val_side = Val(side) An = boundary_normal_area(val_side, grid) From a34c927030851c844657bba259a8bd20e4acd4ff Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 12 Dec 2024 19:34:33 +0000 Subject: [PATCH 16/24] Rename --- src/boundary_normal_mean_velocity.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/boundary_normal_mean_velocity.jl b/src/boundary_normal_mean_velocity.jl index cf85ecb462..1f6795eb67 100644 --- a/src/boundary_normal_mean_velocity.jl +++ b/src/boundary_normal_mean_velocity.jl @@ -5,18 +5,18 @@ using Oceananigans.BoundaryConditions: BoundaryCondition, Open, PerturbationAdve import Adapt: adapt_structure import Oceananigans.BoundaryConditions: update_boundary_condition! -struct BoundaryNormalMeanVelocity{BV} - boundary_velocity :: BV +struct BoundaryMean{V} + value :: V - BoundaryNormalMeanVelocity(; value::BV = Ref(0.0)) where BV = new{BV}(value) + BoundaryMean(; value::BV = Ref(0.0)) where BV = new{BV}(value) end -@inline (value::BoundaryNormalMeanVelocity)(args...) = value.boundary_velocity[] +@inline (value::BoundaryMean)(args...) = value.boundary_velocity[] -Adapt.adapt_structure(to, mo::BoundaryNormalMeanVelocity) = - BoundaryNormalMeanVelocity(adapt(to, mo.mean_outflow_velocity[])) +Adapt.adapt_structure(to, mo::BoundaryMean) = + BoundaryMean(adapt(to, mo.mean_outflow_velocity[])) -const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryNormalMeanVelocity} +const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryMean} @inline boundary_normal_area(::Union{Val{:west}, Val{:east}}, grid) = GridMetricOperation((Face, Center, Center), Ax, grid) @inline boundary_normal_area(::Union{Val{:south}, Val{:north}}, grid) = GridMetricOperation((Center, Face, Center), Ay, grid) @@ -42,5 +42,5 @@ function update_boundary_condition!(bc::MOPABC, val_side, u, model) Ū = sum(u * An; dims) / total_area - bc.condition.boundary_velocity[] = Ū[i, j, k] + bc.condition.value[] = Ū[i, j, k] end \ No newline at end of file From 0df1abeb7628a2a6677f4a3b74bcaddc5d0d77ee Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 12 Dec 2024 19:35:11 +0000 Subject: [PATCH 17/24] finish renaming --- src/Oceananigans.jl | 508 +++++++++--------- ...rmal_mean_velocity.jl => boundary_mean.jl} | 0 2 files changed, 254 insertions(+), 254 deletions(-) rename src/{boundary_normal_mean_velocity.jl => boundary_mean.jl} (100%) diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index c8561680b6..661c9c0aef 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -1,254 +1,254 @@ -""" -Main module for `Oceananigans.jl` -- a Julia software for fast, friendly, flexible, -data-driven, ocean-flavored fluid dynamics on CPUs and GPUs. -""" -module Oceananigans - -export - # Architectures - CPU, GPU, - - # Grids - Center, Face, - Periodic, Bounded, Flat, - RectilinearGrid, LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, - nodes, xnodes, ynodes, znodes, λnodes, φnodes, - xspacings, yspacings, zspacings, λspacings, φspacings, - minimum_xspacing, minimum_yspacing, minimum_zspacing, - - # Pointwise spacing, area, and volume operators - xspacing, yspacing, zspacing, λspacing, φspacing, xarea, yarea, zarea, volume, - - # Immersed boundaries - ImmersedBoundaryGrid, - GridFittedBoundary, GridFittedBottom, PartialCellBottom, - ImmersedBoundaryCondition, - - # Distributed - Distributed, Partition, - - # Advection schemes - Centered, UpwindBiased, WENO, - VectorInvariant, WENOVectorInvariant, FluxFormAdvection, - - # Boundary conditions - BoundaryCondition, - FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, OpenBoundaryCondition, - FieldBoundaryConditions, - - # Fields and field manipulation - Field, CenterField, XFaceField, YFaceField, ZFaceField, - Average, Integral, CumulativeIntegral, Reduction, Accumulation, BackgroundField, - interior, set!, compute!, regrid!, - - # Forcing functions - Forcing, Relaxation, LinearTarget, GaussianMask, AdvectiveForcing, - - # Coriolis forces - FPlane, ConstantCartesianCoriolis, BetaPlane, NonTraditionalBetaPlane, - - # BuoyancyModels and equations of state - Buoyancy, BuoyancyTracer, SeawaterBuoyancy, - LinearEquationOfState, TEOS10, - BuoyancyField, - - # Surface wave Stokes drift via Craik-Leibovich equations - UniformStokesDrift, StokesDrift, - - # Turbulence closures - VerticalScalarDiffusivity, - HorizontalScalarDiffusivity, - ScalarDiffusivity, - VerticalScalarBiharmonicDiffusivity, - HorizontalScalarBiharmonicDiffusivity, - ScalarBiharmonicDiffusivity, - SmagorinskyLilly, - Smagorinsky, - LillyCoefficient, - DynamicCoefficient, - AnisotropicMinimumDissipation, - ConvectiveAdjustmentVerticalDiffusivity, - CATKEVerticalDiffusivity, - RiBasedVerticalDiffusivity, - VerticallyImplicitTimeDiscretization, - viscosity, diffusivity, - - # Lagrangian particle tracking - LagrangianParticles, - - # Models - NonhydrostaticModel, HydrostaticFreeSurfaceModel, ShallowWaterModel, - ConservativeFormulation, VectorInvariantFormulation, - PressureField, fields, - - # Hydrostatic free surface model stuff - VectorInvariant, ExplicitFreeSurface, ImplicitFreeSurface, SplitExplicitFreeSurface, - HydrostaticSphericalCoriolis, PrescribedVelocityFields, - - # Time stepping - Clock, TimeStepWizard, conjure_time_step_wizard!, time_step!, - - # Simulations - Simulation, run!, Callback, add_callback!, iteration, - iteration_limit_exceeded, stop_time_exceeded, wall_time_limit_exceeded, - - # Diagnostics - CFL, AdvectiveCFL, DiffusiveCFL, - - # Output writers - NetCDFOutputWriter, JLD2OutputWriter, Checkpointer, - TimeInterval, IterationInterval, WallTimeInterval, AveragedTimeInterval, - SpecifiedTimes, FileSizeLimit, AndSchedule, OrSchedule, written_names, - - # Output readers - FieldTimeSeries, FieldDataset, InMemory, OnDisk, - - # Abstract operations - ∂x, ∂y, ∂z, @at, KernelFunctionOperation, - - # Utils - prettytime - -using CUDA -using DocStringExtensions -using FFTW - -function __init__() - threads = Threads.nthreads() - if threads > 1 - @info "Oceananigans will use $threads threads" - - # See: https://github.com/CliMA/Oceananigans.jl/issues/1113 - FFTW.set_num_threads(4threads) - end - - if CUDA.has_cuda() - @debug "CUDA-enabled GPU(s) detected:" - for (gpu, dev) in enumerate(CUDA.devices()) - @debug "$dev: $(CUDA.name(dev))" - end - - CUDA.allowscalar(false) - end -end - -##### -##### Abstract types -##### - -""" - AbstractModel - -Abstract supertype for models. -""" -abstract type AbstractModel{TS} end - -""" - AbstractDiagnostic - -Abstract supertype for diagnostics that compute information from the current -model state. -""" -abstract type AbstractDiagnostic end - -""" - AbstractOutputWriter - -Abstract supertype for output writers that write data to disk. -""" -abstract type AbstractOutputWriter end - -# Callsites for Callbacks -struct TimeStepCallsite end -struct TendencyCallsite end -struct UpdateStateCallsite end - -##### -##### Place-holder functions -##### - -function run_diagnostic! end -function write_output! end -function initialize! end # for initializing models, simulations, etc -function location end -function instantiated_location end -function tupleit end -function fields end -function prognostic_fields end -function tracer_tendency_kernel_function end -function boundary_conditions end - -##### -##### Include all the submodules -##### - -# Basics -include("Architectures.jl") -include("Units.jl") -include("Grids/Grids.jl") -include("Utils/Utils.jl") -include("Logger.jl") -include("Operators/Operators.jl") -include("BoundaryConditions/BoundaryConditions.jl") -include("Fields/Fields.jl") -include("AbstractOperations/AbstractOperations.jl") -include("TimeSteppers/TimeSteppers.jl") -include("ImmersedBoundaries/ImmersedBoundaries.jl") -include("Advection/Advection.jl") -include("Solvers/Solvers.jl") -include("OutputReaders/OutputReaders.jl") -include("DistributedComputations/DistributedComputations.jl") - -# TODO: move here -#include("MultiRegion/MultiRegion.jl") - -# Physics, time-stepping, and models -include("Coriolis/Coriolis.jl") -include("BuoyancyModels/BuoyancyModels.jl") -include("StokesDrifts.jl") -include("TurbulenceClosures/TurbulenceClosures.jl") -include("Forcings/Forcings.jl") -include("Biogeochemistry.jl") - -# TODO: move above -include("Models/Models.jl") - -# Output and Physics, time-stepping, and models -include("Diagnostics/Diagnostics.jl") -include("OutputWriters/OutputWriters.jl") -include("Simulations/Simulations.jl") - -# Abstractions for distributed and multi-region models -include("MultiRegion/MultiRegion.jl") - -include("boundary_normal_mean_velocity.jl") - -##### -##### Needed so we can export names from sub-modules at the top-level -##### - -using .Logger -using .Architectures -using .Utils -using .Advection -using .Grids -using .BoundaryConditions -using .Fields -using .Coriolis -using .BuoyancyModels -using .StokesDrifts -using .TurbulenceClosures -using .Solvers -using .OutputReaders -using .Forcings -using .ImmersedBoundaries -using .DistributedComputations -using .Models -using .TimeSteppers -using .Diagnostics -using .OutputWriters -using .Simulations -using .AbstractOperations -using .MultiRegion - -end # module +""" +Main module for `Oceananigans.jl` -- a Julia software for fast, friendly, flexible, +data-driven, ocean-flavored fluid dynamics on CPUs and GPUs. +""" +module Oceananigans + +export + # Architectures + CPU, GPU, + + # Grids + Center, Face, + Periodic, Bounded, Flat, + RectilinearGrid, LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, + nodes, xnodes, ynodes, znodes, λnodes, φnodes, + xspacings, yspacings, zspacings, λspacings, φspacings, + minimum_xspacing, minimum_yspacing, minimum_zspacing, + + # Pointwise spacing, area, and volume operators + xspacing, yspacing, zspacing, λspacing, φspacing, xarea, yarea, zarea, volume, + + # Immersed boundaries + ImmersedBoundaryGrid, + GridFittedBoundary, GridFittedBottom, PartialCellBottom, + ImmersedBoundaryCondition, + + # Distributed + Distributed, Partition, + + # Advection schemes + Centered, UpwindBiased, WENO, + VectorInvariant, WENOVectorInvariant, FluxFormAdvection, + + # Boundary conditions + BoundaryCondition, + FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, OpenBoundaryCondition, + FieldBoundaryConditions, + + # Fields and field manipulation + Field, CenterField, XFaceField, YFaceField, ZFaceField, + Average, Integral, CumulativeIntegral, Reduction, Accumulation, BackgroundField, + interior, set!, compute!, regrid!, + + # Forcing functions + Forcing, Relaxation, LinearTarget, GaussianMask, AdvectiveForcing, + + # Coriolis forces + FPlane, ConstantCartesianCoriolis, BetaPlane, NonTraditionalBetaPlane, + + # BuoyancyModels and equations of state + Buoyancy, BuoyancyTracer, SeawaterBuoyancy, + LinearEquationOfState, TEOS10, + BuoyancyField, + + # Surface wave Stokes drift via Craik-Leibovich equations + UniformStokesDrift, StokesDrift, + + # Turbulence closures + VerticalScalarDiffusivity, + HorizontalScalarDiffusivity, + ScalarDiffusivity, + VerticalScalarBiharmonicDiffusivity, + HorizontalScalarBiharmonicDiffusivity, + ScalarBiharmonicDiffusivity, + SmagorinskyLilly, + Smagorinsky, + LillyCoefficient, + DynamicCoefficient, + AnisotropicMinimumDissipation, + ConvectiveAdjustmentVerticalDiffusivity, + CATKEVerticalDiffusivity, + RiBasedVerticalDiffusivity, + VerticallyImplicitTimeDiscretization, + viscosity, diffusivity, + + # Lagrangian particle tracking + LagrangianParticles, + + # Models + NonhydrostaticModel, HydrostaticFreeSurfaceModel, ShallowWaterModel, + ConservativeFormulation, VectorInvariantFormulation, + PressureField, fields, + + # Hydrostatic free surface model stuff + VectorInvariant, ExplicitFreeSurface, ImplicitFreeSurface, SplitExplicitFreeSurface, + HydrostaticSphericalCoriolis, PrescribedVelocityFields, + + # Time stepping + Clock, TimeStepWizard, conjure_time_step_wizard!, time_step!, + + # Simulations + Simulation, run!, Callback, add_callback!, iteration, + iteration_limit_exceeded, stop_time_exceeded, wall_time_limit_exceeded, + + # Diagnostics + CFL, AdvectiveCFL, DiffusiveCFL, + + # Output writers + NetCDFOutputWriter, JLD2OutputWriter, Checkpointer, + TimeInterval, IterationInterval, WallTimeInterval, AveragedTimeInterval, + SpecifiedTimes, FileSizeLimit, AndSchedule, OrSchedule, written_names, + + # Output readers + FieldTimeSeries, FieldDataset, InMemory, OnDisk, + + # Abstract operations + ∂x, ∂y, ∂z, @at, KernelFunctionOperation, + + # Utils + prettytime + +using CUDA +using DocStringExtensions +using FFTW + +function __init__() + threads = Threads.nthreads() + if threads > 1 + @info "Oceananigans will use $threads threads" + + # See: https://github.com/CliMA/Oceananigans.jl/issues/1113 + FFTW.set_num_threads(4threads) + end + + if CUDA.has_cuda() + @debug "CUDA-enabled GPU(s) detected:" + for (gpu, dev) in enumerate(CUDA.devices()) + @debug "$dev: $(CUDA.name(dev))" + end + + CUDA.allowscalar(false) + end +end + +##### +##### Abstract types +##### + +""" + AbstractModel + +Abstract supertype for models. +""" +abstract type AbstractModel{TS} end + +""" + AbstractDiagnostic + +Abstract supertype for diagnostics that compute information from the current +model state. +""" +abstract type AbstractDiagnostic end + +""" + AbstractOutputWriter + +Abstract supertype for output writers that write data to disk. +""" +abstract type AbstractOutputWriter end + +# Callsites for Callbacks +struct TimeStepCallsite end +struct TendencyCallsite end +struct UpdateStateCallsite end + +##### +##### Place-holder functions +##### + +function run_diagnostic! end +function write_output! end +function initialize! end # for initializing models, simulations, etc +function location end +function instantiated_location end +function tupleit end +function fields end +function prognostic_fields end +function tracer_tendency_kernel_function end +function boundary_conditions end + +##### +##### Include all the submodules +##### + +# Basics +include("Architectures.jl") +include("Units.jl") +include("Grids/Grids.jl") +include("Utils/Utils.jl") +include("Logger.jl") +include("Operators/Operators.jl") +include("BoundaryConditions/BoundaryConditions.jl") +include("Fields/Fields.jl") +include("AbstractOperations/AbstractOperations.jl") +include("TimeSteppers/TimeSteppers.jl") +include("ImmersedBoundaries/ImmersedBoundaries.jl") +include("Advection/Advection.jl") +include("Solvers/Solvers.jl") +include("OutputReaders/OutputReaders.jl") +include("DistributedComputations/DistributedComputations.jl") + +# TODO: move here +#include("MultiRegion/MultiRegion.jl") + +# Physics, time-stepping, and models +include("Coriolis/Coriolis.jl") +include("BuoyancyModels/BuoyancyModels.jl") +include("StokesDrifts.jl") +include("TurbulenceClosures/TurbulenceClosures.jl") +include("Forcings/Forcings.jl") +include("Biogeochemistry.jl") + +# TODO: move above +include("Models/Models.jl") + +# Output and Physics, time-stepping, and models +include("Diagnostics/Diagnostics.jl") +include("OutputWriters/OutputWriters.jl") +include("Simulations/Simulations.jl") + +# Abstractions for distributed and multi-region models +include("MultiRegion/MultiRegion.jl") + +include("boundary_mean.jl") + +##### +##### Needed so we can export names from sub-modules at the top-level +##### + +using .Logger +using .Architectures +using .Utils +using .Advection +using .Grids +using .BoundaryConditions +using .Fields +using .Coriolis +using .BuoyancyModels +using .StokesDrifts +using .TurbulenceClosures +using .Solvers +using .OutputReaders +using .Forcings +using .ImmersedBoundaries +using .DistributedComputations +using .Models +using .TimeSteppers +using .Diagnostics +using .OutputWriters +using .Simulations +using .AbstractOperations +using .MultiRegion + +end # module diff --git a/src/boundary_normal_mean_velocity.jl b/src/boundary_mean.jl similarity index 100% rename from src/boundary_normal_mean_velocity.jl rename to src/boundary_mean.jl From ff0c79ea4c9034b9e6b5af59d3563c24b5d33dde Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 12 Dec 2024 20:37:36 +0000 Subject: [PATCH 18/24] renamed --- src/boundary_mean.jl | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/src/boundary_mean.jl b/src/boundary_mean.jl index 1f6795eb67..54154397b8 100644 --- a/src/boundary_mean.jl +++ b/src/boundary_mean.jl @@ -1,42 +1,54 @@ using Adapt +using Oceananigans: location +using Oceananigans.Fields: Center, Face using Oceananigans.AbstractOperations: GridMetricOperation, Ax, Ay, Az using Oceananigans.BoundaryConditions: BoundaryCondition, Open, PerturbationAdvection import Adapt: adapt_structure import Oceananigans.BoundaryConditions: update_boundary_condition! -struct BoundaryMean{V} +struct BoundaryAdjacentMean{V} value :: V - BoundaryMean(; value::BV = Ref(0.0)) where BV = new{BV}(value) + BoundaryAdjacentMean(; value::BV = Ref(0.0)) where BV = new{BV}(value) end -@inline (value::BoundaryMean)(args...) = value.boundary_velocity[] +@inline (bam::BoundaryAdjacentMean)(args...) = bam.value[] -Adapt.adapt_structure(to, mo::BoundaryMean) = - BoundaryMean(adapt(to, mo.mean_outflow_velocity[])) +Adapt.adapt_structure(to, mo::BoundaryAdjacentMean) = + BoundaryAdjacentMean(adapt(to, mo.mean_outflow_velocity[])) -const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryMean} +const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryAdjacentMean} @inline boundary_normal_area(::Union{Val{:west}, Val{:east}}, grid) = GridMetricOperation((Face, Center, Center), Ax, grid) @inline boundary_normal_area(::Union{Val{:south}, Val{:north}}, grid) = GridMetricOperation((Center, Face, Center), Ay, grid) @inline boundary_normal_area(::Union{Val{:bottom}, Val{:top}}, grid) = GridMetricOperation((Center, Center, Face), Az, grid) -@inline boundary_adjacent_index(::Val{:east}, grid) = (size(grid, 1), 1, 1), (2, 3) -@inline boundary_adjacent_index(::Val{:west}, grid) = (2, 1, 1), (2, 3) +@inline boundary_adjacent_index(::Val{:east}, grid, loc) = (size(grid, 1), 1, 1), (2, 3) +@inline boundary_adjacent_index(::Val{:west}, grid, loc) = (first_interior_index(side, loc, grid), 1, 1), (2, 3) -@inline boundary_adjacent_index(::Val{:north}, grid) = (1, size(grid, 2), 1), (1, 3) -@inline boundary_adjacent_index(::Val{:south}, grid) = (1, 2, 1), (1, 3) +@inline boundary_adjacent_index(::Val{:north}, grid, loc) = (1, size(grid, 2), 1), (1, 3) +@inline boundary_adjacent_index(::Val{:south}, grid, loc) = (1, first_interior_index(side, loc, grid), 1), (1, 3) -@inline boundary_adjacent_index(::Val{:top}, grid) = (1, 1, size(grid, 3)), (2, 3) -@inline boundary_adjacent_index(::Val{:bottom}, grid) = (1, 1, 2), (2, 3) +@inline boundary_adjacent_index(::Val{:top}, grid, loc) = (1, 1, size(grid, 3)), (2, 3) +@inline boundary_adjacent_index(::Val{:bottom}, grid, loc) = (1, 1, first_interior_index(side, loc, grid)), (2, 3) + +@inline first_interior_index(::Union{Val{:west}, Val{:east}}, ::Tuple{<:Center, <:Any, <:Any}, grid) = 1 +@inline first_interior_index(::Union{Val{:west}, Val{:east}}, ::Tuple{Face, <:Any, <:Any}, grid) = 2 + +@inline first_interior_index(::Union{Val{:south}, Val{:north}}, ::Tuple{<:Any, <:Center, <:Any}, grid) = 1 +@inline first_interior_index(::Union{Val{:south}, Val{:north}}, ::Tuple{<:Any, Face, <:Any}, grid) = 2 + +@inline first_interior_index(::Union{Val{:bottom}, Val{:top}}, ::Tuple{<:Any, <:Any, <:Center}, grid) = 1 +@inline first_interior_index(::Union{Val{:bottom}, Val{:top}}, ::Tuple{<:Any, <:Any, Face}, grid) = 2 function update_boundary_condition!(bc::MOPABC, val_side, u, model) grid = model.grid + loc = location(u) An = boundary_normal_area(val_side, grid) - (i, j, k), dims = boundary_adjacent_index(val_side, grid) + (i, j, k), dims = boundary_adjacent_index(val_side, grid, loc) total_area = sum(An; dims)[i, j, k] From 3257eec19a039627d26475727f70d14f713ab823 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 12 Dec 2024 21:26:09 +0000 Subject: [PATCH 19/24] kind of GPU frieldly --- src/boundary_mean.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/boundary_mean.jl b/src/boundary_mean.jl index 54154397b8..7da3921821 100644 --- a/src/boundary_mean.jl +++ b/src/boundary_mean.jl @@ -16,7 +16,7 @@ end @inline (bam::BoundaryAdjacentMean)(args...) = bam.value[] Adapt.adapt_structure(to, mo::BoundaryAdjacentMean) = - BoundaryAdjacentMean(adapt(to, mo.mean_outflow_velocity[])) + BoundaryAdjacentMean(; value = adapt(to, mo.value[])) const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryAdjacentMean} @@ -50,9 +50,11 @@ function update_boundary_condition!(bc::MOPABC, val_side, u, model) (i, j, k), dims = boundary_adjacent_index(val_side, grid, loc) - total_area = sum(An; dims)[i, j, k] + total_area = CUDA.@allowscalar sum(An; dims)[i, j, k] Ū = sum(u * An; dims) / total_area - bc.condition.value[] = Ū[i, j, k] + bc.condition.value[] = CUDA.@allowscalar Ū[i, j, k] + + return nothing end \ No newline at end of file From 07cc4711a3cf24ba61f29f67477e4b2328673b2f Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 13 Dec 2024 17:59:17 +0000 Subject: [PATCH 20/24] typo --- src/boundary_mean.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/boundary_mean.jl b/src/boundary_mean.jl index 7da3921821..0175973e42 100644 --- a/src/boundary_mean.jl +++ b/src/boundary_mean.jl @@ -25,13 +25,13 @@ const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryAdja @inline boundary_normal_area(::Union{Val{:bottom}, Val{:top}}, grid) = GridMetricOperation((Center, Center, Face), Az, grid) @inline boundary_adjacent_index(::Val{:east}, grid, loc) = (size(grid, 1), 1, 1), (2, 3) -@inline boundary_adjacent_index(::Val{:west}, grid, loc) = (first_interior_index(side, loc, grid), 1, 1), (2, 3) +@inline boundary_adjacent_index(side::Val{:west}, grid, loc) = (first_interior_index(side, loc, grid), 1, 1), (2, 3) @inline boundary_adjacent_index(::Val{:north}, grid, loc) = (1, size(grid, 2), 1), (1, 3) -@inline boundary_adjacent_index(::Val{:south}, grid, loc) = (1, first_interior_index(side, loc, grid), 1), (1, 3) +@inline boundary_adjacent_index(side::Val{:south}, grid, loc) = (1, first_interior_index(side, loc, grid), 1), (1, 3) @inline boundary_adjacent_index(::Val{:top}, grid, loc) = (1, 1, size(grid, 3)), (2, 3) -@inline boundary_adjacent_index(::Val{:bottom}, grid, loc) = (1, 1, first_interior_index(side, loc, grid)), (2, 3) +@inline boundary_adjacent_index(side::Val{:bottom}, grid, loc) = (1, 1, first_interior_index(side, loc, grid)), (2, 3) @inline first_interior_index(::Union{Val{:west}, Val{:east}}, ::Tuple{<:Center, <:Any, <:Any}, grid) = 1 @inline first_interior_index(::Union{Val{:west}, Val{:east}}, ::Tuple{Face, <:Any, <:Any}, grid) = 2 From 9ae0b6683d341d8c1830af329f361b47b211149f Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 13 Dec 2024 18:03:06 +0000 Subject: [PATCH 21/24] bug --- src/boundary_mean.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/boundary_mean.jl b/src/boundary_mean.jl index 0175973e42..d94f40606d 100644 --- a/src/boundary_mean.jl +++ b/src/boundary_mean.jl @@ -25,22 +25,22 @@ const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryAdja @inline boundary_normal_area(::Union{Val{:bottom}, Val{:top}}, grid) = GridMetricOperation((Center, Center, Face), Az, grid) @inline boundary_adjacent_index(::Val{:east}, grid, loc) = (size(grid, 1), 1, 1), (2, 3) -@inline boundary_adjacent_index(side::Val{:west}, grid, loc) = (first_interior_index(side, loc, grid), 1, 1), (2, 3) +@inline boundary_adjacent_index(side::Val{:west}, grid, loc) = (first_interior_index(side, loc), 1, 1), (2, 3) @inline boundary_adjacent_index(::Val{:north}, grid, loc) = (1, size(grid, 2), 1), (1, 3) -@inline boundary_adjacent_index(side::Val{:south}, grid, loc) = (1, first_interior_index(side, loc, grid), 1), (1, 3) +@inline boundary_adjacent_index(side::Val{:south}, grid, loc) = (1, first_interior_index(side, loc), 1), (1, 3) -@inline boundary_adjacent_index(::Val{:top}, grid, loc) = (1, 1, size(grid, 3)), (2, 3) -@inline boundary_adjacent_index(side::Val{:bottom}, grid, loc) = (1, 1, first_interior_index(side, loc, grid)), (2, 3) +@inline boundary_adjacent_index(::Val{:top}, grid, loc) = (1, 1, size(grid, 3)), (2, 3) +@inline boundary_adjacent_index(side::Val{:bottom}, grid, loc) = (1, 1, first_interior_index(side, loc)), (2, 3) -@inline first_interior_index(::Union{Val{:west}, Val{:east}}, ::Tuple{<:Center, <:Any, <:Any}, grid) = 1 -@inline first_interior_index(::Union{Val{:west}, Val{:east}}, ::Tuple{Face, <:Any, <:Any}, grid) = 2 +@inline first_interior_index(::Union{Val{:west}, Val{:east}}, ::Tuple{Center, <:Any, <:Any}) = 1 +@inline first_interior_index(::Union{Val{:west}, Val{:east}}, ::Tuple{Face, <:Any, <:Any}) = 2 -@inline first_interior_index(::Union{Val{:south}, Val{:north}}, ::Tuple{<:Any, <:Center, <:Any}, grid) = 1 -@inline first_interior_index(::Union{Val{:south}, Val{:north}}, ::Tuple{<:Any, Face, <:Any}, grid) = 2 +@inline first_interior_index(::Union{Val{:south}, Val{:north}}, ::Tuple{<:Any, Center, <:Any}) = 1 +@inline first_interior_index(::Union{Val{:south}, Val{:north}}, ::Tuple{<:Any, Face, <:Any}) = 2 -@inline first_interior_index(::Union{Val{:bottom}, Val{:top}}, ::Tuple{<:Any, <:Any, <:Center}, grid) = 1 -@inline first_interior_index(::Union{Val{:bottom}, Val{:top}}, ::Tuple{<:Any, <:Any, Face}, grid) = 2 +@inline first_interior_index(::Union{Val{:bottom}, Val{:top}}, ::Tuple{<:Any, <:Any, Center}) = 1 +@inline first_interior_index(::Union{Val{:bottom}, Val{:top}}, ::Tuple{<:Any, <:Any, Face}) = 2 function update_boundary_condition!(bc::MOPABC, val_side, u, model) grid = model.grid From ec9dd88b2fc49c77953a4bd84bc42114e58192a6 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 13 Dec 2024 18:05:05 +0000 Subject: [PATCH 22/24] bug --- src/boundary_mean.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/boundary_mean.jl b/src/boundary_mean.jl index d94f40606d..3c2d4bbd41 100644 --- a/src/boundary_mean.jl +++ b/src/boundary_mean.jl @@ -1,5 +1,5 @@ using Adapt -using Oceananigans: location +using Oceananigans: instantiated_location using Oceananigans.Fields: Center, Face using Oceananigans.AbstractOperations: GridMetricOperation, Ax, Ay, Az using Oceananigans.BoundaryConditions: BoundaryCondition, Open, PerturbationAdvection @@ -25,13 +25,13 @@ const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryAdja @inline boundary_normal_area(::Union{Val{:bottom}, Val{:top}}, grid) = GridMetricOperation((Center, Center, Face), Az, grid) @inline boundary_adjacent_index(::Val{:east}, grid, loc) = (size(grid, 1), 1, 1), (2, 3) -@inline boundary_adjacent_index(side::Val{:west}, grid, loc) = (first_interior_index(side, loc), 1, 1), (2, 3) +@inline boundary_adjacent_index(val_side::Val{:west}, grid, loc) = (first_interior_index(val_side, loc), 1, 1), (2, 3) @inline boundary_adjacent_index(::Val{:north}, grid, loc) = (1, size(grid, 2), 1), (1, 3) -@inline boundary_adjacent_index(side::Val{:south}, grid, loc) = (1, first_interior_index(side, loc), 1), (1, 3) +@inline boundary_adjacent_index(val_side::Val{:south}, grid, loc) = (1, first_interior_index(val_side, loc), 1), (1, 3) @inline boundary_adjacent_index(::Val{:top}, grid, loc) = (1, 1, size(grid, 3)), (2, 3) -@inline boundary_adjacent_index(side::Val{:bottom}, grid, loc) = (1, 1, first_interior_index(side, loc)), (2, 3) +@inline boundary_adjacent_index(val_side::Val{:bottom}, grid, loc) = (1, 1, first_interior_index(val_side, loc)), (2, 3) @inline first_interior_index(::Union{Val{:west}, Val{:east}}, ::Tuple{Center, <:Any, <:Any}) = 1 @inline first_interior_index(::Union{Val{:west}, Val{:east}}, ::Tuple{Face, <:Any, <:Any}) = 2 @@ -44,7 +44,7 @@ const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryAdja function update_boundary_condition!(bc::MOPABC, val_side, u, model) grid = model.grid - loc = location(u) + loc = instantiated_location(u) An = boundary_normal_area(val_side, grid) From 96495b10ab82858428f20e8db406b87415709159 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Sun, 15 Dec 2024 19:26:40 +0000 Subject: [PATCH 23/24] allow forward euler integration --- ...advection_open_boundary_matching_scheme.jl | 71 +++++++++++++++++-- 1 file changed, 65 insertions(+), 6 deletions(-) diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl index 0f12b30c50..6714f4aa6b 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl @@ -11,7 +11,8 @@ see latex document for now TODO: check what the coriolis is doing, and check what happens if U is the mean velocity """ -struct PerturbationAdvection{FT} +struct PerturbationAdvection{VT, FT} + backward_step :: VT inflow_timescale :: FT outflow_timescale :: FT end @@ -21,10 +22,11 @@ Adapt.adapt_structure(to, pe::PerturbationAdvection) = adapt(to, pe.inflow_timescale)) function PerturbationAdvectionOpenBoundaryCondition(val, FT = Float64; + backward_step = true, outflow_timescale = Inf, inflow_timescale = 300.0, kwargs...) - classification = Open(PerturbationAdvection(inflow_timescale, outflow_timescale)) + classification = Open(PerturbationAdvection(Val(backward_step), inflow_timescale, outflow_timescale)) @warn "`PerturbationAdvection` open boundaries matching scheme is experimental and un-tested/validated" @@ -33,7 +35,64 @@ end const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} -@inline function step_right_boundary!(bc, l, m, boundary_indices, boundary_adjacent_indices, +const BPAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection{Val{true}}}} +const FPAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection{Val{false}}}} + +@inline function step_right_boundary!(bc::BPAOBC, l, m, boundary_indices, boundary_adjacent_indices, + grid, u, clock, model_fields, ΔX) + Δt = clock.last_stage_Δt + + Δt = ifelse(isinf(Δt), 0, Δt) + + ūⁿ⁺¹ = getbc(bc, l, m, grid, clock, model_fields) + + uᵢⁿ = @inbounds getindex(u, boundary_indices...) + uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, boundary_adjacent_indices...) + + U = max(0, min(1, Δt / ΔX * ūⁿ⁺¹)) + + pa = bc.classification.matching_scheme + + τ = ifelse(ūⁿ⁺¹ >= 0, pa.outflow_timescale, pa.inflow_timescale) + + τ̃ = Δt / τ + + uᵢⁿ⁺¹ = uᵢⁿ + U * (uᵢ₋₁ⁿ⁺¹ - ūⁿ⁺¹) + + @inbounds setindex!(u, uᵢⁿ⁺¹, boundary_indices...) + + return nothing +end + +@inline function step_left_boundary!(bc::BPAOBC, l, m, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, + grid, u, clock, model_fields, ΔX) + Δt = clock.last_stage_Δt + + Δt = ifelse(isinf(Δt), 0, Δt) + + ūⁿ⁺¹ = getbc(bc, l, m, grid, clock, model_fields) + + uᵢⁿ = @inbounds getindex(u, boundary_secret_storage_indices...) + uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, boundary_adjacent_indices...) + + U = min(0, max(-1, Δt / ΔX * ūⁿ⁺¹)) + + pa = bc.classification.matching_scheme + + τ = ifelse(ūⁿ⁺¹ <= 0, pa.outflow_timescale, pa.inflow_timescale) + + τ̃ = Δt / τ + + u₁ⁿ⁺¹ = uᵢⁿ - U * (uᵢ₋₁ⁿ⁺¹ - ūⁿ⁺¹) + + @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_indices...) + @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_secret_storage_indices...) + + return nothing +end + + +@inline function step_right_boundary!(bc::FPAOBC, l, m, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, ΔX) Δt = clock.last_stage_Δt @@ -52,14 +111,14 @@ const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}} τ̃ = Δt / τ - uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + U) + uᵢⁿ⁺¹ = uᵢⁿ + U * (uᵢ₋₁ⁿ⁺¹ - ūⁿ⁺¹) + (ūⁿ⁺¹ - uᵢⁿ) * τ̃ @inbounds setindex!(u, uᵢⁿ⁺¹, boundary_indices...) return nothing end -@inline function step_left_boundary!(bc, l, m, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, +@inline function step_left_boundary!(bc::FPAOBC, l, m, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, grid, u, clock, model_fields, ΔX) Δt = clock.last_stage_Δt @@ -78,7 +137,7 @@ end τ̃ = Δt / τ - u₁ⁿ⁺¹ = (uᵢⁿ - U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ - U) + u₁ⁿ⁺¹ = uᵢⁿ - U * (uᵢ₋₁ⁿ⁺¹ - ūⁿ⁺¹) + (ūⁿ⁺¹ - uᵢⁿ) * τ̃ @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_indices...) @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_secret_storage_indices...) From c49ed2645d0a91428631a648f246bdd9afcbe539 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Sun, 15 Dec 2024 19:34:43 +0000 Subject: [PATCH 24/24] fix adapt --- .../perturbation_advection_open_boundary_matching_scheme.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl index 6714f4aa6b..aa929eb971 100644 --- a/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl @@ -18,7 +18,8 @@ struct PerturbationAdvection{VT, FT} end Adapt.adapt_structure(to, pe::PerturbationAdvection) = - PerturbationAdvection(adapt(to, pe.outflow_timescale), + PerturbationAdvection(adapt(to, pe.backward_step), + adapt(to, pe.outflow_timescale), adapt(to, pe.inflow_timescale)) function PerturbationAdvectionOpenBoundaryCondition(val, FT = Float64;