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:

• [x] Adds a matching_scheme property to Open boundary classifications to allow different fill_X_halo! to be dispatched
• [x] Introduces update_boundary_condition! to be called before halo fills allowing the mathcing_scheme to have properties which can evolve with the model
• [x] Adds a kwarg to update_boundary_condition! so it can know if the model is currently in the predictor or corrected state (if the pressure correction has occurred)
• [x] Make the existing tests pass

And implements the following matching schemes:

• [x] zero gradient (in the last cell center)
• [ ] bulk mean outflow

And adds the validation cases with:

• [ ] these matching schemes
• [ ] matching schemes + relaxation region

(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:

https://github.com/CliMA/Oceananigans.jl/assets/26657828/1509c8fa-ea24-40a2-84d5-33da424567b1

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 https://github.com/CliMA/Oceananigans.jl/issues/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)

https://github.com/CliMA/Oceananigans.jl/assets/26657828/7d69fe2b-a0dd-4f2e-996a-ebf185dc7e87

[jagoosw]

https://github.com/CliMA/Oceananigans.jl/assets/26657828/11a948b1-d342-4e17-9059-16f1c09366df

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:

https://github.com/CliMA/Oceananigans.jl/assets/26657828/f0326138-a589-4622-9373-ce836d604bc3

[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]

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?

https://github.com/CliMA/Oceananigans.jl/assets/33547697/e1995799-69f3-48f8-989f-6730ff9e9e1b

https://github.com/CliMA/Oceananigans.jl/assets/33547697/a00a23f7-c5e6-4d77-b359-b9ffe898a7a0

[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.

[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

[jagoosw]

Just to clarify, for a Face field in a Bounded direction the first halo points are at 0 and N + 1, or is in N+2 with 1 to N+1 being interior/on the boundary points?

[navidcy]

You are right! For Bounded and Face, the N+1 point is part of the interior.

julia> using Oceananigans.Grids: total_length

help?> total_length
search:

  total_length(loc, topo, N, H=0, ind=Colon())

  Return the total length of a field at location along one dimension of topology with N centered cells and H halo cells. If ind is provided the total_length is
  restricted by length(ind).

julia> N = 10
10

julia> total_length(Face(), Bounded(), N)
11

julia> total_length(Face(), Periodic(), N)
10

julia> total_length(Center(), Bounded(), N)
10

julia> total_length(Face(), LeftConnected(), N)
11

julia> total_length(Face(), RightConnected(), N)
10

[jagoosw]

using Oceananigans.Grids: total_length

This is strange because if I set the Nx + 2 point as the boundary it doesn't end up working (for face bounded fields)

[simone-silvestri]

I think you have to set the velocities at the boundary points to ensure that the pressure solver works correctly. So basically you should set velocities perpendicular to open boundaries at 1, N+1 and tracers at 0, N+1

[glwagner]

It's as @simone-silvestri said, more broadly it seems that we need to think of boundary conditions differently between Face fields and Center fields. Face fields have nodes on the boundary (sometimes we call those "peripheral nodes"), so we simply impose boundary values to satisfy boundary conditions, eg for determining wall-normal velocities. Center fields have nodes within and outside, but not on the boundary. This is the core of the discussion @simone-silvestri and I had about tracer boundary conditions: we can choose either to set the tracer halos to produce a result when tracers are reconstructed on the boundary (that's how Value and Gradient boundary conditions work -- the halos are determined by extrapolation). Or, we can set the tracer halos as if the halo region is part of some other prognostic state (which is what I envisioned Open should do). Finally as @simone-silvestri says the tracer halos are 0, N+1 while wall-normal velocity points on the boundary are 1 and N+1.

I think its fine to pass additional arguments to fill_halo_regions!, and use these for open boundary conditions. In fact, DiscreteBoundaryCondition and ContinuousBoundaryCondition appear to already be designed to accomodate additional arguments that are unused by those user interfaces (note the args...):

https://github.com/CliMA/Oceananigans.jl/blob/643b484e81e0aeb038b3038266912ad051bce9b8/src/BoundaryConditions/discrete_boundary_function.jl#L45

https://github.com/CliMA/Oceananigans.jl/blob/643b484e81e0aeb038b3038266912ad051bce9b8/src/BoundaryConditions/continuous_boundary_function.jl#L136-L137

[glwagner]

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 Δt to the halo fill since it doesn't even make it into the update_state! so any suggestions would be appreciated.

To clarify, the statement "when U_b is negative" means more generally "when U_b implies inflow into the domain". It makes sense to me that inflow is prescribed, whereas outflow is a little tricker, requiring us to smoothly advect information out of the domain.

Having the boundary condition depend on the previous time-step does seem to break assumptions we have made about how information flows. For example, we call update_state! during model construction. How does the first time-step / initialization work with this kind of outflow model? It might help to consider that.

Practically speaking it would not be hard to pass the time-step to update state. But it might be nice to understand the implications before committing, because a different design / redesign could make sense. For example I supposed we would need to have a default "time step" that is 0.

[jagoosw]

From what I've read the reason we might need to not just exactly prescribe inflows is because if the interior solution varies from the prescribed inflow (and possibly outflows) and then suddenly switches it can cause shocks.

As for the first-time step. It looks like the pressure solver deals with this by setting the default timestep to one, and that's whats used when the fields are set etc. Perhaps changing the order in update_state! would make sure that the tendency was correct on the first time step? I'm not sure what other implications that might have and they would probably also need to be recalculated after the halos were filled.