Error when using `OpenBoundaryCondition` with `HydrostaticFreeSurfaceModel`

[tomchor]

When I use a function to set an OpenBoundaryCondition on a HydrostaticFreeSurfaceModel I get an error and I don't understand why. Here's a MWE:

using Oceananigans
grid = RectilinearGrid(topology = (Bounded, Flat, Bounded), size = (4, 4), extent = (1, 1))

u₀ = 1
@inline u_func(z, t) = u₀

u_bcs = FieldBoundaryConditions(east = OpenBoundaryCondition(u_func), west = OpenBoundaryCondition(u_func))

model = HydrostaticFreeSurfaceModel(; grid, boundary_conditions = (; u = u_bcs,))
set!(model, u = u₀)
time_step!(model, 0.1)

This gives me the error:

ERROR: LoadError: TaskFailedException

    nested task error: MethodError: objects of type Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Nothing, Center, Center, Oceananigans.BoundaryConditions.LeftBoundary, typeof(u_func), Nothing, Tuple{}, Tuple{}, Tuple{}} are not callable
    Stacktrace:
     [1] getbc
       @ ~/.julia/packages/Oceananigans/OHYQj/src/BoundaryConditions/boundary_condition.jl:115 [inlined]
     [2] _fill_west_halo!
       @ ~/.julia/packages/Oceananigans/OHYQj/src/BoundaryConditions/fill_halo_regions_open.jl:34 [inlined]
     [3] #25
       @ ~/.julia/packages/Oceananigans/OHYQj/src/BoundaryConditions/fill_halo_regions.jl:260 [inlined]
     [4] ntuple
       @ ./ntuple.jl:50 [inlined]
     [5] cpu__fill_west_and_east_halo!
       @ ~/.julia/packages/KernelAbstractions/HAcqg/src/macros.jl:287 [inlined]
     [6] __thread_run(tid::Int64, len::Int64, rem::Int64, obj::KernelAbstractions.Kernel{…}, ndrange::Nothing, iterspace::KernelAbstractions.NDIteration.NDRange{…}, args::Tuple{…}, dynamic::KernelAbstractions.NDIteration.DynamicCheck)
       @ KernelAbstractions ~/.julia/packages/KernelAbstractions/HAcqg/src/cpu.jl:115
     [7] (::KernelAbstractions.var"#18#21"{…})()
       @ KernelAbstractions ~/.julia/packages/KernelAbstractions/HAcqg/src/cpu.jl:90

which points to this line https://github.com/CliMA/Oceananigans.jl/blob/d4bcc095be66c7b5c98a462106285a6f6d341fe1/src/BoundaryConditions/boundary_condition.jl#L115 which is a fallback method.

If I set-up exactly the same code but using a NonhydrostaticModel things work as expected. Things also work if I set the OpenBoundaryCondition using a constant, such as:

using Oceananigans
grid = RectilinearGrid(topology = (Bounded, Flat, Bounded), size = (4, 4), extent = (1, 1))
u₀ = 1
u_bcs = FieldBoundaryConditions(east = OpenBoundaryCondition(u₀), west = OpenBoundaryCondition(u₀))
model = HydrostaticFreeSurfaceModel(; grid, boundary_conditions = (; u = u_bcs,))
set!(model, u = u₀)
time_step!(model, 0.1)

Am I missing something obvious here, or is something wrong? Also I couldn't find any scripts in the repo that exemplify OpenBoundaryConditions, so I'm happy to set-up a validation script or something once this is figured out if you guys think it's a good addition.

[glwagner]

Non-zero OpenBoundaryConditions probably doesn't work in practice either even if it can be run without an error... @jagoosw is working on developing numerical methods for open boundaries for the nonhydrostatic model. But nobody has shown interest in developing open boundaries for the hydrostatic model yet and I don't think that will work without some modification to the treatment of the free surface.

[glwagner]

@simone-silvestri

[simone-silvestri]

As @glwagner pointed out, the code is not set up to support this feature.

ContinuousBoundaryFunctions and DiscreteBoundaryFunctions are supported only for Flux boundary conditions at the moment, where the boundary contribution is added as an additional tendency rather than filling the halo regions.

If we want to support these more complicated BC for Open, Value and Gradient, we have to make sure that in all places where we use the fill_halo_regions! we also pass as arguments the clock and all the model_fields, which might be doable.

In this case you are hitting this problem in the hydrostatic model because in the ImplicitFreeSurfaceSolver there is a fill_halo_regions! for velocities that does not pass clock and model_fields, making julia complain

https://github.com/CliMA/Oceananigans.jl/blob/d4bcc095be66c7b5c98a462106285a6f6d341fe1/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl#L137

If you use the SplitExplicit the error disappears. This is not to say that the solution will be correct because still the barotropic velocities and barotropic volume fluxes will not account for those open boundary conditions. Also for Nonhydrostatic, even if the script does not error, there might be some subtleties for which the result is incorrect. The correct implementation of open boundary conditions should be in #3842

[tomchor]

Thanks for the clarification. I was under the impression that, even though open boundaries weren't implemented in the sense of radiative/Orlansky/Flather boundary conditions, they still worked in the sense of imposing an inflow/outflow at every boundary point (like it's demonstrated in the first animation here).

This also clarifies for me some of the discussion I've been seeing in #3482.

[glwagner]

Ok but @simone-silvestri that problem with the implicit solver is indeed a bug. It's trivially fixed too.

[glwagner]

Thanks for the clarification. I was under the impression that, even though open boundaries weren't implemented in the sense of radiative/Orlansky/Flather boundary conditions, they still worked in the sense of imposing an inflow/outflow at every boundary point (like it's demonstrated in the first animation here).

This also clarifies for me some of the discussion I've been seeing in #3482.

You may be able to impose inflow and outflow, it's just that any non-trivial simulation will not succeed without further work

[glwagner]

It's fine to keep this open...

[simone-silvestri]

Ok but @simone-silvestri that problem with the implicit solver is indeed a bug. It's trivially fixed too.

I wouldn't consider this a bug but rather a non-implemented feature. It is not the only place where this pattern appears in the code, in a lot of places we write

fill_halo_regions!(field)

without passing clock and model_fields. If we want to make sure that we support complex BCs also for Value, Open, and Gradient we need to ensure that this pattern changes everywhere in the code. Indeed, it is not difficult to do it, @tomchor if you want to can give it a try

[tomchor]

Yeah, sure I can try.

To get a random line as an example, would it be just re-writing

https://github.com/CliMA/Oceananigans.jl/blob/d4bcc095be66c7b5c98a462106285a6f6d341fe1/src/Models/HydrostaticFreeSurfaceModels/update_hydrostatic_free_surface_model_state.jl#L40

as

fill_halo_regions!(model.diffusivity_fields, model.clock, fields(model); only_local_halos = true) 

?

And so on for other instances throughout the code?

[simone-silvestri]

yeah, the inside workings of the code are quite easy, the trickier bits will be coming up with a system that works for the diagnostics and the API.

for example this line https://github.com/CliMA/Oceananigans.jl/blob/d4bcc095be66c7b5c98a462106285a6f6d341fe1/src/AbstractOperations/computed_field.jl#L72

or in general, coming up with a system that makes sure that the API is compatible with the choice of complex BCs For example, in a script people should not be allowed to do

u, v, w = model.velocities
fill_halo_regions!((u, v, w))

[tomchor]

or in general, coming up with a system that makes sure that the API is compatible with the choice of complex BCs For example, in a script people should not be allowed to do

u, v, w = model.velocities
fill_halo_regions!((u, v, w))

I understand the first part, but can you explain why the latter isn't desirable and why it would be hard to come up with something? The way I see it, if we set in stone that fill_halo_regions!() needs to have clock and fields(model) passed, then can't we just remove the method that would make fill_halo_regions!((u, v, w)) possible?

[simone-silvestri]

Sorry, let me rephrase, users shouId pass clock and fields only when boundary conditions require it, so we need some type of warning to make sure the error is well documented.

I think it is nice to have fill_halo_regions!(field), it will be a bit cumbersome to always have to pass a clock and fields also when non using complex BCs

[glwagner]

Yeah, sure I can try.

To get a random line as an example, would it be just re-writing

https://github.com/CliMA/Oceananigans.jl/blob/d4bcc095be66c7b5c98a462106285a6f6d341fe1/src/Models/HydrostaticFreeSurfaceModels/update_hydrostatic_free_surface_model_state.jl#L40

as

fill_halo_regions!(model.diffusivity_fields, model.clock, fields(model); only_local_halos = true) 

?

And so on for other instances throughout the code?

No no, we don't want to do this everywhere. Just for the prognostic fields. I was really referring just to the one line that @simone-silvestri pointed out.

[glwagner]

or in general, coming up with a system that makes sure that the API is compatible with the choice of complex BCs For example, in a script people should not be allowed to do

u, v, w = model.velocities
fill_halo_regions!((u, v, w))

I understand the first part, but can you explain why the latter isn't desirable and why it would be hard to come up with something? The way I see it, if we set in stone that fill_halo_regions!() needs to have clock and fields(model) passed, then can't we just remove the method that would make fill_halo_regions!((u, v, w)) possible?

No, we definitely don't want that. That would make it impossible to use Fields without a model. We want simple / default boundary conditions to work without an entire model state.

One has to understand the additional arguments to fill_halo_regions! as part of a system to support additional features for fields, such as boundary conditions that depend on other parts of the model state. There's a circular depedancy issue: when field boundary conditions depend on other fields, the entire system of fields must be created simultaneously for the computation to work properly. We don't have a general system for this. We just support it within the context of a model, like the hydrostatic or nonhydrostatic model. This works because whenever we fill halo regions within the model time-stepping algorithm, we can indeed incorporate the whole model state.

[simone-silvestri]

Then a better course of action would be to make Value and Gradient operate as Flux so that it would still be possible to write

fill_halo_regions!(u)

and use compute! without having to pass additional arguments

[glwagner]

Sorry, let me rephrase, users shouId pass clock and fields only when boundary conditions require it, so we need some type of warning to make sure the error is well documented.

I think it is nice to have fill_halo_regions!(field), it will be a bit cumbersome to always have to pass a clock and fields also when non using complex BCs

It's not even possible to construct a field with ContinuousForcing outside a model constructor. The model constructor has to do things like compute the index of fields(model) that a field name belongs to. So its more than cumbersome, its not possible. Also we need to prioritize. We are trying to do simulations, not have uber fancy flexible field abstractions (this is just an auxiliary feature that's very nice).

[glwagner]

Then a better course of action would be to make Value and Gradient operate as Flux so that it would still be possible to write

fill_halo_regions!(u)

and use compute! without having to pass additional arguments

It's not Value and Gradient that's the problem, its function boundary conditions that don't work. When you use a function you expect a signature. The additional arguments to the discrete form boundary condition are exactly the arguments that are passed to fill_halo_regions. You can actually use any function you want in your own field system if you pass the arguments to fill_halo_regions.

[glwagner]

In other words if you try to use a boundary condition (discrete form) that's

bc(i, j, grid, a, b, c, d, e)

then you need to call fill halo regions:

fill_halo_regions!(field, a, b, c, d, e)

[glwagner]

And continuous form doesn't work at all, and it can't. You would have to build the continuous form manually, I don't think there's a point.

[glwagner]

Ok but @simone-silvestri that problem with the implicit solver is indeed a bug. It's trivially fixed too.

I wouldn't consider this a bug but rather a non-implemented feature. It is not the only place where this pattern appears in the code, in a lot of places we write

fill_halo_regions!(field)

without passing clock and model_fields. If we want to make sure that we support complex BCs also for Value, Open, and Gradient we need to ensure that this pattern changes everywhere in the code. Indeed, it is not difficult to do it, @tomchor if you want to can give it a try

Doesn't it also mean you can't use function Flux top and bottom boundary conditions with the implicit solver? Like quadratic bottom drag?

[simone-silvestri]

if you do it like that (allow a function bc to be applied by filling the halos rather than as an additional tendency), then a field initialized in a model with a function bc (like a velocity for example) will always expect additional signatures for fill_halo_regions! and computing within an abstract operation would be impossible if not passing the additional arguments.

I am suggesting to treat all complex BCs like we treat Flux so that the API will not need drastic changes but we can incorporate the feature that @tomchor is talking about.

[glwagner]

if you do it like that (allow a function bc to be applied by filling the halos rather than as an additional tendency), then a field initialized in a model with a function bc (like a velocity for example) will always expect additional signatures for fill_halo_regions! and computing within an abstract operation would be impossible if not passing the additional arguments.

I am suggesting to treat all complex BCs like we treat Flux so that the API will not need drastic changes but we can incorporate the feature that @tomchor is talking about.

No that's not true --- that's how it works now. The arguments are propagated via args.... If the bc is not a function, the args are thrown away:

https://github.com/CliMA/Oceananigans.jl/blob/d4bcc095be66c7b5c98a462106285a6f6d341fe1/src/BoundaryConditions/boundary_condition.jl#L121-L122

if you use a function then you get:

https://github.com/CliMA/Oceananigans.jl/blob/d4bcc095be66c7b5c98a462106285a6f6d341fe1/src/BoundaryConditions/boundary_condition.jl#L115

it's actually fairly simple to understand

[simone-silvestri]

Flux is added as an a

Ok but @simone-silvestri that problem with the implicit solver is indeed a bug. It's trivially fixed too.

I wouldn't consider this a bug but rather a non-implemented feature. It is not the only place where this pattern appears in the code, in a lot of places we write

fill_halo_regions!(field)

without passing clock and model_fields. If we want to make sure that we support complex BCs also for Value, Open, and Gradient we need to ensure that this pattern changes everywhere in the code. Indeed, it is not difficult to do it, @tomchor if you want to can give it a try

Doesn't it also mean you can't use function Flux top and bottom boundary conditions with the implicit solver? Like quadratic bottom drag?

Flux is added as a tendency, not by filling halos. I have to think about it a sec but the drag might be already included in the RHS

[glwagner]

So if you call fill_halo_regions!(field) then args = (i, j, grid) are thrown away.

If you call fill_halo_regions!(field, a, b, c) then args = (i, j, grid, a, b, c).

I guess the code actually breaks if you only have one argument because

https://github.com/CliMA/Oceananigans.jl/blob/d4bcc095be66c7b5c98a462106285a6f6d341fe1/src/BoundaryConditions/discrete_boundary_function.jl#L38

but that is easily changed if we wanted to (the code is written that way merely for clarity, it has nothing to do with function)

[glwagner]

Ok but if we use something like ValueBoundaryCondition(0) then this will work with fill_halo_regions!(field).

I don't really understand the point of not passing arguments to fill_halo_regions! for the implicit time stepper. It really is a bug.

[glwagner]

Flux is added as an a

Ok but @simone-silvestri that problem with the implicit solver is indeed a bug. It's trivially fixed too.

I wouldn't consider this a bug but rather a non-implemented feature. It is not the only place where this pattern appears in the code, in a lot of places we write

fill_halo_regions!(field)

without passing clock and model_fields. If we want to make sure that we support complex BCs also for Value, Open, and Gradient we need to ensure that this pattern changes everywhere in the code. Indeed, it is not difficult to do it, @tomchor if you want to can give it a try

Doesn't it also mean you can't use function Flux top and bottom boundary conditions with the implicit solver? Like quadratic bottom drag?

Flux is added as a tendency, not by filling halos. I have to think about it a sec but the drag might be already included in the RHS

Right ok, we can use Flux. Also constant or array Value and Gradient will work without an error.

Function Value, Gradient, Open will fail.

Functions could be implemented by updating the array in a callback. This can also be used for any example. The functions are just a convenience. But we should support it where its easy and obvious, don't you think?

[glwagner]

For example we fill halos correctly during the pressure correction for the nonhydrostatic model:

https://github.com/CliMA/Oceananigans.jl/blob/d4bcc095be66c7b5c98a462106285a6f6d341fe1/src/Models/NonhydrostaticModels/pressure_correction.jl#L13

We just aren't paying as close attention to the hydrostatic model.

[glwagner]

A small step: https://github.com/CliMA/Oceananigans.jl/pull/3629

[simone-silvestri]

can I add some changes in other fill_halo_regions! ?

[glwagner]

You are free to open as many PRs as you like. As long as they are small ;-)

[glwagner]

Here's another for the discrete form boundary condition issue https://github.com/CliMA/Oceananigans.jl/pull/3630

[glwagner]

can I add some changes in other fill_halo_regions! ?

what do you propose.

[simone-silvestri]

Actually no I thought I saw other instances of the same problem somewhere, but the only issue now is the split explicit free surface that does not allow for complex (value, open gradient) BCs in the barotropic velocities and the free surface. Probably we shouldn't worry about that until we have open BCs for momentum in the hydrostatic model though

[glwagner]

Hmm ok so the boundary conditions have to be "slow", is that right?

[simone-silvestri]

yeah, if they are slow, the horizontal BCs for the barotropic velocities are included in the vertically integrated tendencies. Not sure about eta

[ali-ramadhan]

@tomchor's MWE from his original post works for me (as of Oceananigans.jl v0.91.11) so there's no error anymore, although open BCs for hydrostatic models aren't implemented. Is it worth closing this issue?

If you use the SplitExplicit the error disappears. This is not to say that the solution will be correct because still the barotropic velocities and barotropic volume fluxes will not account for those open boundary conditions.

Originally posted by @simone-silvestri in https://github.com/CliMA/Oceananigans.jl/issues/3628#issuecomment-2178789178

Now that PR #3482 has been merged, is this still the thing missing for proper hydrostatic open BCs? I'm guessing that the implementation is different for each free surface solver?

• The implicit free surface solver needs to account for a non-zero surface volume flux $M$? See equation (9) of the Implicit free surface pressure operator docs.
• The split-explicit free surface solver's evolution kernels has difference functions that need to be modified to support open boundary conditions? See this comment.

If it makes sense, I can open a new issue to discuss hydrostatic open BCs.

[simone-silvestri]

Regarding the split explicit solver, indeed, we need to include the boundary conditions in the operators for both the velocities and the free surface, which should ensure volume conservation. These operators are in the split_explicit_free_surface_kernels.jl file but they should not be specific to the free surface. There was an effort to move these operators in the Operators module (see #3268) probably would be good to revive that PR before embarking in an open boundary condition effort.

For the implicit free surface I am not sure how you would include the boundary conditions but I would guess they would need to be included in the RHS of the linear operation.

[ali-ramadhan]

Thanks for confirming @simone-silvestri! Sounds like it's easier to add support for open BCs with the split-explicit solver. Or at least, it's the preferred solver so it would be nice to add support with split-explicit rather than the implicit solver.

I can have a look at PR #3268 to see what needs to be done to move the operators. It doesn't look too stale!

[simone-silvestri]

We probably need to pass the boundary conditions inside the operators. We also need to automatically construct the open boundary conditions for the barotropic velocity and the free surface given the baroclinic u and v boundary conditions. This regularization could go in here https://github.com/CliMA/Oceananigans.jl/blob/b3be950ef6a1f2139c2f65e083f20f1f7958ea55/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface_kernels.jl#L294-L297

where we initialize the free surface values give the initial condition of velocity (it is called only once when initializing the simulation).

[simone-silvestri]

Also, for open boundary conditions, we probably want to have a z-star coordinate and a non-linear free surface to make sure we conserve volume (@jm-c could advise). There is one open PR for ZStar if you want to take a look and test it out (#3411), it is not stale because I am keeping it updated (and should be working), but it's still in a testing phase, and it is still not correct for immersed boundaries.

[ali-ramadhan]

We probably need to pass the boundary conditions inside the operators. We also need to automatically construct the open boundary conditions for the barotropic velocity and the free surface given the baroclinic u and v boundary conditions. This regularization could go in here

https://github.com/CliMA/Oceananigans.jl/blob/b3be950ef6a1f2139c2f65e083f20f1f7958ea55/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface_kernels.jl#L294-L297

where we initialize the free surface values give the initial condition of velocity (it is called only once when initializing the simulation).

Makes sense! Are you saying that the initial free surface height needs to be computed from the initial u and v? But only for open boundary conditions? Or is it just an initial elliptic solve for η to satisfy continuity?

Also would the boundary condition on the barotropic velocities U and V just be the vertical integral of the boundary conditions on baroclinic u and v?

[simone-silvestri]

Indeed, the boundary conditions for U and V should be the integral of boundaries for u and v. I am unsure about the boundaries for η; I think η is typically prescribed outside the domain in regional simulations, but we probably need to find a reference for that.

[tomchor]

Indeed, the boundary conditions for U and V should be the integral of boundaries for u and v. I am unsure about the boundaries for η; I think η is typically prescribed outside the domain in regional simulations, but we probably need to find a reference for that.

For reference, I think ROMS currently uses the algorithm described here. I haven't read it carefully, but it seems to project both U and η to a common point outside of the domain, calculate a variable that depends on both quantities (Eq. (2)), and specify the BCs for both based on that. (Likewise for V.)

[ali-ramadhan]

I hope it's okay that I opened a new issue to discuss tackling open BC support (see #3828) since this issue was originally focused on the @tomchor's MWE erroring, which does not error anymore. I also think the discussion has gotten intermingled with discussions from PR #3268.

If it's okay with everyone maybe we can close this issue and continue discussing in #3828?

[glwagner]

I'm ok with that