A Generalised Complete Flux scheme for anisotropic advection-diffusion equations

by   Hanz Martin Cheng, et al.
TU Eindhoven

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.


page 3

page 18

page 19

page 20

page 21


A second-order accurate semi-Lagrangian method for convection-diffusion equations with interfacial jumps

In this paper, we present a second-order accurate finite-difference meth...

New Numerical Interface Scheme for the Kurganov-Tadmor second-order Method

In this paper, we develop a numerical scheme to handle interfaces across...

Combining the hybrid mimetic mixed method with the Scharfetter-Gummel scheme for magnetised transport in plasmas

In this paper, we propose a numerical scheme for fluid models of magneti...

A fully local hybridised second-order accurate scheme for advection-diffusion equations

In this paper, we present a fully local second-order upwind scheme, appl...

Fourth-Order Anisotropic Diffusion for Inpainting and Image Compression

Edge-enhancing diffusion (EED) can reconstruct a close approximation of ...

A Petrov-Galerkin method for nonlocal convection-dominated diffusion problems

We present a Petrov-Gelerkin (PG) method for a class of nonlocal convect...

1. Introduction

Let be an open connected subset of with boundary , where and are disjoint. We consider the advection-diffusion problem: Find such that


Here, we assume that , the source term , the velocity profile , and the diffusion tensor is a symmetric positive definite field in .

Remark 1.1 (boundary condition).

If , then (1) has a unique solution . For the case , is only unique up to an additive constant. Hence, for uniqueness of , we impose

where is a specified average value of in .

Stationary advection-diffusion problems (1) are a fundamental component of non-stationary advection-diffusion equations of the form

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., . 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 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 .

  • 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 , 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.

Figure 1. Stencils for the discretisation: cells involved in the stencil are shaded gray. Left(cell unknowns only); right(including edge unknowns).

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.

2. 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 the set of cells and faces (edges in 2D) of our mesh, respectively. The set of faces is then written as a disjoint union of interior faces and exterior (boundary) faces , , where . For a cell , we denote by the set of faces (edges) of the cell . 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 , and use Stokes’ formula to obtain the balance of fluxes:

where is the unit outer normal to the face . Now, if is a face shared by two distinct cells and , then

This is known as the conservation of fluxes.

We now denote the approximations to the diffusive and advective fluxes by and , 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 ,

and for each face shared by distinct cells and ,


Here, is a piecewise constant approximation of the source term in cell . 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 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 be a point in cell such that is star-shaped with respect to (typically, a good choice is taking to be the cell barycenter), and for all , define to be the orthogonal distance between and (see Figure 2). We then set to be the centre of mass of . In dimension two, this is simply the midpoint of the edge .

Figure 2. Notations in a generic cell in dimension .

Following [5], if and is the convex hull of and (see Figure 2), we set


Here, is a linearly exact reconstruction of the gradient, i.e., if interpolates an affine function at the edge midpoints, then . The second term in is a stabilisation term, which ensures the coercivity of the numerical scheme. Owing to the fact that for any constant , we may also write


For the HMM method, the definition of the discrete diffusive flux is based on the bilinear form , which stems from the weak formulation of the advection-diffusion equation (1). To be specific, the definition of the fluxes is inspired by the relation


where is the trace operator, with

Under the assumption that is constant in cell with value , we can write

For , the fluxes are then defined by a discrete counterpart of (5),


We note that the fluxes in (6) are uniquely defined. In particular, we can see from (3) and (4) that is uniquely determined by the value of . Hence, for a given edge , can be uniquely determined from (6) by setting, for example, and for all the other edges .

2.2. SG method for advection

We now discuss in this section the discretisation of the advective fluxes . 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 ,


Here, is a characteristic distance, given by

where denotes the smallest distance between and . The quantity is defined depending on whether is an interior edge shared by cells and , or is a boundary edge. Here, we set


are the eigenvalues of

. The purpose of scaling 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 and . The first sum is obtained by taking such that and 0 elsewhere in (6). The latter sum is obtained by taking the sum over in (7). Combining these, we obtain the discrete flux balance equation (2a)

Now, if is an interior edge, we take such that , and 0 for all other components in (6). This gives us the conservativity of the diffusive fluxes

The conservativity of the advective fluxes

is then imposed separately with 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 , if , that is, is a Neumann boundary edge, then we still take such that , and 0 for all other components in (6) to get . 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.

3. 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 , which results in an extension of the stencil onto neighboring elements. In particular, the discrete balance of fluxes is now given by


where are the inhomogeneous fluxes, which come from the cell and its neighboring elements. The inhomogeneous fluxes are defined so that for each cell , the inhomogeneous fluxes are nonzero only on interior edges . 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: Find such that


We assume that the source term , and for simplicity of exposition, assume that the diffusion and velocity field are constants. Proper boundary conditions (Dirichlet or Neumann) are then imposed at and . We then form a partition of the domain . For , the idea behind the complete flux scheme is to solve a local boundary value problem

where . Now, we define

where is a point between and where the flux has to be evaluated. We also introduce a scaled coordinate so that


It can then be shown that an analytical expression for is given by


Integrating (10) from to , we obtain

Substituting the explicit expression of the flux obtained in (12) to the above relation then gives

where . Multiplying by , we obtain the relation


Integrating (13) from to and dividing both sides of the expression by yields


where .

Remark 3.2 (Homogeneous flux).

We note that taking , the one-dimensional homogeneous flux 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 , , and


and the local Péclet number 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 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 , 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 in (11), which leads to the function in (16). In the literature, ; however, in two dimensions or higher, it might occur, especially when the grid is distorted, that . In this case, the definition of and 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 . We use the notation , and for simplicity of exposition, consider Cartesian meshes. We start by considering an edge being shared by cells and . This is described by , (Figure 3, left). We then find two points and in and , 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 ). We then construct rectangular regions and associated to cells and , respectively (Figure 3, right).

Figure 3. Regions involved for the inhomogeneous fluxes.

To obtain an approximation of the complete flux along , we treat the advection-diffusion equation (1) as a quasi-one-dimensional boundary value problem by writing



are the standard basis vectors in

. We note here that for the particular edge under consideration, . Now, to find the value of , 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 . Along the edge shared by cells and , we present two ways of defining a local Péclet number . We denote by the average value of in cell and by the average value of 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 , and the direction of the velocity field is taken into account by taking the dot product with the unit outer normal . 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 which is very large, regardless of the direction of . To see this, we write an orthogonal diagonalisation , where is the diagonal matrix containing the eigenvalues of and

    is the orthogonal matrix containing the eigenvectors of

    . From this, we have

    Note that the eigenvectors of form an orthonormal basis for . As a consequence, we can write


    where are the column vectors of with

    Here, and are the coordinates of and with respect to the basis . We also note that we can determine the values of and in (20). In particular, for ,


    It then follows that and are vectors with entries and , respectively. Denoting by the eigenvalue of that corresponds to the eigenvector , we then see that


    Hence, the term containing the smallest eigenvalue of will dominate (provided that is nonzero), especially if the condition number of is large. Moreover, we see in (20) and (22) that the strength of advection over diffusion is computed with respect to the basis . That is, measures the strength of advection over diffusion in the direction . Afterwards, the Péclet number along is computed by a weighted average, where the weights are determined by writing as a linear combination of . We note that due to computing the strength of advection over diffusion in the direction , a wrong sign may be obtained for the Péclet number along .

    As an example, given a moderate advection, say , and a constant diffusion tensor


    then we would expect the Péclet number along the - and -directions to be moderate. However, the condition number of is very large. This is due to the fact that the eigenvalues of , and , have very large contrast. Also, we observe that the corresponding eigenvectors are given by

    Since , we see that the dominant term in the sum (22) will come from . It can be computed that . The sign of the Péclet number is then determined by the sign of . For Cartesian meshes, the outward normal vectors along the east and north edges are given by and , respectively. This gives and . It is notable here that due to computing the relative strength of advection over diffusion along and , the Péclet number along the -direction has an incorrect sign. Moreover, (19) yields a Péclet number with a very large magnitude in both the - and -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 . 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 is to directly compute a Péclet number that is oriented towards . In particular, we now see that the numerator is the strength of advection along , and the denominator is the strength of diffusion along . Hence, the Péclet number along is computed directly without having to go through the basis vectors . Moreover, upon writing an orthogonal diagonalisation , and denoting by the vector with th entry equal to (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 . Physically, this means that if material is being transported outside (into) cell , which corresponds to being positive (negative), then the Péclet number preserves this property. Moreover, since is an outward unit normal vector, we have

    This shows that the diffusion along is a weighted average of the eigenvalues of . This is a good weighted average in the sense that if is almost orthogonal to , then , which means that the corresponding eigenvalue would only have a minimal contribution to the Péclet number.

    Considering again 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 - and - directions. Moreover, compared to (19), the choice (24) captures properly the strength of advection. In particular, for the choice , we see that the Péclet number along the -direction is twice as large as that along the -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 at 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 along is given by


where is a generalisation of (11), satisfying


and are the average values of the right hand side of (18) along and , respectively, i.e.

Substituting the above expressions for into (26) and taking the integral over , we obtain an approximation to the inhomogeneous component of the flux . In particular, the inhomogeneous flux is given by: