Fixes bug in FlatExtrapolation matching scheme

[tomchor]

No description provided.

[tomchor]

Should we add a test for this? Or at least make all directions testable on the validation script?

[ali-ramadhan]

Nice find! Fix makes sense to me so I'm approving.

If it's non-trivial to design a good test for these open boundary conditions, then to me it make sense to merge this obvious fix and discuss how to more thoroughly test in a separate issue or PR.

[tomchor]

Nice find! Fix makes sense to me so I'm approving.

If it's non-trivial to design a good test for these open boundary conditions, then to me it make sense to merge this obvious fix and discuss how to more thoroughly test in a separate issue or PR.

I think creating a test is non-trivial, but modifying the existing validation script to be more expansive is trivial since I pretty much have the code ready. Should I just do that?

[glwagner]

Does this fix the problems we see on #3848 ?

[tomchor]

Does this fix the problems we see on #3848 ?

No, this bug makes the code fail with an UndefVarError, but only when FlatExtrapolationOBC is used as a bottom BC.

[tomchor]

I changed oscillating_flow.jl to oscillate in two directions at once. It then oscillates first in the xy direction, then the xz direction, therefore testing all directions in the process. This is the animation it produces for the xy direction:

_cylinder_xy_FlatExtrapolation_Nx_40.jld2.mp4

and this is the animation for the xz direction:

_cylinder_xz_FlatExtrapolation_Nx_40.jld2.mp4

@jagoosw, I'm curious to hear your take on the artifacts that appear at the edges of the right x boundary. If I plot v and w those artifacts are also there (also on the "right" side), so I think this a general "issue" with the algorithm, rather than something wrong with the x direction specifically. Do you have any idea of what this is?

[jagoosw]

@jagoosw, I'm curious to hear your take on the artifacts that appear at the edges of the right x boundary. If I plot v and w those artifacts are also there (also on the "right" side), so I think this a general "issue" with the algorithm, rather than something wrong with the x direction specifically. Do you have any idea of what this is?

I'm not sure I'm seeing what you mean? Do you mean the corner points being different?

[tomchor]

@jagoosw, I'm curious to hear your take on the artifacts that appear at the edges of the right x boundary. If I plot v and w those artifacts are also there (also on the "right" side), so I think this a general "issue" with the algorithm, rather than something wrong with the x direction specifically. Do you have any idea of what this is?

I'm not sure I'm seeing what you mean? Do you mean the corner points being different?

Yes that's what I mean. To be even clearer, these:

The zero points in u being the cause of similar weirdness in the calculations of vorticity. Do you know what's going on there?

[glwagner]

You may need to fill the point Nx + 1 , Ny + 1

[tomchor]

You may need to fill the point Nx + 1 , Ny + 1

If I understand your comment, we are doing that:

Oceananigans.jl/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl

Lines 78 to 92 in 0c1f4fd

	@inline function _fill_east_open_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)
	spacing_factor = Δx₁ / (Δx₂ + Δx₃)
	gradient_free_ϕ = @inbounds ϕ[i - 2, j, k] - (ϕ[i - 1, j, k] - ϕ[i - 3, j, k]) * spacing_factor
	@inbounds ϕ[i, j, k] = relax(j, k, grid, gradient_free_ϕ, bc, clock, model_fields)
	return nothing
	end

[glwagner]

You may need to fill the point Nx + 1 , Ny + 1

If I understand your comment, we are doing that:

Oceananigans.jl/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl

Lines 78 to 92 in 0c1f4fd

	@inline function _fill_east_open_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)
	spacing_factor = Δx₁ / (Δx₂ + Δx₃)
	gradient_free_ϕ = @inbounds ϕ[i - 2, j, k] - (ϕ[i - 1, j, k] - ϕ[i - 3, j, k]) * spacing_factor
	@inbounds ϕ[i, j, k] = relax(j, k, grid, gradient_free_ϕ, bc, clock, model_fields)
	return nothing
	end

Are you sure? What is the range of j (and k)?

[jagoosw]

I'm not sure that the corner points are used in the prognostic calculations right? I guess this would be important if you wanted to do something on the boundary that includes the boundary tangent derivatives, but I haven't come across any schemes that use that

[tomchor]

Are you sure? What is the range of j (and k)?

You're right that j and k (or more generally the directions parallel to the boundary that's being filled) only go up to grid.Ny and grid.Nz. I'm not sure how to change that, though. I'd like to play around with things. Do you mind pointing me in the right direction?

I'm not sure that the corner points are used in the prognostic calculations right? I guess this would be important if you wanted to do something on the boundary that includes the boundary tangent derivatives, but I haven't come across any schemes that use that

You may be right. Still I'd like to investigate. The lack of symmetry between left and right sides is a bit suspicious, no?

[glwagner]

I'm not sure that the corner points are used in the prognostic calculations right?

It depends on the physics. You should check for coriolis, VectorInvariant advection, and biharmonic viscosities. I suspect the corners come into play for those.

Do you mind pointing me in the right direction?

Check out Oceananigans.Utils.KernelParameters and use this when launching the kernel that fills the halos.

[tomchor]

It turns out the open boundary part of the code doesn't use KernelParameters. It instead uses fill_halo_size():

Oceananigans.jl/src/BoundaryConditions/fill_halo_regions_open.jl

Line 18 in d66ed9b

fill_size = fill_halo_size(field, regular_fill, indices, boundary_conditions, loc, grid)

Changing parameters there in order to go from 0 to grid.Ny+1 gets rid of those edges:

_cylinder_xy_FlatExtrapolation_Nx_40.jld2.mp4

I believe this is the preferred way to do things, but it would be nice to hear from other people. If we agree to change things I can open a PR.

[glwagner]

It turns out the open boundary part of the code doesn't use KernelParameters. It instead uses fill_halo_size():

For sure, you just have to change it