Potential bug in the reconstruction of MPAS fields from edges to cell centers

[cbegeman]

The reconstruction from edges to cell centers as written in MPAS-Framework assumes that the values at all boundary edges are zero. There are two MPAS-O fields that use this reconstruction, sshGradient and surfaceVelocity (and their constituents sshGradientZonal, sshGradientMeridional, velocityZonal, and velocityMeridional).

The zero normalVelocity boundary condition at edges is desirable. It is less clear that zero ssh gradient boundary condition is desirable, particularly for sea ice dynamics.

I have run idealized MPAS-Ocean standalone tests to expose the issue, which are detailed in a comment below and show a reduction in the amplitude of the reconstructed field on the order of 10-50% at the boundary cells, with some spatial variability that is likely related to the orientation of the current relative to boundary edges.

MALI also calls the MPAS-Framework routine mpas_reconstruct but I am unsure about whether the zero-valued Dirichlet boundary condition is desirable for those fields.

[cbegeman]

I have set up idealized MPAS-Ocean standalone cases that have a velocity field with an analytic solution and use the config option config_prescribe_velocity along with disabling the tendency terms. This results in the normalVelocity field being reset to the initial condition right before it is reconstructed from edges to cell centers. I have verified that this is the case by outputting the normalVelocity field after 1 time step. We can then compare the reconstructed velocity field on cell centers with the analytic solution.

Fractional reduction in current speed on cell centers in a planar case with velocity parallel to the non-periodic boundaries at the top and bottom of the domain with period boundaries at the left and right:

Fractional reduction in current speed on cell centers in a spherical case with land cells culled around an ice shelf domain (isomip_plus configuration):

[cbegeman]

@proteanplanet I wanted to make you aware of this because it could affect sea ice motion near land (but not land ice).

[xylar]

I have a hunch about what could be causing this.

First, I investigated coeffs_reconstruct and they are not any different at cells with boundary edges than those without.

Second, I can reproduce the problem using the MPAS-Tools vector reconstruction tool: https://mpas-dev.github.io/MPAS-Tools/stable/vector.html#vector-reconstructions I get a good reconstruction of a perfectly zonal velocity except at boundaries.

My hunch is that angleEdge is wrong for about 50% of boundary edges after culling the mesh. If my hunch is correct, this happens because angelEdge simply gets copied from the base mesh: https://github.com/MPAS-Dev/MPAS-Tools/blob/master/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp#L1444 but the new cellsOnEdge variable always gets set up so the first cell at a boundary edge is valid and the second one is invalid: https://github.com/MPAS-Dev/MPAS-Tools/blob/master/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp#L1272-L1280 In cases where the removed cell on a boundary edge used to be the first entry in cellsOnEdge, angleEdge is off by 180 degrees. The more boundary edges on a cell, the higher the odds of getting a flipped angleEdge.

In practice in E3SM simulations, normalVelocity is always zero at boundary edges. I haven't verified it but I believe gradSSH must also be because we are missing one of the 2 cells that would be needed to compute it. So I don't think this has any impact currently except in idealized tests with prescribed normal velocity fields (including at boundary edges).

[xylar]

Still, this should definitely be fixed, presumably in the cell culler.

[cbegeman]

@xylar I confirmed that the solution is identical whether the velocity on the boundary edges is zero or non-zero. So whatever the issue is, it would not affect the velocity reconstruction. I'm updating the header.

[xylar]

@cbegeman, I don't see anywhere in the velocity reconstruction that boundary edges get zeroed out.

In my offline test, I get very different results when I zero out the normal velocity at boundary edges than when I don't. So this confirms my earlier finding that coeffs_reconstruct are not responsible.

I believe you said on Slack the normal velocity you are using is being created right before the call to mpas_reconstruct()? Can you confirm that? To me, the most likely place where normal velocities at edges are getting zeroed out is elsewhere in MPAS-Ocean (e.g. between getting read in from an initial condition and when reconstruction gets called) so I wanted to rule that out.

In the meantime, I'll keep looking into it on my end.

[xylar]

Regarding the computation of gradSSH, the boundary edges are getting zeroed out when we multiply by edgeMask: https://github.com/E3SM-Project/E3SM/blob/f50d2639c68c1607a521407b420af040080c1380/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F#L2988 This is important because otherwise we would be trying to use the value of pressureAdjustedSSH from a dummy cell in the gradient computation at boundary edges.

It appears to me that the vector reconstruction code would handle a nonzero gradSSH value on boundary edges if we were to come up with a way to produce one. But I don't currently have time to play with modifying MPAS-Ocean itself to be sure that this is the case.

[cbegeman]

@xylar Thanks for looking into this.

I think I found where velocity is being set to 0 on boundary edges on init: https://github.com/E3SM-Project/E3SM/blob/f50d2639c68c1607a521407b420af040080c1380/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F#L3188-L3194

The velocity is reset here: https://github.com/E3SM-Project/E3SM/blob/f50d2639c68c1607a521407b420af040080c1380/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F#L2798-L2814

Then mpas_ocn_diagnostics is called with time level = 2, but I don't see normalVelocity getting zeroed out there.

Then mpas_reconstruct is called.

[xylar]

@cbegeman, okay so here is my summary of my understand of things right now:

• there is likely a bug in angleEdge at boundary edges that comes from the cell culler and we should probably fix this.
• in MPAS-Ocean, we seem to use mpas_reconstruct only on normalVelcoity and gradSSH, both of which are explicitly set to zero at boundary edges, so the angleEdge bug probably doesn't affect anything currently in MPAS-Ocean
• if we wanted to prescript nonzero normal velocity or grad ssh at boundaries, it seems like there is nothing in the mpas_reconstruct() that would prevent this, it just might be difficult to come up with sensible values.

[cbegeman]

@xylar I think that's a good summary. I'll assume you'll take on the first bullet unless you say otherwise. With respect to the third bullet, I'll take responsibility for chatting with @proteanplanet and others about the design and priority level of setting a different boundary condition for gradSSH.

[xylar]

@cbegeman, yes, exactly. I'll take on the first bullet. It shouldn't be too hard.

[proteanplanet]

@cbegeman and I just talked about this, and I think that changing the current method is a low priority. So long as we know what's being done, and it can be expressed neatly as a boundary condition - which it can - I'm fine with the current method given there is no evidence it's adversely affecting sea ice in the V3 code base.

[xylar]

I'm fixing the issue with angleEdge in culled meshes in https://github.com/MPAS-Dev/MPAS-Tools/pull/580.

A preliminary test shows that it flips angles on a planar mesh as expected. I don't have an easy way to get the vector reconstruction coefficients for that mesh so I'm rerunning the compass ISOMIP+ test with the spherical mesh to make sure the problem with reconstruction at boundary edges is, in fact, fixed.