Can we impose boundary conditions on functions?

[francispoulin]

Can we define a function and impose say Value (Dirichlet) or Derivative (Neumann) boundary conditions on the function?

To be more precise, suppose we have a 1D domain and that A is a tracer defined at a Center and we want to define dA to be it's derivative on Face.

We can define the function dA with the following:

A = model.tracers.A
dA = $\partial_x (A) $

This will allow us to compute dA on the two boundaries but can we then fill in the first halo regions based on whether we want one of the two above mentioned boundary conditions?

[simone-silvestri]

you could do something like this:

julia> grid = RectilinearGrid(size = (Nx, Ny),
                                   x=(-Lx, Lx), y=(-Lx, Ly),
                                   topology = (Bounded, Bounded, Flat))
10×10×1 RectilinearGrid{Float64, Bounded, Bounded, Flat} on CPU with 3×3×0 halo
├── Bounded  x ∈ [-1.0, 1.0]      regularly spaced with Δx=0.2
├── Bounded  y ∈ [-1.0, 1.0]      regularly spaced with Δy=0.2
└── Flat z

julia> A = CenterField(grid)
10×10×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 10×10×1 RectilinearGrid{Float64, Bounded, Bounded, Flat} on CPU with 3×3×0 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: ZeroFlux, east: ZeroFlux, south: ZeroFlux, north: ZeroFlux, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 16×16×1 OffsetArray(::Array{Float64, 3}, -2:13, -2:13, 1:1) with eltype Float64 with indices -2:13×-2:13×1:1
    └── max=0.0, min=0.0, mean=0.0

julia> boundary_conditions = FieldBoundaryConditions(north = FluxBoundaryCondition(1), south = ValueBoundaryCondition(2), west = nothing, east = nothing, top = nothing, bottom = nothing)

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: Nothing
├── east: Nothing
├── south: ValueBoundaryCondition: 2
├── north: FluxBoundaryCondition: 1
├── bottom: Nothing
├── top: Nothing
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

julia> dA = Field(∂x(A); boundary_conditions)
11×10×1 Field{Face, Center, Center} on RectilinearGrid on CPU
├── grid: 10×10×1 RectilinearGrid{Float64, Bounded, Bounded, Flat} on CPU with 3×3×0 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Value, north: Flux, bottom: Nothing, top: Nothing, immersed: Default
├── operand: Derivative at (Face, Center, Center)
├── status: time=0.0
└── data: 17×16×1 OffsetArray(::Array{Float64, 3}, -2:14, -2:13, 1:1) with eltype Float64 with indices -2:14×-2:13×1:1
    └── max=0.0, min=0.0, mean=0.0

we have problems with fields defined on Faces though (in this case if you want to have an east or west boundary condition)

[francispoulin]

Thank you @simone-silvestri ! That is very helpful.