1. Models and assumptions
1.1. Introduction
In this paper, we introduce an algorithm for implementing characteristicbased schemes for a pure advection model
(1) 
Here, , is an open bounded domain of (), the porosity , the source term , and the velocity are given, with on . The unknown represents the amount of material (a fraction) present at . Note that the boundary is noncharacteristic due to the assumption on , and thus no boundary conditions need to be enforced in (1).
The need to solve advection equations of the form (1) usually forms part of an operator splitting technique used to solve an advectiondiffusion equation
(2) 
where the diffusion tensor
is given. Advectiondiffusion equations of this type are usually encountered in mathematical models for porous media flow (e.g. reservoir simulation, nuclear waste storage) [15, 20], and computational fluid dynamics (e.g. NavierStokes equations) [18], and are usually advection dominated. The diffusive component of the model is discretised separately, by mixed finite elements (MFEM), finite volumes, or other schemes that fit in the framework of the gradient discretisation method (GDM) [13], and will not be detailed in this paper. Here, we only focus on the implementation of characteristicbased schemes in (1), such as the Eulerian Lagrangian Localised Adjoint Method (ELLAM) and the Modified Method of Characteristics (MMOC). The advantage of these schemes lies on the fact that they are based on characteristic methods, and thus capture the advective component of the PDE better than standard upwind schemes.Several variants of the ELLAM, some of which are the finite element (FE) ELLAM [4] and the finite volume (FV) ELLAM [16], as well as a summary of their properties, have been presented in [21]. One of the major issues faced when implementing characteristicbased schemes is the conservation of mass (both local and global). In order to achieve global mass conservation, some adjustments were performed on the MMOC, leading to the development of MMOC with adjusted advection (MMOCAA) [12]. Although the MMOCAA achieves global mass conservation, it does not achieve local mass conservation. On the other hand, from its formulation, ELLAM satisfies global mass conservation; more recent variants of the ELLAM, such as the volume corrected characteristics mixed method (VCCMM) [1, 2, 3], achieve local volume conservation by adjusting the points tracked through the characteristics. These points may also be adjusted by following the algorithm proposed in [10]. Another way to achieve local volume conservation for characteristicbased schemes has been proposed in [8]. This is particularly useful for schemes with piecewise constant approximations, such as hybrid and mixed finite volume type schemes [14]. As an example, in [8], it was used to perform adjustments to make the HMM–ELLAM schemes in [7] locally mass conserving. More details about the convergence analysis and implementation of GDM characteristicbased schemes for (2) and its applications to flows in porous media, can be found in [8, 9]. For schemes with piecewise constant approximations, evaluating the integrals arising from the discretisation of (1) boils down to computing intersections between polytopal regions (polygons in 2D and polyhedra in 3D). Although several algorithms are available for taking the intersection of polygons in 2D, they are quite expensive to implement in practice. Moreover, even though these methods can theoretically be extended to 3D, the main difficulty for a 3D implementation would come from taking intersections between polyhedra. Most of the polyhedral intersection algorithms in 3D are able compute the intersection between two convex polyhedra efficiently, as in [5, 6, 17, 19]. However, even though the cells are initially convex, the tracked cell may not be convex. To our knowledge, the intersection of a convex polyhedron with a general polyhedron has only been dealt with in [11], and even here, the computation of the intersection is not trivial or easy to implement.
The purpose of this paper is to develop a feasible method to implement characteristicbased schemes in 2D and 3D, whilst preserving the important properties of local and global mass conservation. The novelty of this paper is the idea of approximating the polytopal regions by balls (circles in 2D, spheres in 3D). By doing so, we convert the problem of computing polytopal intersections into that of computing intersections of balls, which is trivial to implement and has an essentially zero computational cost. Naively doing so will lead to a loss of mass, and hence we propose an adjustment algorithm which will help reduce the errors induced by this loss of mass. We then design to solve an optimisation problem, with both global and local mass conservation as constraints. We call this process the Ballapproximated characteristics (Bchar) method. Due to its formulation, the Bchar method will yield a scheme that is both locally and globally mass conserving.
The paper is organised as follows. We start by giving some details on the assumptions on the data for the advection equation (1). After which, we give a short summary of the ELLAM scheme used to discretise this equation in Section 2. We also enumerate some of its mass conservation properties, and give a physical interpretation of the scheme. We then give a brief summary of how the ELLAM type schemes were implemented in the literature. The Bchar method will then be introduced in Section 3. In Section 4, numerical tests will first be performed in 2D in order to compare the performance of the Bchar method with the ELLAM scheme obtained from polygonal intersections, with volume adjustments as described in [8]. In these tests, we will see that the Bchar method yields very similar results to the polygonal intersections, with a much cheaper computational cost. Numerical tests are performed to show the applicability of the Bchar ELLAM in 3D.
1.2. Assumptions on the data, and numerical setting
We start by forming a mesh, i.e. a partition of into polygonal (in 2D) or polyhedral (in 3D) sets. Following the notations in [13, Definition 7.2], we then denote to be the set of cells and faces (edges in 2D) of our mesh, respectively. We also use to denote the volume (area in 2D) of a cell . Throughout the article we assume the following properties:
(3)  
Assumption (3) simply states that the initial concentration inside the medium and the porosity of the medium are bounded, which is natural in physical applications. Also, the assumption that is piecewise smooth on is not restrictive, since in practice, is usually taken to be piecewise constant on each cell . As in [8], we describe the numerical method in a general setting, to ensure that our algorithm applies at once to various possible spatial discretisations for the diffusion terms in (2). These can be dealt with using the GDM as shown in [8, 9], and will not be discussed in further detail for this paper. We replace, in the weak formulation of the model, the continuous (infinitedimensional) spaces and corresponding operators by a discrete (finitedimensional) space and function reconstructions. We then define a spacetime discretisation , where

is a finitedimensional real space, describing the unknowns of the chosen scheme,

is a linear operator that reconstructs a piecewise constant function on the mesh from the unknowns,

are the time steps, and we let .
Different choices of lead to different schemes (e.g. finite volume based methods, including hybrid ones with face unknowns like HMM [14], or masslumped finite element methods [22]).
Finally, we assume that
(4) 
and that is approximated on each time interval by a function
(5) 
Remark 1.1 (Approximation of the velocity field).
Although is given in (1), we use an approximation for the velocity field in order to include the more general case where comes from solving a PDE coupled to (2). For example, for flows in porous media, usually comes from Darcy’s law: Given a source term should satisfy the PDE on , with suitable boundary conditions.
In the rest of the paper, the variables are only made explicit in the integrals when there is a risk of confusion. Otherwise we simply write, e.g., .
2. ELLAM scheme for the advection–reaction equation
For simplicity, we assume that there is no source term, i.e. in (1). We then start by multiplying (1) with a sufficiently smooth function , and performing integrations by parts. Using the identity
(1) gives, for any time interval ,
To simplify the second term on the left hand side of the above equation, the ELLAM requires that test functions satisfy
(6) 
with given. The advection equation (1) then leads to the relation
(7) 
We now write the ELLAM scheme, which consists of writing (7) in the discrete context, in which trial and test functions are replaced by reconstructions
applied to trial and test vectors in
.Definition 2.1 (ELLAM scheme).
Given a spacetime discretisation , the ELLAM scheme for (1) reads as: find such that and, for all , satisfies
(8) 
where is the solution to
(9) 
Define the flow such that, for a.e. ,
(10) 
Under Assumptions (3) and (5), the existence of this flow is proved in [9, Lemma 5.1]. The solution to (9) is then understood in the sense: for and a.e. , . In particular,
(11) 
Remark 2.2 (Source term).
For a nonzero source term , we simply replace the right hand side of (8) by an approximation for the spacetime integration of the function , e.g. a trapezoid rule
The construction of the Bchar method in Section 3 will draw inspiration from a physical interpretation of the ELLAM, which uses the fact that is a piecewiseconstant reconstruction on a given mesh . For each cell , we assume that there is such that , where is the function that has a value of 1 in , and 0 elsewhere. Writing and taking as test function, (8) and (11) give
which reduces to
(12) 
where is the available porous volume in a set . The term on the right hand side of (12) tells us that the amount of material present in a particular cell at time is obtained by intersecting a tracked cell and a residing cell . This intersection can be interpreted as locating where the material in cell comes from, hence backtracking the cell to , measuring which fraction of the material is taken from each (by taking their intersection), and deposing this fraction into the cell .
2.1. Global mass conservation
Since the advection equation (1) usually comes from solving a model in computational fluid dynamics or engineering, we would want our numerical scheme to conserve global mass. Essentially, we would want an equation which tells us that the change in is dictated by the amount of inflow/outflow and by the source term. In this case, due to the noflow boundary conditions and the absence of a source term, this simply means that the amount of substance present at time should be the same as the amount of substance present at time . The desired equation is thus given by
(13) 
It can easily be checked that the ELLAM scheme satisfies this property. Indeed, taking the sum over all in (12) yields
which is the discrete form of (13).
Remark 2.3 (Achieving global mass conservation).
We note here that the ELLAM scheme achieves global mass conservation due to
(14) 
for all . An analogue of this identity will be needed to ensure that the Bchar method in Section 3 also achieves global mass conservation.
2.2. Local mass conservation
One of the main difficulties of implementing an ELLAM type scheme is the evaluation of the integral for each cell , where . In general, the region (see Figure 1, left) cannot be exactly described and hence, in the literature, it was approximated by polygons obtained from backtracking the vertices, together with a number of points along the edges of the cell . Figure 1 (right) gives an illustration of the approximate traceback region obtained by tracking the vertices, together with the edge midpoints of the cell .
Remark 2.4 (Reconstruction of polytopes).
In 2D, most of the time, we can reconstruct the polygons approximating the traceback region by following the tracked points in the same order as the original points, since it gives a welldefined polygon. However, in 3D, a face that is tracked may no longer be planar, and the original polyhedron faces need to actually be triangulated to ensure that a polyhedron is created after tracking.
In general, . However, the equality of these volumes is essential, otherwise the numerical scheme will not be able to preserve even a constant solution. Consider, for example, the simple case of a divergence free velocity field in (1), with and . In this test case, the exact solution is given by . In theory, upon implementing an ELLAM scheme with piecewise constant approximations for the unknown , we should have the following simplified form of (12) at the first time step:
However, due to the approximation of the traceback region, we only have
and thus
This example shows that an inaccurate approximation of the volume of the tracked cell renders the numerical scheme unable to recover constant solutions. Hence, we need to perform some adjustments on the polygonal region in order to yield , which we shall define as the local volume constraint for . Several adjustment strategies which would lead to local mass conservation have been studied, as in [1, 8, 10]. In particular, for local mass conservation to be achieved, we should have, for all ,
(15) 
3. Bchar method
In this section, we present the idea of approximating the cells by balls, instead of the usual approximation using polygons. We will call this type of approximation the Ballapproximated Characteristics (BChar). For simplicity of exposition, we consider solenoidal fields, so that and . For each cell , we choose points in its interior. We then assume that each of these points represents centers of disjoint balls which are strictly inside cell . The idea now is to distribute the porous volume in each cell over the balls ; hence, we compute an adjusted porosity so that . The quantities may now be interpreted as the porous volume inside the ball . The main interest of approximating the cells by balls is the fact that computing the intersection of balls is trivial compared to intersecting polytopes. As a consequence, the computational cost is greatly reduced. Moreover, this idea is easily applicable in both 2D and 3D.
Upon working on the assumption that each ball, when tracked, remains as balls, the points are then tracked by solving (10) to obtain , which will be treated as the center of the tracked ball (see Figure 2). Of course, this assumption is not true in general, but gives a good approximation of the volumes, especially if the initial balls are small. Since we work with solenoidal fields, the radii and measure of each of the tracked balls are unchanged, i.e.
As a consequence, the adjusted porosity remains unchanged, , and the porous volume inside each of the tracked cells is given by
(16) 
3.1. Initial approximation for the volume of intersecting regions
We now describe the process for obtaining an initial approximation for the volume of the intersecting regions .
We start by recalling that is interpreted as the amount of material in a cell that comes from a residing cell . In the context of approximation by balls, this reads: each ball contains an amount of material from some residing balls transported by the flow, given by .
Remark 3.1.
The choice of (instead of ) in the approximation comes from the interpretation that the amount of material present in is obtained by measuring how much of the material is taken from each , and deposing this material into the ball .
An initial approach for approximating would then involve taking the sum of the masses of the balls in a residing cell , intersected with the tracked balls that originated from cell , that is . However, since there are gaps between the residing balls, . This will lead to a loss in volume, which will in turn lead to a loss of mass conservation and a poor approximation. Instead, we use this to compute
(17) 
which represents the fraction of the mass in that comes from . From this, we then see that
is the actual amount of mass in the tracked ball that comes from . The quantity is then computed by taking the sum over all tracked balls and residing balls , given by
(18) 
3.2. Mass conservation for the Bchar method
Since are approximations to , in order to achieve local mass conservation, we should have an analogue of (15), given by . We can easily check that in (18) satisfies this relation by using (16). Hence, the approximation of by in (18) leads to a scheme that is locally mass conserving. However, , which means that (14) is not satisfied, and thus global mass conservation is not achieved. This indicates that some adjustments need to be performed on . Upon indexing each tracked cell and each residing cell , , let be the matrix with entries , where . In short, each entry of the matrix gives an approximation of the volume . In terms of the matrix , we would thus need
(19) 
in order to achieve (15), which will lead to local mass conservation. Now, in order to have a globally mass conserving scheme, we would need to satisfy (14), which in this context is equivalent to
(20) 
To build a matrix which satisfies (19) and (20), we start with the assumption that for . This means that at least one of the balls of a residing cell intersects with at least one of the balls that has been tracked from . The only time holds is when the flow traces everything into voids. For example, in Figure 3, none of the balls tracked from , nor did any of the other balls that were tracked from time , intersect with the residing balls in , and so for the corresponding to the particular cell . This can easily be resolved by increasing the number of balls, and making sure that each residing cell is tightly packed with balls. Another instance when is when there is a very strong inflow or outflow at the boundary of the domain, which is not the case for (1).
Remark 3.2 (Presence of inflow or outflow).
In the presence of inflows or outflows, we may avoid by one of the following options:

Take a smaller time step so that the region does not get emptied out.

Create ghost cells at the boundary of the domain so that the cells that get emptied out are the ghost cells.
Now, start by setting
which yields
Next, we set
(21) 
so that
Essentially, the adjustments perform the following twostep process:

Firstly, we redistribute, according to a proportion evaluated by , the mass along each of the intersecting regions so that global mass conservation (20) is achieved.

Next, we redistribute the mass so that local mass conservation (19) is achieved.
Intuitively, we can see that these adjustments involve redistributing the errors and hence scaling them down in each iterate. A naive approach would involve iterating this process, stopping only when the error in global and local mass conservations are less than a certain tolerance value. However, there is no guarantee that such a result is achievable. A more efficient approach would involve, after taking iterations and arriving at the matrix , solving a minimisation problem. In practice, we found that taking to be such that the maximum error in mass conservation is at most is sufficient to give a good initial approximation.
We then assign an unknown corresponding each of the entries of , which gives us unknowns. For such that , the corresponding unknown is fixed to 0. This tells us that if no intersection has been detected between a tracked cell and a residing cell , then our adjustment algorithm should not introduce any volume into these regions. Hence, the number of unknowns in our new system is equal to the number of nonzero entries in . From the nonzero entries of construct a matrix in the following manner: Write
and denote by () the row vector of size obtained by removing the zero entries in . The first rows of are then formed by the matrix
That is, we stagger the vectors so that the coefficients of start at the column after the last coefficient of
, and we pad each row with zeros to ensure we obtain an
matrix.The latter rows of the matrix are then formed in the following manner: for the th row, we look for the nonzero entries () in column of . For each corresponding , we then find the column corresponding to where the coefficient resides in the first rows of . We then set . As an example, if
then would be (with the line separating the first rows from the last rows of )
In practice, the last rows of are assembled simultaneously with its first rows. In the example above, after setting as the first, second and third entries of the first row, are simultaneously set to be the first, second and third entries of the and th row, respectively. After which, when and were set to be the fourth and fifth entry of the second row, they were also set to be the fourth and fifth entries of the and the th row. In this manner, the whole matrix is rather easy to assemble.
The essential property of is that when it is multiplied by an vector consisting of all ones, we have
The sums in the righthand side correspond to the quantities that must be fixed to certain values in order to achieve local and global mass balance, see (19) and (20). To obtain local and global mass conservation, we therefore solve the system
(22) 
where and are the vectors containing the local and global mass constraints given by the right hand side of (19) and (20), respectively. Letting for all we can view as an adjustment (in terms of scaling) of the approximations of .
In general, a tracked cell intersects more than one residing cell, and so , and hence the system is underdetermined and we have to select one of its solutions. This is done through the following minimisation problem: minimise subject to the local and global mass constraints (22). Moreover, since these quantities represent volumes, we also impose the constraint that each coefficient in is positive. Finally, we impose that each entry of the vector should be less than or equal to 2 (so that the maximum change is doubling a given volume). Essentially, this tells us that we want to achieve global and local mass conservation with minimal adjustments on the computed/approximated volumes, which makes sense since we assume that these intersections have been wellapproximated. In terms of computational cost, solving the minimisation problem is not too expensive since a tracked cell usually only intersects a few residing cells (as long as the velocity field and the mesh are not too irregular), and hence the matrix is usually sparse.
Remark 3.3 (Nonsolenoidal fields).
The Bchar method may also be applied on nonsolenoidal fields. In general, given a velocity field , we denote by the measure of the ball at time . We then solve for or by using a generalised Liouville’s formula [9]
Based on the increase or decrease in the volume , the radii of the tracked balls are then scaled accordingly. The adjusted porosity is then computed so that . These can then be used in (18) for the initial approximation of .
3.3. Summary of the Bchar method
To summarise, the Bchar method consists of approximating each cell by a collection of balls and their traceback regions by tracking back the centres of the balls to obtain the centres of the tracked balls . We then perform the following:
4. Numerical tests
In this section, we perform numerical tests on Cartesian type meshes for the advectionreaction equation (1). We start by performing tests in 2D, for which the numerical results presented are obtained using two methods:

polygonal ELLAM, obtained by approximating the cells and their traceback regions with polygons (see, e.g. Figure 1, right), with mass conservation achieved approximately, with a relative error less than , by performing volume adjustments as in [8]. Here, the polygonal intersections are computed using a general polygon clipper (GPC) library, obtained from http://www.cs.man.ac.uk/~toby/gpc/.

Bchar ELLAM, as described in Section 3. For this method, an has to be chosen to stop the iterations (21) before solving the optimisation problem (22). Figure 4 shows the relative error on the mass balances against for a typical test case, and indicates that a reasonable choice is (reducing the errors to about 5%). Further reduction does not bring much improvement.
For each of the test cases, we seek the concentration at time , i.e. . We also assume that the velocity field is provided. Numerical tests in 3D are then performed using the Bchar ELLAM.
The relative errors will be measured in the and norm, by providing for the quantities
where denote the norm in .
4.1. Numerical tests in 2D
The test cases in 2D are performed over the domain .
4.1.1. Test case 1
For the first test case, we consider a simple velocity field that simulates a translation along the axis, . Although the noflow boundary conditions are not satisfied, the final time is small enough so that no relevant characteristic traces outside the domain. The initial condition is set to be
Based on this initial condition and the given velocity, we expect the square block initially on the lower left corner of the domain to be transferred to the lower right corner of the domain, i.e.
We now compute and compare the approximate solutions using the polygonal ELLAM scheme and the Bchar ELLAM. For the Bchar ELLAM, 4 balls are used to approximate each cell. This will be performed starting on a grid with a time step of , and refined for 2 levels in space and time, leading to a test on a grid with a time step of . Upon looking at Figures 5 to 7, we see that the concentration profiles obtained from the polygonal ELLAM and the Bchar ELLAM are quite similar, with the Bchar ELLAM producing maximum concentrations which are slightly closer to 1, as compared to the polygonal ELLAM.
Now, we compare these methods in more detail by looking at Tables 1 and 2. As can be seen, the polygonal ELLAM and the Bchar ELLAM produce results that are quite close to one another. Moreover, upon measuring the CPU runtime (in seconds) for the total process of tracking, computing intersections, and performing volume adjustments, we see that the Bchar ELLAM gets to perform the simulations much faster compared to the polygonal ELLAM. This is mainly due to the fact that ball intersections are much cheaper to compute as compared to polygonal intersections.
Mesh  CPU time  

0.8  0.5175  4.8271e01  3.7277e01  
0.4  6.4640  3.4911e01  3.1673e01  
0.2  97.3994  2.4956e01  2.6898e01 
Mesh  CPU time  

0.8  0.1141  4.7637e01  3.8273e01  
0.4  0.4321  3.4889e01  3.3183e01  
0.2  3.5188  2.5558e01  2.9220e01 
Remark 4.1 (CPU runtime).
The CPU times are only used as an indication to show the advantage of the Bchar method over the polygonal ELLAM. The codes may not be fully optimised, but are implemented in a similar manner for both methods, by taking advantage of the vectorial capacities of MATLAB.
4.1.2. Test case 2
The second test case considers a velocity field . Here, is a divergencefree velocity field which simulates a rotation with some stretching, and the centre of this rotation is located at (see Figure 8).
The initial condition is set to be
Essentially, this assumes that we have a substance near the topleft corner of our domain (see Figure 9, left), being rotated, and somehow stretched for time units. Unlike the first test case, an exact solution is not available. Hence, we compare our results with a benchmark solution, obtained by solving (10) using an Euler method over a very fine grid (to be particular, 2 levels of refinement compared to the mesh being considered), with a very small time step . This is then projected onto the mesh being consideredin the case of Figure 9, right, a mesh consisting of squares.
As with the first test case, we start by comparing the concentration profiles obtained by solving (1) using the polygonal ELLAM and the Bchar ELLAM, with 4 balls being used to approximate each cell for the Bchar ELLAM. Based on Figures 10 to 12, we see that the concentration profile obtained from the Bchar ELLAM is very similar to those obtained from the polygonal ELLAM. Also, as with the first test case, the maximum concentrations for the Bchar ELLAM are closer to 1, compared to the polygonal ELLAM.
Upon performing a more rigorous comparison by looking at Tables 3 and 4, we note that the errors in both the and norm for both methods are quite close to each other. Also, the Bchar ELLAM performs much faster than the polygonal ELLAM.
Mesh  CPU time  

0.8  2.7441  7.3431e01  5.1027e01  
0.4  43.0472  6.3375e01  4.2258e01  
0.2  700.9942  4.9580e01  3.6537e01 
Mesh  CPU time  

0.8  0.1865  7.3138e01  5.0673e01  
0.4  1.3095  6.1391e01  4.1428e01  
0.2  14.5061  4.7916e01  3.5931e01 
4.1.3. Test case 3
Finally, we present a test case using the same velocity field as the second test, but now with a smooth initial condition, given by . Similar to the second test case, this assumes that we have majority of our substance near the topleft corner of our domain (see Figure 13, left). The benchmark solution is also obtained in the same way as the second test case, and is then projected onto a mesh consisting of squares (Figure 13, right).