CliMA / Oceananigans.jl

🌊 Julia software for fast, friendly, flexible, ocean-flavored fluid dynamics on CPUs and GPUs
https://clima.github.io/OceananigansDocumentation/stable
MIT License
1k stars 196 forks source link

Open halo filling #3810

Open jagoosw opened 1 month ago

jagoosw commented 1 month ago

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.

https://github.com/CliMA/Oceananigans.jl/blob/3ea2545331d9910d8b467dd8eb31074fb426af5b/src/BoundaryConditions/fill_halo_regions_open.jl#L86-L91

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 commented 1 month ago

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.

https://github.com/user-attachments/assets/b711f8aa-1848-4e9a-9546-3567f4816a1c

https://github.com/user-attachments/assets/04f8e2bb-ecd7-4af4-8f47-6633e099721e

glwagner commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

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

It happens in this halo fill:

https://github.com/CliMA/Oceananigans.jl/blob/45838a57dd5ebc1153c2c827f83cb848d20e4c92/src/Models/NonhydrostaticModels/pressure_correction.jl#L8-L20

after the tendency integration but before the pressure correction

jagoosw commented 1 month ago

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 commented 1 month ago

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:

https://github.com/CliMA/Oceananigans.jl/blob/45838a57dd5ebc1153c2c827f83cb848d20e4c92/src/Models/NonhydrostaticModels/pressure_correction.jl#L8-L20

after the tendency integration but before the pressure correction

Where is the halo filling code?

glwagner commented 1 month ago

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 commented 1 month ago

The halo filling code is unbelievably gnarly

glwagner commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

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

glwagner commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

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

glwagner commented 1 month ago

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 commented 1 month ago

@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 https://github.com/CliMA/Oceananigans.jl/pull/3854#issuecomment-2435113186) and checking out your fixed code would help.

cc @ali-ramadhan

jagoosw commented 1 month ago

@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.