A Generalised Complete Flux scheme for anisotropic advection-diffusion equations

In this paper, we consider separating the discretisation of the diffusive and advective fluxes in the complete flux scheme. This allows the combination of several discretisation methods for the homogeneous flux with the complete flux (CF) method. In particular, we explore the combination of the hybrid mimetic mixed (HMM) method and the CF method, in order to utilize the advantages of each of these methods. The usage of HMM allows us to handle anisotropic diffusion tensors on generic polygonal (polytopal) grids; whereas the CF method provides a framework for the construction of a uniformly second order method, even when the problem is advection dominated.


Introduction
Let Ω be an open connected subset of R d (d = 2, 3) with boundary ∂Ω = Γ D ∪ Γ N , where Γ D and Γ N are disjoint. We consider the advection-diffusion problem: Find c ∈ H 1 (Ω) such that ∇ • (cV − Λ∇c) = s on Ω, c = g on Γ D , Here, we assume that g, h ∈ L 2 (∂Ω), the source term s ∈ L 2 (Ω), the velocity profile V ∈ L 2 (Ω) d , and the diffusion tensor Λ is a symmetric positive definite field in L 2 (Ω) d×d . which are usually encountered in computational fluid dynamics. To be specific, these types of problems are encountered in applications to plasma physics [12] and porous media flow (e.g. reservoir engineering [7,14] and groundwater flow [6]). In particular, for flows in porous media, heterogeneous and highly anisotropic diffusion tensors are involved, due to varying rock properties (such as permeability and porosity) in the domain [9].
The uniformly second order complete flux (CF) scheme [19] was originally developed on Cartesian grids for (1) with scalar diffusion, i.e., Λ = D(x)I d . Some other uniformly second (or higher order) schemes for the advection-diffusion problem (1) with scalar diffusion, applicable on generic meshes, can be found in [2,13]. On the other hand, some recent approaches to deal with anisotropic diffusion tensors while maintaining high order accuracy involve the use of discontinuous skeletal or hybrid high order methods [15,16], multipoint flux approximations (MPFA) [20], or the introduction of nonlinear fluxes [11]. The aim of this paper is to extend the complete flux schemes in [19] in order to handle anisotropic diffusion, whilst maintaining its uniformly second order convergence properties.
In particular, one of the main advantages in using the complete flux scheme lies in the fact that it is linear; hence, it is much cheaper to implement than nonlinear methods. Also, the generalised complete flux scheme we propose in this work only requires one unknown inside each cell, and another one on each interface. In comparison, hybrid high order methods require piecewise linear approximations to the solution c in order to achieve second order accuracy. In two dimensions, this corresponds to four unknowns inside a cell, and two unknowns on each interface. Hence, in general, the generalised complete flux scheme is computationally cheaper. On the other hand, the stencils for hybrid high order methods are more compact, i.e., they do not depend on values from neighboring cells, whereas the stencil for the generalised complete flux scheme depends on values from the neighboring cells.
The novelties of this work are the following: • Splitting the diffusive and advective components of the flux, which allows for a combination of different numerical discretisations with the CF method.
In particular, we describe the combination of the hybrid mimetic mixed (HMM) and the CF method. • A generalisation of the CF method on nonuniform meshes in one dimension, which paves a way to formulating the CF method on generic meshes in dimension d > 1. • The formulation of a generalised local Péclet number, which allows us to use the CF method for highly anisotropic problems. • An alternative derivation of the CF method in two-dimensional Cartesian meshes, which can be straightforwardly extended into three dimensions. • The introduction of unknowns along the faces (edges), which allows us to have a more compact stencil for the CF method. In particular, the stencil will only involve the direct neighbors of a cell K, i.e., in the case of Figure  1, introduction of edge unknowns will yield a stencil that involves only five cells, whereas without the edge unknowns, the stencil involves nine cells.
The outline of the paper is as follows: In Section 2, we present a finite volume (homogeneous flux) discretisation for (1), and give details on how we discretise the diffusive and advective fluxes. Following this, we describe in Section 3 the introduction of inhomogeneous fluxes, which are fundamental in the construction of the complete flux schemes. Numerical tests are then provided in Section 4 to illustrate the second order accuracy of the combined HMM-CF scheme. Finally, Section 5 provides a summary, as well as possible directions for future research.

Finite volume methods for the advection-diffusion problem
For the discretisation, we start by defining a mesh as in [15, Section 3.1], i.e., a partition of Ω into polyhedral (polygonal) cells. We then denote by T = (M, E) the set of cells K and faces (edges in 2D) σ of our mesh, respectively. The set of faces is then written as a disjoint union of interior faces E int and exterior (boundary) faces E ext , E = E int ∪ E ext , where E int ∩ ∂Ω = ∅. For a cell K ∈ M, we denote by E K ⊂ E the set of faces (edges) of the cell K. Our discretisation is piecewise constant, and hence, involves one unknown on each cell and another one on each of the faces (edges), and we denote the space of unknowns by To write a finite volume scheme, we start by taking the integral of (1) over a control volume K, and use Stokes' formula to obtain the balance of fluxes: where n K,σ is the unit outer normal to the face σ ∈ E K . Now, if σ is a face shared by two distinct cells K and L, then This is known as the conservation of fluxes. We now denote the approximations to the diffusive and advective fluxes by F D K,σ and F A K,σ , respectively, i.e.
Key to the formulation of a finite volume scheme is the definition of discrete fluxes (diffusive and advective) so that the fluxes satisfy a discrete version of the balance and conservation of fluxes [4]. That is, for all K ∈ M, and for each face σ shared by distinct cells K and L, Here, s K is a piecewise constant approximation of the source term s in cell K.
For this paper, we will use the HMM method [5] for the diffusive fluxes, and the Scharfetter-Gummel (SG) method [17] for the advective fluxes.
2.1. HMM method for diffusion. In this section, we discuss the discretisation of the diffusive fluxes F D K,σ by using the HMM method. The choice of using the HMM method is due to the fact that it can handle anisotropic diffusion tensors. To construct the fluxes, we start by introducing a piecewise constant gradient, which is defined on a sub-triangulation of cells. Let x K ∈ K be a point in cell K such that K is star-shaped with respect to x K (typically, a good choice is taking x K to be the cell barycenter), and for all σ ∈ E K , define d K,σ > 0 to be the orthogonal distance between x K and σ (see Figure 2). We then set x σ to be the centre of mass of σ. In dimension two, this is simply the midpoint of the edge σ.  x σ x K Following [5], if K ∈ M and (D K,σ ) σ∈E K is the convex hull of σ and x K (see Figure 2), we set ∀w ∈ X D , ∀x ∈ D K,σ , Here, ∇ K w is a linearly exact reconstruction of the gradient, i.e., if (w σ ) σ∈E K interpolates an affine function A at the edge midpoints, then ∇ K w = ∇A. The second term in ∇ D w(x) is a stabilisation term, which ensures the coercivity of the numerical scheme. Owing to the fact that σ∈E K |σ|w K n K,σ = 0 for any constant w K , we may also write For the HMM method, the definition of the discrete diffusive flux is based on the bilinear form a(u, w) := Ω Λ∇u • ∇w, which stems from the weak formulation of the advection-diffusion equation (1). To be specific, the definition of the fluxes (F D K,σ ) K∈M, σ∈E K is inspired by the relation where γ : Under the assumption that w is constant in cell K with value w K , we can write For c ∈ X D , the fluxes (F D K,σ ) K∈M, σ∈E K are then defined by a discrete counterpart of (5), We note that the fluxes F D K,σ in (6) are uniquely defined. In particular, we can see from (3) and (4) that ∇ D v is uniquely determined by the value of (v σ −v K ). Hence, for a given edge σ ∈ E K , F K,σ can be uniquely determined from (6) by setting, for example, v σ − v K = 1 and v σ − v K = 0 for all the other edges σ ∈ E K . 2.2. SG method for advection. We now discuss in this section the discretisation of the advective fluxes F A K,σ . The original SG method gives a discretisation for both the diffusive and advective fluxes simultaneously. However, this original formulation of the SG method does not work in cases where diffusion is anisotropic. In this section, we use a modification of the SG method, which only gives the discretisation of the advective fluxes, introduced in [3]. This is done by setting Following the ideas in [1,3], we use a hybridised discretisation for the advective flux, using edge unknowns instead of unknowns from neighboring cells. The hybridised SG method for advective fluxes then reads: For each K ∈ M, σ ∈ E K , Here, h σ is a characteristic distance, given by where d(x K , σ) denotes the smallest distance between x K and σ. The quantity λ σ is defined depending on whether σ is an interior edge shared by cells K and L, or σ is a boundary edge. Here, we set where spec(Λ K ) are the eigenvalues of Λ K . The purpose of scaling h σ V K,σ by the minimum eigenvalue of Λ is to bring enough numerical diffusion to ensure a better stability for advection dominated problems [3]. In some instances, however, this definition of λ σ might introduce excessive numerical diffusion. This will be illustrated in Section 3.2; an improvement in the definition of λ σ , which captures directional diffusion better, will then be introduced in equation (25) of Section 3.2.
2.3. Finite volume (homogeneous) fluxes. In this section, we define the finite volume (homogeneous) fluxes, given by For the discrete flux balance equation (2a), we need to compute both σ∈E K F D K,σ and σ∈E K F A K,σ . The first sum is obtained by taking v ∈ X D such that v K = 1 and 0 elsewhere in (6). The latter sum is obtained by taking the sum over σ ∈ E K in (7). Combining these, we obtain the discrete flux balance equation (2a) Now, if σ ∈ E K , K ∈ M is an interior edge, we take v ∈ X D such that v σ = 1, and 0 for all other components in (6). This gives us the conservativity of the diffusive fluxes F D K,σ + F D L,σ = 0. The conservativity of the advective fluxes F A K,σ + F A L,σ = 0, is then imposed separately with F A K,σ defined as in (7). Together, these give us the conservativity of fluxes (2b), which is equivalent to Finally, we describe the fluxes along the boundary edges. Given a cell K ∈ M, if σ ∈ E K ∩ Γ N , that is, σ is a Neumann boundary edge, then we still take v ∈ X D such that v σ = 1, and 0 for all other components in (6) to get F D K,σ . We then impose the Neumann boundary conditions by setting On the other hand, if σ is a Dirichlet boundary edge, we impose At this stage, we note that other methods may be used for discretising the homogeneous diffusive and advective fluxes. To maintain the second order accuracy of the scheme, a natural choice for the discretisation is such that the approximations to both the diffusive and advective fluxes are second order.

Inhomogeneous fluxes and the complete flux scheme
In this section, we discuss the complete flux (CF) scheme. The main idea behind the CF scheme is the introduction of inhomogeneous fluxes F I K,σ , which results in an extension of the stencil onto neighboring elements. In particular, the discrete balance of fluxes is now given by where F I K,σ are the inhomogeneous fluxes, which come from the cell K and its neighboring elements. The inhomogeneous fluxes are defined so that for each cell K ∈ M, the inhomogeneous fluxes are nonzero only on interior edges σ ∈ E K ∩ E int . Inhomogeneous fluxes are not needed for imposing boundary conditions. Remark 3.1 (Localised stencil). With the hybridised discretisation of the advective fluxes in (7), the expression (9) gives a more localised stencil compared to the original formulation of the complete flux scheme [19] (see Figure 1, left and right).

3.1.
Complete flux scheme in one dimension. To better understand the formulation of the complete flux scheme, we start by recalling its derivation in one dimension, following the ideas in [19]. In one dimension, (1) reads: We assume that the source term s ∈ L 2 (Ω), and for simplicity of exposition, assume that the diffusion and velocity field Λ, V ∈ R are constants. Proper boundary conditions (Dirichlet or Neumann) are then imposed at x 0 and x N . We then form a partition of the domain x 0 < x 1 < · · · < x N −1 < x N . For x ∈ (x j , x j+1 ) j=0,1,...,N −1 , the idea behind the complete flux scheme is to solve a local boundary value problem where x σj is a point between x j and x j+1 where the flux has to be evaluated. We also introduce a scaled coordinate α σj ∈ (0, 1) so that It can then be shown that an analytical expression for f is given by Substituting the explicit expression of the flux obtained in (12) to the above relation then gives where f σj = f (x σj ). Multiplying by (Λ −1 e −P ), we obtain the relation Integrating (13) from x j to x j+1 and dividing both sides of the expression by xj+1 xj whereP j =P (x j ).
Remark 3.2 (Homogeneous flux). We note that taking x σj = x j+ 1 2 , the onedimensional homogeneous flux f H σj in (14) corresponds to the original (second order accurate) Scharfetter-Gummel flux [17], where diffusion and advection are discretised simultaneously. The approximation of the homogeneous flux in (14) can be replaced by other discretisations in diffusion and advection. A natural choice for maintaining the accuracy of the numerical scheme is to use second order schemes for discretising both the diffusive and advective fluxes.
We now focus on the inhomogeneous flux. Since diffusion and velocity are assumed to be constant, and we use a piecewise constant approximation for the source term, an explicit expression for the inhomogeneous flux is given by where s j = s(x j ), ∆x j = x j+1 − x j , and and the local Péclet number P σj is defined to be Remark 3.3 (Péclet number). We note that in one dimension, the definition of the (localised) Péclet number (17) is quite straightforward, due to the fact that • No complication arises in taking the quotient in (17), since both V and Λ are scalars. • The direction of the velocity is only either towards the left or right, and hence the sign of the Péclet number is either negative or positive, respectively. This was automatically taken care of in (15). In order to extend to dimension d > 1, a generalised Péclet number would have to be introduced, to ensure that the direction of the velocity and the relative local strength of advection over diffusion are captured properly.
The main difference in this formulation of the one-dimensional complete flux scheme on nonuniform grids compared to [8,18] is the introduction of α σj in (11), which leads to the function Z(P, α) in (16). In the literature, α σj = 0.5; however, in two dimensions or higher, it might occur, especially when the grid is distorted, that α σj = 0.5. In this case, the definition of α σj and Z(P, α) as in (11) and (16) would be required.
3.2. Complete flux scheme in higher dimensions. In this section, we discuss the formulation of the complete flux scheme in higher dimensions, starting with dimension d = 2. We use the notation x = (x, y), and for simplicity of exposition, consider Cartesian meshes. We start by considering an edge σ ∈ E K being shared by cells K and L. This is described by x = x σ , y S < y < y N (Figure 3, left). We then find two points x K and x L in K and L, respectively, and construct a segment orthogonal to σ that passes through these points (in Figure 3, left, this pertains to the segment that lies on the line y = y σ ). We then construct rectangular regions  x To obtain an approximation of the complete flux F K,σ ≈ σ (−Λ∇c + cV) • n K,σ along σ, we treat the advection-diffusion equation (1) as a quasi-one-dimensional boundary value problem by writing where e x = (1, 0), e y = (0, 1) are the standard basis vectors in R 2 . We note here that for the particular edge under consideration, n K,σ = e x . Now, to find the value of (−Λ∇c + cV) • e x , we need to solve the quasi-one-dimensional problem (18). In order to do so, we need to define an extension of the local Péclet number (17), so that it is applicable in dimension d > 1. Along the edge σ shared by cells K and L, we present two ways of defining a local Péclet number P K,σ . We denote by Λ K the average value of Λ in cell K and by V σ the average value of V on σ.
• Firstly, we may define the local Péclet number as This is a straightforward modification of (17), where the quotient is obtained by taking the matrix inverse of Λ K , and the direction of the velocity field is taken into account by taking the dot product with the unit outer normal n K,σ . However, this is highly dependent on the eigenvalues of Λ. If the eigenvalues of Λ have high contrast, then in general, this would always yield a Péclet number P K,σ which is very large, regardless of the direction of n K,σ . To see this, we write an orthogonal diagonalisation where D K is the diagonal matrix containing the eigenvalues of Λ K and U K is the orthogonal matrix containing the eigenvectors of Λ K . From this, we have . Note that the eigenvectors of Λ K form an orthonormal basis for R d . As a consequence, we can write where u K,i are the column vectors of U K with u T K,i u K,j = δ ij . Here, γ i and β i are the coordinates of V σ and n K,σ with respect to the basis {u K,i } i=1,...,d . We also note that we can determine the values of γ i and β i in (20). In particular, for i = 1, . . . , d, It then follows that U T K V σ and U T K n K,σ are d×1 vectors with entries γ i and β i , respectively. Denoting by λ K,i the eigenvalue of Λ K that corresponds to the eigenvector u K,i , we then see that Hence, the term containing the smallest eigenvalue λ K,j of Λ K will dominate (provided that γ j β j is nonzero), especially if the condition number κ(Λ K ) of Λ K is large. Moreover, we see in (20) and (22) that the strength of advection over diffusion ( γi λ K,i ) is computed with respect to the basis {u K,i } i=1,...,d . That is, ( γi λ K,i ) measures the strength of advection over diffusion in the direction u K,i . Afterwards, the Péclet number along n K,σ is computed by a weighted average, where the weights are determined by writing n K,σ as a linear combination of {u K,i } i=1,...,d . We note that due to computing the strength of advection over diffusion in the direction u K,i , a wrong sign may be obtained for the Péclet number along n K,σ .
As an example, given a moderate advection, say V = (1, 2) T , and a constant diffusion tensor then we would expect the Péclet number along the x-and y-directions to be moderate. However, the condition number κ(Λ) = 10 8 of Λ is very large. This is due to the fact that the eigenvalues of Λ, λ 1 = 1 and λ 2 = 10 −8 , have very large contrast. Also, we observe that the corresponding eigenvectors are given by Since 1 λ2 = 10 8 , we see that the dominant term in the sum (22) will come from u 2 . It can be computed that V • u 2 < 0. The sign of the Péclet number is then determined by the sign of n K,σ • u 2 . For Cartesian meshes, the outward normal vectors along the east and north edges are given by e x and e y , respectively. This gives (Λ −1 It is notable here that due to computing the relative strength of advection over diffusion along u 1 and u 2 , the Péclet number along the xdirection has an incorrect sign. Moreover, (19) yields a Péclet number with a very large magnitude in both the x-and y-directions.
At this stage, we also recall that the scaling factor λ σ in (7) was for the purpose of stabilising the SG method for advection dominated regimes. However, in this case, where diffusion and advection are both moderate, λ σ would always be 10 −8 . This results in introducing too much numerical diffusion to the scheme. We will show in the numerical tests in Section 4 that these cause the accuracy of the scheme to degrade to first order, which is not desirable.
• Another option would involve taking The rationale behind the usage of n T K,σ Λ K n K,σ is to directly compute a Péclet number that is oriented towards n K,σ . In particular, we now see that the numerator V σ • n K,σ is the strength of advection along n K,σ , and the denominator n T K,σ Λ K n K,σ is the strength of diffusion along n K,σ . Hence, the Péclet number along n K,σ is computed directly without having to go through the basis vectors {u K,i } i=1,...,d . Moreover, upon writing an orthogonal diagonalisation Λ K = U K D K U T K , and denoting by B K the d × 1 vector with ith entry equal to β i (see (21)), we obtain From this, we see that the quantity on the denominator is always positive; hence it does not affect the sign of V σ • n K,σ . Physically, this means that if material is being transported outside (into) cell K, which corresponds to V σ • n K,σ being positive (negative), then the Péclet number preserves this property. Moreover, since n K,σ is an outward unit normal vector, we have This shows that the diffusion along n K,σ is a weighted average of the eigenvalues of Λ K . This is a good weighted average in the sense that if u K,i is almost orthogonal to n K,σ , then β i = n T K,σ u K,i ≈ 0, which means that the corresponding eigenvalue λ K,i would only have a minimal contribution to the Péclet number.
Considering again V = (1, 2) T and Λ as in (23), we have Here, we now see that the Péclet numbers both have correct signs and are moderate along both the x-and y-directions. Moreover, compared to (19), the choice (24) captures properly the strength of advection. In particular, for the choice V = (1, 2) T , we see that the Péclet number along the y-direction is twice as large as that along the x-direction. Drawing inspiration from the definition (24) of the local Péclet number, we redefine the scaling factor λ σ in the modified Scharfetter-Gummel flux (7) to be These choices mitigate the introduction of excessive numerical diffusion; hence allowing us to preserve the second order accuracy expected from the complete flux scheme, as will be demonstrated in the numerical tests in Section 4. Upon solving the quasi-one-dimensional problem (18), the value of (−Λ∇c + cV) • e x at x σ will of course consist of both a homogeneous and an inhomogeneous component. Since we already have a discretised homogeneous flux in Section 2, we only focus on the inhomogeneous component. By a generalisation of (15), the inhomogeneous component of (−Λ∇c + cV) • e x along x = x σ is given by where α K,σ is a generalisation of (11), satisfying s K andŝ L are the average values of the right hand side of (18) along x K < x < x σ and x σ < x < x L , respectively, i.e.
Substituting the above expressions forŝ K ,ŝ L into (26) and taking the integral over σ, we obtain an approximation to the inhomogeneous component F I K,σ of the flux σ (−Λ∇c + cV) • n K,σ . In particular, the inhomogeneous flux F I K,σ is given by: wheres Remark 3.4 (Extension into 3D). The formulation of the inhomogeneous flux F I K,σ by taking the integral of (26) over σ can straightforwardly be applied to obtain inhomogeneous fluxes in 3D. The only modification to (28) would be the definition ofs K ands L , due to an additional term that would come from the partial derivative with respect to z in (18).
We now detail how to computes K , which consists of two terms. The first term involves the integral of the source term, which, under the assumption that s is piecewise constant with value s K on cell K, can be written as On the other hand, we recognise that the second term can be written as a sum of fluxes, namely where F K σ ,σ N and F K σ ,σ S are the fluxes along the northern and southern edges, respectively, of the rectangular region K σ (see Figure 4, left). Since the outward x normal vectors at σ N and σ S are both orthogonal to the normal vector at σ, we call these the cross-fluxes of σ. Assuming that the fluxes along the northern and southern edges of K are constant and approximated by the homogeneous fluxes F H K,σ N and F H K,σ S respectively, we have A similar approximation holds for F K σ ,σ S . The terms ins L are then computed using a similar argument. This completes the definition of the inhomogeneous flux F I K,σ along the interior edge σ ∈ E K . A similar process is used to obtain the inhomogeneous fluxes along the other edges of K.
To summarise, the inhomogeneous flux F I K,σ along the edge σ of cell K is given by On square or rectangular meshes, we can write the cross fluxes in terms of the homogeneous fluxes as in (29) to obtain Lemma 3.5 (conservativity of the inhomogeneous fluxes). The definition of the inhomogeneous fluxes (28) is conservative. That is, if σ is an edge shared by cells K and L, then We start by taking note that P K,σ = −P L,σ and α L,σ = 1 − α K,σ . Substituting these expressions into (28) then gives for σ ∈ E K ∩ E int do 3: Find the corresponding neighbor L that shares σ with K.

4:
Find points x K , x L in K and L so that the segment passing through x K , x L is orthogonal to σ.

8:
Obtain an approximation for the cross-fluxes of σ in K σ and L σ .

9:
Obtain an approximation for the value of the source term s at K σ and L σ .

10:
Substitute the obtained quantities into (30) to obtain F I K,σ .

11:
end for 12: end for Algorithm 1 presents a short summary of how to compute the inhomogeneous fluxes. Although the expression for the inhomogeneous fluxes F I K,σ were derived on two-dimensional Cartesian meshes, Algorithm 1 is also applicable for threedimensional Cartesian meshes (cf. Remark 3.4). Moreover, it also provides a framework for computing inhomogeneous fluxes on irregular meshes. As an illustration, Figure 4, right, shows the cells K σ and L σ and the cross-fluxes of σ on a triangular mesh. The only caveat here is that on non-rectangular meshes, the cross-fluxes no longer lie on cell faces. Hence, we need to perform interpolations and a change of coordinates (similar to those described in [18]) in order to obtain approximations for the cross-fluxes and source terms in lines 8 and 9 of Algorithm 1. This is not straightforward, and will be discussed in a separate paper.

Numerical tests
In this section, we present some numerical tests to demonstrate the second order accuracy of the generalised complete flux scheme for the advection-diffusion equation (1). These will be performed on Cartesian meshes with square cells over the domain Ω = (0, 1) × (0, 1). In particular, for the tests presented below, the homogeneous fluxes will be computed via HMM (6) for diffusion, and SG (7) with λ σ in (25) for advection.

Convergence tests.
We start with three test cases, where the velocity field V and the diffusion tensor Λ are constant over the domain Ω. For these tests, we solve (1) with a given exact solution c(x, y) = sin(πx) sin(πy). Denoting by π D c the piecewise constant function reconstructed from the discrete unknowns, we measure the relative errors in the L 1 -norm Here, we fix V = (1, 2) T , and take different Λ, corresponding to strong or mild anisotropy, and also to whether the problem is advection dominated or not. Here, the strength of anisotropy is quantified by the condition number κ(Λ) of Λ; we say that anisotropy is strong if κ(Λ) ≥ 10 4 .
• Test case 1: constant, homogeneous diffusion, advection dominated. We start with an advection dominated test case, where the diffusion is a scalar, i.e., Λ = 10 −8 I d . • Test case 2: strong anisotropy (not aligned with the mesh), moderate advection. Here, advection is said to be moderate in the sense that the Péclet number defined in (24) P K,σ < 1 for all K ∈ M, σ ∈ E K . For this test case, we consider the diffusion tensor The purpose of this choice for the diffusion tensor is twofold: firstly, it shows that the generalised Péclet number defined in (24) is better than the one in (19). Secondly, this shows that the combination of the HMM and CF yields a second order scheme for problems with strong anisotropy (even if the anisotropy is not aligned with the mesh). In this case, the condition number of Λ is κ(Λ) = 10 8 . • Test case 3: strong anisotropy (almost aligned with the mesh), advectiondominated. For our third test, we consider a diffusion tensor which features strong anisotropy that is almost aligned with the mesh, where diffusion is very weak in the y-direction, but moderate in the x-direction. This is done by taking Λ = 1.5 10 −4 10 −4 10 −8 .
Here, the condition number is κ(Λ) ≈ 4.5 × 10 8 , and the eigenvalues are given by λ 1 ≈ 3.33 × 10 −9 , λ 2 ≈ 1.5, with corresponding eigenvectors This aims to show that even for a strongly anisotropic diffusion tensor and a relatively strong advection, by making a proper choice for the Péclet number, the combined HMM-CF method still has second order accuracy. For the first test case, we note that the Péclet numbers (19) and (24) are identical. In particular, since the diffusion is a scalar, the generalised complete flux scheme returns to the formulation presented in [19], which we expect to achieve second order accuracy. This is confirmed by the numerical tests presented in Table 1. Now, we see in Tables 2 and 3 that the Péclet number (24) maintains the second order accuracy of the complete flux scheme. On the other hand, using the Péclet number in (19) reduces the complete flux scheme into a first order scheme. This is due to the fact that (19) causes the numerical scheme to decide whether the problem is advection dominated based on the eigenvalues of Λ (i.e., the Péclet number is always large when the eigenvalues of Λ have high contrast), which is due to computing the strength of advection over diffusion along the directions of the eigenvectors of Λ, as discussed in Section 3.2. This phenomenon is also illustrated in the third test case, for which diffusion is only weak in the y-direction. By using (19), the Péclet number has a magnitude of 3.99 × 10 4 in the x-direction and a magnitude of 5.99 × 10 8 in the y-direction; hence, diffusion is interpreted to be weak in all directions. On the other hand, by using (24), the Péclet number has a magnitude of 6.66 × 10 −1 and 2 × 10 8 in the x-and y-directions, respectively, so diffusion is weak only in the y-direction.
Having illustrated the second order convergence of the scheme, we now proceed to test the limits of the scheme by performing two extreme tests. In these cases, an analytical solution is not available, and hence, we analyse the qualitative aspects of the numerical solutions.
The velocity field considered is V = (40x(2y − 1)(x − 1), −40y(2x − 1)(y − 1)) T , which simulates a counterclockwise rotation (see Figure 5). The source term is a ring positioned at a distance of 0.35 from the center of the domain, i.e., s(x, y) = 10 −2 exp(−(r − 0.35) 2 /0.005), where r 2 = (x − 0.5) 2 + (y − 0.5) 2 (see Figure 6). We now observe that the source is evenly distributed, and the velocity field is oriented along the direction of increasing diffusivity. Consider now the interface shared by Ω 4 and Ω 1 . The velocity field causes the solution profile to carry the source in Ω 4 towards the interface it shares with Ω 1 . Due to the low diffusivity along the y-direction in Ω 4 , we expect that the value of the solution to be small near the interface. On the other hand, due to the strong diffusivity along the y-direction in Ω 1 , we expect the solution to be large near the interface it shares with Ω 4 . Hence, we expect the exact solution to form internal layers near the interfaces that separate the subdomains. This can be observed in the solution profiles presented in Figure 7. A better visualisation is given by 3D plots in Figure  9. Moreover, we observe that the numerical solution does not yield non-physical negative values (minimum value is non-negative). We also note that in the eyeball norm, the numerical solution at the coarse mesh with 60 × 60 cells is already very close to what we obtain on the very fine mesh of 480 × 480 cells, which is expected from the complete flux scheme. Also, the maximum value of 7.3 × 10 −4 observed here is also very close to the maximum value ranging from 6.7 × 10 −4 − 6.9 × 10 −4 observed in [3,6].
Another point of comparison can be made by reversing the flow,that is, by taking −V, whilst retaining all the other parameters used for the test case. In this case, the velocity field is now oriented along the direction of decreasing diffusivity; hence, we expect the solution to be continuous near the interfaces. This behavior can be seen in Figure 8 and Figure 9, right. As with the initial test, we observe that the numerical solution at the coarse mesh with 60 × 60 cells is already very close to  what we obtain on the very fine mesh of 480 × 480 cells. Also, the maximum value of 7.9 × 10 −4 observed here is also very close to the maximum value of 7.3 × 10 −4 observed in [6]. We note however that negative values are present in the numerical solutions, but the effect is not very significant, since they are only very small in magnitude: 10 −6 on the coarse mesh, and 10 −10 on the fine mesh.

Monotonicity test.
Finally, we present a test from [11], which we refer to as test case 5. This is an extreme test in the sense that many linear methods result in a significant violation of the discrete maximum principle, and consequently produce numerical solutions that exhibit non-physical oscillations. We consider Ω = (0, 1) 2 \ [4/9, 5/9] 2 ; that is, the domain is a square with a hole punched in the middle, resulting to a boundary consisting of two disjoint parts Γ 1 and Γ 2 (see Figure 10). Here, the source term is set to be s = 0, and Dirichlet boundary conditions are imposed, with c = 0 on Γ 1 and c = 2 on Γ 2 .   The velocity field is given by V = (700, 700) T , and the diffusion tensor reads Λ = R(−π/6) 1000 0 0 1 R(π/6), R(θ) = cos θ sin θ − sin θ cos θ .
For this test case, the exact solution is expected to be bounded by 0 and 2. Here, we see in Figure 11 that the maximum value is close to what is expected, with slight overshoots occurring at the top right corner of the boundary Γ 2 . On the other hand, we see that the complete flux scheme does not guarantee the non-negativity of the numerical solution, with undershoots occurring near the lower left corner of Γ 2 . Non-physical oscillations are also detected in these regions. Although we observe these non-physical qualities, it still offers an improvement over some classical methods, e.g. the lowest-order RaviartThomas mixed finite element method, for which the numerical solution is negative over almost half of the computational domain [10]. In particular, the overshoots and undershoots for the generalised complete flux scheme are only observed locally near the top right and lower left regions of Γ 2 and only cover 16% and 6% of the domain for the coarse and fine meshes, respectively.

Summary and outlook
In this paper, we were able to present a generalised complete flux scheme for anisotropic advection-diffusion problems. The main novelty comes from splitting the diffusive and advective components of the flux, which allows for a combination of different numerical discretisations with the CF method. This resulted in a scheme which can be applied to problems with anisotropic diffusion tensors. Related to this, we introduced a generalised local Péclet number (24), which allows us to capture properly the local strength of advection over diffusion. Another important contribution is an alternative presentation of the CF method in two dimensions, which can straightforwardly be extended into three dimensions. This also gave us a framework of how to extend the complete flux scheme onto more generic meshes; which would be a topic for future research. The numerical tests presented in Section 4.1 illustrate the second order accuracy of the generalised CF method, even for strong anisotropy and advection dominated problems. Moreover, the tests performed in Section 4.2 showcase the ability of the generalised CF method to handle strongly anisotropic heterogeneous diffusion tensors. We however notice in Section 4.3 the presence of non-physical oscillations in the numerical solution. This was not unexpected as it has already been remarked in [11] that these non-physical oscillations are present in many linear methods, and future work will involve working on how to mitigate these non-physical oscillations. Another interesting aspect for future work will involve the study of non-stationary advection-diffusion problems.