Adds a pertubation advection open boundary matching scheme

[jagoosw]

This PR adds a PertubationAdvection open boundary condition which assumes that there is some mean flow (i.e. prescribed in a channel or from a coarser resolution model) which is homogeneous on the boundary but can vary in time.

On the boundary, we approximate (for an x-boundary) to:

$$\frac{\partial u}{\partial t} \approx -U\frac{\partial u\prime}{\partial x},$$

where $u=u\prime+U$
I have chosen to take an explicit/backwards Euler step:

$$\frac{u\prime(t^{n+1}) - u\prime(t^n)}{\Delta t} + \frac{U\prime(t^{n+1}) - U\prime(t^n)}{\Delta t} \approx -U\frac{u\prime(x_b, t^{n+1}) - u\prime(x_{b-1}, t^{n+1})}{\Delta x},$$

because then we don't have to store any other information (i.e. if we did a forward step we would need $u\prime(x_{b-1}, t^n)$ ). This results in:

$$u(t^{n+1}) = \frac{u^n+\tilde{U}u(x_{b-1}, t^{n+1})}{1+\tilde{U}},$$

where $\tilde{U}=\min\left(1, U\Delta t / \Delta x\right)$. I have also added restoring to $U$ (i.e. damping $u\prime\to0$) in some timescale $\tau$ which can be different for inflow and outflow to allow the scheme to work when there is inflow, which results in:

$$u(t^{n+1}) = \frac{u^n+\tilde{U}u(x_{b-1}, t^{n+1}) + U\tilde{\tau}}{1+\tilde{U}+\tilde{\tau}},$$

where $\tilde{\tau}=\Delta t / \tau$.

This has the limitation that if we want to extend to have wall tangent mean velocity ($\partial_t u \approx -(\vec{U}\cdot\nabla) u$) then we either have to solve the whole boundary at once since every point is going to depend on its boundary neighbours, or we need to take a forward Euler step in which case we need the first interior point at the previous time step.

I think we probably will need to address this for real nesting cases at some point.

u.mp4

For now, it seems to be working okay (left plot vorticity, right plot x-velocity):

(old plot)

https://github.com/user-attachments/assets/21ae3eb2-5d3b-4c08-86fb-bb9c222e8e6a

and the drag coefficient is ~1.4 at Re = 100 and I'm running a case at Re=1000 now.

I have also checked the Strouhal number which matches exactly to the expected values (~0.17).

[jagoosw]

@tomchor, do you think something like this:

Oceananigans.jl/validation/open_boundaries/cylinder.jl

Lines 18 to 76 in 32cd354

	"""
	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

could belong in Oceanostics?

[jagoosw]

I have also run the same case with different domain lengths (6m, 12m - shown, 18m and 30m) and they all produce very similar Cd and St

[tomchor]

@tomchor, do you think something like this:

Oceananigans.jl/validation/open_boundaries/cylinder.jl

Lines 18 to 76 in 32cd354

	"""
	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

could belong in Oceanostics?

Sure! Why not? :)

[tomchor]

I think we could implement everything using getindex() so that we implement it once (or at max twice; one for left and other for right BC) and can use it in three dimensions.

For example this

    uᵢⁿ     = @inbounds u[i, j, k]
    uᵢ₋₁ⁿ⁺¹ = @inbounds u[i - 1, j, k]

would become

    uᵢⁿ     = @inbounds getindex(u, boundary_indices...)
    uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, one_off_boundary_indices...)

This would avoid errors when transcribing to other dimensions (I remember catching a couple for the flat extrapolation matching scheme).

[jagoosw]

yeah, I was trying to think of a clean way to do this

[jagoosw]

Attempted this now (but I've still not written out the other directions in case we want to change things first)

[glwagner]

There's a typo in your top equation right? I think it should be $\partial_t (u + U) = U \partial_x u$.

[glwagner]

A problem I think this scheme might have is that as the flow changes direction the assumption that the perturbation is small falls apart, which is maybe why it's been a bit unstable so far in an oscillating case.

What do you mean? What is the consequence of this approximation?

[jagoosw]

What do you mean? What is the consequence of this approximation?

I was thinking about this more and really to remain valid then in :

$$\partial_t (u\prime+U)\approx u\prime\partial_xu\prime + U\partial_xu\prime + (U - u\prime)/\tau$$

as $U\to0$ then we need $u\prime\partial_xu\prime\ll u\prime/\tau$ to still be able to throw away the first term, so I think its actually okay as long as you have the relaxation.

I found a typo in the west boundary which might be why the oscillating case wasn't working, rerunning now

[jagoosw]

There's a typo in your top equation right? I think it should be ∂ t ( u + U ) = U ∂ x u .
Yeah thanks for catching!

[glwagner]

What do you mean? What is the consequence of this approximation?

I was thinking about this more and really to remain valid then in :

as U → 0 then we need u ′ ∂ x u ′ ≪ u ′ / τ to still be able to throw away the first term, so I think its actually okay as long as you have the relaxation.

I found a typo in the west boundary which might be why the oscillating case wasn't working, rerunning now

Okay, I agree with you. You can also analyze the update formula / backward Euler step in the limit U -> 0. It doesn't look like it needs to be regularized in any way in that case though, it looks ok.

[tomchor]

Is there a reference for this method? If there is, please include it on the methods docstring. If not, it'd be good to include a brief explanation like the one at the top comment in the dosctring. I don't think this radiation method is very well known.

[jagoosw]

Is this not rendered correctly too @tomchor ?

[tomchor]

If I use FiraCode on chrome then some of the stuff in the first line of code doesn't render correctly:

On firefox it all renders fine:

And when I manually switch to JuliaMono it all renders fine.

But after finding out that this is all browser/machine dependent, I do think we shouldn't spend too much time trying to "fix" this. I say we move on from tweaking the unicode and address the rest of the comments so that we can merge.

[jagoosw]

I can no longer see the buildkite logs (its says I have to be part of the buildkite organisation) so could someone check what the problem with those tests was please?

[tomchor]

A couple of them are failing with an Out of GPU memory error, but there's also this error:

[tomchor]

We still have a comment that hasn't been addressed, right? About what to do with the boundary_mean.jl file. After that's dealt with, I can make tests pass if you want.

[jagoosw]

We still have a comment that hasn't been addressed, right? About what to do with the boundary_mean.jl file. After that's dealt with, I can make tests pass if you want.

yeah that's the only problem, whats wrong with the tests?

[tomchor]

We still have a comment that hasn't been addressed, right? About what to do with the boundary_mean.jl file. After that's dealt with, I can make tests pass if you want.

yeah that's the only problem, whats wrong with the tests?

No idea what's up with GPU errors. @glwagner mentioned at the some point that that server is pretty full, so I'm assuming some cleanup will be done soon. (I don't have access to that server btw, so I only know what people tell me.)

The error I posted in #3977 (comment) seems to be the only actual error rn though.

[jagoosw]

On the boundary mean issue, perhaps we should make a new module called something like HelperFunctions since I imagine there's other things that are easiest to define with all of the other infrastructure (I.e. abstract operators, etc.) so could be user defined, but are also nice to have in the source code?

@glwagner any thoughts?

[glwagner]

On the boundary mean issue, perhaps we should make a new module called something like HelperFunctions since I imagine there's other things that are easiest to define with all of the other infrastructure (I.e. abstract operators, etc.) so could be user defined, but are also nice to have in the source code?

@glwagner any thoughts?

Can you elaborate? What you're suggesting might make sense but I need more background, motivation, design proposal, etc.

[jagoosw]

Can you elaborate? What you're suggesting might make sense but I need more background, motivation, design proposal, etc.

Basically, the best and most general way to implement a boundary condition that uses the boundary mean velocity is by using AbstractOperations, but the problem is that they're defined after boundary conditions so it can't go in the boundary conditions module.

I think there is kind of a grey area of things like this that aren't part of the core functionality of the code and could be implemented by a user, but are important and useful for us to support and maintain.

[glwagner]

Can you elaborate? What you're suggesting might make sense but I need more background, motivation, design proposal, etc.

Basically, the best and most general way to implement a boundary condition that uses the boundary mean velocity is by using AbstractOperations, but the problem is that they're defined after boundary conditions so it can't go in the boundary conditions module.

I think there is kind of a grey area of things like this that aren't part of the core functionality of the code and could be implemented by a user, but are important and useful for us to support and maintain.

Possibly, this belongs in Models then. Are we assuming something specific about the presence of a velocity field vs tracer field? In that case, this is not a "finite volume engine" capability, but rather a "fluid dynamics model" capability, and it does not belong in BoundaryConditions (which has primitive utilities that should be independent of fluid dynamics).

[jagoosw]

Possibly, this belongs in Models then. Are we assuming something specific about the presence of a velocity field vs tracer field? In that case, this is not a "finite volume engine" capability, but rather a "fluid dynamics model" capability, and it does not belong in BoundaryConditions (which has primitive utilities that should be independent of fluid dynamics) so isn't fluid specific.

That's a good idea, it doesn't make any assumptions about what it is and is just computing the mean of whatever field it is associated with on a boundary:

Oceananigans.jl/src/boundary_mean.jl

Lines 45 to 56 in 045a6b3

	function update_boundary_condition!(bc::MOPABC, val_side, u, model)
	grid = model.grid
	loc = instantiated_location(u)
	An = boundary_normal_area(val_side, grid)
	(i, j, k), dims = boundary_adjacent_indices(val_side, grid, loc)
	total_area = CUDA.@allowscalar sum(An; dims)[i, j, k]
	Ū = sum(u * An; dims) / total_area
	bc.condition.value[] = CUDA.@allowscalar Ū[i, j, k]
	return nothing
	end

[glwagner]

shouldn't you define a kernel for sum? the way implemented seems like its going to be slow

[glwagner]

Also this function takes in model which to me clearly says it belongs in Models. This must only work with specific models right? It will not work as a generic boundary condition for any Field / with any fill_halo_regions!.

[jagoosw]

Also this function takes in model which to me clearly says it belongs in Models. This must only work with specific models right? It will not work as a generic boundary condition for any Field / with any fill_halo_regions!.

Model doesn't actually do anything (if we wrote u.grid rather than model.grid) so maybe I should rewrite this like:

function (bam::BoundaryAdjacentMean)(val_side, u)
   ... the above function
   return bam.value
end

update_boundary_condition!(bc::MOPABC, val_side, u, model) = (bc.condition(val_side, u); nothing)

[jagoosw]

shouldn't you define a kernel for sum? the way implemented seems like its going to be slow

Maybe? Doesn't sum go back to this reduction:

Oceananigans.jl/src/Fields/field.jl

Lines 695 to 716 in b67ad56

	function Base.$(reduction)(f::Function,
	c::AbstractField;
	condition = nothing,
	mask = get_neutral_mask(Base.$(reduction!)),
	dims = :)
	T = filltype(Base.$(reduction!), c)
	loc = reduced_location(location(c); dims)
	r = Field(loc, c.grid, T; indices=indices(c))
	conditioned_c = condition_operand(f, c, condition, mask)
	initialize_reduced_field!(Base.$(reduction!), identity, r, conditioned_c)
	Base.$(reduction!)(identity, r, conditioned_c, init=false)
	if dims isa Colon
	return CUDA.@allowscalar first(r)
	else
	return r
	end
	end
	Base.$(reduction)(c::AbstractField; kwargs...) = Base.$(reduction)(identity, c; kwargs...)
	end

I guess its more efficient to preallocate and use the preallocating version? We could put a field like:

Oceananigans.jl/src/Fields/field.jl

Line 703 in b67ad56

r = Field(loc, c.grid, T; indices=indices(c))

in the boundary adjacent mean struct.

I guess summing is always inefficient since it has to be serialized can't be straightforwardly parallelised without making a race condition though right?

[glwagner]

Where is update_boundary_condition! defined? In my opinion, update_boundary_condition! is defining a model-specific interface.

[glwagner]

shouldn't you define a kernel for sum? the way implemented seems like its going to be slow

Maybe? Doesn't sum go back to this reduction:

Oceananigans.jl/src/Fields/field.jl

Lines 695 to 716 in b67ad56

	function Base.$(reduction)(f::Function,
	c::AbstractField;
	condition = nothing,
	mask = get_neutral_mask(Base.$(reduction!)),
	dims = :)
	T = filltype(Base.$(reduction!), c)
	loc = reduced_location(location(c); dims)
	r = Field(loc, c.grid, T; indices=indices(c))
	conditioned_c = condition_operand(f, c, condition, mask)
	initialize_reduced_field!(Base.$(reduction!), identity, r, conditioned_c)
	Base.$(reduction!)(identity, r, conditioned_c, init=false)
	if dims isa Colon
	return CUDA.@allowscalar first(r)
	else
	return r
	end
	end
	Base.$(reduction)(c::AbstractField; kwargs...) = Base.$(reduction)(identity, c; kwargs...)
	end

I guess its more efficient to preallocate and use the preallocating version? We could put a field like:

Oceananigans.jl/src/Fields/field.jl

Line 703 in b67ad56

r = Field(loc, c.grid, T; indices=indices(c))

in the boundary adjacent mean struct.
I guess summing is always inefficient since it has to be serialized can't be straightforwardly parallelised without making a race condition though right?

I'm confused. Is this working inside a kernel? You cannot launch a kernel (to compute sum) from within a kernel. Does this work on GPU? Furthermore, @allowscalar should never be used to fetch a generic index (eg i, j, k). The only acceptable usage is to fetch a single known index, typically first or 1, 1, 1.

[jagoosw]

Where is update_boundary_condition! defined? In my opinion, update_boundary_condition! is defining a model-specific interface.

The fallback is defined in boundary conditions, but this is the first time we've actually implemented on in the source code. Happy for it to go in the model modules.

[jagoosw]

shouldn't you define a kernel for sum? the way implemented seems like its going to be slow

Maybe? Doesn't sum go back to this reduction:

Oceananigans.jl/src/Fields/field.jl

Lines 695 to 716 in b67ad56

	function Base.$(reduction)(f::Function,
	c::AbstractField;
	condition = nothing,
	mask = get_neutral_mask(Base.$(reduction!)),
	dims = :)
	T = filltype(Base.$(reduction!), c)
	loc = reduced_location(location(c); dims)
	r = Field(loc, c.grid, T; indices=indices(c))
	conditioned_c = condition_operand(f, c, condition, mask)
	initialize_reduced_field!(Base.$(reduction!), identity, r, conditioned_c)
	Base.$(reduction!)(identity, r, conditioned_c, init=false)
	if dims isa Colon
	return CUDA.@allowscalar first(r)
	else
	return r
	end
	end
	Base.$(reduction)(c::AbstractField; kwargs...) = Base.$(reduction)(identity, c; kwargs...)
	end

I guess its more efficient to preallocate and use the preallocating version? We could put a field like:

Oceananigans.jl/src/Fields/field.jl

Line 703 in b67ad56

r = Field(loc, c.grid, T; indices=indices(c))

in the boundary adjacent mean struct.
I guess summing is always inefficient since it has to be serialized can't be straightforwardly parallelised without making a race condition though right?

I'm confused. Is this working inside a kernel? You cannot launch a kernel (to compute sum) from within a kernel. Does this work on GPU? Furthermore, @allowscalar should never be used to fetch a generic index (eg i, j, k). The only acceptable usage is to fetch a single known index, typically first or 1, 1, 1.

Not from within a kernel, and yeah I used in on GPU. We only want a single value where 2 of i, j, k and k will be 1 but the other will either be 1 or the end. I don't know how else we can pull the single value out?

[glwagner]

Maybe I was confused then. You're saying that i, j, k don't vary and take just a single value per boundary condition?