journals.ametsoc.org

Mimetic Interpolation of Vector Fields on Arakawa C/D Grids

️Wolfgang Hayek
️Tue Jan 01 2019

1. Introduction

Interpolation is ubiquitous in climate and weather models. It is used to initialize simulations, in postprocessing applications where output data need to be remapped, to imprint atmospheric fluxes onto an ocean grid and vice versa, to advect tracers, to compare model output, and in visualization tasks.

A number of software libraries have been developed in the past decades to address these needs, for instance, the ESMF library (Hill et al. 2004). These libraries typically implement linear/bilinear/multilinear and conservative interpolation (Jones 1999). Conservative interpolation, also sometimes called area-weighted interpolation, is often applied to water or energy, quantities that need to be conserved globally over the entire planet. Conservative interpolation is more computationally intensive than linear interpolation, as one needs to compute the overlaps of target cell areas with source cells. In spite of this, conservative interpolation has gained significant traction and is probably the most widely used interpolation method for regridding purposes.

A common misconception, in our view, is that the interpolation method (linear or conservative) can be chosen freely. Our approach follows Pletzer and Fillmore (2015), where it is argued that the field staggering determines the interpolation. Vector fields on Arakawa C/D grids require interpolation methods that are neither linear nor conservative, but borrow aspects from each. One outcome of this paper is to clarify when linear interpolation, conservative interpolation, or their vector variants should be used.

Arakawa and Lamb (1977) described the different placement of fields within cells; the staggerings of interest here are referred to as Arakawa C and D. In Fig. 1, the Arakawa C/D vector field (u, υ, w) components are shown on a two-dimensional, horizontal grid cell. This picture is adequate for the grids originally considered by Arakawa and Lamb, typically uniform latitude–longitude-based grids. Such grids cause numeric problems near the poles, and most climate/weather atmospheric and ocean models have since moved or are in the process of moving to grids that have no geometric singularity over the domain of interest. Such grids are characterized by cells that are curvilinear in latitude–longitude coordinates, with sometimes nonorthogonal angles between cell edges. To extend the meaning of Arakawa staggering to such grids, we define Arakawa C and D as a discretized vector field where the data are attached to cell faces and edges, respectively.

In addition to edges and faces, fields can also be attached to nodes and cells. For readers familiar with exterior calculus (Frankel 2012), the above nodal, edge, face, and cell fields correspond to the 0-, 1-, 2-, and 3-forms of differential geometry.

A 0-form is a scalar function of space, whereas a 1-form can be thought of as a vector field. A 2-form has as many components as a vector field in three dimensions; however, it behaves differently under a mirror reflection, prompting some to refer to it as a pseudovector or axial vector. A pseudovector arises when taking the cross product of two vectors. Taking the dot product of a pseudovector with a vector yields a pseudoscalar, a single component quantity that behaves in many respects like a scalar but changes sign under mirror reflection. Pseudoscalars are typically volume densities, and these map to 3-forms in differential geometry.

Table 1 summarizes the correspondence among field types, field staggering/placement, and interpolation methods. Software is widely available for nodal/bilinear interpolation, and a selected list of packages exists for conservative/cell-centered interpolation (e.g., SCRIP and ESMF). The authors are not aware of any openly available software package capable of conservatively interpolating vector fields with staggered components (rows 2 and 3), hence the focus of this article.

Table 1.

Different interpolation methods are required for different types of fields (scalar, vector, pseudovector, and pseudoscalar) in 3D. Each type of field has its own staggering/field placement within grid cells.

Table 1.

The aim of this paper is to develop interpolation methods that are suitable for vector fields attached to edges and faces. The ideas behind vector interpolation go back to Whitney (1957), who introduced what are now known as Whitney basis forms on triangular grids. Similar basis functions for vector fields have been independently developed multiple times and are generally known as and mixed finite elements (Nédélec 1980; van Welij 1985). Generalizations of those basis functions to n-dimensional hexahedral cells have been used by Pletzer and Fillmore (2015). Mixed finite elements were designed to avoid numerical instabilities/pollution and other spurious numerical artifacts arising in the discretization of operators (Ferziger and Perić 2002), and this led to the emergence of “mimetic” discretization methods such as discrete exterior calculus (Hiptmair 2001; Hirani 2003; Bossavit 2005) and finite element exterior calculus (Cotter and Shipton 2012; Cotter and Thuburn 2014). Mimetic methods aim to preserve the mathematical identities .

The novelty of this paper is to apply vector interpolation methods to Earth science grids and demonstrate how these methods satisfy mimetic properties on a cell-by-cell basis and for fractional cells. A by-product of these methods is that the total vorticity/flux budget over lines and surfaces is conserved. Apart from conservation, the methods naturally handle grid distortion and data masking; they support fields that are partially valid over edges and faces.

Our focus here is on two-dimensional grids at the surface of the Earth. Owing to the geometry of planets, most three-dimensional Earth grids are constructed from two-dimensional surface grids that are extruded vertically. When the target lines are at constant elevation or the target surfaces are lateral (i.e., tangential to the vertical columns), three-dimensional interpolation then reduces to two-dimensional interpolation with a separate vertical interpolation step. In other words, horizontal and vertical interpolation weights are computed independently. The split between horizontal grid and vertical axis can be exploited by software to reduce memory footprint and execution time, and this warrants a careful treatment of the two-dimensional case.

This paper is organized as follows. In section 2, the interpolation weights for two-dimensional Arakawa C/D grids are derived. In section 3, the mimetic and conservation properties are proven. In section 4, edge-/face-centered interpolation is tested on a variety of grids and conditions that are relevant to Earth science, and in section 5, we summarize our results.

2. Interpolation generalized to all types of field staggering

Before proceeding, some terminology needs to be clarified. Interpolation is often thought of as approximating a function on a set of target points. A more useful definition is to think of interpolation as a projection onto a set of n-dimensional target objects. Depending on the type of interpolation, the target objects can be points, lines, surfaces, volumes, or hypervolumes in higher-dimensional space. This definition of interpolation agrees with bilinear interpolation when the targets are points and with area-weighted (conservative) interpolation when targets are areas in two dimensions or volumes in three dimensions.

a. A single interpolation formula

All interpolation methods can be written as a linear combination

of field values

, with interpolation weights

as coefficients. Here,

is the interpolated field, and σ tags the field elements that contribute to the interpolation. In the case of bilinear interpolation, for example, σ denotes one of the grid nodes in the neighborhood of the target point (Fig. 2, left). For conservative interpolation, σ would be one of the grid cells that overlaps with the target area/volume (Fig. 2, right).

Note that for conservative interpolation, it is often more convenient to have

and

represent integrated values rather than density fields, so

and

should be thought of as total amount of mass, energy, and so forth (i.e., the quantities we would like to conserve). Equation (1) can be made to work for vector field staggerings in a similar way by stating that

and

are line/surface integrals of the field. For a scalar field f, vector field F, pseudovector field G, and pseudoscalar field ρ, (1) becomes

where

, and

are integrals of the interpolated fields

, and

on a line, surface, and volume, respectively. In the following, integrals of the interpolated fields are indicated with an overbar symbol, and these are always scalar quantities. In contrast, the (nonintegrated) interpolated fields are denoted with a tilde, and these can be either scalars or vectors, depending on the field. (The bar and the tilde quantities do not share the same units since the former are integrated over lines, surfaces, etc.) We have a strong preference for using the bar quantities whenever possible; the reason will become apparent later.

The field-integrated values,

on each cell object σ are assumed to be given—they are obtained by integrating the field over the various objects σ that make up the cell (nodes, edges, faces, etc.). For nodal fields, the integral is over a zero-dimensional object and thus simply amounts to evaluating the function on that node.

Naturally, we cannot assume the interpolation weights

to be the same for all types of field α = 0, 1, 2, and 3. We should expect the coefficients

to depend on the choice of basis functions

. By analogy, with bilinear and conservative interpolation, we write

where T is the target: a point

, line, surface, or volume. Expressions (2)–(4) are very general; they apply to all field staggerings in 3D, irrespective of the order of the interpolation scheme.

Let us look at specific examples to see how these expressions work in practice. In the case of nodal interpolation, σ is a grid vertex, and the target is, as mentioned before, . The interpolation weight is then the linear basis function evaluated at the target point, which is the well-known formula .

For conservative interpolation, the sum is over all the source cells σ for which gives a nonzero contribution when integrated over target T. It happens that for the lowest-order conservative interpolation, the basis function is piecewise constant over cells and thus can be taken out of the integral when integrating over a cell. However, it is always present and should be included for higher-order conservative interpolation. Comparing linear with conservative interpolation, it becomes clear that the integral is also present in linear interpolation, but in this case, the integration is over a zero-dimensional point—in other words, a “no operation.” This should make the expressions for and in (4) quite plausible.

To compute the interpolation weights, we need basis functions. Recall that bilinear interpolation basis functions are attached to grid nodes, and these have the property to take the value of one on one of the cell nodes and zero on all other nodes. (High-order basis functions satisfy a similar property on cell vertices and some internal nodes.) This condition is critical since it allows the interpolated field to take the exact value as the target point approaches a grid vertex. Similar conditions are desirable for edge, face, and cell interpolation. The main difference is that the orthogonality conditions now involve integrals:

where

is the Kronecker symbol (

). Equations (5) state that the integrals of the edge, face, or cell basis function

attached to cell object σ evaluate to one when integrated over this object and to zero when integrated on any other cell edge, face, or cell.

b. Basis functions for two-dimensional vector fields

So far, general interpolation formulas for nodal-, edge-, face-, and cell-centered fields were derived in three space dimensions. In the following, only two-dimensional, horizontal geometry will be considered, with emphasis on edge and face interpolation. In the case where the field has a dependence on elevation, the interpolation will involve an additional separate calculation of vertical interpolation coefficients, which will not be described here. The full 3D interpolation weights are then the product of horizontal and vertical interpolation weights. The separation between horizontal and vertical dimensions is warranted if the target runs parallel to the horizontal grid. (A case where the current split between horizontal and vertical dimensions would not work is when the source grid uses pressure levels for its vertical axis and the target grid a terrain-following coordinate, or vice versa.)

Let , be some coordinates that are aligned to cell faces. These can be chosen to be the index coordinates of a structured grid or, more generally, the barycentric coordinates of the cell. In addition, there is a vertical coordinate z, but because our main interest is in two dimensions, all surface integrals will implicitly integrate over the vertical z dimension. Thus, only line integrals in the – plane will survive for face interpolation. (Faces and edges are the same thing in two dimensions, although the basis functions are different.)

Since each edge/face of the cell will hold a basis function, we need a scheme to index edges and faces. Figure 3 shows how a combination of symbols—“0,” “1,” and “*”—can be used to uniquely identify any cell object (node, edge, face, etc.). (For high-order interpolation, additional symbols are required as the cell objects split into subcell objects; e.g., “0,” “1,” “2,” and “*” for second order.) Symbol “0” corresponds to the “low” end of the cell, “1” means high end, and the wild-card symbol “*” indicates that the element varies along that direction. Cell nodes are labeled “00” (lower-left node), “10” (lower right), “11” (upper right), and “01” (upper left). The cell itself is labeled with “**.” As before, superscript “(1)” stands for edge and superscript “(2)” for face interpolation.

The following vector low-order, basis functions can be shown to satisfy (5):

where

is the Jacobian of the cell and

the unit vector pointing out of the page. In the case of orthogonal cells, the vector edge basis functions are tangential to the supporting edge and linearly decaying toward the opposite edge. The vector face basis functions are perpendicular to the face and vanishing on the opposite face. Under more general conditions, Fig. 4 shows the edge and face basis functions for a curvilinear cell with nonorthogonally intersecting edges; see Pletzer and Fillmore (2015) for a representation of the edge and face basis functions in 3D. Notice that in this case, the edge basis function is no longer tangential to the supporting east cell edge. More generally, the direction of the vector edge basis function will turn and adjust so as to be perpendicular to adjacent edges at the cell vertices. This ensures that the line integral of an edge basis function over adjacent edges remains zero; that is, the interpolation picks up line integral contributions from that edge only. The integral over the opposite edge is zero because of the linear dependency of basis function along the perpendicular direction to the edge. Also observe how the vector face basis functions rotate, this time in order to ensure that the flux integrals of the basis function over adjacent faces are zero.

c. Computation of the two-dimensional edge- and face-centered weights

To derive the interpolation weights

in (4), the line target T must be expressed in some form. A possibility is to write the target as a polyline with each segment satisfying

where

and

are the starting point and endpoint of the segment, respectively (see Fig. 5). Equation (8) and

can then be plugged into (4) to get

(see appendix).

First, note that the weights are the same for edge and face interpolation in two dimensions. This means lateral fluxes can be computed using the edge interpolation method. In the remainder, we will use interchangeably the terms edge and face interpolation and line and flux integrals, depending on the physical context, with the understanding that the interpolation method works the same way in both cases.

The weights are the product of two terms. The first term represents the projection of the segment ξ_ib–ξ_ia onto the cell edge. This term is zero if the segment is perpendicular to the edge. The second term is a weighted average in the direction perpendicular to the edge, which gives zero if the segment is on the opposite edge. Together, these two terms ensure that the weight becomes one when the segment occupies the entire edge and zero when the segment is on any other edge.

In (10), we assumed the segments to reside entirely within the cell. Care should be taken to break the target line into subsegments if the line intersects multiple cells. The numerical procedure then involves computing all the intersection points of the target line with the grid, creating cell-contained subsegments between each pair of intersection points.

3. Mimetic and conservation properties of edge and face interpolation

In this section, the previously derived basis functions and weights will be shown to satisfy the conservation property

for the interpolated vector field

given in (2). We start by proving Stokes’ theorem for the cases shown in Fig. 6 where the path integrals (dashed lines) are inside a cell C (solid lines). In the first case (Fig. 6a), the integration path follows the four edges

where

are the projected field values on the edges obtained using (3). By convention, the cell edges’ orientation is from left to right and bottom to top—the minus sign in two of the above terms arises when the edge runs in the opposite direction to the anticlockwise integration path. For the surface integral, we get

where the basis functions

are given by (6):

and

Thus,

evaluates to ±1, and this proves the interpolated version of Stokes’ theorem at the level of a single cell.

In the second case (Fig. 6b), we cut the cell at some

location,

, with

resulting in the full cell. Denoting the smaller cell by

(shaded area), the line integral in Stokes’ theorem becomes

In the above, the horizontal edges

and

have been scaled by λ, and the contribution from

was replaced by a λ-weighted average. The first case is recovered by setting

, while

gives zero as expected. Effectively, the surface integral only covers the fraction of the cell whose area is equal to λ:

so that Stokes’ theorem also holds for a split cell. The third case (Fig. 6c) is equivalent to the second case. In the fourth case (Fig. 6d), the cell is cut diagonally:

where the contribution along the diagonal line, obtained using (10), consists of averages over opposite edges.

We now have all the building blocks in place for constructing more complex paths inside a cell. Consider Fig. 7a, where the integration path (dashed line) cuts across the lower-right corner of the cell. Using the linearity of the edge integrals, the cells can be broken into smaller subcells (dotted lines) where edge integrals and are replaced by and , respectively. This results in a path integral that is equivalent to Fig. 6d.

Next, consider a case with two line segments inside the same cell (Fig. 7b). Using the same arguments as before, the cell is cut such that each line segment lies in a separate subcell with edge integrals and adjusted accordingly. The additional integration paths (dot–dashed lines) along the edges of subcells cancel out in the total integral, enabling a split into two instances of Fig. 6d (or Figs. 6b or 6c if the line segment happened to be horizontal or vertical, respectively) and one instance of Fig 6a. The case of multiple cells (Fig. 7c) can be shown to work similarly.

We now show how the previous results apply to the special case of a vector field that derives from a scalar potential. Let

and

be the projected value of F on a line target with parametric representation

, where

are the starting/endpoints of the line from a to b in parametric space. One can use (1) to approximate this projection while noting that the edge values

are just the difference of potential between the edge’s endpoints in this case (i.e.,

). Denoting by

, and

the scalar potential values at the lower-left, lower-right, upper-right, and upper-left corners of the cell, respectively, we find the integrated projection

after some straightforward but tedious manipulations, where

are the bilinear interpolation weights. Thus, (17) amounts to a difference of potential values linearly interpolated at the endpoints of the line. An important observation is that projection

will be exact if the interpolation of the potential function is exact, and this will be the case whenever the endpoints a and b fall onto grid nodes. Furthermore,

will be exact if a and b are the same points since the interpolation errors will then exactly cancel out to give

. Closed loop integrals of interpolated vector fields deriving from a gradient are zero, consistent with our earlier discussion of Stokes’ theorem.

In summary, Stokes’ theorem holds for interpolated fields integrated along arbitrary collections of straight lines. It should be noted that the integrals along full edges are exact if the are exact, and numerical errors only arise when integrating along partial edges, over partial cells, and in oblique directions. Furthermore, the proof does not assume a particular shape for the cells; it applies to any distorted cells for which a parametric coordinate system , can be defined.

It is instructive to compare the above results with the case where the field F is expanded in bilinear basis functions instead of the mimetic bases. One can show in this case that taken along a straight line segment reduces to , where is the average value of the bilinear interpolation of along the segment. This quantity equals to the difference only when the potential varies linearly (constant ), and this demonstrates that a bilinear interpolation of a vector field will, in general, not be mimetic.

4. Applications

To test mimetic interpolation in two dimensions, we manufactured analytic fields that have zero divergence everywhere, except possibly at one location. In such cases, open path integrals will depend on the endpoints only, and closed loop integrals will be zero if the contour does not contain any field singularity. The grids are either uniform or curvilinear. All results are computed using 64-bit precision software.

a. Divergence-free vector field in Cartesian coordinates

Let

be a divergence-free vector field with streamfunction ψ,

where

are the Cartesian coordinates, and

is the unit vector pointing in the vertical direction. Figure 8 shows vector field v and corresponding contours of the streamfunction ψ.

The flux (e.g., velocity flow) across any path a → b can easily be estimated. Writing , we get . In other words, the flux is the difference of streamfunction values between the endpoints of the path. (The streamfunction plays here a role similar to the scalar potential of section 3.) We use differences of streamfunction values to set the exact flux on cell edges. Any observed interpolation error can then be entirely attributed to the interpolation rather than to the source field.

Families of target paths are generated by discretizing the parametric equations

over n linear segments, where

is the parametric coordinate of the curve (Fig. 8). Path A has a single, straight segment and was chosen to verify that the interpolation method preserves symmetry; the upward flux exactly cancels the downward flux. Path B is a closed loop and hence should return a zero flux, regardless of the number of straight segments (4 or 32). Path C begins and finishes at grid nodes to yield a nonzero, exact flux of 1/20. Path D (not shown) is similar to C but is displaced along the x axis by half a grid cell. From the analysis in section 3, path D is the only curve that is expected to give any numerical error.

The results of the numerical flux integrations on a 20 × 10 grid are shown in Table 2. For A, B, and C, we observe that the numerical integrals are exact within machine accuracy. To estimate the numerical error of the projection of the vector field onto path D, which starts at (0.25, 0.2) and ends at (1.85, −0.3), only the distance of the starting and endpoints of the path to the nearest cell edges needs to be considered—in this case, 0.05 for both the starting and endpoints. The error ε is due to linear interpolation of the streamfunction and is given by the standard linear interpolation error estimate

where h is the cell size, and

is the maximum value of the second derivative of the streamfunction within the cell. We have

, and for path D, whose starting and endpoints are displaced in the x direction from the nearest grid nodes,

at the starting point and

at the endpoint. With these numbers, the upper bound of the error is found to be ≤2.3 × 10⁻³, in good agreement with the observed value of 2.3 × 10⁻³.

Table 2.

Numerical results obtained for the flux integration of paths A, B, C, and D for a vector field that derives from a streamfunction.

Table 2.

b. Vector field with a singularity

To demonstrate that mimetic interpolation works for a singular field, let

be a radial field with singularity located at (0, 0) (see Fig. 9). The 11 × 11 grid includes the point (0, 0). Note that while the vector field’s magnitude

is infinite at the singularity, all the edge integrals are finite, provided the singularity falls inside a cell. By Gauss’ theorem, the closed contour flux integral of

gives

To verify (27), the total flux is computed across a sequence of closed contours

which are broken into 16 equal-length, straight segments. The contours were chosen to be away from the singularity (

and

), to fully encircle the singularity (

, and

), or to cut through the grid cell that contains the singularity (

, and

Figure 10 confirms the mimetic interpolation property—the flux integrals are zero for contours and and one for , , , and , with numerical error comparable to machine accuracy (~10⁻¹¹). Contours , , , and produce fluxes with intermediate values—naturally, the method cannot resolve subgrid features.

Also shown in Fig. 10 are the flux values (dotted line) obtained by applying a trapezoidal quadrature method. Bilinear interpolation was used to evaluate the vector field on the midsegment locations of the trapezoidal integration rule. In contrast to the mimetic interpolation, flux integrals do not amount to zero and one, respectively. A residual error of up to 3% can be observed for contours , , , and . In general, the numerical error of bilinear interpolation will be largest near the singularity.

c. Cubed-sphere grid

To show that the interpolation algorithm works for nonorthogonal cells in spherical geometry, the velocity field is set to

with streamfunction

where λ is the longitude (°), θ is the latitude (°), and r is the radial distance from the center of the coordinate system. A cubed-sphere grid similar to that used by atmospheric dynamical cores (Lin and Rood 1996) was generated using LFRic mesh generation tools (Melvin et al. 2017). As in the previous tests, the exact edge values were set using the difference in streamfunction values between the two endpoints of the integration path.

Figure 11 shows the grid laid out in latitude–longitude space for a low-resolution case with the two flux integration paths:

and

for

. The orientation of the paths is from lower left to upper right for both A and B. Two points were used to define path A and 11 points for path B. Note that the paths cross multiple panels of the cubed-sphere grid and that many cells are highly distorted in latitude–longitude space.

In the computer program, the grid was represented as a collection of individual, quadrilateral cells with each cell having its own set of latitude–longitude vertices. The cells are disjoint; there is no connectivity between cells. A benefit of this approach is that the number of edges and vertices relates simply to the number of cells. Another advantage is that neighboring cell vertices may take different values, and this can be used to add and subtract a period to ensure that all cells have a well-defined, positive orientation.

The octree-based class vtkCellLocator from the open-source Visualization Toolkit (Schroeder et al. 2006) was used to quickly locate cells that are intersected by edges. When a target segment runs along two neighboring cells, the intersection between the segment and cell edges produces a singular system; we detect this collinearity and assign one or more intersection points to record the overlap between segment and edge.

It is found that the open path flux integration error, shown in Fig. 12, decreases quadratically with the average gridcell size , where N is the number of grid cells. The error depends on the proximity of the start and endpoints from a grid vertex, and this can cause the convergence of the error to fluctuate depending on where the start and endpoints of the path fall within grid cells. Such behavior can also be observed with bilinear interpolation where the error depends on the distance between the target point and the nearest vertex.

Global flux conservation depends on the ability of the interpolation to conserve fluxes at the level of each target grid cell. Figure 13 shows the maximum loop-integrated flux error obtained after regridding vector field (29) from a uniform, latitude–longitude grid to a cubed-sphere grid. In the process, each edge-integrated field on the source grid is projected onto the target gridcell edges.

Apart from the accumulation of floating point roundoff errors, mimetic regridding produces closed flux conservation errors <10⁻¹³ for all source grid resolutions. In contrast, the dashed line shows the error from bilinear regridding; that is, treating the edge field as if it were nodal, which amounts to bilinearly interpolating each edge value at the starting point of the edge. Bilinear interpolation yields an error proportional to the gridcell size . As the source grid resolution increases, the bilinear regridding error decreases but cannot match the mimetic error by many orders of magnitude for the considered source grid resolutions (down to 4 km on a global grid). It would take ~10¹⁵ source grid cells or a resolution of a meter for bilinear interpolation to achieve comparable flux conservation accuracy to the mimetic method.

5. Summary

Edge- and face-centered interpolation was presented in two dimensions with segmented lines as targets. Although the basis functions for edge and face interpolation are different, the interpolation weights turn out to be identical in two dimensions, and hence only one type of interpolation needs to be implemented.

In contrast to bilinear interpolation, the proposed mimetic interpolation method is suitable for vector fields. Mimetic interpolation satisfies Stokes’ theorem for arbitrary, piecewise linear paths, and the method produces exact line/flux integrals when the vector field is divergence free and the start and endpoints fall onto grid nodes. The method is immune to pole-like singularities, and line/flux integrals automatically take care of partially valid cells—there is no need to introduce masking arrays as the assigned fluxes take split edges into account. Finally, the method handles curvilinear coordinates with no modification.

In three dimensions, a face-centered field target becomes a surface, and the implementation becomes more complex. One then needs to break the area target into triangles that are fully contained inside the source grid cells. This requires collecting the intersection points of source grid edges with the target face and the intersection points of target edges with the source grid faces and finding the source grid cells that contain the target nodes. From all these points, a triangulation can be assembled, which has the property that the triangles are entirely contained within source cells. Formulas similar to (10) can be derived for the computation of the interpolation weights and then applied to each cell-contained triangle.

It remains to be seen how edge/face interpolation can be generalized to higher-order elements. There are some indications that future atmospheric and oceanic dynamical cores will employ pseudospectral elements. It is too early to assess the impact of high-order elements on the presented method, although this will hopefully just affect the integrals to be computed. Provided the basis functions are polynomials of the parametric cell variables, the integrals can be computed exactly, a property that is critical for numerical conservation.

It is our hope that the presented method will find new applications in several areas. Because climate change rests on small discrepancies (on the order of 0.1%) in the radiative energy budget, climate scientists will be among the first to benefit from mimetic edge and face interpolation. Second, with the increased need to share datasets across institutions, edge/face interpolation might provide independent measures of the accuracy of numerical data (e.g., testing the conservation of vortex tubes in Rossby waves). Third, the method can be used to compute fluxes across irregular, political boundaries to provide global estimates of how much water and energy enters and leaves a region in a given interval of time and how such quantities might change under different greenhouse gas emission scenarios. Finally, edge/face interpolation could provide a basis for advecting vector fields based on (pseudo-) Lagrangian approach.

Software implementing the methods described in this paper can be downloaded from https://github.com/pletzer/mint.

Acknowledgments

The authors thank Jorge Bornemann and the anonymous reviewers for their comments and for suggesting improvements to the manuscript, Samantha Adams and Mike Hobson (Met Office) for supplying grids and a cubed-sphere mesh generation tool, and the LFRic project for supporting this work. The authors also wish to acknowledge the contribution of NeSI high-performance computing facilities to the results of this research. New Zealand’s national facilities are provided by the New Zealand eScience Infrastructure and funded jointly by NeSI’s collaborator institutions and through the Ministry of Business, Innovation and Employment’s Research Infrastructure program (https://www.nesi.org.nz).

APPENDIX

Derivation of Some Interpolation Weights

We start by computing weight

on the west edge of the cell. From (4) and (6), we have

where the integration is along the target line (8) passing through the points

and

. First, note that

. Further replace

with

and likewise for

to get

which reduces to the first expression of (10).

The weight on the east edge

can be computed similarly to yield

To compute the weight associated with the east face

we write

for the east face and note that

after using (7). Hence, face weight

has the same expression as edge weight

REFERENCES

Ferziger, J. H., and M. Perić, 2002: Computational Methods for Fluid Dynamics. Springer, 389 pp.
- Crossref
- Export Citation
Frankel, T., 2012: The Geometry of Physics: An Introduction. 3rd ed. Cambridge University Press, 684 pp.
Hiptmair, R., 2001: Higher order Whitney forms. Geometrical Methods for Computational Electromagnetics, F. Teixeira, Ed., EMW Publishing, 271–299.
- Crossref
- Export Citation
Schroeder, W., K. Martin, and B. Lorensen, 2006: The Visualization Toolkit: An Object-Oriented Approach to 3-D Graphics. 4th ed. Kitware, 528 pp.
- Crossref
- Export Citation