Unable to use continuous function with ValueBoundaryCondition

[bgroenks96]

Hi guys,

Thanks for all of your outstanding work on Oceananigans, really happy to be working with you :)

For unit testing in our land model, I would like to be able to apply a simple periodic upper boundary condition to the temperature field. I figured this would be possible with ValueBoundaryCondition(f) where f is the function. This works fine for FluxBoudnaryCondition, as indicated in the docs. However, with ValueBoundaryCondition, I get an error. Here is the MWE:

using Oceananigans
using Oceananigans.BoundaryConditions: FieldBoundaryConditions, ValueBoundaryCondition, NoFluxBoundaryCondition, fill_halo_regions!

f(x, z, t) = sin(2π*t)
bc = ValueBoundaryCondition(f)
grid = RectilinearGrid(CPU(), size=(1,10), x=(0,1), z=(0,1), topology=(Periodic,Flat,Bounded))
field_bc = FieldBoundaryConditions(grid, (Center, Center, Nothing), top=bc, bottom=NoFluxBoundaryCondition())
field = CenterField(grid, boundary_conditions=field_bc)
fill_halo_regions!(field)

Error:

ERROR: MethodError: objects of type Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Nothing, Nothing, Nothing, Nothing, typeof(f), Nothing, Tuple{}, Nothing, Nothing} are not callable
Stacktrace:
  [1] getbc
    @ ~/.julia/packages/Oceananigans/bKiml/src/BoundaryConditions/boundary_condition.jl:113 [inlined]
  [2] getbc
    @ ~/.julia/packages/Oceananigans/bKiml/src/BoundaryConditions/boundary_condition.jl:112 [inlined]
  [3] right_gradient
    @ ~/.julia/packages/Oceananigans/bKiml/src/BoundaryConditions/fill_halo_regions_value_gradient.jl:13 [inlined]

@glwagner mentioned that regularize_field_boundary_conditions might fix this, but it didn't seem to work when I tried it.

[simone-silvestri]

You cannot fill_halo_regions! of a field when the boundary conditions has possible dependencies outside the field itself (in this case, model time). For function value boundary conditions, you need to provide the dependencies to the fill_halo_regions! function.
For example, if you have a model, the correct way to fill the halo regions is:

fill_halo_regions!(field, model.clock, fields(model))

(In general, the signature should be fill_halo_regions!(field, clock, other_fields))

A couple of syntax notes on the MWE:

• If you want to use a function as a BC for a field that has a Flat topology in one direction (in your case, y), you have to exclude that dimension from the function signature. The correct way to define f would be f(x, t) here.
• Since recently, you need to pass instantiated locations to the FieldBoundaryConditions constructor.

A MWE that would work here is

using Oceananigans
using Oceananigans.BoundaryConditions: FieldBoundaryConditions, ValueBoundaryCondition, NoFluxBoundaryCondition, fill_halo_regions!
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions

grid = RectilinearGrid(CPU(), size=(1,10), x=(0,1), z=(0,1), topology=(Periodic,Flat,Bounded))

f(x, t) = sin(2π*t) + 10
bc = ValueBoundaryCondition(f)
bcs = FieldBoundaryConditions(grid, (Center(), Center(), nothing), top=bc)
bcs = regularize_field_boundary_conditions(bcs, grid, :c)

field = CenterField(grid, boundary_conditions=bcs)
clock = Clock(; time = 10)

fill_halo_regions!(field, clock, field)

@show field.data[1, 1, 11]

[simone-silvestri]

The reason why it works for FluxBoundaryCondition and not for ValueBoundaryCondition is that they are different in how they work.

FluxBoundaryCondition adds a term to the tendency fields during timestepping through the compute_flux_bcs! function.
The fill_halo_regions! function always acts as a FluxBoundaryCondition is a NoFluxBoundaryCondition, so you do not need to pass dependencies to the function, as these are passed to compute_flux_bcs! where the flux is computed. (You'll see that after a fill_halo_regions! of a field with a FluxBoundaryCondition, the interior value is just copied in the halos)

A ValueBoundaryCondition acts directly on the field through fill_halo_regions! by modifying the halos, so you need to pass those extra terms to fill_halo_regions!

We have had a lot of discussions on how to unify the treatment of boundary conditions, but we never got around to doing it.

[bgroenks96]

Thanks Simone! Should we close this issue then?

[simone-silvestri]

Sure, let's close it

[bgroenks96]

Ok two quick quesitons, though:

Does the :c in the call the regularize_field_boundary_conditions have any significance?

And is there any way for me to get rid of the warning about immersed boundary conditions?