Boundary conditions cannot be enforced with closure = nothing

[glwagner]

In test_dynamics.jl, an experiment is set up with inhomogeneous boundary conditions on temperature:

Oceananigans.jl/test/test_dynamics.jl

Lines 310 to 312 in f83fce8

	y_bc = GradientBoundaryCondition(∂T∂z * g̃[2])
	z_bc = GradientBoundaryCondition(∂T∂z * g̃[3])
	T_bcs = TracerBoundaryConditions(grid, bottom=z_bc, top=z_bc, south=y_bc, north=y_bc)

while using closure = nothing:

Oceananigans.jl/test/test_dynamics.jl

Lines 314 to 321 in f83fce8

	model = IncompressibleModel(
	architecture = arch,
	grid = grid,
	buoyancy = buoyancy,
	tracers = (:T, :S),
	closure = nothing,
	boundary_conditions = (T=T_bcs,)
	)

This is at best misleading --- with closure = nothing, boundary conditions are not enforced. It may also be innocuous though, because there is no diffusive flux anywhere -- not only across boundaries.

It may in fact make sense to throw a warning when using closure = nothing with non-default boundary conditions ("Boundary conditions are not enforced when isnothing(closure)").

[francispoulin]

Good point, but no-normal flow boundary conditions should still be enforced, even without any closure scheme, right?

[glwagner]

That's a good point ! Yes, we still have "first order" boundary conditions on wall-normal velocity components. Technically we could "enforce" these by "refusing" to update those velocities (but this isn't how the algorithm works, so we do need fill_halo_regions! for wall-normal velocity bcs).

[francispoulin]

Glad we agree. If we have a warning maybe we can specify so the user knows no-normal flow is still enforced?

[glwagner]

I think the warning should read

Flux, Value, and Gradient boundary conditions are not enforced when isnothing(closure).

[francispoulin]

That is very clear and accurate. Agreed!

[tomchor]

I understand why FluxBoundaryConditions can't be enforced with closure=nothing since they depend on a diffusivity. But why can't Gradient and Value BCs be enforced?

[francispoulin]

When you turn off closure you reduce the order of the differential equations and therefore no longer have the freedom to impose anything else but no-normal flow. It is due entirely to the order of the differential equations and not do to whether the boundary conditiosn are flux, gradient or value. Does that make sense?

[tomchor]

Ah, that makes sense. I'd never thought of that...

I agree with the suggestion of adding a warning btw

[glwagner]

@francispoulin explained the situation well for Gradient and Value boundary conditions. I'd like to restate what he said in the context of the finite volume method and clarify the situation for Flux boundary conditions because I said something wrong above...

Oceananigans uses a weak formulation of PDEs via the finite volume method. We evolve the average value of fields, integrated over a cell volume. In this integral formulation there are two contributions to the tendency of a field: fluxes across cell interfaces (advective and diffusive usually, but also possibly others via user-defined forcing), and interior "sources" (such as pressure gradients, Coriolis forces, gravitational forces).

When users invoke Gradient or Value boundary conditions, then the same model is used for fluxes across interior cell interfaces as for "exterior" interfaces (domain boundaries). Thus if users set closure=nothing and zero out diffusive fluxes, boundary fluxes are always zero regardless of what Gradient or Value are specified. More specifically, we enforce Gradient or Value boundary conditions by filling halo regions and then calculating exterior fluxes using the same kernel that's evaluated in the interior. As @francispoulin points out this is not some quirk of our algorithm or numerics; this follows from the properties of the underlying PDE.

Something different happens when Flux boundary conditions are specified and this is where I was wrong above. With Flux boundary conditions, users are really specifying two different models for fluxes: one across "interior" cell interfaces, and another for fluxes across domain boundaries. Setting closure=nothing in this case only nullifies interior diffusive fluxes. But since boundary fluxes are explicitly specified in this case, they still do have an impact on interior tendency values.

A further subtlety is that Value and Gradient boundary conditions do actually have an impact if there is non-zero NormalFlow advecting the field across a boundary; in this case the flux is advective (first-order) and thus mathematically consistent with the underlying PDE.

So we should probably amend the warning to read:

Value and Gradient boundary conditions are not enforced across Impenetrable boundaries when isnothing(closure).

[francispoulin]

Thanks @glwagner for the detailed explanation and putting it context of the finite volume formulation, which is critical to Oceanangians. I'm happy with the new suggestion.

It occurs to me that there is some mention of finite volume in the docs, but not really a derivation. If we wanted the user to be able to understand all the details, would it be worth while having the derivation of the finite volume method? If yes I'd be happy to do some of this as I do think it would add some clarity. But only if there is a consensus.

[tomchor]

If yes I'd be happy to do some of this as I do think it would add some clarity. But only if there is a consensus.

I definitely agree with that, but if you're looking for consensus, it might be best to create a separate issue to vote this out.

[glwagner]

I'm closing this issue because I'm judging that it's not of current, timely relevance to Oceananigans development. If you would like to make it a higher priority or if you think the issue was closed in error please feel free to re-open.