Adds a pertubation advection open boundary matching scheme

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

src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl

@@ -0,0 +1,90 @@
+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
+
+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)
+
+
+    τ̃ = Δt / τ
+
+    uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + 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)
+
+    τ̃ = Δ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))
+end

validation/open_boundaries/cylinder.jl

@@ -1,127 +1,214 @@
-# 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.Models.NonhydrostaticModels: ConjugateGradientPoissonSolver
+using Oceananigans.Solvers: DiagonallyDominantPreconditioner
+using Oceananigans.Operators: ℑxyᶠᶜᵃ, ℑxyᶜᶠᵃ
+using Oceananigans.Solvers: FFTBasedPoissonSolver
+using Printf
+using CUDA
 
-using Oceananigans, CairoMakie
-using Oceananigans.BoundaryConditions: FlatExtrapolationOpenBoundaryCondition
+using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition
 
+u∞ = 1
+r = 1/2
+arch = GPU()
+stop_time = 200
 
-@kwdef struct Cylinder{FT}
-    D :: FT = 1.0
-   x₀ :: FT = 0.0
-   y₀ :: FT = 0.0
-end
+cylinder(x, y) = (x^2 + y^2) ≤ r^2
 
-@inline (cylinder::Cylinder)(x, y) = ifelse((x - cylinder.x₀)^2 + (y - cylinder.y₀)^2 < (cylinder.D/2)^2, 1, 0)
+"""
+    drag(modell; bounding_box = (-1, 3, -2, 2), ν = 1e-3)
 
-architecture = GPU()
+Returns the drag within the `bounding_box` computed by:
 
-# model parameters
-Re = 200
-U = 1
-D = 1.
-resolution = D / 40
+```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)
 
-# add extra downstream distance to see if the solution near the cylinder changes
-extra_downstream = 0
+    u, v, _ = model.velocities
 
-cylinder = Cylinder(; D)
+    uᶜ = Field(@at (Center, Center, Center) u)
+    vᶜ = Field(@at (Center, Center, Center) v)
 
-x = (-5, 5 + extra_downstream) .* D
-y = (-5, 5) .* D
+    xc, yc, _ = nodes(uᶜ)
 
-Ny = Int(10 / resolution)
-Nx = Ny + Int(extra_downstream / resolution)
+    i₁ = findfirst(xc .> bounding_box[1])
+    i₂ = findlast(xc .< bounding_box[2])
 
-ν = U * D / Re
+    j₁ = findfirst(yc .> bounding_box[3])
+    j₂ = findlast(yc .< bounding_box[4])
 
-closure = ScalarDiffusivity(;ν, κ = ν)
+    uₗ² = Field(uᶜ^2, indices = (i₁, j₁:j₂, 1))
+    uᵣ² = Field(uᶜ^2, indices = (i₂, j₁:j₂, 1))
 
-grid = RectilinearGrid(architecture; topology = (Bounded, Periodic, Flat), size = (Nx, Ny), x, y)
+    uvₗ = Field(uᶜ*vᶜ, indices = (i₁:i₂, j₁, 1))
+    uvᵣ = Field(uᶜ*vᶜ, indices = (i₁:i₂, j₂, 1))
 
-@inline u(y, t, U) = U * (1 + 0.01 * randn())
+    ∂₁uₗ = Field(∂x(uᶜ), indices = (i₁, j₁:j₂, 1))
+    ∂₁uᵣ = Field(∂x(uᶜ), indices = (i₂, j₁:j₂, 1))
 
-u_boundaries = FieldBoundaryConditions(east = FlatExtrapolationOpenBoundaryCondition(),
-                                       west = OpenBoundaryCondition(u, parameters = U))
+    ∂ₜuᶜ = Field(@at (Center, Center, Center) model.timestepper.Gⁿ.u)
 
-v_boundaries = FieldBoundaryConditions(east = GradientBoundaryCondition(0),
-                                       west = GradientBoundaryCondition(0))
+    ∂ₜu = Field(∂ₜuᶜ, indices = (i₁:i₂, j₁:j₂, 1))
 
-Δt = .3 * resolution / U
+    p = model.pressures.pNHS
 
-u_forcing = Relaxation(; rate = 1 / (2 * Δt), mask = cylinder)
-v_forcing = Relaxation(; rate = 1 / (2 * Δt), mask = cylinder) 
+    ∫∂ₓp = Field(∂x(p), indices = (i₁:i₂, j₁:j₂, 1))
 
-model = NonhydrostaticModel(; grid, 
-                              closure, 
-                              forcing = (u = u_forcing, v = v_forcing),
-                              boundary_conditions = (u = u_boundaries, v = v_boundaries))
+    a_local = Field(Integral(∂ₜu))
 
-@info "Constructed model"
+    a_flux = Field(Integral(uᵣ²)) - Field(Integral(uₗ²)) + Field(Integral(uvᵣ)) - Field(Integral(uvₗ))
 
-# initial noise to induce turbulance faster
-set!(model, u = U, v = (x, y) -> randn() * U * 0.01)
+    a_viscous_stress = 2ν * (Field(Integral(∂₁uᵣ)) - Field(Integral(∂₁uₗ)))
 
-@info "Set initial conditions"
+    a_pressure = Field(Integral(∫∂ₓp))
 
-simulation = Simulation(model; Δt = Δt, stop_time = 300)
+    return a_local + a_flux + a_pressure - a_viscous_stress
+end
 
-wizard = TimeStepWizard(cfl = 0.3)
+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
+=#
 
-simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(100))
+Ny = 2048 
+Nx = Ny
 
-progress(sim) = @info "$(time(sim)) with Δt = $(prettytime(sim.Δt)) in $(prettytime(sim.run_wall_time))"
+prefix = "flow_around_cylinder_Re$(Re)_Ny$(Ny)"
 
-simulation.callbacks[:progress] = Callback(progress, IterationInterval(1000))
+ϵ = 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...)
 
-simulation.output_writers[:velocity] = JLD2OutputWriter(model, model.velocities,
-                                                        overwrite_existing = true, 
-                                                        filename = "cylinder_$(extra_downstream)_Re_$Re.jld2", 
-                                                        schedule = TimeInterval(1),
-                                                        with_halos = true)
+grid = ImmersedBoundaryGrid(grid, GridFittedBoundary(cylinder))
 
-run!(simulation)
+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
 
-# load the results 
+uᵢ(x, y) = 1e-2 * randn()
+vᵢ(x, y) = 1e-2 * randn()
+set!(model, u=uᵢ, v=vᵢ)
 
-u_ts = FieldTimeSeries("cylinder_$(extra_downstream)_Re_$Re.jld2", "u")
-v_ts = FieldTimeSeries("cylinder_$(extra_downstream)_Re_$Re.jld2", "v")
+Δx = minimum_xspacing(grid)
+Δt = max_Δt = 0.002#0.2 * Δx^2 * Re
 
-u′, v′, w′ = Oceananigans.Fields.VelocityFields(u_ts.grid)
+simulation = Simulation(model; Δt, stop_time)
+#conjure_time_step_wizard!(simulation, cfl=1.0, IterationInterval(3); max_Δt)
 
-ζ = Field((@at (Center, Center, Center) ∂x(v′)) - (@at (Center, Center, Center) ∂y(u′)))
+u, v, w = model.velocities
+#d = ∂x(u) + ∂y(v)
 
-# there is probably a more memory efficient way todo this
+# Drag computation
+drag_force = drag(model; ν=1/Re)
+compute!(drag_force)
 
-ζ_ts = zeros(size(grid, 1), size(grid, 2), length(u_ts.times)) # u_ts.grid so its always on cpu
+wall_time = Ref(time_ns())
 
-for n in 1:length(u_ts.times)
-    set!(u′, u_ts[n])
-    set!(v′, v_ts[n])
-    compute!(ζ)
-    ζ_ts[:, :, n] = interior(ζ, :, :, 1)
+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
 
-@info "Loaded results"
+add_callback!(simulation, progress, IterationInterval(100))
 
-# plot the results
+ζ = ∂x(v) - ∂y(u)
 
-fig = Figure(size = (600, 600))
+p = model.pressures.pNHS
 
-ax = Axis(fig[1, 1], aspect = DataAspect())
+outputs = (; u, v, p, ζ)
 
-xc, yc, zc = nodes(ζ)
+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)
 
-n = Observable(1)
+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))
+
+run!(simulation)
 
-ζ_plt = @lift ζ_ts[:, :, $n]
 
-contour!(ax, xc, yc, ζ_plt, levels = [-2, 2], colorrange = (-2, 2), colormap = :roma)
+# 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
 
-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
+# Re = 1000
+# with 12m ~ HZ and Cd ~