Can we impose a surface boundary condition that depends on the surface average of a field?

[tomchor]

Basically I'd like to know if it's currently possible to impose a surface boundary condition for a field (say u) that depends on the surface average of that field (in this case an xy average of u at the point closest to the surface)?

I know that the answer to this question used to be a no (I'm guessing because the kernels only access local information?), but with the recent improvements to reductions, callbacks, etc., I figured there could be a way to achieve that?

That said, every attempt that I've tried today failed completely... heh

If it's not indeed possible yet, maybe we can discuss ideas of how to implement that? The easiest to me seems to be changing the boundary condition of the given field with a callback at every iteration, but I get an error every time I try that:

ERROR: MethodError: Cannot `convert` an object of type 
  BoundaryCondition{Oceananigans.BoundaryConditions.Flux{},Float64} to an object of type 
  BoundaryCondition{Oceananigans.BoundaryConditions.Flux{},Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Center, Center, Nothing, 64, typeof(buoyancy_flux), NamedTuple{(:initial_buoyancy_flux, :shut_off_time), Tuple{Float64, Float64}}, Tuple{Symbol}, Tuple{Int64}, Tuple{typeof(Oceananigans.Operators.ℑxᶜᵃᵃ)}}}
Closest candidates are:
  convert(::Type{T}, ::T) where T at essentials.jl:205
  BoundaryCondition{C, T}(::Any, ::Any) where {C<:Oceananigans.BoundaryConditions.AbstractBoundaryConditionClassification, T} at /home/tomas/repos/Oceananigans.jl/src/BoundaryConditions/boundary_condition.jl:9
Stacktrace:
 [1] setproperty!(x::FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Gradient, Float64}, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Center, Center, Nothing, 64, typeof(buoyancy_flux), NamedTuple{(:initial_buoyancy_flux, :shut_off_time), Tuple{Float64, Float64}}, Tuple{Symbol}, Tuple{Int64}, Tuple{typeof(Oceananigans.Operators.ℑxᶜᵃᵃ)}}}, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, f::Symbol, v::BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Float64})
   @ Base ./Base.jl:34
 [2] top-level scope
   @ REPL[49]:1

so my guess is that currently the BCs are implemented with immutable structs or something along those lines...

[glwagner]

There are two primary ways I can think of (but probably many more):

1. Use callbacks:

• Build the boundary condition with an array (for example u_top_bc = FluxBoundaryCondition(zeros(grid, Nx, Ny))).
• Compute model.velocities.u.boundary_conditions.top.condition (for example) at every time-step with a function called update_boundary_condition! and simulation.callbacks[:bc] = Callback(update_boundary_condition!).
• Note: if you are using RungeKutta3TimeStepper, this will not update the boundary condition during time-stepper substages (only at the end of every time-step). Usually this is "ok" (maybe not if your boundary conditions fluctuates very fast, nearly as fast as a time-step). Some source code changes would be required to have callbacks within time-stepper stages.

2. Pre-build velocities = VelocityFields(grid).

• Then define, eg U = Field(Average(velocities.u, dims=(1, 2)))
• Then build the boundary conditions in terms of U, eg

u_flux(x, y, t, U) = # something
u_top_bc = FluxBoundaryCondition(u_flux, parameters=U)

• Then add U an an auxiliary field, eg model = NonhydrostaticModel(; grid, velocities, auxiliary_fields = (; U)).
• Because auxiliary fields are computed every time-stepper stage, this will indeed update the boundary condition every time-stepper stage:

Oceananigans.jl/src/Models/NonhydrostaticModels/update_nonhydrostatic_model_state.jl

Lines 24 to 27 in c74dfc6

	# Compute auxiliary fields
	for aux_field in model.auxiliary_fields
	compute!(aux_field)
	end

[tomchor]

Thanks for the help @glwagner.

I think I got the option 1 working, but I can't make option 2 work. I'm interested in making that one work to explore the auxiliary_fields functionality (is that documented anywhere? That's pretty useful) and because it seems to be easier to adapt for GPUs.

If I understand correctly, this is what the second option would look like (notice that this example makes no physical sense; it's just something numerically simple that tests the functionality):

using Oceananigans
using Oceananigans.Fields: VelocityFields, TracerFields

Nx = 64
Nz = 32

grid = RectilinearGrid(topology=(Periodic, Flat, Bounded),
                       size=(Nx, Nz),
                       x=(0, 1), z=(0, 1),
                       halo=(3,3),
                       )

velocities = VelocityFields(grid)
U = Field(Average(velocities.u, dims=(1, 2)), indices=(:,:,Nz:Nz)) # Trying to get the surface here but for some reason ` indices=(:,:,Nz)` doesn't work


u_flux(x, y, t, U) = U
u_top_bc = FluxBoundaryCondition(u_flux, parameters=U)
u_bc = FieldBoundaryConditions(top=u_top_bc)
                           
model = NonhydrostaticModel(; grid,
                            velocities=velocities,
                            auxiliary_fields = (; U),
                            boundary_conditions=(u=u_bc,),
                           )
u_ic(x, y, z) = randn()
set!(model, u=u_ic)

u, v, w = model.velocities

using Oceananigans.Grids: min_Δz
simulation = Simulation(model, Δt=0.1min_Δz(grid)/maximum(u), stop_time=10)

run!(simulation)

But any time I try to run it I get a segfault error which doesn't give me much to go from (line 21 is the call to NonhydrostaticModel):

julia> include("op2.jl")
[ Info: Oceananigans will use 4 threads

signal (11): Segmentation fault
in expression starting at /home/tomas/Dropbox/tests/bcdep/op2.jl:21
zsh: segmentation fault (core dumped)  ~/bin/julia-1.6.3/bin/julia --project -t 4

Any ideas about what's causing this?

[glwagner]

This line:

u_flux(x, y, t, U) = U

tries to set every x, y location to the field U, which fails. Instead you need to use the discrete_form:

@inbounds u_flux(i, j, grid, clock, fields, U) = @inbounds U[i, j, 1]

It doesn't really make sense to set the flux of "u" to the average U does it? What is the real example being used here? Is there a computation that must be done on U?

[tomchor]

Oh yeah, this makes no physical sense. It's just the simplest computational form to check.

I guess I should make that clear with the code. I'll make an edit

Thanks for the suggestion though. It makes sense.

[tomchor]

Just for reference this is how I ended up implementing a MWE of an adaptive top BC based on @glwagner's option 1 (with callbacks):

using Oceananigans
using Oceananigans.Fields: VelocityFields, TracerFields

Nx = 64
Nz = 32

grid = RectilinearGrid(topology=(Periodic, Flat, Bounded),
                       size=(Nx, Nz),
                       x=(0, 1), z=(0, 1),
                       halo=(3,3),
                       )

u_top_bc = FluxBoundaryCondition(zeros(grid, grid.Nx, grid.Ny))
u_bc = FieldBoundaryConditions(top=u_top_bc)

model = NonhydrostaticModel(; grid,
                            boundary_conditions=(u=u_bc,),
                            closure = ScalarDiffusivity(ν=1e-3),
                           )
u_ic(x, y, z) = randn()
set!(model, u=u_ic)

u, v, w = model.velocities

using Oceananigans.Grids: min_Δz
simulation = Simulation(model, Δt=0.1min_Δz(grid)/maximum(u), stop_time=10)

using Statistics: mean
function update_boundary_condition!(simulation)
    u, v, w = simulation.model.velocities
    U = mean(interior(u)[:,:,grid.Nz]) # this and the following line probably require adaptation to run on GPUs
    model.velocities.u.boundary_conditions.top.condition .= U
    return nothing
end
simulation.callbacks[:bc] = Callback(update_boundary_condition!)


wizard = TimeStepWizard(max_change=1.1, cfl=0.7)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(2))


output_interval = 5
fields_to_output = merge(model.velocities, model.tracers,)
simulation.output_writers[:fields] = NetCDFOutputWriter(model, fields_to_output,
                     schedule = TimeInterval(output_interval),
                     filepath = "op1.nc",
                     mode = "c")

run!(simulation)

This turned out to be pretty close to what I had originally tried, with the exception that I was trying to use callbacks to change model.velocities.u.boundary_conditions.top directly, which doesn't work. The trick is to change model.velocities.u.boundary_conditions.condition, which is an array.

Thanks again to @glwagner for the help.

[tomchor]

So just to make sure I understand you're saying this is easier for GPUs:

    u, v, w = simulation.model.velocities
    U = Field(Average(u, dims=(1,2)))
    model.velocities.u.boundary_conditions.top.condition .= interior(U, :,:,1)
    return nothing
end

(I've just tested this on a CPU and I think it works, but haven't done so for a GPU.)

[glwagner]

Close. I would pre-allocate U:

U = Field(Average(u, dims=(1, 2)))

function update_boundary_condition!(sim)
    compute!(U)
    sim.model.velocities.u.boundary_conditions.top.condition .= interior(U, :, :, 1)
    return nothing
end

Also I tested it and you do need to drop the trailing z dimension, unless you build the boundary condition array as zeros(Nx, Ny, 1).