Modern interpolation methods for unstructured grids.

[lewisn-met]

Is your feature request related to a problem? Please describe.

While the existing finite element interpolation method for unstructured grids implemented in Atlas is sufficiently accurate for most use cases, the algorithm chosen has a number of downsides, namely,

• It is prone to C1 discontinuity: The approach to interpolating on quads involves dividing the quad into two triangles, however, the result of the interpolation is generally not smooth and of lower accuracy along the dividing line.
• It is not invariant under cycling of polygon vertices: if the mesh is read in a different order, the result of the interpolation may change, generating a reproducibility issue. Furthermore, the orientation of the C1 discontinuity varies under this cycling, so the difference can be quite stark (see below).
• It is not well adapted to the sphere: the current method does not account for the curvature of the sphere in the interpolation and interpolates bilinearly along straight lines in the plane, as opposed to great circle geodesics.

Describe the solution you'd like

I propose that the current method be supplemented by a more modern alternative developed for interpolation on the sphere, such as the spherical mean value interpolation described by M. Floater. This algorithm is

• Universal across convex spherical polygons: The same formulae can be applied to both the quad and triangle (and even pentagons, if that is ever an issue which needs addressing), without having to divide the polygon in such a way that C1 discontinuities might be introduced.
• Is invariant under the cycling of polygon vertices: The same result is obtained regardless of the order in which the mesh is read.
• Is well adapted to the sphere: In a preliminary test suite interpolating real spherical harmonics up to 4th order on the cubed sphere, spherical mean value coordinates were found to be 33% more accurate than the current interpolation method, since they account for the curvature of the underlying surface.

Note that, while spherical mean value coordinates are more accurate at interpolating varying functions on the sphere, they are not suitable for interpolating uniform fields, where values at each of the cell vertices are equal. This is due to the weights not satisfying the partition of unity constraint, causing the curvature to dominate and the interpolant to increase near the centre of the cell. As such, the weights would need to be normalised or the existing method used in these cases.

Describe alternatives you've considered

Other interpolation methods, such as (spherical/planar) Wachspress coordinates have been considered, however if the goal is increased flexibility and accuracy, spherical mean value coordinates seem to be the most suitable.

Additional context

Images from preliminary testing.

The spherical mean value (smv, right) interpolation method has visibly improved accuracy over the currently used (barycentric, left) method when applied to spherical harmonics (such as f=xy) on the cubed sphere. Yellow peaks in error occur on the lines along which the quads making up the cubed sphere mesh are divided into triangles.

Illustration of the C1 discontinuity arising in the interpolation and error when using the existing interpolation method.

Note how the C1 discontinuity rotates when the vertices of the polygon are cycled before being fed into the algorithm.

Organisation

Met Office

[wdeconinck]

Hi @lewisn-met this looks like a very nice feature and I'd be happy to see this be developed!