Higher-order conservative interpolation between control-volume meshes: Application to advection and multiphase ﬂow problems with dynamic mesh adaptivity

A general, higher-order, conservative and bounded interpolation for the dynamic and adaptive meshing of control-volume ﬁelds dual to continuous and discontinuous ﬁnite element representations is presented. Existing techniques such as node-wise interpolation are not conservative and do not readily generalise to discontinuous ﬁelds, whilst conservative methods such as Grandy interpolation are often too diffusive. The new method uses control-volume Galerkin projection to interpolate between control-volume ﬁelds. Bounded solutions are ensured by using a post-interpolation diffusive correction. Example applications of the method to interface capturing during advection and also to the modelling of multiphase porous media ﬂow are presented to demonstrate the generality and robustness of the approach. of an immiscible displacement using mesh adaptivity to accurately resolve the shock front and the growth of viscous ﬁngers. The simulation was carried out using a higher-order CVFE simulator The simulator solves a modiﬁed form of the equation system above where a new velocity ( u α = q α / S α ), distinct from the Darcy velocity, is solved for. Velocity is eliminated from the equation system and the resulting pressure-saturation equation is then solved. The simulations were run with a P2DG-P1DG element (second order discontinuous ﬁnite element representation for velocity, ﬁrst order discontinuous ﬁnite element representation for pressure). The numerical solution is interpolated between meshes after each adapt. the ﬁrst whilst the value 7000 corresponds to the approximate number of elements at the end of the simulation. The number of elements in the initial mesh is approximately 20000 whilst the domain ﬁxed mesh we compare the adaptive results against has approximately 40000 elements.


Introduction
Dynamic, adaptive meshing is often used during numerical simulations of transient fluid flows to improve accuracy [1,2]. The mesh is refined in regions where properties are changing rapidly in space and (usually) coarsened in regions where the properties change more slowly in order to improve accuracy whilst minimising computational effort. The mesh may change every time-step depending on the error metrics used to control refinement and coarsening. This inevitably means that data must be mapped from one mesh to another using an interpolation algorithm [3]. The interpolation must be done efficiently, to minimise computational overheads, and in such a way that the field being re-meshed remains bounded. The method must also be conservative when dealing with fields constrained by a conservation law (e.g. mass).
It is advantageous to model transient flow problems using a Control-Volume Finite Element (CVFEM) formulation as this helps enforce mass conservation. Mesh adaptivity for control-volume fields may be used to better capture shocks and interfaces between fluids as well as material property boundaries in both inertial and porous media flow [4][5][6][7][8]. In the latter example, discontinuous representations of control-volume fields (that is to say when the dual 3 finite element representation is discontinuous) may further help to capture abrupt changes in material properties (e.g. rock permeability) [9,10].
Whilst efficient techniques already exist in the literature for the interpolation of finite element fields [11,12], as will be explicitly demonstrated in this paper, existing methods of control-volume interpolation produce poor results when applied to discontinuous fields. The interpolation introduces a large error into the solution which can become irreparably distorted, highlighting the need for a new method without these deficiencies.
The paper is structured as follows. In Section 1, existing methods of mesh to mesh interpolation are reviewed outlining their strengths and weaknesses. Having motivated the need for a new method of control-volume interpolation, Section 2 discusses the theory underlying the new method, including the implementation of a bounded solution correction. Section 3 provides three distinct applications of the method to interface capturing problems, finishing with an application to viscous fingering in porous media flow, in order to demonstrate the robustness of the method. The examples demonstrate scenarios where other methods of control-volume interpolation can fail to give acceptable solutions. Finally, some key conclusions drawn from the work are presented in Section 4.

Node-wise interpolation
The simplest approach to interpolation, referred to here as node-wise interpolation proceeds by setting the values of a field at the nodes on the new (target) mesh to be equal to function values at those physical positions on the old (donor) mesh [13]. This method of interpolation is sometimes known as collocation interpolation [13] or consistent interpolation [11] (1) where N old j are basis functions (finite element or control-volume depending on the nature of the field ψ old ) defined on the old mesh. An illustration of the node-wise interpolation procedure for a 1D (linear) finite element field is shown in Fig. 1.
Node-wise interpolation has several drawbacks. In particular, it is not conservative in the sense that integrals (volume or surface) of the resultant field are not in general the same as integrals of the original field when taken over the same region of physical space. In addition, it has been demonstrated that node-wise interpolation does not in general preserve the extrema of the interpolant [12,14] and moreover is unsuitable for application to discontinuous fields that are not pointwise defined without substantial modification to the procedure [15,16]. Finally, for higher-order discretisations, node-wise interpolation does not guarantee bounded solutions.

Galerkin projection/interpolation
Galerkin interpolation is defined as the interpolation method that is optimal in the L 2 norm [11]. In addition to being conservative, the Galerkin method circumvents many of the above issues that afflict node-wise interpolation and in particular is applicable to both continuous and discontinuous interpolants [3,17]. Consider a field ψ old defined on the old mesh, and ψ , some interpolation of this field that is defined on the new mesh. The requirement of maximal accuracy in the L 2 norm is that ∂ ∂ψ i ψ old − ψ L 2 = 0 (where the ψ i are the coefficients of ψ when expanded in target mesh basis functions) or equivalently: (2) Expanding ψ = i N i ψ i then leads to the defining equation of a Galerkin projection: The issue that arises in the evaluation of this integral is that the donor mesh basis functions are only guaranteed to be (continuous) polynomials over each element of the donor mesh. Over elements of the target mesh, the basis functions may be discontinuous or piece-wise defined. If the integral of the product of the basis functions is evaluated at Gauss quadrature points on the target mesh, the integral will not be exact, leading to a loss of accuracy and conservativeness in the method.
The solution to this proceeds by constructing a geometry known as the supermesh defined as the intersection of the donor and target meshes. On the elements of the supermesh, the donor mesh basis functions are, by construction, continuous polynomials (since the supermesh is an intersection of donor and target meshes and the donor basis functions are continuous polynomials on the donor mesh), allowing high accuracy evaluation of this integral [11].
In the case of interpolating a field that is defined on a general control-volume mesh, direct construction of a supermesh is challenging. Heuristically this is because in the (3D) finite element case, the supermesh is usually given by the intersection of two tetrahedral meshes which is itself tetrahedral, whilst in the control-volume case it involves the intersection of more complex polyhedra, the classification of which is mathematically challenging. The approach must therefore be adapted in order to easily apply it in the control-volume case.

Grandy interpolation
The next method of interpolation we consider is Grandy interpolation [18]. It is a conservative approach that is often effective but can be highly diffusive due to its low-order nature. Grandy interpolation proceeds by mapping from donor to target mesh by calculating the volume of intersection of overlapping polyhedra on the donor and target meshes. The method assumes the interpolant to be constant across a donor element and it is this feature that accounts for its diffusive nature. The calculation of the volume of intersection of overlapping polyhedra in this case is an example of the construction that must be done in order to build a control-volume supermesh. The Grandy method may therefore be thought of as an example of performing Galerkin projection to a control-volume field ψ old in Eq. (3) but where the donor field is constant across donor elements. Since the method can be regarded as a special case of a Galerkin projection it is conservative.

Control-volume Galerkin interpolation
In this section, the new control-volume interpolation method proposed in this paper is explained in detail. Section 2.1 presents the method as a three step procedure that begins by projecting the control-volume field onto a finite element representation on the donor mesh. The resulting finite element field is then projected to a finite element field on the target mesh by assembling a supermesh. Finally, this field is transformed back into a control-volume representation on the target mesh. Having described the interpolation method, Section 2.2 describes how one can obtain a bounded solution after the interpolation is complete.

Three projections
Fields defined on control-volume meshes may be interpolated using a three step procedure that is now discussed in detail.

Step one -projection to a finite element representation on the old mesh
First the control-volume (CV) field on the old mesh is mapped via Galerkin projection onto a finite element (FEM) representation within an element which may be either a continuous (CG) or discontinuous (DG) Galerkin representation: A graphical representation of the projection defined by Eq. (4) is shown in 1D in Fig. 2 and in 2D in Fig. 3.

2.1.2.
Step two -projection from old to new mesh by building a supermesh The resulting FEM representation ψ old FEM interpolating the control-volume field ψ old CV can then be mapped onto the new mesh via a strictly FEM Galerkin projection, greatly simplifying the implementation of this mapping: Solving Eq. (6) for ψ new FEM requires computation of the mixed finite element mass matrix (MM) new-old ij = N new FEM i N old FEM j dV which is a complex problem but can be tackled by constructing a supermesh from the intersection of the old and new meshes. This procedure was outlined in Section 2.1, but the reader is referred to [12] for a detailed account of global supermeshing, and [11] for a discussion of local supermeshing and also to [19] for a more detailed account of supermeshing as a whole. The mapping is greatly simplified as it now only requires calculation of intersections of finite element meshes as opposed to control-volume meshes or a mixture of the two.

Step three -projection back to CV representation on the new mesh
This concludes the three step interpolation procedure. Note that there are advantages to using a DG FEM representation within the mapping rather than a continuous FEM representation as it is more accurate and does not propagate Gibbs oscillations beyond an element as there is no communication between elements. Note also that the complete mapping (Eq. (4)-Eq. (7)) is conservative and does not change the interpolant (that is ψ new CV = ψ old CV ) if the old and the new mesh are the same. This statement follows by considering (Eq. (4)-Eq. (7)) when the new mesh and old mesh are the same: and thus M new , k is an element of the set of CVDG's (discontinuous sub-control-volumes) that are contained inside control-volume CV i .

Implementation of bounded solution corrections
Despite having many advantages over other interpolation methods, interpolation of a control volume field by way of Galerkin projection onto a finite element basis will not necessarily result in a solution that is bounded. This issue also arises in node-wise interpolation when the order of the finite element basis functions is quadratic or higher and so is not unique to the Galerkin method. We now outline a procedure whereby given a solution to the above interpolation, ψ new CV , one can obtain a new solution that is bounded.
In order to measure the deviation from boundedness, begin by defining a field ψ new CV devk control-volume wise such that at each control-volume k: where boundedness is achieved when ψ new within the solution bounds (ψ new CV mink , ψ new CV maxk ), the control-volume does not need mass adjustment in order to obtain a bounded solution. We describe such a control-volume as having 'absorptive capacity'. The strategy to obtain a bounded solution is then to introduce diffusion to the deviation field in order to spread non-zero deviation into neighbouring control-volumes. Iterating this procedure should then lead to a state where all control-volumes are at absorptive capacity. This process must be done in such a way that conservativeness is maintained. This diffusion algorithm in now described in greater detail.
This operation is trivial since by construction, the mass matrix M L is diagonal (see below). Since the distributed mass matrix The interpolant vector ψ new CV is then modified using: . (13) Note that the above operations have no effect on the integral of ψ new CV and hence conservativeness. The mass matrices in Eq. (11) are defined as: and where m CV i is the volume of control-volume i and m SURi contains the sum of the volumes of the surrounding controlvolumes: and neig{i} labels the control-volumes that share a face with control-volume i, including itself. One concern with the iterative scheme we have presented above is that it substantially adjusts control volumes that are already bounded. This problem can be lessened by the following modification. When iterating over control-volumes, for a given control-volume k, we identify the neighbouring control-volumes to k with the maximum positive and the minimum negative defect adjustment needed and mass is exchanged exclusively between these until no further mass adjustment can be made for k. One then proceeds to control-volume k + 1. More precisely: 2) For control-volumes i and j set: 3 The mass adjustment associated with step 3 above takes the control-volume with the minimal needed mass adjustment, alt j m CV j and control-volume i otherwise, and adds its mass deficiency (negative or positive mass) to the other control-volume. The effect of this is to eliminate the donor control-volume's mass deficiency and reduce the mass deficiency of the receiving control-volume so that either ψ new alt i ← 0 or ψ new alt j ← 0.
In practice, the algorithm above (Eq. (17)-Eq. (20)) may not be able to eliminate all mass deficiencies and it is therefore used together with the mass distribution approach (Eq. (11)-Eq. (13)). By iteratively solving Eq. (17)-Eq. (20) until no further adjustment can be made and then solving Eq. (11)-Eq. (13), and repeating this two-stage process until the maximum control-volume field adjustment is small (e.g. 10 −6 in the examples in this paper), one can obtain a bounded solution after interpolation.

Interpolation during mesh refinement
Having discussed the theoretical background underlying the proposed new method of control-volume interpolation, a series of applications of the method are now considered.
The first three examples are to advection problems. We consider the dynamics of a scalar field ψ that is passively transported with the (prescribed) bulk velocity u of the underlying fluid according to the advection equation: Table 1 Parameters used in the static scalar simulation described in the main text. The variable t is physical time. All measures of distance in this paper may be regarded as dimensionless (hence the absence of units in the specification of domain size and minimum/maximum edge-length).

Size
Initial ∂ψ ∂t The mesh adapts to the curvature of the scalar field ψ which is represented control-volume wise and the behaviour of different methods of interpolation used to map this scalar field between meshes are compared. The scalar field ψ is defined on a control-volume mesh that is dual to an n'th order finite element pressure mesh that may be discontinuous or continuous. The pressure representation will be referred to using the notation PnDG in the former case and PnCG in the latter case. The value of n used, as well as the mesh continuity will be stated in each example in turn. The finite element pressure representation is used to calculate the higher-order fluxes for advection as detailed in [4,6].
As a first example, a case where the bulk fluid velocity u = 0 is examined. Consider a 2D scalar field ψ(x, y) that is non-zero on a square region and zero elsewhere as shown in Fig. 4. Despite the lack of any fluid motion, the mesh is made to change over the course of the simulation by decreasing the minimum element edge-length for mesh adaptivity with time (see Table 1 for its functional form). The mesh therefore gradually refines along the square interfaces of the problem and the scalar field must be interpolated at every time-step. With no change in minimum element edge-length there would of course be essentially no interpolation to do in this example. Note also that this simulation involves no advective transport and hence constitutes a pure test of the interpolation method itself and in particular of the conservativeness of the method. The results after adapting the mesh 100 times using node-wise interpolation and the Grandy method are compared with those obtained by the control-volume Galerkin method presented in this paper. The numerical parameters used in the simulation are specified in Table 1. The resolution specifies the approximate number of elements in the initial mesh as well as in the final mesh. Note that here and in all examples, the final resolutions differ for the three interpolation methods and a single representative value is quoted. The minimum/maximum element edge-lengths bound the element sizes that enter the metric tensor in adaptivity. They control the smallest and largest element sizes permitted in a mesh after an adapt and hence effectively control the mesh resolution for a fixed interpolation error bound and gradation settings (the scaling ratio between adjacent elements). A concise overview of pertinent features of mesh adaptivity can be found in [19]. The minimum/maximum edge-lengths are fixed to be the same for the three interpolation methods. As discussed previously, the minimum element edge-length decreases as a function of time to a minimum value of 0.001.
Using the notation introduced at the start of this section, a P1DG pressure discretisation was used in this example. A diffusive correction as outlined in Section 2.2 was also applied after interpolation to keep the solution fields within physical bounds. The same correction is applied to both node-wise and control-volume Galerkin methods, but is not required in the Grandy case because this is automatically bounded.
The initial state is shown in Fig. 4 and the results of the three interpolation methods after 100 mesh adapts are shown in Fig. 5. It is clear that node-wise interpolation no longer sharply captures the interfaces in the problem and inhomogeneities have begun to appear along the boundary of the square. This is to be contrasted with Grandy and Galerkin interpolation which are able to capture the boundaries very closely. Qualitatively, the Grandy and Galerkin results are similar in this example. Fig. 6(a) shows the integral of the scalar field over the domain for each interpolation method, and Fig. 6(b) shows the error in the three interpolation methods. The integral fluctuates in the node-wise case every time the mesh adapts and conservativeness is lost. The fluctuations in the control-volume Galerkin and Grandy cases on the other hand are negligible -of the order machine precision, implying that the methods are conservative. In all cases, the L 2 error decreases with time as the mesh resolution increases (due to the decreasing minimum element size) as expected but the error is smallest in the control-volume Galerkin case. The error is reduced compared to Grandy interpolation due to a reduction in inhomogeneities  node-wise interpolation (circles) and Grandy interpolation (triangles). (b) L 2 error in the static scalar field plotted against the number of mesh adapts for control-volume Galerkin projection (squares), node-wise interpolation (circles) and Grandy interpolation (triangles).

Table 2
Parameters used in the Gaussian scalar simulation discussed in the main text.

Size
Target along the interface boundary. Note as mentioned in Section 1, node-wise interpolation can be problematic for discontinuous interpolants. Whilst not so pronounced in this example, it will become so in what follows.

Higher-order convergence
As discussed in Section 2.1, the control-volume Galerkin method begins by projecting a field defined control-volume wise onto a finite element representation. In contrast to Grandy interpolation, this yields (in general) a higher-order polynomial representation as opposed to a flat function. It is expected that this higher-order representation in finite element space should lead to higher-order numerical convergence and this is indeed the case as is now demonstrated by example. We also further motivate this higher-order convergence in Appendix A.
In order to study the convergence of the interpolation methods, we consider a sequence of simulations beginning on a high resolution fixed mesh and mapping once (through mesh adaptivity) to a new mesh with different minimum edgelength (and hence resolution). The initial fixed mesh resolution is kept the same whilst the minimum element edge-length on the target mesh is varied over the range specified in Table 2. The three interpolation methods are used in turn to map the scalar fields from the old mesh onto the new mesh and the L 2 error in the interpolation results is plotted as a function of minimum element edge-length on the target mesh. The same P1DG pressure discretisation as in the previous example is adopted.
In this example, the scalar field ψ(x, y) that is interpolated from the old mesh to the new mesh is taken to be a 2D Gaussian function of the form. The initial scalar distribution together with an example target mesh are shown in Fig. 7 and the numerical parameters used in the simulation are shown in Table 2. The target resolution specifies the approximate number of elements in the target mesh and the minimum/maximum element edge-lengths bound the element sizes that enter the metric tensor in adaptivity. The initial mesh that the scalar field is initially defined on before being mapped onto the target mesh has 20000 elements (Fig. 7b). This mesh is by construction very fine so as to represent the Gaussian initial state with a very small error. A sequence of 11 simulations are performed (for each interpolation method in turn) varying the minimum element edge-length between 0.06 and 0.25.
A log-log plot of discrete L 2 error against minimum element edge-length for the three interpolation methods is shown in Fig. 8. It is clear that node-wise interpolation and Grandy interpolation exhibit slower than linear convergence, whilst the control-volume Galerkin method shows close to second order convergence and is hence a higher-order method. The order m of the trend line through the control-volume Galerkin data is m ∼ 1.85, whilst m ∼ 0.77 for node-wise interpolation and m ∼ 0.96 for Grandy interpolation.

Advection of a square scalar field
The next step in complexity is to introduce a non-zero bulk velocity into the interpolation problem. We consider again initial data where the scalar is zero everywhere except on a (rotated) square. The system is initialised with a positive velocity in the x direction and the mesh is adapted to the curvature of the scalar. Since there is now advective transport of the scalar, the mesh changes with time automatically and there is no need to vary the minimum element edge-length to engineer this. Once again the results of interpolating via node-wise and Grandy interpolation are compared with those obtained by the control-volume Galerkin method presented in this paper.
The numerical parameters used in the simulation are specified in Table 3. The CFL number reported here corresponds to the largest value reached after approximately 140 adapts. The velocity is that of the initial square, and the size refers to the whole domain dimensions. Note the total domain length is purposely chosen to be much larger than the total distance travelled by the square during the simulation to minimise the influence of errors due to material advection through the outlet boundary on conservation integrals. Thus, the behaviour of the interpolation methods can be studied in a way that is decoupled as much as possible from the properties of the advection method. The resolution specifies the approximate number of elements in the initial mesh (this remains approximately constant throughout the simulation in this case for all interpolation methods except in the node-wise case, where the solution field is gradually destroyed) and the minimum/maximum edgelengths bound the element sizes that are permitted during adaptivity. Note that now the minimum element edge-length is fixed unlike in the previous example so there is no 'global' mesh refinement. The same P1DG pressure element as in the previous example is adopted together with the compressive advection scheme outlined at the start of Section 3. Table 3 Parameters used in the square advection simulation described in the main text. The initial state is shown in Fig. 9 and the results of the three interpolation methods after approximately 140 mesh adapts are shown in Fig. 10. Fig. 11(a) shows the integral of the scalar field over the domain for all three methods and Fig. 11(b) shows the L 2 error in the solutions. With advection included, there are now major qualitative differences in the results of the three interpolation methods at late times. In the case of node-wise interpolation, the solution field has been severely distorted. In the Grandy case, the initial state is well preserved in the interior of the square but the interfaces have become spread out. This behaviour can be attributed to the diffusivity associated with Grandy interpolation and can be problematic in certain classes of problems as will be illustrated in the final example discussed in this paper.
Of all the methods, control-volume Galerkin interpolation most closely retains the structure of the initial square, keeping both the uniform field values internally and the sharp interface boundaries. There is some overall distortion to the square-tip geometry common to all three methods but this can be attributed to numerical errors in advection and not the interpolation. Fig. 11(a) verifies conservation in the Galerkin and Grandy cases and lack thereof in the case of node-wise interpolation. This example demonstrates that conservativeness is maintained even in scenarios with advection. Finally, Fig. 11(b) demonstrates that as expected, the error in the case of node-wise interpolation increases rapidly compared to the other two methods. The error is again smallest in the case of the control-volume Galerkin method, remaining approximately constant over the course of the simulation.
It is also instructive to compare the computational costs of the three interpolation methods in this example. The inclusion of advection, makes the problem non-trivial and the costs are representative of the general behaviour observed in other examples. The control-volume Galerkin method runs approximately 12% slower than the Grandy method, but correspondingly leads to an approximate 45% reduction in the L 2 error at the end of the simulation. Node-wise interpolation is in general very similar in cost to Grandy interpolation (only 1%-5% faster than the Grandy method) except in situations such as this example where the interpolant is discontinuous and the solution field is damaged by the interpolation. In these cases, the node-wise method is in fact much more costly (running approximately 54% slower than the Grandy method) as solving the equation system becomes progressively harder with time, ultimately leading to a very large error in the final state. These simulation times were obtained running in serial with a 2.6 GHz dual Intel Xeon processor and 128 GB of RAM. In summary, the control-volume Galerkin interpolation method is able to significantly reduce the interpolation error during mesh adaptivity compared to node-wise and Grandy interpolation whilst incurring little extra computational cost.

Advection of a circular scalar field
We now consider interpolation for a final pure advection problem, this time with a continuous interpolant. In this case, the 2D scalar field ψ(x, y) is initialised to unity on a disc which is then set rotating in the anticlockwise direction for a quarter of a turn after which the sign of the velocity is changed and the disc is rotated back to its initial state. Since the flow is reversible, the final scalar field distribution should be identical to the initial state. We compare how closely node-wise    Fig. 12. Initial circular scalar field distribution (a) and initial mesh (b). The disk has radius r = π 5 and is centered at π 2 , π +1 5 . and control-volume Galerkin interpolation are able to capture this reversibility by analysing the scalar fields in the final state.
The initial state of the system is shown in Fig. 12 and the numerical settings used in the simulation are shown in Table 4. The Courant number reported here corresponds to the largest value reached. The initial resolution refers to the approximate number of elements in the initial/final mesh whilst the inter(mediate) resolution specifies the number of elements after a quarter turn before the velocity is reversed. The minimum/maximum edgelength control the element sizes in adaptivity. A P2CG pressure element was used in the simulation.
We show the results after a quarter turn in Fig. 13 and after returning to the initial state in Fig. 14. After this time 400 mesh adapts have taken place. In this example, there is again a qualitative difference between the results of the two interpolation methods. In particular, it is apparent from Fig. 14 that the control-volume Galerkin results more closely capture the original scalar field distribution. In particular, the node-wise interpolation results have developed oscillatory modes along the boundary of the disc that are not present in the initial state.

Viscous instabilities in two phase immiscible displacement in a porous medium
As a final and more complex application of the new interpolation method presented in this paper, we consider simulation of an immiscible two-phase displacement in a porous medium. This process is of great importance occurring in the water-flooding of oil reservoirs [20]. In such scenarios where a more viscous fluid is displaced by a less viscous fluid, it can be shown that the interface between the two phases is unstable to infinitesimal transverse perturbations. This instability may be viewed as analogous to the Saffman-Taylor instability observed in Hele-Shaw cell experiments but now in the context of two-phase porous media flow [21]. These perturbations can grow, becoming macroscopic viscous fingers. Simulation of viscous fingering is known to be challenging requiring a combination of higher-order numerical methods and very fine resolution to minimise the influence of numerical diffusion [22][23][24]. Mesh adaptivity is advantageous in such situations, allowing very fine mesh resolution to be maintained in the vicinity of the fingers as they grow with coarser resolution elsewhere [25].
When modelling unstable fluid phenomena with mesh adaptivity it is important that the interpolation method adopted is conservative since any fluctuations in field values introduced by interpolation can lead to artificial, unphysical perturbations that may grow and affect the solution field. In this example, mesh adaptivity is used to simulate viscous fingering at early-times, using a discontinuous representation of water saturation, and contrasting the results of the new interpolation method presented in the paper with those of Grandy interpolation. A discontinuous representation is used to maximise the effective number of control-volumes to capture the instability. Note we do not use node-wise interpolation since, as shown by the previous examples, it is unreliable in applications to discontinuous fields. This example also serves as a demonstration of when the Grandy method, which has hitherto performed well, can yield poor results due to its highly diffusive nature.
The equations governing two phase flow in a porous medium are the generalised Darcy law for the volumetric fluid flux q α : together with the continuity equation: and finally the constraint that the volume fractions occupied by each of the two fluids sum to unity [26]: In the above, K is the permeability tensor, α indexes the two fluids, k r,α is the relative permeability of each fluid, p α the pressure, S α the phase volume fraction (saturation) and φ the porosity of the porous medium. Note that capillary pressure is neglected in the discussion here.
In the case of 1D symmetric solutions of the form S(x, y) = S(x), Eq. (23) and Eq. (24) reduce to the Buckley-Leverett equation [20]: where M is the viscosity ratio of the two fluids and q w is the fluid inflow velocity. Eq. (26) admits shock solutions composed of an advancing compression wave together with a rarefaction. In order to simulate viscous fingering we are interested in the behaviour of pressure and saturation perturbations about solutions to the Buckley-Leverett equation (26). Solution of the linearised equations can be shown to reduce to an eigenvalue problem that relates the growth rate of perturbations to their frequency. Further theoretical discussions of viscous fingering can be found in for example [22,23,[26][27][28][29].
We consider simulation of an immiscible displacement using mesh adaptivity to accurately resolve the shock front and the growth of viscous fingers. The simulation was carried out using a higher-order CVFE simulator developed at Imperial College detailed in [7,10]. The simulator solves a modified form of the equation system above (Eq. (23)-Eq. (25)) where a new velocity (u α = q α /S α ), distinct from the Darcy velocity, is solved for. Velocity is eliminated from the equation system and the resulting pressure-saturation equation is then solved. The simulations were run with a P2DG-P1DG element (second order discontinuous finite element representation for velocity, first order discontinuous finite element representation for pressure). The numerical solution is interpolated between meshes after each adapt. Galerkin interpolation (c, d). The initial state is overlaid on these images in white to compare how closely it is retained by the interpolation methods.

Table 5
Multiphase flow parameters used in the immiscible displacement simulation described in the main text. The last two columns specify the irreducible saturations S wr , S or of the two phases used in the relative permeability functions.

Viscosity ratio
Permeability Porosity Irred. sat. phase 1 Irred. sat. phase 2 μ = 303 As can be seen from Eq. (26), closure of the system of equations requires specification of the functional form of the relative permeabilities. A Corey-type correlation [30] given by [24] is used: (27) where S wr and S or are the immobile fractions of the two fluids (the fractions of the two fluids that cannot be displaced). The remaining multiphase and numerical parameters used in the simulation are specified in Table 5 and Table 6 respectively. The displacing phase (water) is injected from the left with the specified velocity. The resolution specifies the approximate number of elements in the adaptive mesh. The maximum value 10500 is the number of elements after the first adapt, whilst the value 7000 corresponds to the approximate number of elements at the end of the simulation. The number of elements in the initial mesh is approximately 20000 whilst the whole domain fixed mesh we compare the adaptive results against has approximately 40000 elements. Viscous fingering can be triggered numerically by introducing a perturbation to the saturation, pressure or permeability fields. It is triggered here by a saturation perturbation across the inlet boundary at the first time-step (shown in Fig. 15(c)). The form of the perturbation used is the same as that considered in [24]. This perturbation is allowed to grow on a fixed mesh until the fingers reach just under half-way across the domain at which point the mesh adapts to the water saturation field. The mesh is initially kept fixed so as to allow the fingering pattern time to form. Adapting too early before the fingering pattern has fully formed could damage the solution field. The saturation profile and mesh immediately before the first adapt are shown in Fig. 15 whilst the results after the fingers have nearly reached the outlet (after 680 mesh adapts) are shown in Fig. 16 for both methods of interpolation together with results on a fixed mesh (with edge-length equal to the minimum allowed in the adaptive simulations) at the same point in time.
It is clear from these results that the control-volume Galerkin method is able to preserve the structure of the viscous fingers after many adapts whilst in contrast the diffusivity of the Grandy method smears the fingering pattern, capturing only the large-scale structure of the front and diffusing away any small scale details. This behaviour can be attributed to the lower order nature of Grandy interpolation -smaller variations in field values are eroded during the interpolation process which leads to a runaway adaptive coarsening of the mesh.
Further evidence to support the conclusion that the control-volume Galerkin method is able to capture viscous fingering with comparable accuracy to results obtained on a fixed mesh can be seen in saturation cross sections through the fingering pattern for the three simulations. From Fig. 17(a) it is clear that mesh adaptivity with control-volume Galerkin interpolation produces very similar average saturations to those found in the fixed mesh simulations. In contrast, Grandy interpolation shows deviations from the fixed mesh in the region where the viscous fingers are forming. A water saturation cross section as a function of vertical spatial coordinate along a transverse slice (x = 0.04) that intersects the fingering pattern is shown in Fig. 17(b). The control-volume Galerkin results more closely match the pattern of maxima and minima in saturation seen on the fixed mesh when compared to the Grandy method which tends to average out nearby extrema. Similarly, the saturation cross section along a longitudinal slice ( y = 0.022) is shown in Fig. 17(c). The control-volume Galerkin results again more closely match the fixed mesh results. The growth rate in particular appears to be under-predicted in the Grandy case with the fingers being too short. Finally, a plot of the water flux across the outlet boundary as a function of dimensionless simulation time is shown in Fig. 17(d). The water flux is initially zero until the viscous fingers reach the outlet boundary and cause water 'breakthrough'. The control-volume Galerkin interpolation results show a very similar breakthrough time to that seen on the fixed mesh whilst the Grandy results show a later breakthrough time due to the diffused fingering pattern.

Conclusions
We have presented a general higher order, conservative and bounded method for interpolating fields between controlvolume meshes. This is achieved by first mapping the control-volume field into a finite element representation on the donor mesh by Galerkin projection, then interpolating onto the target mesh by constructing a finite element supermesh from the intersection of the donor and target meshes and finally projecting back to a control-volume representation on the target mesh. Boundedness is achieved by using a simple diffusion algorithm on the resulting interpolated solution on the target mesh.
A series of test cases have demonstrated that the method is superior to both node-wise and Grandy interpolation for control-volume fields with a discontinuous dual finite element representation. Node-wise interpolation was shown to result in significant distortion in the solution field when using a discontinuous (higher order) interpolant. Grandy interpolation did not result in such distortion but was found to be much more diffusive than the Galerkin projection method. For the test case investigating the simulation of immiscible viscous fingering in porous media this meant that the Grandy interpolation smeared out the structure of the viscous fingers after a large number of mesh adapts. Overall the Galerkin projection method has a lower error than both the node-wise and Grandy interpolation methods, when compared with a fine fixed mesh, and causes less numerical diffusion.
The error in the Galerkin projection method grows much more slowly than in both the node-wise and Grandy interpolation so would be the method of choice for simulations requiring frequent mesh adapts. It is slightly more CPU intensive that the Grandy method (12% slower in our test case) but the L 2 error is reduced at the end of the simulation (45% in this test case). The node-wise interpolation was even slower than the Galerkin projection method for discontinuous interpolants because it distorts the solution field making the equation system harder to solve as the simulation progresses. For lower order interpolants the performance of the node-wise interpolation was similar to that of the Grandy interpolation.
The proposed control-volume Galerkin projection method is currently the only conservative, higher order interpolation technique available for control-volume fields. The method was applied to a CVFE model in this paper but it should be possible to generalise this to other interpolations of control-volume fields.

Appendix A. Higher-order convergence of the CV interpolation method
In this appendix, we provide arguments to support the claim that the control-volume interpolation introduced in this paper exhibits higher-order convergence. As discussed in Section 2.1, the method is composed of three steps: A map from a control-volume to finite element representation on the old mesh, a map between finite element representations on the old and new meshes and finally a mapping from a finite element representation back into control-volumes on the new mesh. The overall convergence of the method can be studied by bounding the interpolation errors in each of these three projections.
A.1. Error analysis for projections (1) and (3) We begin by considering the error in the projection between control-volume and finite element meshes. In the case of a finite element representation, the errors in Eq. (4) and Eq. (7) are identically zero. This follows from the definition of finite element interpolation from or onto a control-volume mesh and can be seen clearly in Fig. 2 by regarding the error as the difference between the finite element interpolation and the initial control-volume functions.

A.2. Error analysis for projection (2)
We now consider the error in Eq. (6), the FEM Galerkin projection from old to new mesh. Let an element E k be a tetrahedron over a 3D domain , I new and I old are the interpolation operators over the new and old meshes respectively.
where E k ∈ (E 1 , . . . , E K ) (K is the total number of elements over the domain ), d = 3 is the dimensional size of the domain, h i (E k ) are suitably defined element sizes on the element E k [31].
The error bound in Eq. (A.1) can be further estimated [32] as: where A E k is the edge shape matrix, σ min is the minimal singular value of A E k , κ is the order of approximation of the finite element degrees of freedom and V E k is the volume of the element E k . Since we have argued that the convergence behaviour of the interpolation method is governed by the second step, we can conclude from the above that the method is higher-order. The error bound Eq. (A.2) is consistent with the higher-order numerical example in the main text. In that case, a first order discretisation (κ = 1) for pressure was used, and thus quadratic h 2 (E k ) convergence is expected, as seen in Fig. 8.