Add free-slip, partial-slip and no-slip momentum BCs

[dengwirda]

This PR updates the way horizontal momentum boundary conditions are computed --- highlighting that the current MPAS discretisation applies no-slip conditions, introducing new free-slip and partial-slip conditions, and cleaning-up the way cell-vertex remapping is handled across partially masked cells/duals in general.

No-slip and free-slip BCs can be expressed in terms of $\nabla \times \mathbf{u}$ on the boundaries --- no-slip conditions (the current implementation) correspond to a 'full' relative vorticity computed on boundary vertices, while free-slip conditions are imposed by setting $\nabla \times \mathbf{u} = 0$ on boundary vertices (corresponding to no induced torque when a free-stream velocity flows tangential to a wall). I've also introduced so-called 'partial-slip' conditions (see e.g. NEMO) as an interpolation between the no-slip and free-slip cases, leading to a reduced vorticity on boundaries. See e.g. Roullet and Gaillard for confirmation of how momentum BCs should be applied in a vector-invariant formulation.

A new lateral_walls namelist section for BC settings is introduced, where config_wall_slip_factor controls the amount of slip at horizontal walls (both coastlines at the surface and masked boundaries against bathymetry in general). config_wall_slip_factor = 0.0 is a no-slip condition, config_wall_slip_factor = 1.0 is a free-slip condition, and values between 0.0 and 1.0 correspond to partial-slip cases.

No-slip (left) vs free-slip (right) behaviour can be seen in the channel case below (developed with @hetland) in which the initially uniform zonal flow generates strong vorticity sheets at the north/south boundaries in the no-slip case, while the flow remains quiescent in the free-slip case.

This PR also updates the way vorticity and layer-thickness is computed and remapped at partially masked cells and duals at boundaries, ensuring that the numerical stencils reflect the masking. In the case of computing potential vorticity $\frac{1}{h}(\nabla \times \mathbf{u} + f)$ it's necessary to remap layer thickness from cells to duals (triangles) using the intersecting cell-dual 'kite' areas. Currently, I feel we are not accounting for masked cells/duals correctly. Consider a case at a boundary in which one of the three cells surrounding a boundary vertex is masked, and where a spatially uniform flow with layer-thickness $h$ exists inside the domain.

In the current implementation the thicknesses from the two valid cells is remapped along with a 0 value for the masked cell. The result is divided by the full area of the dual (triangle) resulting in a vertex layer-thickness of only $\frac{2}{3} h$, inconsistent with the spatially uniform interior thickness. An alternative approach implemented here is to accumulate only the unmasked area of the dual (triangle) and use this in the remapping, which will results the expected vertex thickness of $h$.

$$h_{v_{old}} = \frac{A_1 h_{c_1} + A_2 h_{c_2} + 0}{A_{tria}} = \frac{A_{tria}}{A_{tria}} (\frac{1}{3}h_{c_1} + \frac{1}{3}h_{c_2}) = \frac{2}{3}h$$

$$h_{v_{new}} = \frac{A_1 h_{c_1} + A_2 h_{c_2}}{A_1 + A_2} = \frac{A_{tria}}{A_{tria}} (\frac{\frac{1}{3}h_{c_1} + \frac{1}{3}h_{c_2}}{\frac{1}{3} + \frac{1}{3}}) = h$$

$$A_{1,2,3} = \frac{1}{3}A_{tria}$$

I've tested both in channel and global standalone spin-ups. Results are non-BFB, with somewhat increased kinetic energy developing in global spin-ups with free-slip BCs due to the reduced dissipation compared to no-slip lateral boundaries.

• Roullet, G. and Gaillard, T., 2022. A fast monotone discretization of the rotating shallow water equations. Journal of Advances in Modeling Earth Systems, 14(2), p.e2021MS002663.

• https://www.nemo-ocean.eu/doc/node58.html

[xylar]

@dengwirda, same questions as for #48: Is this ready to review? If so, who would you like to review it before it goes to E3SM?

[dengwirda]

@xylar, let me know how you'd like the 'discussion' aspect of this to proceed --- this branch is ready to run, but I'd say that additional testing is required before review.

[xylar]

@dengwirda, typically these PRs on E3SM-Ocean-Discussion only proceed if you ping specific people to indicate that you want their feedback. You can do that by adding them as a reviewer or just by @ing them and saying what kind of feedback you're looking for. Neither of those typically requires the code to be tested, though it could make things a bit clearer if you just ping people to initiate discussion and then flag them for a review when things are tested. The review here is typically cursory -- are things far enough along that we don't think we'll have a lengthy discussion (both in time and number of posts) on E3SM?

[dengwirda]

@xylar please feel free to add extra reviewers and/or review yourself here.
The slip BCs are hopefully not controversial, though the implications of using them in different config.'s is something to discuss.
The updates to layerThicknessVertex at boundaries will/may change behaviour in general.

[xylar]

Is boundaryVertex also needed in the present() arguments?

[xylar]

@dengwirda, your list of reviewers looks good to me.

[xylar]

@dengwirda, I suggest you make one PR to E3SM to fix dualArea as you have done here. I think we just want to do that once and for all, and it will be a non-BFB PR for sure. Then, we can do a second PR with the boundary conditions. That one can leave no-slip as the default for now but we can test the other options and consider switching to them before the v3 code freeze.

[mark-petersen]

@dengwirda this looks like a good addition. When you compare the current code with this PR using config_wall_slip_factor = 0.0, are the differences reasonably small? It that case, it looks like order-of-operations differences, so should be <1e-10 for the first few time steps. If that is the case, I think we could proceed to a PR to E3SM-Project with the default set to config_wall_slip_factor = 0.0.

[xylar]

@mark-petersen, I could be wrong but my understanding was that, because the area for vertices was being computed with full triangles rather than from the sum of the kites that actually exist, the differences will not be small and we will need to treat this as a (hopefully) non-climate-changing PR. That is why I was suggesting that we make the bug fix for the vertex area first and then the rest.

[mark-petersen]

@xylar I see what you are saying. This is the change where invAreaTri1 includes all three kites, while areaDual only includes the kites from active cells. This would be answer changing, and the changes in the loop above could go into a single PR.

[mark-petersen]

The changes in this loop are what is needed to include wall_slip_factor. For this section, it is true that wall_slip_factor = 0.0 should only change the solution due to order of operations. For a second PR for just the addition of wall_slip_factor and these lines, it will be easier to show that only minor changes occurred.

[sbrus89]

@dengwirda, is there a reason to exchange loop order here? It seems like the order of operations issue could be avoided by applying the (1.0_RKIND - wall_slip_factor * boundaryVertex(k, iVertex)) factor in a separate loop k loop following the original do i = 1, vertexDegree loop.

[dengwirda]

@mark-petersen, @xylar, yes I think you're both right here --- the wall_slip_factor change by itself should hopefully be non-climate-changing (though perhaps not BFB due to order-of-op type differences) but the boundary vertex thickness issue will change the way vorticity is computed in general --- at all coastlines for the surface layer, and at any cells/duals where layers are maksed against bathymetry in the deep ocean.

One reason I've grouped this work together here is that implicitly both changes relate to the way momentum boundary conditions are imposed. The wall_slip_factor should allow users to interpolate between no-slip and free-slip conditions explicitly, but the vertex thickness change will alter the way PV is computed on boundaries in all cases. In the example I have above, it looks like we may currently be computing too small a boundary vertex thickness (due to the masked kites), which would give too large a boundary PV value. Considering that free-slip has PV=0 on boundaries, and that no-slip uses a full PV value, this may mean that some kind of more-than-no-slip condition is being used currently. What do you think?

[xylar]

@dengwirda, I see the logic in including them together.

I think it will still be simpler in terms of how to frame these changes to the broader project if we first introduce a bug fix to the area weighting of layer thickness on vertices, knowing that the changes will not be small but that it's a bug. Then, we can still point to that PR when we introduce the other changes here but we would point out that the change in behavior after the bug fix is small as long as you are using the no-slip boundary condition. We would further point out that the "stealth" option to have free-slip or partial slip boundary conditions would be expected to introduce changes of a similar magnitude a and in a similar direction to the previous bug fix. If testing is sufficient at that point, we could also advocate for making something other than no-slip the default at that point.

But I think it really will be necessary to introduce these changes in those 2 stages, rather than in one PR.

[mark-petersen]

I'm happy for these changes to move to E3SM-Project. Thanks.