Open halo filling

[jagoosw]

Hi all,

I think there's been a mistake in the open boundary filling that's only becoming a problem now that we're trying to fill non-zero value.

Oceananigans.jl/src/BoundaryConditions/fill_halo_regions_open.jl

Lines 86 to 91 in 3ea2545

	@inline _fill_west_open_halo!(j, k, grid, c, bc::OBC, loc, args...) = @inbounds c[1, j, k] = getbc(bc, j, k, grid, args...)
	@inline _fill_east_open_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_open_halo!(i, k, grid, c, bc::OBC, loc, args...) = @inbounds c[i, 1, k] = getbc(bc, i, k, grid, args...)
	@inline _fill_north_open_halo!(i, k, grid, c, bc::OBC, loc, args...) = @inbounds c[i, grid.Ny + 1, k] = getbc(bc, i, k, grid, args...)
	@inline _fill_bottom_open_halo!(i, j, grid, c, bc::OBC, loc, args...) = @inbounds c[i, j, 1] = getbc(bc, i, j, grid, args...)
	@inline _fill_top_open_halo!(i, j, grid, c, bc::OBC, loc, args...) = @inbounds c[i, j, grid.Nz + 1] = getbc(bc, i, j, grid, args...)

The open fill has always set point at index 1 on the right hand side and grid.N+1 on the right hand side, but 1 is part of the prognostic domain and halo points we need are just for computing gradients at the face point, which should be at 0.

I came across this because I've only been testing open boundaries on the right side, but was checking it worked in the other directions and realised it always failed when I just switched the direction and sides for a simple case.

Am I missing something here?

[jagoosw]

For reference, when I fill u[0, j, k] instead of u[1, j, k] I can just reverse the flow in this cylinder case by changing the sign and swapping the boundary conditions round, but before it was failing instantly.

2d_cylinder.mp4

2d_cylinder_reversed.mp4

[glwagner]

For face indexing convention, 1 is the interface forming left boundary of the domain and N+1 is the interface forming the right boundary.

For center indexing, the first cell on the left is 1 and the last cell on the right is N.

I'm not 100% sure what you are asking but this is the definition of the indices.

I think if you set N+1 for right-sided open boundaries, you should set 1 for left-sided open boundaries. If you set N+2 for the right side, then you would set 0 for the left side.

Maybe there is a bug somewhere else?

[jagoosw]

I've worked out where my problem is coming from. For the wall-normal velocity: first, we compute and apply the tendencies from 1:N face points, then compute the pressure correction at 1:N center points, then fill the boundary points at 1 and N+1, and apply it at 1:N face points (except the gradient is zero across the 1 face point so this doesn't do anything to the boundary. The N+1 boundary point is fine because we can just set it to anything, or time integrate something at the point since nothing else effects its value.

The same is true if we prescribe a value at the 1 face point because (even though we redundantly integrate the tendencies there) it just gets reset to whatever we want. The problem is if we try to integrate something like a radiation condition there then we actually end up with $u(1, j, k) = \int (G_u + B_u) dt$ where $B_u$ is whatever integration we're trying to do at the boundary.

On bounded topology I don't think we ever want to integrate the tendency right? But it might be more complicated to do that than to think of a different way todo it.

[glwagner]

That's right. Purely for simplicity we launch all the tendency kernels from 1:N, though for Face-fields in Bounded directions, we only require 2:N.

In fact using tendencies only from 2:N could allow an optimization where we don't need to "enforce" no-penetration boundary conditions. It'd be hard to achieve this optimization though because users can write things like parent(u) .= 1 so I'm not sure we can get away with this being guaranteed correct.

This has never been a problem because we simply overwrite the boundary velocity and therefore simply discard the tendency at index 1.

The problem is if we try to integrate something like a radiation condition

Can you point me to where in the code this goes down?

On bounded topology I don't think we ever want to integrate the tendency right? But it might be more complicated to do that

I think that's right that we don't need the tendency. This has been part of the algorithm since time immemorial and back in the mists of time it was indeed more complicated than worthwhile. The complication is that KernelAbstractions assumes indices start at 1...

However, we now have a way of offsetting indices in kernels via our KernelParameters. So it's not very hard to do this anymore. I can give it a start. If we make this change, we also want to take a step back and look at all the kernels we are launching currently to make sure everything makes sense.

For example, here is a question: while we don't want to integrate the velocity tendencies on boundaries, what do we do about diagnostics? Do we want to compute vorticity on the boundary, for example, if we are computing a vorticity diagnostic? It seems simpler if we don't, that way we don't have special cases...

[glwagner]

I just looked over the source code, and while I think this will be easy for the serial case, the code for distributed models is... fun...

I have wanted to refactor the distributed stuff to make it understandable for a while anyways though. So I think we can come out on top.

[jagoosw]

Can you point me to where in the code this goes down?

It happens in this halo fill:

Oceananigans.jl/src/Models/NonhydrostaticModels/pressure_correction.jl

Lines 8 to 20 in 45838a5

	function calculate_pressure_correction!(model::NonhydrostaticModel, Δt)
	# Mask immersed velocities
	foreach(mask_immersed_field!, model.velocities)
	fill_halo_regions!(model.velocities, model.clock, fields(model))
	solve_for_pressure!(model.pressures.pNHS, model.pressure_solver, Δt, model.velocities)
	fill_halo_regions!(model.pressures.pNHS)
	return nothing
	end

after the tendency integration but before the pressure correction

[jagoosw]

Do you mean when we're computing a diagnostic like vorticity should the kernel include the boundary point? It does seem the simplest would be that anything on a bounded face only launches over 2:N

[glwagner]

Do you mean when we're computing a diagnostic like vorticity should the kernel include the boundary point?

Can you point me to where in the code this goes down?

It happens in this halo fill:

Oceananigans.jl/src/Models/NonhydrostaticModels/pressure_correction.jl

Lines 8 to 20 in 45838a5

	function calculate_pressure_correction!(model::NonhydrostaticModel, Δt)
	# Mask immersed velocities
	foreach(mask_immersed_field!, model.velocities)
	fill_halo_regions!(model.velocities, model.clock, fields(model))
	solve_for_pressure!(model.pressures.pNHS, model.pressure_solver, Δt, model.velocities)
	fill_halo_regions!(model.pressures.pNHS)
	return nothing
	end

after the tendency integration but before the pressure correction

Where is the halo filling code?

[glwagner]

Do you mean when we're computing a diagnostic like vorticity should the kernel include the boundary point? It does seem the simplest would be that anything on a bounded face only launches over 2:N

That's right.

[glwagner]

The halo filling code is unbelievably gnarly

[glwagner]

#3814 is one approach.

Another approach which is a lot simpler actually would be to multiply the tendencies by !peripheral_node(i, j, k, grid, loc...). Let's see how #3814 turns out and we can try that instead if we want.

[jagoosw]

Where is the halo filling code?

Something like:

@inline function _fill_west_open_halo!(j, k, grid, u, bc::PAOBC, loc, clock, model_fields)
    Δt = clock.last_stage_Δt

    Δt = ifelse(isinf(Δt), 0, Δt)

    Δx = xspacing(1, j, k, grid, Face(), Center(), Center())

    C  = getbc(bc, j, k, grid, clock model_fields)

    u₀ⁿ   = @inbounds u[0, j, k]
    u₁ⁿ⁺¹ = @inbounds u[2, j, k]

    C = min(0, max(1, Δt / Δx * C))

    @inbounds u[1, j, k] = (u₀ⁿ - C * u₁ⁿ⁺¹) / (1 - C)
end

This is a simplification (won't work if the flow is inwards) but I'll link when I put it on GH.

[glwagner]

Oh this isn't implemented? How do we test it?

[glwagner]

I think in that case let's wait until this code is implemented, then we'll be able to test that whatever fix we devise is working

[jagoosw]

We only put the no normal gradient matching scheme in the source code which just overwrites the boundary point so this wasn't a problem.

I though we weren't going to put lots of matching schemes in the source code since its not clear what is the best/correct. We could put in a simple scheme that integrates the boundary point if that would be better?

[glwagner]

Why don't we want to put the numerics in the source code? It certainly seems appropriate as a fairly low-level feature.

[glwagner]

I remember discussing a strategy for working on the design of open boundary conditions, and for that I advocated for finding a simple scheme to implement and focusing on the overall design. The purpose of that is to allow us to think clearly and logically about the software design without getting tangled up in numerics.

Once we have a good design (I'm not sure that we do unfortunately...) then the door is open to work on numerics, hopefully without being hindered too much (the point of a good design). Then we can make rapid progress. But this sort of strategy to focus on "one thing at a time" is not a comment about whether we should put numerics in the source code or not. It's a strategy for software development, not a comment about package organization.

[tomchor]

@jagoosw can you share the branch in which you fix this issue? I'm trying to understand it better as I'm seeing some artifacts that may be due to it (see #3854 (comment)) and checking out your fixed code would help.

cc @ali-ramadhan

[jagoosw]

@tomchor, the "fix" I made was just a hack to keep the previous value in the next unused halo point:
https://github.com/jagoosw/OpenBoundaries.jl/blob/e2d041b0418dc015ea6a85a6d408d0c47d483b8e/src/perturbation_advection.jl#L94-L95
But @glwagner has made changes to the kernel launching so it should be possible todo this properly now

I think that this bug is also only relevant for boundary methods that have to timestep the boundary value, so not for the flat extrapolation one etc where it is diagnostically set.