An efficient implementation of mass conserving characteristic-based schemes in 2D and 3D

TL;DR

This paper introduces the ball-approximated characteristics (B-char) method, a computationally efficient algorithm for implementing mass-conserving characteristic-based schemes in 2D and 3D, using ball intersections instead of polytopes.

Contribution

The paper presents a novel ball-based approximation approach for characteristic schemes that reduces computational complexity while maintaining mass conservation.

Findings

Achieves similar accuracy to polytope-based methods

Significantly reduces computational cost in 3D implementations

Ensures local and global mass conservation through optimization

Abstract

In this paper, we develop the ball-approximated characteristics (B-char) method, which is an algorithm for efficiently implementing characteristic-based schemes in 2D and 3D. Core to the implementation of numerical schemes is the evaluation of integrals, which in the context of characteristic-based schemes with piecewise constant approximations boils down to computing the intersections between two regions. In the literature, these regions are approximated by polytopes (polygons in 2D and polyhedra in 3D) and, due to this, the implementation in 3D is nontrivial. The main novelty in this paper is the approximation of the regions by balls, whose intersections are much cheaper to compute than those of polytopes. Of course, balls cannot fully tessellate a region, and hence some mass may be lost. We perform some adjustments, and also solve an optimisation problem, in order to yield a scheme that is both locally and globally mass conserving. This algorithm can achieve results that are similar to those obtained from an implementation which uses polytopal intersections, with a much cheaper computational cost.

Tables

Table 1. Table 1. CPU runtime and errors in the concentration profiles, test case 1, polygonal ELLAM, final time T = 8 𝑇 8 T=8 .

Mesh	$\delta t$	CPU time	$E_1$	$E_2$
		(one time step)
$16\times 16$	0.8	0.5175	4.8271e-01	3.7277e-01
$32\times 32$	0.4	6.4640	3.4911e-01	3.1673e-01
$64\times 64$	0.2	97.3994	2.4956e-01	2.6898e-01

Table 2. Table 2. CPU runtime and errors in the concentration profiles, test case 1, B-char, final time T = 8 𝑇 8 T=8 .

Mesh	$\delta t$	CPU time	$E_1$	$E_2$
		(one time step)
$16\times 16$	0.8	0.1141	4.7637e-01	3.8273e-01
$32\times 32$	0.4	0.4321	3.4889e-01	3.3183e-01
$64\times 64$	0.2	3.5188	2.5558e-01	2.9220e-01

Table 3. Table 3. CPU runtime and errors in the concentration profiles, test case 2, polygonal ELLAM, final time T = 8 𝑇 8 T=8 .

Mesh	$\delta t$	CPU time	$E_1$	$E_2$
		(one time step)
$16\times 16$	0.8	2.7441	7.3431e-01	5.1027e-01
$32\times 32$	0.4	43.0472	6.3375e-01	4.2258e-01
$64\times 64$	0.2	700.9942	4.9580e-01	3.6537e-01

Table 4. Table 4. CPU runtime and errors in the concentration profiles, test case 2, B-char, final time T = 8 𝑇 8 T=8 .

Mesh	$\delta t$	CPU time	$E_1$	$E_2$
		(one time step)
$16\times 16$	0.8	0.1865	7.3138e-01	5.0673e-01
$32\times 32$	0.4	1.3095	6.1391e-01	4.1428e-01
$64\times 64$	0.2	14.5061	4.7916e-01	3.5931e-01

Table 5. Table 5. CPU runtime and errors in the concentration profiles, test case 3, polygonal ELLAM, final time T = 8 𝑇 8 T=8 .

Mesh	$\delta t$	CPU time	$E_1$	$E_2$
		(one time step)
$16\times 16$	0.8	2.7276	1.4141e-01	1.4244e-01
$32\times 32$	0.4	42.6126	9.3292e-02	9.9508e-02
$64\times 64$	0.2	724.9937	5.6584e-02	6.7150e-02

Table 6. Table 6. CPU runtime and errors in the concentration profiles, test case 3, B-char, final time T = 8 𝑇 8 T=8 .

Mesh	$\delta t$	CPU time	$E_1$	$E_2$
		(one time step)
$16\times 16$	0.8	0.2146	1.4961e-01	1.5055e-01
$32\times 32$	0.4	1.3015	9.1979e-02	9.8428e-02
$64\times 64$	0.2	16.0974	5.6735e-02	6.7733e-02

Table 7. Table 7. CPU runtime and errors in the concentration profile in 3D, B-char, final time T = 8 𝑇 8 T=8 .

Test case	$\delta t$	CPU time	$E_1$	$E_2$
		(one time step)
1	0.8	37.2347	4.8130e-01	4.0692e-01
2	0.8	63.5058	9.6106e-01	6.2141e-01
3	0.8	63.2127	2.3673e-01	2.4150e-01

Equations

$$\{\begin{array}{ll}\varphi \frac{\partial c}{\partial t}+\mathrm{div}(\mathbf{u}c)=0& \text{mbox}on{Q}_T:=\Omega \times (0,T),\\ c(\cdot ,0)=c_{\mathrm{ini}}& \text{mbox}on\Omega .\end{array}$$

$$\left\{\begin{array}[]{ll}\phi\dfrac{\partial c}{\partial t}+{\mathop{\rm div}}(\mathbf{u}c-\Lambda\nabla c)=f(c)&\qquad\mbox{ on }Q_{T}:=\Omega\times(0,T),\\ \Lambda\nabla c\cdot\mathbf{n}=0&\qquad\mbox{ on }\partial\Omega\times(0,T),\\ c(\cdot,0)=c_{\rm ini}&\qquad\mbox{ on }\Omega,\end{array}\right.$$

$$c_{\mathrm{ini}}\in L^{\infty }(\Omega ),$$

$$\varphi \in L^{\infty }(\Omega )\text{mbox}ispiecewiseconstanton\mathcal{M},\text{mbox}and$$

$$\displaystyle\mbox{there exists $\phi_{*}>0$ s.t. }\phi\geq\phi_{*}\mbox{ a.e.~{}on $\Omega$}.$$

$$\mathbf{u}\in L^{\infty }(0,T;L^2(\Omega {)}^d)\text{mbox}and\mathrm{div}\mathbf{u}\in L^{\infty }(Q_T),$$

$${\mathbf{u}}^{(n+1)}\in H(\mathrm{div},\Omega )\text{mbox}whichispiecewisepolynomialon\mathcal{M}.$$

$$\psi \frac{\partial c}{\partial t}=\frac{\partial (c\psi )}{\partial t}-c\frac{\partial \psi }{\partial t},$$

$${\int }_{t^{(n)}}^{t^{(n+1)}}$$

$$\displaystyle-\int_{t^{(n)}}^{t^{(n+1)}}\int_{\Omega}c({\boldsymbol{x}},t)\bigg{[}\phi({\boldsymbol{x}})\dfrac{\partial\psi}{\partial t}({\boldsymbol{x}},t)+\mathbf{u}({\boldsymbol{x}},t)\cdot\nabla\psi({\boldsymbol{x}},t)\bigg{]}d{\boldsymbol{x}}dt=0.$$

$$\varphi \frac{\partial \psi }{\partial t}+\mathbf{u}\cdot \nabla \psi =0\text{ on }\Omega \times (t^{(n)},t^{(n+1)}),$$

$${\int }_{\Omega }\varphi (\mathbf{x})(c\psi )(\mathbf{x},t^{(n+1)})d\mathbf{x}-{\int }_{\Omega }\varphi (\mathbf{x})(c\psi )(\mathbf{x},t^{(n)})d\mathbf{x}=0.$$

$${\int }_{\Omega }\varphi {\Pi }_{\mathcal{C}}{c}^{(n+1)}{\Pi }_{\mathcal{C}}z-{\int }_{\Omega }\varphi {\Pi }_{\mathcal{C}}{c}^{(n)}{\psi }_z(t^{(n)})=0\forall z\in X_{\mathcal{C}},$$

$$\phi\partial_{t}\psi_{z}+\mathbf{u}^{(n+1)}\cdot\nabla\psi_{z}=0\quad\text{ on }(t^{(n)},t^{(n+1)})\,,\mbox{ with $\psi_{z}(\cdot,t^{(n+1)})=\Pi_{{\mathcal{C}}}z$ }.$$

$$\dfrac{dF_{t}({\boldsymbol{x}})}{dt}=\dfrac{\mathbf{u}^{(n+1)}(F_{t}({\boldsymbol{x}}))}{\phi(F_{t}({\boldsymbol{x}}))}\quad\mbox{ for $t\in[-T,T]$},\qquad F_{0}({\boldsymbol{x}})={\boldsymbol{x}}.$$

$${\psi }_z(\cdot ,t^{(n)})={\Pi }_{\mathcal{C}}z(F_{{\delta t}^{(n+\frac{1}{2})}}(\cdot )).$$

$${\int }_K\varphi {\Pi }_{\mathcal{C}}{c}^{(n+1)}d\mathbf{x}={\int }_{\Omega }\varphi \sum _{M\in \mathcal{M}}{c}_M^{(n)}\text{mathbbm}{1}_M(\mathbf{x})\text{mathbbm}{1}_K(F_{{\delta t}^{(n+\frac{1}{2})}}(\mathbf{x}))d\mathbf{x},$$

$$|K{|}_{\varphi }{c}_K^{(n+1)}=\sum _{M\in \mathcal{M}}|M\cap F_{-{\delta t}^{(n+\frac{1}{2})}}(K){|}_{\varphi }{c}_M^{(n)},$$

$${\int }_{\Omega }\varphi (\mathbf{x})c(\mathbf{x},t^{(n+1)})d\mathbf{x}={\int }_{\Omega }\varphi (\mathbf{x})c(\mathbf{x},t^{(n)})d\mathbf{x}.$$

$$\sum _{K\in \mathcal{M}}{c}_K^{(n+1)}|K{|}_{\varphi }=\sum _{M\in \mathcal{M}}{c}_M^{(n)}|M{|}_{\varphi },$$

$$\sum _{K\in \mathcal{M}}|M\cap F_{-{\delta t}^{(n+\frac{1}{2})}}(K){|}_{\varphi }=|M{|}_{\varphi },$$

$$|K|c_K^{(1)}$$

$$=|F_{-\frac{1}{2}}(K)|$$

$$\displaystyle=|K|\qquad\mbox{(since $\mathbf{u}$ is divergence free)}.$$

$$|K|c_K^{(1)}=\sum _{M\in \mathcal{M}}|M\cap \tilde{K}|\times 1=|\tilde{K}|$$

$$c_K^{(1)}=\frac{|\tilde{K}|}{|K|}\ne 1.$$

$$\sum _{M\in \mathcal{M}}|M\cap \tilde{K}{|}_{\varphi }=|F_{-{\delta t}^{(n+\frac{1}{2})}}(K){|}_{\varphi }.$$

$$\frac{d}{dt}{\int }_{F_t(A)}\varphi (\mathbf{y})d\mathbf{y}={\int }_{F_t(A)}\mathrm{div}\mathbf{u}(\mathbf{y})d\mathbf{y}.$$

$$|F_{-{\delta t}^{(n+\frac{1}{2})}}(K){|}_{\varphi }=|K{|}_{\varphi }.$$

$${\rho }_K\sum _{s=1}^{n_K}|B_{K,s}{|}_{\varphi }=|K{|}_{\varphi }.$$

$${\rho }_K|{\hat{B}}_{K,s}{|}_{\varphi }={\rho }_K{\hat{\varphi }}_{K,s}|{\hat{B}}_{K,s}|.$$

$${\varphi }_K|B_{K,s}|={\hat{\varphi }}_{K,s}|{\hat{B}}_{K,s}|.$$

$${\hat{r}}_{K,s}=r_{K,s},|{\hat{B}}_{K,s}|=|B_{K,s}|.$$

Full text

1. Models and assumptions

1.1. Introduction

In this paper, we introduce an algorithm for implementing characteristic-based schemes for a pure advection model

[TABLE]

Here, $T>0$, $\Omega$ is a polytopal domain in $\mathbb{R}^{d}$ ($d\geq 1$), the porosity $\phi$, and the velocity $\mathbf{u}$ are given, with $\mathbf{u}\cdot\mathbf{n}=0$ on $\partial\Omega$. The unknown $c({\boldsymbol{x}},t)$ represents the amount of material (a fraction) present at $({\boldsymbol{x}},t)$. Note that the boundary is non-characteristic due to the assumption $\mathbf{u}\cdot\mathbf{n}=0$ on $\partial\Omega$, 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 advection-diffusion equation

[TABLE]

where the source term $f$ and the diffusion tensor $\Lambda$ are given. Although (1) is presented with $f=0$, it is not difficult to handle equations with nonzero reaction/source terms. This can be done, e.g., by a splitting approach: solve the pure advection equation first (using the method presented here), then solve the reaction equation.

Advection-diffusion equations of the form (2) are usually encountered in mathematical models for porous media flow (e.g. reservoir simulation, nuclear waste storage) [17, 26], and computational fluid dynamics (e.g. Navier-Stokes equations) [20], 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) [14], and will not be detailed in this paper. Here, we only focus on the implementation of characteristic-based schemes in (1), such as the Eulerian Lagrangian Localised Adjoint Method (ELLAM) and the Modified Method of Characteristics (MMOC). The advantage of these schemes stems from 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 [5] and the finite volume (FV) ELLAM [18], as well as a summary of their properties, have been presented in [27]. One of the major issues faced when implementing characteristic-based 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) [13]. 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 [11]. Another way to achieve local volume conservation for characteristic-based schemes has been proposed in [9]. This is particularly useful for schemes with piecewise constant approximations, such as hybrid and mixed finite volume type schemes [15]. As an example, in [9], it was used to perform adjustments to make the HMM–ELLAM schemes in [8] locally mass conserving. More details about the convergence analysis and implementation of GDM characteristic-based schemes for (2) and its applications to flows in porous media, can be found in [9, 10]. 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 to compute the intersection between two convex polyhedra efficiently, as in [6, 7, 19, 24]. 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 [12], 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 characteristic-based 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 conservations as constraints. We call this process the Ball-approximated characteristics (B-char) method. Due to its formulation, the B-char 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 B-char method is introduced in Section 3. In Section 4, numerical tests are first performed in 2D in order to compare the performance of the B-char method with the ELLAM scheme obtained from polygonal intersections, with volume adjustments as described in [9]. In these tests, we see that the B-char method yields very similar results to the polygonal intersections, with a much cheaper computational cost. The B-char method is applied on some benchmark test cases; here too the results demonstrate the accuracy of the method. Finally, numerical tests are performed to show the applicability of the B-char ELLAM in 3D.

1.2. Assumptions on the data, and numerical setting

We start by forming a mesh, i.e. a partition of $\Omega$ into polygonal (in 2D) or polyhedral (in 3D) sets. Following the notations in [14, Definition 7.2], we then denote $\mathcal{T}=({\mathcal{M}},{\mathcal{E}})$ to be the set of cells $K$ and faces (edges in 2D) $\sigma$ of our mesh, respectively. We also use $|K|$ to denote the volume (area in 2D) of a cell $K$. Throughout the article we assume the following properties:

[TABLE]

Assumption (3) simply states that the initial concentration inside the medium is bounded and that the porosity $\phi$ of the medium does not vanish, which is natural in physical applications. The piecewise constant assumption on $\phi$ is also satisfied in practical applications; the value of $\phi$ on a cell $K$ will be denoted by $\phi_{K}$. As in [9], 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 [9, 10], and will not be discussed in further detail for this paper. We replace, in the weak formulation of the model, the continuous (infinite-dimensional) spaces and corresponding operators by a discrete (finite-dimensional) space and function reconstructions. We then define a space-time discretisation ${\mathcal{C}}=(X_{\mathcal{C}},\Pi_{\mathcal{C}},\mathcal{I}_{{\mathcal{C}}},(t^{(n)})_{n=0,\dots,N})$, where

• •

$X_{{\mathcal{C}}}$ is a finite-dimensional real space, describing the unknowns of the chosen scheme,

• •

$\Pi_{\mathcal{C}}:X_{\mathcal{C}}\to\mathbb{P}^{0}({\mathcal{M}})$ is a linear operator that reconstructs a piecewise constant function on the mesh ${\mathcal{M}}$ from the unknowns,

• •

$\mathcal{I}_{\mathcal{C}}$ is a rule to map $c_{\rm ini}$ onto an element $\mathcal{I}_{\mathcal{C}}c_{\rm ini}\in X_{{\mathcal{C}}}$,

• •

$0=t^{(0)}<t^{(1)}<\dots<t^{(N)}=T$ are the time steps, and we let ${\delta\!t}^{(n+\frac{1}{2})}=t^{(n+1)}-t^{(n)}$.

Different choices of ${\mathcal{C}}$ lead to different schemes (e.g. finite volume based methods, including hybrid ones with face unknowns like HMM [15], or mass-lumped finite element methods [28]).

Finally, we assume that

[TABLE]

and that $\mathbf{u}$ is approximated on each time interval $(t^{(n)},t^{(n+1)})$ by a function

[TABLE]

Remark 1.1** (Approximation of the velocity field).**

Although $\mathbf{u}$ is given in (1), we use an approximation for the velocity field $\mathbf{u}$ in order to include the more general case where $\mathbf{u}$ comes from solving a PDE coupled to (2). For example, for flows in porous media, $\mathbf{u}$ usually comes from Darcy’s law: Given a source term $g,\mathbf{u}$ should satisfy the PDE $-{\mathop{\rm div}}(\mathbf{u})=g$ on $\Omega$, 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., $\int_{\Omega}\phi d{\boldsymbol{x}}$.

2. ELLAM scheme for the advection–reaction equation

We multiply (1) with a sufficiently smooth function $\psi$, and perform integrations by parts. Using the identity

[TABLE]

(1) gives, for any time interval $(t^{(n)},t^{(n+1)})$,

[TABLE]

To simplify the second term on the left hand side of the above equation, the ELLAM requires that test functions $\psi$ satisfy

[TABLE]

with $\psi(\cdot,t^{(n+1)})$ given. The advection equation (1) then leads to the relation

[TABLE]

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 $\Pi_{\mathcal{C}}$ applied to trial and test vectors in $X_{\mathcal{C}}$.

Definition 2.1** (ELLAM scheme).**

Given a space-time discretisation ${\mathcal{C}}$, the ELLAM scheme for (1) reads as: find $(c^{(n)})_{n=0,\ldots,N}\in X_{{\mathcal{C}}}^{N+1}$ such that $c^{(0)}=\mathcal{I}_{{\mathcal{C}}}c_{\rm ini}$ and, for all $n=0,\ldots,N-1$, $c^{(n+1)}$ satisfies

[TABLE]

where $\psi_{z}$ is the solution to

[TABLE]

Define the flow $F_{t}:\Omega\to\Omega$ such that, for a.e. ${\boldsymbol{x}}\in\Omega$,

[TABLE]

Under Assumptions (3) and (5), the existence of this flow is proved in [10, Lemma 5.1]. The solution to (9) is then understood in the sense: for $t\in(t^{(n)},t^{(n+1)}]$ and a.e. ${\boldsymbol{x}}\in\Omega$, $\psi_{z}({\boldsymbol{x}},t)=\Pi_{\mathcal{C}}z(F_{t^{(n+1)}-t}({\boldsymbol{x}}))$. In particular,

[TABLE]

The construction of the B-char method in Section 3 will draw inspiration from a physical interpretation of the ELLAM, which uses the fact that $\Pi_{\mathcal{C}}$ is a piecewise-constant reconstruction on a given mesh ${\mathcal{M}}$. For each cell $K\in{\mathcal{M}}$, we assume that there is $z_{K}\in X_{\mathcal{C}}$ such that $\Pi_{\mathcal{C}}z_{K}=\mathbbm{1}_{K}$, where $\mathbbm{1}_{K}$ is the function that has a value of 1 in $K$, and 0 elsewhere. Writing $\Pi_{{\mathcal{C}}}c^{(k)}=\sum_{M\in{\mathcal{M}}}c_{M}^{(k)}\mathbbm{1}_{M}$ and taking $z_{K}$ as test function, (8) and (11) give

[TABLE]

which reduces to

[TABLE]

where $|E|_{\phi}=\int_{E}\phi$ is the available porous volume in a set $E\subset\mathbb{R}^{d}$. The term on the right hand side of (12) tells us that the amount of material $c_{K}^{(n+1)}$ present in a particular cell $K\in{\mathcal{M}}$ at time $t^{(n+1)}$ is obtained by intersecting a tracked cell $F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)$ and a residing cell $M$. This intersection can be interpreted as locating where the material in cell $K$ comes from, hence back-tracking the cell $K$ to $F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)$, measuring which fraction of the material $c^{(n)}_{M}$ is taken from each $M\in{\mathcal{M}}$ (by taking their intersection), and depositing this fraction into the cell $K$.

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 $c$ is dictated by the amount of inflow/outflow and by the source term. In this case, due to the no-flow boundary conditions and the absence of a source term, this simply means that the amount of substance present at time $t^{(n+1)}$ should be the same as the amount of substance present at time $t^{(n)}$. The desired equation is thus given by

[TABLE]

It can easily be checked that the ELLAM scheme satisfies this property. Indeed, taking the sum over all $K\in{\mathcal{M}}$ in (12) yields

[TABLE]

which is the discrete form of (13).

Remark 2.2** (Achieving global mass conservation).**

We note here that the ELLAM scheme achieves global mass conservation due to

[TABLE]

for all $M\in{\mathcal{M}}$. An analogue of this identity will be needed to ensure that the B-char 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 $\int_{\widehat{K}}\phi\Pi_{{\mathcal{C}}}c^{(n)}d{\boldsymbol{x}}$ for each cell $K$, where $\widehat{K}=F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)$. In general, the region $\widehat{K}$ (see Figure 1, left) cannot be exactly described and hence, in the literature, it was approximated by polygons obtained from back-tracking the vertices, together with a number of points along the edges of the cell $K$. Figure 1 (right) gives an illustration of the approximate trace-back region $\widetilde{K}$ obtained by tracking the vertices, together with the edge midpoints of the cell $K$.

Remark 2.3** (Reconstruction of polytopes).**

In 2D, most of the time, we can reconstruct the polygons approximating the trace-back region by following the tracked points in the same order as the original points, since it gives a well-defined 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, $|\widetilde{K}|_{\phi}\neq|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)|_{\phi}$. 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 $\phi=1$ and $c_{\rm{ini}}=1$. In this test case, the exact solution is given by $c({\boldsymbol{x}},t)=1$. In theory, upon implementing an ELLAM scheme with piecewise constant approximations for the unknown $c$, we should have the following simplified form of (12) at the first time step:

[TABLE]

However, due to the approximation of the trace-back region, we only have

[TABLE]

and thus

[TABLE]

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 $\widetilde{K}$ in order to yield $|\widetilde{K}|_{\phi}=|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)|_{\phi}$, which we shall define as the local volume constraint for $K$. Several adjustment strategies which would lead to local mass conservation have been studied, as in [1, 9, 11]. In particular, for local mass conservation to be achieved, we should have, for all $K\in{\mathcal{M}}$,

[TABLE]

For simplicity of exposition, we consider solenoidal fields, so that ${\mathop{\rm div}}\mathbf{u}=0$. The generalised Liouville’s formula [10, equation (26)] determines the evolution of available porous volume in a given domain: for any measurable set $A\subset\Omega$,

[TABLE]

Under the assumption that ${\mathop{\rm div}}\mathbf{u}=0$, (16) with $A=K$ gives

[TABLE]

3. B-char 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 Ball-approximated Characteristics (B-Char). For each cell $K$, we choose $n_{K}$ points $C_{K,s},(s=1,\dots,n_{K})$ in its interior. We then assume that each of these points represents centers of disjoint balls $B_{K,s}$, with radius $r_{K,s}$, which are strictly inside cell $K$. The idea now is to distribute the porous volume in each cell $K$ over the balls $B_{K,s}$. To do so, we introduce a porous density $\rho_{K}$ so that

[TABLE]

The porous density $\rho_{K}$ is fixed and does not change throughout the tracking. Here, the quantity $\rho_{K}|B_{K,s}|_{\phi}$ may be interpreted as an equivalent porous volume inside the ball $B_{K,s}$. 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 a ball, the points $C_{K,s}$ are then tracked by solving (10) to obtain $\widehat{C}_{K,s}$, which will be treated as the center of the tracked ball $\widehat{B}_{K,s}$ (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 $B_{K,s}$ are not too large. Typically, a good enough approximation is obtained if the maximum radius of the balls $B_{K,s}$ is one-fourth or less of the diameter of the smallest face in the cell $K\in{\mathcal{M}}$.

Since $\rho_{K}$ does not change throughout the tracking, the equivalent porous volume inside the tracked ball $\widehat{B}_{K,s}$ is given by $\rho_{K}|\widehat{B}_{K,s}|_{\phi}$. Denoting by $\widehat{\phi}_{K,s}$ the (unknown, at this stage) average porosity over the region covered by $\widehat{B}_{K,s}$, we may write

[TABLE]

Now, upon applying the generalised Liouville formula (16) (recalling that ${\mathop{\rm div}}\mathbf{u}=0$) with $A=B_{K,s}$, we have $|B_{K,s}|_{\phi}=|\widehat{B}_{K,s}|_{\phi}$, which implies

[TABLE]

We now describe how to find the radius $\widehat{r}_{K,s}$ (and thus volume) of the tracked ball $\widehat{B}_{K,s}$, considering separately the case where $\phi$ is constant and where $\phi$ varies in the domain.

• •

$\phi$ constant in $\Omega$. In this case, $\widehat{r}_{K,s}$ can be exactly computed. (20) implies that the volume, and thus the radius, of the tracked ball is unchanged, i.e.

[TABLE]

• •

$\phi$ non-constant. The radius $\widehat{r}_{K,s}$ cannot be exactly computed in general, only approximated with some additional computational cost. In this situation, aside from the center $C_{K,s}$, we also track points on the circumference of $B_{K,s}$. Considering the two-dimensional case for illustration, letting $C_{K,s}=(x_{K,s},y_{K,s})$, we also track the four cardinal points $C_{K,s_{1}}=(x_{K,s},y_{K,s}+r_{K,s})$, $C_{K,s_{2}}=(x_{K,s},y_{K,s}-r_{K,s})$, $C_{K,s_{3}}=(x_{K,s}+r_{K,s},y_{K,s})$, $C_{K,s_{4}}=(x_{K,s}-r_{K,s},y_{K,s})$. The radius $\widehat{r}_{K,s}$ of the tracked ball $\widehat{B}_{K,s}$ is then approximated by

[TABLE]

where $(\widehat{r}_{K,s_{j}})_{j=1,\dots,4}$ denote the distances between the tracked center $\widehat{C}_{K,s}$ and the tracked cardinal point $(\widehat{C}_{K,s_{j}})_{j=1,\ldots,4}$. Of course, we may track more than 4 points on the circumference of the ball in order to get a better approximation of $\widehat{r}_{K,s}$. Finding the optimal number of points to be tracked in order to get a good approximation for $\widehat{r}_{K,s}$ depends on how strongly we expect the velocity field to distort the tracked region, and will be a topic for future research. The volume $|\widehat{B}_{K,s}|$ can then be obtained from $\widehat{r}_{K,s}$ defined by (22).

After computing $|\widehat{B}_{K,s}|$, we use (20) to determine the value of $\widehat{\phi}_{K,s}$. Similarly to (18), we can then write the porous volume inside each tracked cell in the following form:

[TABLE]

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 $|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)\cap M|_{\phi}\approx V_{\widehat{K},M}$.

We start by recalling that $|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)\cap M|_{\phi}$ is interpreted as the amount of material in a cell $K$ that comes from a residing cell $M$. In the context of approximation by balls, this reads: each ball $B_{K,s}$ contains an amount of material from some residing balls $B_{M,m}$ transported by the flow, given by $\rho_{M}\phi_{M}|\widehat{B}_{K,s}\cap B_{M,m}|$.

Remark 3.1**.**

The choice of $\rho_{M}\phi_{M}$ (instead of $\rho_{K}\phi_{K}$) in the initial approximation $\rho_{M}\phi_{M}|\widehat{B}_{K,s}\cap B_{M,m}|$ comes from the interpretation that the amount of material present in $B_{K,s}$ is obtained by measuring how much of the material is taken from each $B_{M,m}$, and by depositing this material into the ball $B_{K,s}$.

An initial approach for approximating $V_{\widehat{K},M}$ would then involve taking the sum of the masses of the balls in a residing cell $M$, intersected with the tracked balls that originated from cell $K$, that is $\sum_{s=1}^{n_{K}}\sum_{m=1}^{n_{M}}\rho_{M}\phi_{M}|\widehat{B}_{K,s}\cap B_{M,m}|$. However, since there are gaps between the residing balls, $\sum_{M\in{\mathcal{M}}}\sum_{m=1}^{n_{M}}|\widehat{B}_{K,s}\cap B_{M,m}|\neq|\widehat{B}_{K,s}|$. 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

[TABLE]

which represents the fraction of the mass in $\widehat{B}_{K,s}$ that comes from $B_{M,m}$. From this, we then see that

[TABLE]

is the actual amount of mass in the tracked ball $\widehat{B}_{K,s}$ that comes from $B_{M,m}$. The quantity $V_{\widehat{K},M}$ is then computed by taking the sum over all tracked balls $\widehat{B}_{K,s}$ and residing balls $B_{M,m}$, given by

[TABLE]

3.2. Mass conservation for the B-char method

Since $V_{\widehat{K},M}$ are approximations to $|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)\cap M|_{\phi}$, in order to achieve local mass conservation, we should have an analogue of (15), given by $\sum_{M\in{\mathcal{M}}}V_{\widehat{K},M}=|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)|_{\phi}$. We can easily check that $V_{\widehat{K},M}$ in (25) satisfies this relation by using (23). Hence, the approximation of $|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)\cap M|$ by $V_{\widehat{K},M}$ in (25) leads to a scheme that is locally mass conserving. However, $\sum_{K\in{\mathcal{M}}}V_{\widehat{K},M}\neq|M|_{\phi}$, 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 $V_{\widehat{K},M}$.

Remark 3.2** (Overlapping balls).**

Two given tracked balls, say $\widehat{B}_{K,s_{i}}$ and $\widehat{B}_{L,s_{j}}$, may overlap. In this case, the mass in the overlapping region $\widehat{B}_{K,s_{i}}\cap\widehat{B}_{L,s_{j}}$ would be allocated twice (once each for the contribution to $B_{K,s_{i}}$ and $B_{L,s_{j}}$). Recalling Remark 3.1, this would imply that excess mass is deposited into cell $K$ or $L$ (or both), since these are the cells where the tracked balls $\widehat{B}_{K,s_{i}}$ and $\widehat{B}_{L,s_{j}}$ come from, respectively. This corresponds to a loss in mass that should be deposited in another cell, say $M$. This is implicitly reflected by the failure of the initial approximation $V_{\widehat{K},M}$ defined in (25) to achieve global mass conservation. The algorithm provided below to adjust $V_{\widehat{K},M}$, and recover both global and local mass conservations, also corrects the excessive and missing mass deposits induced by overlapping balls.

Upon indexing each tracked cell $K_{i}$ and each residing cell $M_{j}$, $i,j=1,\dots n_{c}$, let $A^{(n)}$ be the matrix with entries $a_{ij}^{(n)}$, where $a_{ij}^{(0)}=V_{\widehat{K}_{i},M_{j}}$. In short, each entry $a_{ij}^{(n)}$ of the matrix $A^{(n)}$ gives an approximation of the volume $|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K_{i})\cap M_{j}|_{\phi}$. In terms of the matrix $A^{(n)}$, we would thus need

[TABLE]

in order to achieve (15), which will lead to local mass conservation. Owing to (17), the right hand side $|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K_{i})|_{\phi}$ of (26) is equal to $|K_{i}|_{\phi}$, which can be easily and exactly computed since $\phi$ is constant in each cell. Now, in order to have a globally mass conserving scheme, we would need to satisfy (14), which in this context is equivalent to

[TABLE]

To build a matrix which satisfies (26) and (27), we start with the assumption that $\sum_{i=1}^{n_{c}}a_{ij}^{(0)}>0$ for $j=1,\dots,n_{c}$. This means that at least one of the balls $B_{M_{j},\ell}$ of a residing cell $M_{j}$ intersects with at least one of the balls $B_{\widehat{K}_{i},s}$ that has been tracked from $K_{i}$. The only time $\sum_{i=1}^{n_{c}}a_{ij}^{(0)}=0$ holds is when the flow traces everything into voids. For example, in Figure 3, none of the balls $\widehat{B}_{K,\ell}$ tracked from $K$, nor did any of the other balls that were tracked from time $t^{(n+1)}$, intersect with the residing balls $B_{M,\ell}$ in $M$, and so $\sum_{i=1}^{n_{c}}a_{ij}^{(0)}=0$ for the $j$ corresponding to the particular cell $M$. 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 $\sum_{i=1}^{n_{c}}a_{ij}^{(0)}=0$ 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.3** (Presence of inflow or outflow).**

In the presence of inflows or outflows, we may avoid $\sum_{i=1}^{n_{c}}a_{ij}^{(0)}=0$ by one of the following options:

1. a.

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

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

[TABLE]

which yields

[TABLE]

Next, we set

[TABLE]

so that

[TABLE]

Essentially, the adjustments perform the following two-step process:

• •

Firstly, we redistribute, according to a proportion evaluated by $a_{ij}^{(n)}$, the mass along each intersecting region so that global mass conservation (27) is achieved.

• •

Next, we redistribute the mass so that local mass conservation (26) is achieved.

Intuitively, we can see that these adjustments involve re-distributing 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 $N$ iterations and arriving at the matrix $A^{(N)}$, solving a minimisation problem. In practice, we found that taking $N$ to be such that the maximum error in mass conservation is at most $5\%$ is sufficient to give a good initial approximation.

We then assign an unknown corresponding to each entry of $A^{(N)}$, which gives us $n_{c}\times n_{c}$ unknowns. For $i,j$ such that $a_{ij}^{(N)}=0$, the corresponding unknown is fixed to 0. This tells us that if no intersection has been detected between a tracked cell $\widehat{K}_{i}$ and a residing cell $M_{j}$, 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 $n_{z}$ in $A^{(N)}$. From the nonzero entries of $A^{(N)}$ construct a $2n_{c}\times n_{z}$ matrix $\hat{A}$ in the following manner: Write

[TABLE]

and denote by $\hat{\mathbf{r}}_{i}$ ($i=1,\dots,n_{c}$) the row vector of size $\leq n_{c}$ obtained by removing the zero entries in $\mathbf{r}_{i}$. The first $n_{c}$ rows of $\hat{A}$ are then formed by the $n_{c}\times n_{z}$ matrix

[TABLE]

That is, we stagger the vectors $(\hat{\mathbf{r}}_{i})_{i=1,\ldots,n_{c}}$ so that the coefficients of $\hat{\mathbf{r}}_{j+1}$ start at the column after the last coefficient of $\hat{\mathbf{r}}_{j}$, and we pad each row with zeros to ensure we obtain an $n_{c}\times n_{z}$ matrix.

The latter $n_{c}$ rows of the matrix $\hat{A}$ are then formed in the following manner: for the $n_{c}+m$th row, we look for the nonzero entries $a_{j,m}$ ($j=1,\dots n_{c}$) in column $m$ of $A$. For each corresponding $j$, we then find the column $\hat{\jmath}$ corresponding to where the coefficient $a_{j,m}$ resides in the first $n_{c}$ rows of $\hat{A}$. We then set $\hat{a}_{n_{c}+m,\hat{\jmath}}=a_{j,m}$. As an example, if

[TABLE]

then $\hat{A}$ would be (with the line separating the first $n_{c}$ rows from the last $n_{c}$ rows of $\hat{A}$)

[TABLE]

In practice, the last $n_{c}$ rows of $\hat{A}$ are assembled simultaneously with its first $n_{c}$ rows. In the example above, after setting $a,b,c$ as the first, second and third entries of the first row, $a,b,c$ are simultaneously set to be the first, second and third entries of the $n_{c}+1,n_{c}+2,$ and $n_{c}+3$th row, respectively. After which, when $d$ and $e$ 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 $n_{c}+1$ and the $n_{c}+3$th row. In this manner, the whole matrix $\hat{A}$ is rather easy to assemble.

The essential property of $\hat{A}$ is that when it is multiplied by an $n_{z}\times 1$ vector $\mathbf{1}$ consisting of all ones, we have

[TABLE]

The sums in the right-hand side correspond to the quantities that must be fixed to certain values in order to achieve local and global mass balance, see (26) and (27). To obtain local and global mass conservation, we therefore solve the system

[TABLE]

where $\mathbf{b}_{\mathrm{loc}}$ and $\mathbf{b}_{\mathrm{glob}}$ are the vectors containing the local and global mass constraints given by the right hand side of (26) and (27), respectively. Letting $\mathbf{x}=(x_{j})_{j=1,\ldots,n_{c}}$ for all $j=1,\ldots,n_{c}$ we can view $1+x_{j}$ as an adjustment (in terms of scaling) of the approximations $a_{ij}^{(N)}$ of $|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K_{i})\cap M_{j}|_{\phi}$.

In general, a tracked cell intersects more than one residing cell, and so $n_{z}\geq 2n_{c}$, and hence the system is underdetermined and we have to select one of its solutions. This is done through the following minimisation problem: minimise $\mathbf{x}^{T}\mathbf{x}$ subject to the local and global mass constraints (29). Moreover, since these quantities represent volumes, we also impose the constraint that each coefficient in $(\mathbf{1}+\mathbf{x})$ is positive. Finally, we impose that each entry of the vector $(\mathbf{1}+\mathbf{x})$ 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 well-approximated. 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 $\hat{A}$ is usually sparse.

Remark 3.4** (Non-solenoidal fields).**

The B-char method may also be applied with non-solenoidal fields. In general, given a velocity field $\mathbf{u}$, the idea is to approximate the radius $\widehat{r}_{K,s}$, and hence the volume $|\widehat{B}_{K,s}|$, by (22). The quantity $\widehat{\phi}_{K,s}$ is then obtained from an analogue of (20), which comes from applying the generalised Liouville formula (16). The quantity $\rho_{K}\widehat{\phi}_{K,s}$ and the tracked balls $\widehat{B}_{K,s}$ with radius $\widehat{r}_{K,s}$ can then be used in (25) for the initial approximation of $|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)\cap M|_{\phi}$.

3.3. Summary of the B-char method

To summarise, the B-char method consists of approximating each cell $K\in{\mathcal{M}}$ by a collection of balls $B_{K,s}$ and their trace-back regions $F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)$ by tracking back the centres of the balls $B_{K,s}$ to obtain the centres of the tracked balls $\widehat{B}_{K,s}$. We then perform the steps outlined in Algorithm 1:

4. Numerical tests

In this section, we perform numerical tests on Cartesian type meshes for the pure advection equation (1). We start by performing tests in 2D, for which the numerical results presented are obtained using two methods:

1. (1)

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

B-char ELLAM, as described in Section 3. For this method, an $N$ has to be chosen to stop the iterations (28) before solving the optimisation problem (29). Figure 4 shows the relative error on the mass balances against $N$ for a typical test case, and indicates that a reasonable choice is $N=10$ (reducing the errors to about 5%). Further reduction does not bring much improvement.

For the first three test cases, we seek the concentration at time $T=8$, i.e. $c({\boldsymbol{x}},8)$. We also assume that the velocity field $\mathbf{u}$ is provided. These simple test cases aim at showing that the B-char ELLAM can achieve, with a much cheaper computational cost, numerical results that are essentially the same as those obtained from the polygonal ELLAM. Following these, we then use the B-char ELLAM for some 2D benchmark test cases, to show the robustness of the numerical scheme. Simple numerical tests (similar to the first three test cases in 2D) are then performed to show the applicability and efficiency of the B-char ELLAM in 3D.

The relative errors will be measured in the $L^{1}$ and $L^{2}$ norm, by providing for $p=1,2$ the quantities

[TABLE]

where $\left\|\cdot\right\|_{p}$ denote the norm in $L^{p}(\Omega)$.

4.1. Numerical tests in 2D

The test cases in 2D are performed over the domain $\Omega=(0,1)\times(0,1)$.

4.1.1. Test case 1

For the first test case, we consider a simple velocity field that simulates a translation along the $x$ axis, $\mathbf{u}=(\frac{1}{16},0)$. Although the no-flow boundary conditions are not satisfied, the final time $T=8$ is small enough so that no relevant characteristic traces outside the domain. The initial condition is set to be

[TABLE]

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.

[TABLE]

We now compute and compare the approximate solutions using the polygonal ELLAM scheme and the B-char ELLAM. For the B-char ELLAM, 4 balls are used to approximate each cell. This will be performed starting on a $16\times 16$ grid with a time step of $\delta t=0.8$, and refined for 2 levels in space and time, leading to a test on a $64\times 64$ grid with a time step of $\delta t=0.2$. Upon looking at Figures 5 to 7, we see that the concentration profiles obtained from the polygonal ELLAM and the B-char ELLAM are quite similar, with the B-char 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 B-char ELLAM produce results that are quite close to one another. Moreover, upon measuring the CPU runtime (in seconds) for one time step for the total process of tracking, computing intersections, and performing volume adjustments, we see that the B-char 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.

Remark 4.1** (CPU runtime).**

The CPU times are only used as an indication to show the advantage of the B-char 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 $\mathbf{u}=((1-2y)(x-x^{2}),-(1-2x)(y-y^{2}))$. Here, $\mathbf{u}$ is a divergence-free velocity field which simulates a rotation with some stretching, and the centre of this rotation is located at $(0.5,0.5)$ (see Figure 8).

The initial condition is set to be

[TABLE]

Essentially, this assumes that we have a substance near the top-left corner of our domain (see Figure 9, left), being rotated, and somehow stretched for $t=8$ 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 $\delta t=0.001$. This is then projected onto the mesh being considered$-$in the case of Figure 9, right, a mesh consisting of $16\times 16$ squares.

As with the first test case, we start by comparing the concentration profiles obtained by solving (1) using the polygonal ELLAM and the B-char ELLAM, with 4 balls being used to approximate each cell for the B-char ELLAM. Based on Figures 10 to 12, we see that the concentration profile obtained from the B-char ELLAM is very similar to those obtained from the polygonal ELLAM. Also, as with the first test case, the maximum concentrations for the B-char 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 $L^{1}$ and $L^{2}$ norm for both methods are quite close to each other. Also, the B-char ELLAM performs much faster than the polygonal ELLAM.

4.1.3. Test case 3

Now, we present a test case using the same velocity field as the second test, but now with a smooth initial condition, given by $c({\boldsymbol{x}},0)=\exp(-10((x-\frac{1}{4})^{2}+(y-\frac{3}{4})^{2}))$. Similar to the second test case, this assumes that we have majority of our substance near the top-left 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 $16\times 16$ squares (Figure 13, right).

From Figures 14 to 16 and Tables 5 and 6, a similar observation can be made as with the first two test cases: the concentration profiles obtained from the B-char ELLAM are very close to those obtained from the polygonal ELLAM, and the B-char ELLAM also performs much faster than the polygonal ELLAM. It is also notable that due to the continuous initial condition, the errors for this test case are smaller than those obtained from the second test case (by more than a factor of 5).

In general, we see that for both the polygonal ELLAM scheme and B-char ELLAM, numerical diffusion is most prominent on the coarse $16\times 16$ grid. This numerical diffusion becomes less prominent as the grid and the time step are refined. Moreover, for all test cases, using the B-char ELLAM achieves the same level of accuracy as the polygonal ELLAM, while at the same time, reducing the required computational cost (up to 40 times faster on the finest mesh in our tests).

4.1.4. Solid body rotation test case

Now that we have established that the B-char ELLAM achieves the same level of accuracy as the polygonal ELLAM with a cheaper computational cost, we perform a numerical simulation on a benchmark test case, the solid body rotation problem proposed in [23]. Instead of the domain $\Omega=(0,1)\times(0.5,1.5)$ considered therein, we retain our domain at $\Omega=(0,1)\times(0,1)$, and modify the initial conditions accordingly, as done in [4]. To be specific, the initial condition is made up of three components: a bump described by

[TABLE]

together with a cone and a slotted cylinder of radius $0.15$ and height $1$, centered at $(0.5,0.25)$ and $(0.5,0.75)$, respectively. The slot in the cylinder is created by removing the region $[0.475,0.525]\times[0.6,0.85]$ (see Figure 17, left). For this test case, the velocity field is given by $\mathbf{u}=((0.5-y),(x-0.5))$; we note that it does not satisfy the no-flow boundary conditions $\mathbf{u}\cdot\mathbf{n}=0$ on the boundary of the domain, but this has no impact on the B-char method since none of the balls that intersect non-zero values of the solution are tracked near the boundary. This velocity field simulates a rotation, and completes a revolution at time $T=2\pi$. Hence, we expect the solution profile at time $T=2\pi$ to be identical to the initial condition. Numerical simulations are performed using the B-char ELLAM with 4 balls in each cell. As with the test cases provided in [4, 21], we take a mesh with $128\times 128$ cells, but to take advantage of the ELLAM, we use a much larger time step of $\delta t=2\pi/10$, compared to $\delta t=2\pi/810$ or smaller in the literature.

As can be seen in Figure 17, the numerical solution obtained after one cycle of revolution is almost the same as the initial condition. In particular, the computed relative errors are $E_{1}=1.3630$e-$01$ and $E_{2}=2.1178$e-$01$. The absolute errors obtained in [4] (using a different scheme) in $L^{1}$ and $L^{2}$ norms are $1.7$e-$02$ and $7.5$e-$02$, respectively. These absolute errors correspond to relative errors of $1.4560$e-$01$ and $2.4020$e-$01$, respectively. These show that, despite its cheap computational cost, the B-char ELLAM method produces solutions with an accuracy that is comparable to those obtained in the literature. Moreover, the numerical solution preserves the slot in the cylinder, which indicates that the numerical diffusion introduced by the B-char ELLAM is minimal.

4.1.5. Deformational flow test

Finally, we present a test case for which the velocity field gives a very strong deformation. In [22], such a velocity field was given in terms of spherical coordinates, but [25] provided an analogue in Cartesian coordinates:

[TABLE]

This velocity field simulates a deformational flow, which reaches zero and then changes direction at half time $T/2$, resulting in a reversal of the flow field. We expect the solution profile to be the most highly deformed at $T/2$, and to be identical to the initial condition at time $T$. Following the test proposed in [16, Section 5.7.4], we set $T=5$ and take our initial condition to be a cosine bell function:

[TABLE]

We perform numerical tests on the following meshes: a mesh with $64\times 64$ cells and a refinement which consists of $128\times 128$ cells. In order to take advantage of ELLAM, we take larger time steps of $\delta t=0.5$ and $\delta t=0.25$, respectively, compared to $\delta t=0.01$ in the literature. As expected, in our numerical simulations, we obtain the most deformed solution profile at time $T=2.5$ (see Figure 18). As noticed in Figures 19 and 20, right, the strong distortion at time $T=2.5$ entails a distortion of the numerical solution at final time $T=5$. This distortion is more prominent for the coarse mesh. Aside from this, we notice the presence of some numerical diffusion, since the maximum amplitude of the solution has decreased, by around $27\%$ for the coarse mesh, and $18\%$ for the refined mesh. These issues were also encountered in [16] (with similar levels of distortion and diffusion), and can be resolved by further refinement of the mesh.

4.2. Numerical tests in 3D

In 3D, we perform the numerical tests over the domain $\Omega=(0,1)\times(0,1)\times(0,1)$. These are only performed using the B-char ELLAM, by approximating each cell with 8 balls. The polygonal ELLAM was not used for the numerical tests here as an easy/efficient implementation of intersecting a convex and non-convex polyhedron is not readily available. The numerical tests are performed over a mesh with $16\times 16\times 16$ cubes, with a time step $\delta t=0.8$.

As with the 2D case, the first test case involves a translation about the $x$ axis by considering the velocity field $\mathbf{u}=(\frac{1}{16},0,0)$. The concentration is initially set to be at the cubic block on the lower left corner of the cube, given by

[TABLE]

For the second test case, we consider a translation along the $z$ axis, accompanied by some rotation and stretching along the $x$ and $y$ directions by considering the velocity field $\mathbf{u}=((1-2y)(x-x^{2}),-(1-2x)(y-y^{2}),\frac{1}{16})$. Initially, the concentration is located on a cylinder described by

[TABLE]

Finally, the third test case also uses the velocity field $\mathbf{u}=((1-2y)(x-x^{2}),-(1-2x)(y-y^{2}),\frac{1}{16})$, with a continuous initial condition

[TABLE]

for which majority of the substance is concentrated around a neighborhood of the point $(\frac{1}{4},\frac{3}{4},\frac{3}{16})$.

We note that for all test cases, the no-flow boundary conditions are not satisfied, but the final time $T=8$ is such that none of the relevant characteristics trace outside the domain. For the first test case, an exact solution is available, given by

[TABLE]

For the second and third test cases, benchmark solutions are computed using a similar method described in the second test case in Section 4.1.

Here, we note that although the simulations were only performed over a mesh with $16\times 16\times 16$ cubes, with a time step $\delta t=0.8$, the errors obtained in the numerical simulations in 3D is of a similar magnitude as those obtained in 2D. As with the tests in 2D, by using a continuous initial condition in the third test case, the error dropped by at least a factor 2.5 compared to that of the second test case. Moreover, the computational times, which ranged from 37 to 64 seconds for the three test cases, indicate that through the B-char method, characteristic-based schemes are not too costly to implement in 3D.

These tests demonstrate that the B-char approach is a cost-effective and accurate way of implementing characteristic-based schemes such as the ELLAM, in both 2D and 3D.

Remark 4.2** (Number of balls used for B-char).**

Upon increasing the number of balls used for the B-char method, the initial approximation of the volumes becomes better. For the tests presented above, this allows us to achieve results with a similar quality even if we take $N<10$ iterations of (28) before solving the optimisation problem (29). The choice here is between two options: increasing the number of balls and ball intersections, or increasing the number of iterations of (28); we found that, for the tests presented above, the choice of 4 balls (in 2D) or 8 balls (in 3D) was producing a good accuracy/cost ratio.

5. Conclusion and possible outlooks

In this paper, we have developed the B-char method, which is a cheap and efficient way to implement characteristic-based schemes, whilst preserving the important properties of local and global mass conservation, applicable in both 2D and 3D. Numerical tests were provided, and showed that the quality of the solutions obtained from B-char ELLAM is similar to the quality obtained by tracking and computing intersections of polygonal cells. The computational cost of the B-char method, however, is considerably lower than the computational cost of implementing the polygonal ELLAM. The tests run on benchmarks from the literature also showed that, in most situations, the B-char ELLAM produces rather accurate solutions, comparable to those obtained via other numerical schemes and much larger time steps.

Extension of the B-char method onto generic meshes can also be achieved, provided that we can design an algorithm to pack the balls inside the cells. Although the numerical tests were only performed on ELLAM type schemes, it would be interesting to explore the extension of the B-char method onto other characteristic-based schemes which are globally mass conserving, such as the MMOCAA. Future work would involve extending the B-char method onto non divergence-free velocity fields, and also those that are approximated by numerical schemes, e.g. piecewise polynomial velocity fields. The problem here lies with the fact that the divergence of the velocity field may be different from cell to cell, which makes it difficult to obtain a good approximation for the measure of the trace-back balls $|\widehat{B}_{K,s}|$. This requires further development and study of the idea proposed for approximating $\widehat{r}_{K,s}$ in (22) for non-constant $\phi$, so that (25) gives a good enough initial approximation for the volume of the intersecting region $|F_{-{\delta\!t}^{(n+\frac{1}{2})}}(K)\cap M|_{\phi}|$.