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
• 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)
• Make the existing tests pass

And implements the following matching schemes:

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

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 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)

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

[jagoosw]

didn't mean to push this file yet

[johnryantaylor]

In this example the sponge on the inflow probably isn't necessary, but for more complicated inflows, it probably is needed.  Probably not a bad thing to keep it in the example, but I don’t have a strong view either way. On Mar 28, 2024, at 12:01 PM, Jago Strong-Wright ***@***.***> wrote: @jagoosw sorry to be dropping into this discussion a bit late, but I don't quite understand why the sponge is needed at the inflow in the example above. Shouldn't the prescribed inflow BC be enough? You are probably right it shouldn't be necessary here! —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>

[jagoosw]

Now that you have a basic example that illustrates an open boundary condition implementation with no matching scheme, you are in a position to implement a non-trivial matching scheme, and demonstrate its benefit. If the matching scheme has some benefit, then we are motivated to support it by adding source code and tests...

Okay I started doing this and realised that the other source code change we need is some way to update the boundary conditions as discussed a while ago. Are you happy with how I've done this? It's a line in update_state! just before fill_halo_regions! but I guess it could be moved into fill_halo_regions!. I'll also add it to the hydrostatic model for future use.

[jagoosw]

I also think that we should consider also another example problem like a flow around a cylinder because the internal wave problem seems to be hard to tolerate without large damping regions which makes it quite hard to tell how well the boundary condition is doing.

For example, in this flow around a cylinder it "works" by just setting the boundary values:

cylinder.mp4

But then if I run it with a 4x longer domain (not shown) you can see visually that the interior solution has been modified by the boundary and the vortexs are getting squashed as they approach the boundary:

cylinder_long.mp4

Here a damping region would similarly modify the internal solution, but if upwind-implicit Euler the boundary point with the boundary normal mean velocity we get:

cylinder_bulk_outflow_explicit.mp4

which looks closer to the long solution and the vortices appear to be leaving the domain with less deformation. We can also quantify the effect of the boundary condition with e.g. the Strouhal and the drag coefficient (which I'll work out for these examples in a minute).

[jagoosw]

I've been thinking about open boundaries more and have some things I want to test, but they're more questions about the numerics and I don't think this is the place for it.

From the last examples I sent + other experiments I've done with it I think matching schemes are justified because even in the nested case we don't want the boundary value to modify the outflow to modify the internal solution as in the cylinder example. For example, in a nested case, this would prevent higher resolution eddies from exiting the domain without un-physically modifying the upstream solution.

To resolve this PR I could tidy up a simple matching scheme where we compute the mean outflow on the boundary and do a 1D advection for the boundary point, or relax to the external state if there is inflow. I think this shows how to use all parts of the new infrastructure, and is a satisfactory boundary condition for some cases. Would this be okay?

[glwagner]

I also think that we should consider also another example problem like a flow around a cylinder

Please consider as many different cases as possible or as you need. I think a real nesting case would be even better --- for example, generate a turbulent solution using Oceananigans, and then illustrate embedding a high resolution domain inside that solution. A two-dimensional turbulent case would be sufficient.

To resolve this PR I could tidy up a simple matching scheme where we compute the mean outflow on the boundary and do a 1D advection for the boundary point, or relax to the external state if there is inflow. I think this shows how to use all parts of the new infrastructure, and is a satisfactory boundary condition for some cases. Would this be okay?

A simple matching scheme is a great place to start to justify changing the Open classification.

For example, in a nested case, this would prevent higher resolution eddies from exiting the domain without un-physically modifying the upstream solution.

I apologize if this is just a problem with communication style, but I don't think this statement can possibly be correct. The above example seems to demonstrate, for example, that the cylinder wake can be represented well by either making the domain large enough (with sponge layers), or by using a better numerical method with a smaller domain.

Basically this allows us to exclude statements like "we can't model X flow". The question is not whether or not we can represent the physics of a flow. The question is rather about representing the flow in a domain that is as small as possible / with optimal computational expense. We're pursuing performance optimization. Our language has to take the form "this scheme allows us to reduce the sponge layer thickness by XX%", or "this scheme reduces the width of the near-boundary region where the interior physics are modified". The statement "this scheme does not unphysically modify the solution" makes no sense to me. The presence of the open boundary is unphysical by definition...

I personally don't need much convincing that a matching scheme of some kind is needed. I'm just confused about how all of the different pieces that seem to enter into nesting interact with one another --- sponge layers, matching scheme, and also the particular physics of the flow in question (idealized in and outflows vs true nesting). With so many degrees of freedom and complexity I think a systematic approach is absolutely essential.

I really hope that the outcome of this work are a set of solid recommendations for nesting (better yet, a recipe that does not have to be changed --- a recommended matching scheme, or a combination of matching scheme and sponge layer) which is supported by unequivocal and rationally presented evidence.

[jagoosw]

I apologize if this is just a problem with communication style, but I don't think this statement can possibly be correct. The above example seems to demonstrate, for example, that the cylinder wake can be represented well by either making the domain large enough (with sponge layers), or by using a better numerical method with a smaller domain.

Yeah sorry, that was very imprecise of me. I more mean that that case was a good demonstration of a situation where we need a very large domain to not effect the interior solution in the area we are interested in, or have a matching scheme.

I really hope that the outcome of this work are a set of solid recommendations for nesting (better yet, a recipe that does not have to be changed --- a recommended matching scheme, or a combination of matching scheme and sponge layer) which is supported by unequivocal and rationally presented evidence.

That sounds good in the long term. Do you want that all for this PR or can we merge this part after I've implemented a simple matching scheme as a demonstration and then work on the rest elsewhere?

[glwagner]

Do you want that all for this PR or can we merge this part after I've implemented a simple matching scheme as a demonstration and then work on the rest elsewhere?

That's up to you. Smaller PRs can be easier because you will have less risk of merge conflicts. However you should make sure that the code in any individual PR is motivated and tested (ie if you implement a new type then it'd be best to have a use case for it, plus a test).

[johnryantaylor]

Thanks @jagoosw and @glwagner for your thoughts on this. @glwagner, the approach to nesting is going to be highly dependent on the particular case. The nesting strategy for a simulation with mesoscale eddies passing across the boundary will be very different from the strategy to nest an LES inside a low resolution simulation. I think that this PR should provide different types of boundary conditions (the 'ingredients') and the user will need to decide how to use them for their particular case (the 'recipe' which will rely on different combinations of the 'ingredients'). Does that make sense?

[glwagner]

Thanks @jagoosw and @glwagner for your thoughts on this. @glwagner, the approach to nesting is going to be highly dependent on the particular case. The nesting strategy for a simulation with mesoscale eddies passing across the boundary will be very different from the strategy to nest an LES inside a low resolution simulation. I think that this PR should provide different types of boundary conditions (the 'ingredients') and the user will need to decide how to use them for their particular case (the 'recipe' which will rely on different combinations of the 'ingredients'). Does that make sense?

Yes, that makes sense! I do think we need to have a way to test the different methods... ideally, the tests are non-trivial enough to illustrate that the numerics "works" for each ingredient that we implement (either independently or when combined with other ingredients to make a recipe). But the tests don't have to be exhaustive; I imagine it will take some time to develop our best practices.

[jagoosw]

In my last commit, I changed how boundary conditions are internally initialised as before we were passing the name of the classification so that we could catch functions and turn them into Discrete/ContinuousBoundaryFunctions, but this was causing a problem when I wanted the classification to also have properties (i.e. the matching scheme in this case).

Instead, I've changed it so we pass an instance of the classification (so e.g. we can pass Open(SomeMatchingScheme())) and this seems to work fine still.

The boundary conditions setup tests now pass locally and I will wait to see how the rest of the tests do, but is there any problem with this change otherwise @glwagner? This is also not an API change as it only affects how boundary conditions are initialised internally.

[jagoosw]

I'm not sure how the adapt_structures in boundary_condition.jl were working before but the method definition in boundary_condition_classifcation.jl fails if this isn't here

[jagoosw]

unclear to me what change would make this now fail (scalar indexing error), or why this was written as mean(interior( before?

[glwagner]

mean has not always been defined for Field and this is old code...

As for an out of bounds error, isn't that a bug?

[jagoosw]

Ah I see, it was a scalar indexing error here which I can reproduce:

julia> using Oceananigans
julia> grid = RectilinearGrid(GPU();size=(8, 8, 8), extent = (1, 1, 1))
julia> c = CenterField(grid)
julia> mean(interior(c))
ERROR: Scalar indexing is disallowed.
...

but I don't think that's from a bug here its from trying to use the base mean on a CuArray.

[glwagner]

Its on a view of a CuArray

[glwagner]

The problem might be related to the fact that interior(c) is not a contiguous view into parent(c) (because it omits halo regions in x, y). Possibly (but I'm not sure), such mean is not supported in CUDA. But I am a little surprised and it could also be a bug / misidrected method issue

[jagoosw]

I've been thinking and we definitely need a way to know if the halos are being filled before or after the pressure correction. The main reason is because if we are explicitly stepping the boundary point (e.g. in the bulk outflow case) then we will step twice when we only meant to step once. It does leave the question of what the correct thing todo to the boundary point after the pressure correction, but that is a different matter.

I'm going to update the top comment with a list of changes made here.

[glwagner]

I've been thinking and we definitely need a way to know if the halos are being filled before or after the pressure correction. The main reason is because if we are explicitly stepping the boundary point (e.g. in the bulk outflow case) then we will step twice when we only meant to step once. It does leave the question of what the correct thing todo to the boundary point after the pressure correction, but that is a different matter.

I'm going to update the top comment with a list of changes made here.

The way our algorithm works, we must enforce the penetration / no-penetration boundary conditions on wall-normal velocities prior to the pressure correction. Calling fill_halo_regions! does a bit more; we actually don't need to fill the other halos.

Then we fill_halo_regions! after the time step is complete.