`Solvers.initialize_matrix` does not infer boundary conditions from `template_field`

[navidcy]

I noticed that similar(::Field) does not transfer boundary conditions

https://github.com/CliMA/Oceananigans.jl/blob/a86ef32581192d82a655e1b2ff4411e79917a379/src/Fields/field.jl#L217-L227

so we might need to update

https://github.com/CliMA/Oceananigans.jl/blob/c8a65a4fdeff25722104ced3a0c74d2d921ae1cf/src/Solvers/multigrid_solver.jl#L225-L244

and in particular

https://github.com/CliMA/Oceananigans.jl/blob/c8a65a4fdeff25722104ced3a0c74d2d921ae1cf/src/Solvers/multigrid_solver.jl#L230

[navidcy]

Does anybody know what is the rationale for similar(::Field) to drop the boundary conditions?

[simone-silvestri]

Hmmm, not sure about it. Maybe you don't want to transfer field-specific complex boundary functions? (for example drag or restoration BCs). In the end the "default" BCs (Periodic, No-flux, Non-penetration, Communication) should be transferred

[glwagner]

Does anybody know what is the rationale for similar(::Field) to drop the boundary conditions?

I think it depends on how you view "boundary conditions". If we adopt a "weak formulation philosophy" then boundary conditions are part of the definition of an equation set and don't belong to fields at all. This is how ImmersedBoundaryGrid works; inhomogeneous boundary conditions can only be enforced by adding boundary fluxes to a tendency. There's no such thing as a "boundary condition" outside the context of time stepping. Unfortunately, Oceananigans is not consistent in how this philosophy is applied and when we are on "non-immersed" grids we have more of a "strong formulation" philosophy.

I believe the rationale for not transferring boundary conditions is the expectation that we will eventually adopt a "weak formulation philosophy" consistently throughout the code. In that case, non-default boundary conditions are meaningless on anything but prognostic fields that are evolved during time-stepping.

[glwagner]

Maybe we should create a function similar(bcs::FieldBoundaryConditions) that has the behavior we want. We can inspect bc.condition. If its Nothing, Number, or AbstractArray then I think copying makes sense. Otherwise we should not.

Note too that because immersed boundary conditions are always "weak" conditions I think we shouldn't copy them... ?

So

# Utility for boundary conditions
TransferrableBoundaryConditions = Union{Nothing, Number, AbstractArray}
similar(bc::BoundaryCondition) = bc.condition isa TransferrableBoundaryCondition ? BoundaryConditions(bc.classification, bc.condition) : NoFluxBoundaryCondition(nothing)
similar(bc::PBC) = PeriodicBoundaryCondition()

# Utility for FieldBoundaryCondition
similar(fbc::FieldBoundaryConditions) = FieldBoundaryConditions(similar(fbc.top), etc...)

Then we can use this utility for "similar" boundary conditions in similar(::Field).