(0.91.7) Open boundary conditions for NonhydrostaticModel

[jagoosw]

This PR implements the infrastructure for open boundary conditions in the NonhydrostaticModel as well as a few simple methods.

This PR:

• Adds a matching_scheme property to Open boundary classifications to allow different fill_X_halo! to be dispatched
• Introduces update_boundary_condition! to be called before halo fills allowing the mathcing_scheme to have properties which can evolve with the model
• Make the existing tests pass
• and implements a zero gradient matching scheme

(Others please feel free to update this comment as necessary.)

---

Hi all,

Following discussion with @glwagner, @simone-silvestri, and @jm-c this is a first attempt at implementing open boundary conditions. First I will try to get it working for the non-hydrostatic model which seems to be relatively straightforward.

As a first step, I have implemented east/west boundaries which allow a flow to be prescribed or to travel out of the domain (i.e. if you set OpenBoundaryCondition(nothing) then my code is assuming the flow will be travelling into that boundary). The outflow condition is equivalent to having a nondimensional phase speed of 1 (sec 3.1 of https://doi.org/10.1016/S0924-7963(97)00023-7) which seems to work fine.

When I vary the inflow velocity I do see waves in the velocity field reflecting from the downstream boundary. I gather that we expect this with any outflowing boundary and would remedy this with a sponge layer, but maybe this is where we need to add something to the pressure solver.

With a constant inflow (0.1 m/s, 0.625m resolution, AMD and SeawaterBuoyancy), and a turbulent outflow, it seems to work okay:

nonhydrostatic_channel.mp4

I am going to try and implement the test cases in the paper mentioned below next.

I have a few things I think we should think about:

• Does it make more sense for Open to be a topology because otherwise, users have to manually specify boundary conditions on everything? Then we can automatically set OpenBoundaryCondition(nothing) unless a user sets something else.
• Should we have an automatic way to setup a sponge layer?
• Do we think it is correct not to add a term to the pressure for time-varying inflow? I think we do not (see below).

How I understand the maths:

Starting from the momentum equation:

$\partial_t\vec u = \vec{G_u} - \nabla P$

We split the time step by defining:

$\vec u^\star - \vec u^n = \int_{t_n}^{t_{n+1}}\vec{G_u}dt$

so

$\vec u^{n+1} - \vec u^\star = -\int_{t_n}^{t_{n+1}}\nabla P_{non}dt \approx -\Delta t\nabla P_{non}$

we take the divergence of this (and requiring $\nabla \cdot \vec u^{n+1}$) to give:

$\nabla^2P_{non} = \frac{\nabla \cdot \vec u^\star}{\Delta t}$

To enforce boundary conditions on $u^{n+1}$ they are instead enforced on $u^\star$. In the case where we prescribe an inflow (e.g. $u(x = 0) = 1$) this results in the $\nabla \cdot \vec u^\star$ term having the term u_{2jk} - u_{1jk} at the i=1 point (in the central difference case, I don't know where spatial operation approximations are defined in the code) change from u_{2jk} - 0 in the no penetration case to `u_{2jk} - u_{open boundary}$.

In the case where the interior velocity is already the same as the specified velocity (and everything else is uniform), this means that at i=1 $\nabla^2 P_{non}$ changes from being positive to zero, so we go from having a pressure gradient at the wall counteracting the flow away from the wall (or I suppose a reconfiguration to u=0 everywhere to enforce compressibility) to having no pressure gradient across the wall.

Thinking about it this is exactly what happens for periodic domains where we are essentially specifying the flow outside of the domain. That makes me think that we don't need to do anything else for time-varying inflows.

As for making the outflows more correct, I think we should be able to extend to the method for calculating ht phase velocity by the method in 3.3 of https://doi.org/10.1016/0021-9991(83)90127-4 (linked in #833) which doesn't depend on previous time steps as some other Orlanski methods do, but it does depend on the time difference of the interior solution with the next step which currently does not get passed to the boundary conditions. Perhaps it might be most straightforward to evaluate $c=-\frac{\partial\phi/\partial t}{\partial\phi/\partial x}$ by passing the tendencies and using $\partial\phi/\partial t = G_n$ (although this isn't quite correct for the velocity so maybe $-\nabla P$.

I think passing the tendencies automatically is going to require some materialization step when the model is setup to pass $G^n$ into the boundary conditions but I know we're trying to avoid doing that so any other suggestions would be useful. I've started testing this by initialising the timestepper first but it is a bit clumsy.

[jagoosw]

Bounded, Bounded, Bounded domain with initial velocity then the boundary velocity 0.1tanh(t/10)

wave.mp4

[jagoosw]

cylinder.mp4

Qualitatively this looks like it is generally okay, but I am having problems when the right boundary ends up with an inward flow as it seems to kind of stick to it (as seen early on in this video) and sometimes grows to form a constantly inflowing region which then doesn't disappear. Perhaps this is a more general problem with Orlanski boundaries?

Example of it going wrong:

cylinder_new.mp4

[navidcy]

Suggested change

	@inline function _fill_west_halo!(j, k, grid, ϕ, bc::OBC, loc, args...)
	ϕ[1, j, k] = getbc(bc, j, k, grid, args...)
	end
	@inline function _fill_west_halo!(j, k, grid, ϕ, bc::OBC, loc, args...)
	@inbounds ϕ[1, j, k] = getbc(bc, j, k, grid, args...)
	return nothing
	end

[navidcy]

why not use @inbounds here?

also, what's eps(0.0)? is this intentional?

[glwagner]

Probably that should be a parameter of the boundary condition.

We'll have to figure out the logic of the code as well. Perhaps we can add appropriate tendencies within regularize_boundary_conditions in the model constructor. Another option is to pass the tendencies into fill_halo_regions! which could make sense...

[jagoosw]

I had to have the eps(0.0) for when ϕ[i, grid.Ny, k] - ϕ[i, grid.Ny - 1, k] = 0 for example if you set the field to a constant value

[glwagner]

Why do we have to pass the tendency here? Isn't there a wave speed that can be set?

[jagoosw]

I don't think so, it seems like the wave speed is usually estimated as $-\frac{\partial_t\phi}{\partial_x\phi}$ where $\phi$ is whatever is going through the boundary

[glwagner]

Ah I see. Sounds like @simone-silvestri has a different method though that determines the wave speed differently?

[glwagner]

If we only fill one point into the halo regions, then we can keep the current advection scheme logic where we limit to second order advection on the boundary. Alternatively though, it seems that we could fill more points and then do ordinary advection. In that case though, we may need a new topology for open boundaries. Not sure.

[glwagner]

Can you do a validation case that has implicit vertical diffusion? Does it work in that case?

[glwagner]

Should we have an automatic way to setup a sponge layer?

We do have Relaxation, but do you mean something more than that? Code example?

[glwagner]

Does it make more sense for Open to be a topology because otherwise, users have to manually specify boundary conditions on everything? Then we can automatically set OpenBoundaryCondition(nothing) unless a user sets something else.

This choice should not merely be a question about user interface / convenience but also about how the code internals work.

One problem is that the topology refers to both sides. We want to support domains that are, for example, bounded on the west but open on the east.

We do have an abstraction called RightConnected for distributed cases. Possibly, we can implement topologies that represent doubly open and single-sided open.

But to motivate such an abstraction, I think this needs to have implications on the grid level --- not just a way to generate boundary conditions conveniently. I think there are other solutions for generating boundary conditions conveniently.

[glwagner]

I think passing the tendencies automatically is going to require some materialization step when the model is setup to pass into the boundary conditions but I know we're trying to avoid doing that so any other suggestions would be useful. I've started testing this by initialising the timestepper first but it is a bit clumsy.

Well, we already do materialize boundary conditions, so possibly this isn't such a big deal. Another possibility is to pass the tendencies through to fill_halo_regions!, but that has wider implications that we'd have to think about (for example should the tendencies also be passed on to discrete-form boundary condition functions?)

[simone-silvestri]

for tracers you need to specify ValueBoundaryCondition. Open will change the value at 1 but you want to change the value outside the domain for temperature

[glwagner]

Hm no I think these are different. ValueBoundaryCondition sets the tracer value on the boundary. Whereas for an open boundary condition, we are specifying the tracer values in a cell that lies outside the boundary.

[simone-silvestri]

OpenBoundaryCondition sets the value at 1 which for the tracer is inside the domain, not a boundary. Value should be used to set the boundary for tracers (another option is to redefine Open for centered fields).

[glwagner]

Value doesn't set the tracer value --- it computes the halo point by extrapolating from the interior point through the boundary point.

[glwagner]

OpenBoundaryCondition is supposed to fill halo regions. If it doesn't that's a bug.

[simone-silvestri]

That is what I was saying, to set a boundary for a tracer, Value is what we need. Open is for velocities. If we use Open for both T and u we end up with boundary conditions set at different points.

[glwagner]

Then we have to fix Open to be "location-aware".

Value determines the halo value by extrapolating across the boundary. This is not always what we want for an open boundary condition. I'm thinking of an open boundary as connecting to another domain that is adjacent to the prognostic domain. For this we want to specify the halo state directly. The reconstructed value of a tracer should be the average of the specified halo state and the interior state. This is different from Value.

Possibly there are cases where we want to combine Value on tracers with Open on velocities (this would have the feature of fixing the tracer fluxes), but I think for realistic regional modeling applications we will instead want to fix the value of the tracer within the halo region. This can be thought of as fixing the prognostic state in cells outside the domain.

Maybe @jm-c has more to say.

[simone-silvestri]

I think for nonhydrostatic flows instead of specifying a phase speed, it is more physically sound to set an outlet velocity.
For example: $$\frac{\partial \phi}{\partial t} + U_b \frac{\partial \phi}{\partial n} = 0$$ where $U_b$ is the bulk velocity perpendicular to the outlet.
These are examples of developing (nonhydrostatic) heated channel flows calculated using the above outflow boundary condition and no changes to the pressure solver (the Ub is calculated at the last interior grid point).
The code is a staggered c-grid so the same algorithm should work well.
However, there is no buoyancy term, so that might create some problems?

out.mp4

outW.mp4

[glwagner]

I think for nonhydrostatic flows instead of specifying a phase speed, it is more physically sound to set an outlet velocity. For example: ∂ϕ∂t+Ub∂ϕ∂n=0 where Ub is the bulk velocity perpendicular to the outlet. These are examples of developing (nonhydrostatic) heated channel flows calculated using the above outflow boundary condition and no changes to the pressure solver (the Ub is calculated at the last interior grid point). The code is a staggered c-grid so the same algorithm should work well. However, there is no buoyancy term, so that might create some problems?

out.mp4
outW.mp4

Isn't this the same thing? We're just calling "Ub" a phase speed, since it isn't a physical velocity.

[simone-silvestri]

$U_b$ changes with the flow so it allows information to exit the domain at the same speed of the flow evolution

[glwagner]

I gather that we expect this with any outflowing boundary and would remedy this with a sponge layer, but maybe this is where we need to add something to the pressure solver.

I don't think so --- I think we are free to add a sponge how we want. The sponge would restore velocities and tracers to the externally imposed state?

[glwagner]

Ub changes with the flow so it allows information to exit the domain at the same speed of the flow evolution

Okay but just to clarify, you are proposing to dynamically estimate the phase speed somehow rather than specifying it, right? Either way, we have an outlet phase speed, somehow.

[simone-silvestri]

Yep, I would call it bulk velocity though, instead of phase speed, and change the name from Orlanski to something more descriptive like AdvectiveOutflow.

[glwagner]

Yep, I would call it bulk velocity though, instead of phase speed, and change the name from Orlanski to something more descriptive like AdvectiveOutflow.

But these are mathematically identical right? Orlanski called is a "phase speed", but "outflow velocity" is equally valid and refers to exactly the same mathematical object. The reference you posted says

The test results also confirm that this type of boundary condition, which was originally designed by Orlanski primarily for equations which are hyperbolic in nature, also performs well for parabolic problems.

I think we can keep the name "Orlanski" and provide a generic interface for specifying the outflow speed (whatever you want to call it). It can be user-specific, dynamically computed, etc. The code can be extensible, so if users want to experiment with different methods for computing the outflow speed, this is possible.

[glwagner]

Probably OrlanskiCondition is the right term here, since this is intended to be used in the condition property.

Btw another possibility is to simply use "Orlanski" as a synonym for "Open", and add properties to the Open classification. Not sure what's best.

[jagoosw]

I don't think "Open" always means "Orlanski" though since Orlanski is just one way of approximating the boundary point but e.g. specifying the value (as I have done for in flows) is another valid "Open" boundary?

[glwagner]

I see what you're saying. In that case, perhaps just using either Orlanski or OrlanskiCondition is the improvement we need here.

[jagoosw]

I think for nonhydrostatic flows instead of specifying a phase speed, it is more physically sound to set an outlet velocity. For example: ∂ϕ∂t+Ub∂ϕ∂n=0 where Ub is the bulk velocity perpendicular to the outlet. These are examples of developing (nonhydrostatic) heated channel flows calculated using the above outflow boundary condition and no changes to the pressure solver (the Ub is calculated at the last interior grid point). The code is a staggered c-grid so the same algorithm should work well. However, there is no buoyancy term, so that might create some problems?
out.mp4
outW.mp4

I think this is what I'm doing because I'm estimating the "phase speed" as $-\frac{\partial_t\phi}{\partial_x\phi}$ which is the same as $U_b$ here, is that what you mean?

[jagoosw]

Does it make more sense for Open to be a topology because otherwise, users have to manually specify boundary conditions on everything? Then we can automatically set OpenBoundaryCondition(nothing) unless a user sets something else.

This choice should not merely be a question about user interface / convenience but also about how the code internals work.

One problem is that the topology refers to both sides. We want to support domains that are, for example, bounded on the west but open on the east.

We do have an abstraction called RightConnected for distributed cases. Possibly, we can implement topologies that represent doubly open and single-sided open.

But to motivate such an abstraction, I think this needs to have implications on the grid level --- not just a way to generate boundary conditions conveniently. I think there are other solutions for generating boundary conditions conveniently.

I see, I think having an open boundary does necessarily have grid level implications because every tracer needs to have some open boundary specified if the grid has an open boundary right? I suppose the other way todo this is to check if some tracer has an open boundary specified and add them to all of the others. It would be a bit confusing to keep having this as a boundary condition applied to a Bounded direction though since the grid isn't bounded in the same way anymore.

[jagoosw]

If we only fill one point into the halo regions, then we can keep the current advection scheme logic where we limit to second order advection on the boundary. Alternatively though, it seems that we could fill more points and then do ordinary advection. In that case though, we may need a new topology for open boundaries. Not sure.

I guess this might be useful for outflows but we would have the problem of having insufficient information if the flow is inflowing, or if we are specifying the values as it might be that the specified inflow is from a much coarser grid and we only have one value to give?

[glwagner]

I see, I think having an open boundary does necessarily have grid level implications because every tracer needs to have some open boundary specified if the grid has an open boundary right?

What are those implications?

For example, Periodic implies that there is one fewer interface for Face fields in that direction. As a result, Periodic direction changes the way that a grid needs to be constructed. This is what I mean by "grid-level implication". What's the difference between Open and Bounded in terms of data layout? It seems to me that Open is a boundary condition rather than dictating a different data layout. For example, would the different topology change the way that we construct the grid?

[glwagner]

I suppose the other way todo this is to check if some tracer has an open boundary specified and add them to all of the others. It would be a bit confusing to keep having this as a boundary condition applied to a Bounded direction though since the grid isn't bounded in the same way anymore.

I would approach the problem with a utility function, something like

open_bcs = open_boundary_conditions(:south, :north, :west)
u_bcs = FieldBoundaryConditions(top=wind_stress, bottom=u_drag; open_bcs...)

for example. We could also add a kwarg to the model constructor like

model = NonhydrostaticModel(; grid, boundary_conditions = (; u=u_bcs),
                              open_boundaries = (:south, :north, :west))

which could automatically edit the boundary conditions for every prognostic field accordingly. There's a lot of other possibilities too. The point is that convenience features are easy to design. We shouldn't motivate superfluous topologies with a convenience argument, when it has no material effect on the grid itself.

[jagoosw]

Ah yes, I see your point. I suppose having an argument for in the boundary condition regularization would be the easiest way todo this.

I'll try clean this up a bit in the next few days.

[jagoosw]

Playing around with an internal wave test case I think we actually need todo something more like the adaptive boundary described in section 4.1 of this paper https://doi.org/10.1016/S1463-5003(00)00013-5 as I have come across two problems: when the flow is directed out of the domain on a prescribed interface (e.g. u = cos(pi/h(z+h)) then information can't get out, and on the "Orlanski" side where information is travelling into the domain I am getting instability as it is just keeping the boundary value constant which by default is zero.

This might present some more user interface issues as it is going to require us to set a "known" value on every open boundary unless we're confident that the flow will only be leaving.

[jagoosw]

I've just realised as implemented now I've got the bulk velocity wrong by a factor of $\Delta t$.

Boundary conditions don't get given $\Delta t$ and as far as I can tell we either need the previous step values (which I suppose we could store in the boundary condition), or we need tendencies * $\Delta t$ as a guess at them. This would not really work with the separation of model and simulation though.

[jagoosw]

I've played more with this now and think the way we can do this with the minimum code changes is as follows. The only other way I can see is for the boundary to store the previous timestep of the boundary adjacent points but since Fields depends on BoundaryConditions that isn't really possible.

The key problem really is that previous solutions to this problem have been written in codes that store at least two time levels which we don't have.

So, if we consider for now just a 1D problem with boundary point $\phi_b$ and interior points at $\phi_{b-1}$ etc. When we go to update the boundary point we have $\phi_b^n$ and $\phi_{b-1}^{n+1}$ and we want $\phi_b^{n+1}$. As per previous work we assume that the bulk speed is the same at both $b$ and $b-1$ but we don't have both the spatial and time derivatives of the $\phi$ interior at the same step so we first need to approximate the previous step as:

$\phi_{b-i}^n = \phi_{b-i}^{n+1}-\Delta t G^n_{b+1}$

We can then find the bulk velocity at timestep n as:

$U_b^n = -\frac{2\Delta x_{b-1}G^n_{b-1}}{\phi^n_b - \phi^{n+1}_{b-2} - \Delta t G^n _{b-2}}$

We then have all the information to step $\phi_b$ to:

$\phi_b^{n+1}=\phi_b^n-\frac{\Delta t}{\Delta x_b}U_b(\phi_b^n-\phi_{b-1}^{n+1} - \Delta t G^n_{b-1})$.

This will require us to give the boundary condition both the tendencies and $\Delta t$, but this seems to be the easiest thing to change. I also think this is the only way we can get a physically sensible bulk speed where all of the components are calculated at the same timestep.

I have also realised that we need to have an exterior value for every open boundary for when the flow spontaneously becomes an inflow so I think it would make sense to have every open boundary be the same and just be OpenBoundaryCondition(external_value), and then put the tendencies etc in as arguments to _fill_X_halo!. Then when U_b is negative we either set the value or do nudging like in ROMS to prevent shocks (but I think this is a question for further down the line). We might also want to consider oblique waves but that shouldn't be too hard to extend to we might just have to calculate a lot of surrounding points which might motivate some other way to do it.

Does this make sense to everyone? I also can't see an obvious way to get $\Delta t$ to the halo fill since it doesn't even make it into the update_state! so any suggestions would be appreciated.

Update:

I'm not sure if

[glwagner]

the lessons here about coding style will also be important for OceanBioME if you ever want to support automatic differentiation there.

[glwagner]

Try this

@inline boundary_conditions(f::Field) = f.boundary_conditions
@inline boundary_conditions(w::WindowedField) = FieldBoundaryConditions(w.indices, w.boundary_conditions)

The @inline sometimes encourage type inference. WindowedField is defined just above boundary_conditions in field.jl.

[jagoosw]

It looks like that worked for the enzyme tests!

[jagoosw]

Are doctest failures where the output writer is reporting the file size differently to what the doctrine from a different PR and fixed now?

[glwagner]

The doctests should pass on main so you should be to able to rely on that. The annoying bit is if your PR also changes the doctests, then you need to take some care.

[jagoosw]

Happy for me to merge this now @glwagner?

[glwagner]

merge away!

[jagoosw]

Before I merge, I don't think this is a breaking change any more because I added a method:

BoundaryCondition(Classification::DataType, args...) = BoundaryCondition(Classification(), args...)

so the "old" syntax of e.g. BoundaryCondition(Value, 1) (v.s. BoundaryCondition(Value(), 1)) still works, do we want to support this legacy way of creating a boundary condition (and change the version number to just a patch increase), or should I remove this method and make this a breaking change?

[simone-silvestri]

I would keep the "old" way to define boundaries. It does not seem damaging to have and it is a convenience function that still stands.

[jagoosw]

Okay I'll change the version number and then merge