An entropy stable high-order discontinuous Galerkin method for cross-diffusion gradient flow systems

As an extension of our previous work in Sun et.al (2018) [41], we develop a discontinuous Galerkin method for solving cross-diffusion systems with a formal gradient flow structure. These systems are associated with non-increasing entropy functionals. For a class of problems, the positivity (non-negativity) of solutions is also expected, which is implied by the physical model and is crucial to the entropy structure. The semi-discrete numerical scheme we propose is entropy stable. Furthermore, the scheme is also compatible with the positivity-preserving procedure in Zhang (2017) [42] in many scenarios. Hence the resulting fully discrete scheme is able to produce non-negative solutions. The method can be applied to both one-dimensional problems and two-dimensional problems on Cartesian meshes. Numerical examples are given to examine the performance of the method.


Introduction
Cross-diffusion systems are widely used to model multi-species interactions in various applications, such as population dynamics in biological systems [39], chemotactic cell migration [32], spread of surfactant on the membrane [26], interacting particle systems with volume exclusion [5,2], and pedestrian dynamics [24]. In many situations, the systems are associated with neither symmetric nor positive-definite diffusion matrices, which not only complicate the mathematical analysis but also hinder the development of numerical methods. Recently, progress has been made in analyzing a class of cross-diffusion systems having the structure of a gradient flow associated with a dissipative entropy (or free energy) functional, see [28,30,31] and references therein. We will exploit this gradient flow structure to develop a high-order stable discontinuous Galerkin (DG) method for these cross-diffusion systems.
We also assume F to be an m × m positive-semidefinite ρ-dependent matrix, in the sense that z · F z ≥ 0. Note F can be non-symmetric. Equation (1.2) The system can be rewritten as ∂ t ρ = ∇ · F (ρ)∇ δE δρ . One can see at least for classical solutions that 3) with the usual notation for the matrix product of square matrices A : B := i,j a ij b ij . The Liapunov functional (1.2) due to (1.3) indicates certain well-posedness of the initial value problem (1.1). However, in applications with ρ representing non-negative physics quantities, such as species densities, mass fractions and water heights, the well-posedness usually relies heavily on the positivity of the solution. For example, with F (ρ) = diag(ρ), F (ρ) is semi-positive definite only when ρ is non-negative. Violating the non-negativity may result in a non-decreasing entropy in (1.3). Furthermore, in problems that a logarithm entropy E = Ω m l=1 ρ l (log ρ l − 1)dx is considered, the entropy may not even be well-defined when ρ admits negative values.
The entropy structure is crucial to understand various theoretical properties of crossdiffusion systems, such as existence, regularity and long time asymptotics of weak solutions, see [28,30,31]. It is desirable to design numerical schemes that preserve the entropy decay in (1.3), and would also be positivity-preserving if the solution to the continuum equation is non-negative, due to the concern we have mentioned. Various efforts have been spent on numerical methods for scalar problems with a similar entropy structure, including the mixed finite element method [4], the finite volume method [3,6], the direct DG method [33,34], optimal mass transportation based methods [11,27,12,10], the particle method [9] and the blob method [21,7]. Some methods can also be generalized to systems, for example the Poisson-Nerst-Plack system [35] and system of interacting species with cross-diffusion [8].
In this paper, we extend the DG method in [41] for scalar gradient flows to cross-diffusion systems. The DG method is a class of finite element methods utilizing discontinuous piecewise polynomial spaces, which was originally designed for solving transport equations [38] and was then developed for nonlinear hyperbolic conservation laws [18,17,16,15,20]. The method has also been generalized for problems involving diffusion and higher order derivatives, for example, the local DG method [1,19], the ultra-weak DG method [14] and the direct DG method [36]. The method preserves local conservation, achieves high-order accuracy, is able to handle complex geometry and features with good h-p adaptivity and high parallel efficiency. As an effort to incorporate these potential advantages in gradient flow simulations, in [41], we adopted the technique from [13] to combine the local DG methods with Gauss-Lobatto quadrature rule, and designed a DG method for scalar gradient flows that is entropy stable on the semi-discrete level. This method also features weak positivity with Euler forward time stepping method. Hence after applying a scaling limiter and the strong-stability-preserving Runge-Kutta (SSP-RK) time discretization [23], the fully discrete scheme produces non-negative solutions. This positivity-preserving procedure is established in [43,42,40].
The main idea for handling cross-diffusion systems is to formally rewrite the problem into decoupled equations and apply our previous numerical strategy in [41] for scalar gradient flows to the unknown density vector component-wise. In fact, the system (1.1) can be rewritten as In the particular case of one dimension, we are reduced to apply our previous scheme in [41] to ∂ t ρ l = ∂ x (ρ l v l ), for l = 1, . . . , m. The numerical fluxes are properly chosen to ensure the entropy decay and the weak positivity. In this approach, the boundedness of diag(ρ) −1 F (ρ) is required for v to be well-defined. This assumption does hold in various applications, see examples in Section 4 and Section 5.
The rest of the paper is organized as follows. In Section 2 and Section 3, we propose our DG schemes for one-dimensional and two-dimensional gradient flow systems respectively.
Each section is composed of three parts: notations, semi-discrete scheme and entropy inequality, and fully discrete scheme as well as positivity-preserving techniques for producing non-negative solutions (with central and Lax-Friedrichs fluxes). In Section 4 and Section 5, we show several numerical examples to examine the performance of the numerical schemes.
Finally, concluding remarks are given in Section 6.
2 Numerical method: one-dimensional case

Notations
In this section, we consider the one-dimensional problem For simplicity, the compact support or periodic boundary condition is assumed, while it can also be extended to the zero-flux boundary condition.
) and I = ∪ N i=1 I i be a regular partition of the domain Ω. The length of the mesh cell is denoted by for some fixed constant c mesh . The numerical solutions are defined in the tensor product space of discontinuous piecewise polynomials for the jumps. Let {x r i } k+1 r=1 be the k + 1 Gauss-Lobatto quadrature points on I i and {w r i } k+1 r=1 be the corresponding weights associated with the normalized interval [−1, 1]. We introduce the following notations for the quadrature rule.
Here I is a component-wise interpolation operator. Namely, Iψ : R → R m , and the l-th component of Iψ is the k-th order interpolation polynomial of ψ l at Gauss-Lobatto points.
We also define

Semi-discrete scheme and entropy stability
We firstly rewrite (2.1) into a first-order system.
On each mesh cell I i , we multiply with a test function and integrate by parts to get The numerical scheme is obtain by taking trial and test functions from the finite element space, replacing cell interface values with numerical fluxes and applying the quadrature rule.
More precisely, we seek We have the following choices for F u and ξ.
(1) Central and Lax-Friedrichs fluxes: where and | · | ∞ is the l ∞ norm on R m .
(2) Alternating fluxes: which is formally reduced to the scalar case discussed in [41].
We define the discrete entropy As is stated in Theorem 2.1, the numerical scheme has a decaying entropy as that for the continuum system. E h is the associated discrete entropy defined in (2.5). Then (1) for central and Lax-Friedrichs fluxes, for alternating fluxes,

Proof. A direct computation yields
Here we have used the fact that I(F h u h ) · ∂ x I(ξ h ) is a polynomial of degree 2k − 1 and that the Gauss-Lobatto quadrature rule is exact. Then again with integrating by parts and the exactness of the quadrature, one can get where in the last identity we used the scheme (2.2b) with ψ h = F h u h . After summing over the index i and using the periodicity, we obtain The proof is completed by substituting numerical fluxes in (2.3) and (2.4). Note that This is referred to as the weak positivity. Then we can apply a scaling limiter, squashing the solution polynomials towards cell averages to avoid inadmissible negative values. It is shown in [42] that the limiter preserves the high-order spatial accuracy.
then the preliminary solution ρ n+1,pre h,l obtained through the Euler forward time stepping has non-negative cell averages.

Let
Then ρ n+1 h,l (x r i ) is non-negative. Furthermore, the interpolation polynomial ρ n+1 h,l (x) maintains spatial accuracy in the sense that where C k is a constant depending only on the polynomial degree k.
We now invoke Lemma 2.1 in [41] for the scalar case implying (ρ h,l ) n+1 i ≥ 0. The nonnegativity of the solution is ensured through the definition of θ l,i . For the accuracy result we refer to Theorem 4 in [42]. for most scenarios, since the scheme is defined only through these points. One can also ensure pointwise non-negativity of solution polynomials on the whole interval by taking

Remark 2.5 (Time steps).
As has been analyzed for the scalar case in Remark 2.2 in [41], we expect the time step to be τ ≤ µh 2 for some constant µ. In practice, we assume a diffusion number µ. If a negative cell average emerges, we halve the time step and redo the computation. Theorem 2.2 guarantees that the halving procedure will end after finite loops.
Remark 2.6 (The scaling parameter). For robustness of the algorithm, especially for dealing with log-type entropy functionals, we introduce a parameter ε l,i = min 10 −13 , (ρ h,l ) n+1,pre i and use θ ε l,i = min we scale the polynomial to require it takes values not smaller than 10 −13 at Gauss-Lobatto points; otherwise, we set the solution to be the constant (ρ h,l ) n+1,pre i in the interval. Note as long asρ h,l (x) ≥ 10 −13 , the accuracy could still be maintained.

Remark 2.7 (Other bounds).
In general, it would be difficult to preserve other bounds besides positivity though this procedure. Similar difficulty has also been encountered in [40].

High-order time discretization
We adopt SSP-RK methods for high-order time discretizations. Since the time step will be chosen as τ = µh 2 , the first-order Euler forward time stepping would be sufficiently accurate for k = 1 to achieve the second-order accuracy. For k = 2 and k = 3, we apply the second-order SSP-RK method For k = 4 and k = 5, a third-order time discretization method should be used As one can see, SSP-RK methods can be rewritten as convex combinations of Euler forward steps. Hence if the time step is chosen to be sufficiently small, and the scaling limiter is applied at each Euler forward stage, then the solution will remain non-negative after one full time step.
3 Numerical method: two-dimensional case

Notations
In this section, we generalize the previous ideas to two-dimensional systems of the form The problem is set on a rectangular domain Ω = I × J with the compact support or periodic boundary condition. We consider a regular Cartesian mesh on Ω, with Then the solution is sought in the following finite element space. .
Same notations will be used in the scalar case v h,l ∈ V h . Finally, for the quadrature rule, we denote by

Semi-discrete scheme and entropy inequality
The semi-discrete scheme is given as follows.
, y) dy As that in the one-dimensional case, two choices of numerical fluxes can be used.
(1) Central and Lax-Friedrichs flux (2) Alternating fluxes The numerical scheme mimics a similar entropy decay behavior as the continuum equation.
We define the discrete entropy We state the following entropy decay property for the semi-discrete scheme.
Theorem 3.1. Let ρ h , u x h and u y h be obtained from the semi-discrete scheme for twodimensional problems. E h defined in (3.1) is the discrete entropy.
(1) Suppose central and Lax-Friedrichs fluxes are used, then (2) Suppose alternating fluxes are used, then Proof. Following the blueprint of the one dimensional case, we get by a direct computation Then using the scheme, we deduce For each fixed y, each component of -th order polynomial with respect to y for each fixed x. Hence the Gauss-Lobatto quadrature with k + 1 node is exact. We then replace the quadrature rule with the exact integral, and integrate by parts to get , y) dy Again by changing back to the quadrature rule and applying the scheme, we obtain , y) dy Since we have assumed periodicity, all cell interface terms cancel out with alternating fluxes.
For central and Lax-Friedrichs fluxes, there remain additional penalty terms.
since F h is positive-semidefinite and e is convex as in the one dimensional case.

Fully discrete scheme and preservation of positivity
.
Then the preliminary solution (ρ h,l ) n+1,pre i,j obtained through the Euler forward time discretization has non-negative cell averages. Define ρ n+1 h,l so that

Let
Then ρ n+1 h,l (x r i , y s j ) is non-negative. Furthermore, ρ n+1 h,l would maintain the spatial accuracy achieved by ρ n+1,pre h,l .
Proof. Using the scheme we obtain Note that the numerical fluxes can be rewritten as Then we have which formally reduces to the scalar case. One can follow the proof of Theorem 3.1 in [41] to show thatρ n+1,pre h,l ≥ 0, if the prescribed time step restriction is satisfied. The non-negativity of ρ n+1 h,l (x r i , y s j ) can be justified using the definition of θ l,i,j . Remark 3.1. As before, pointwise non-negative solutions can be obtained by taking ρ min l,i,j = inf (x,y)∈I i ×J j ρ n+1,pre h,l (x, y). The time step τ = µh 2 will be used for time marching. As negative cell averages emerge, we halve the time step. We will also use θ ε l,i,j = min instead of θ l,i,j when applying the scaling limiter. Here ε l,i,j = min 10 −13 , (ρ h,l ) n+1,pre i,j .
We compute up to t = 0.002, with time step set as τ = 0.001h 2 . Positivity-preserving limiter is not activated. We use the Lax-Friedrichs flux and alternating fluxes respectively for computation. As one can see from Table 4.1 and Table 4.2, optimal order of convergence is achieved with the alternating fluxes. The order of accuracy seems to be degenerated when we use central and Lax-Friedrichs fluxes with odd-order polynomials. This may related with the fact that the problem is smooth hence the jump term α 2 [ρ h ] is too small, hence the Lax-Friedrichs flux F u is close to the central flux when the mesh is not well refined, which would cause order reduction. In Table 4.3, we document the error with different Lax-Friedrichs constants, with α replaced byα = 0, 2α, 10α. As one can see, the optimal order is retrieved asα becomes large. We also remark, choosing a larger Lax-Friedrichs constant does not affect the compatibility with positivity-preserving procedure, while it does make the time step more restrictive.  and Teramoto [39] for the second accuracy test. All the cross-diffusion and self-diffusion coefficients are taken as 1. The system is written as follows.
The numerical solution at the next mesh level is set as a reference to evaluate the error. The

System (4.2) comes from the tumor encapsulation model proposed by Jackson and Byrne
in [25], which describes the formation of a dense, fibrous connective tissue surrounding the benign neoplastic mass. In the system, ρ 1 corresponds to the concentration of tumors, and ρ 2 is the concentration of the surrounding tissue. In [29], the authors pointed out the system can be formulated as a gradient flow with the entropy Then ξ = (log ρ 1 1−ρ 1 −ρ 2 , log ρ 2 1−ρ 1 −ρ 2 ) T and the corresponding F is semipositive definite with the prescribed parameters. In the numerical test, we firstly consider the case β = 0.0075 and γ = 10. The initial condition is set as Example 4.4 (Surfactant spreading). We perform numerical simulations of a system modelling the surfactant spreading on a thin viscous film, which can be used for analyzing the delivery of aerosol for curing the respiratory distress syndrome. The system was derived by Jensen and Grotberg in [26] and was then analyzed by Escher et al. in [22]. In this system, ρ 1 represents the film thickness and ρ 2 is the concentration of the surfactant. The parameter g corresponds to a gravitational force. (4.3a) (4.3b) Following the derivation in [22], the system is associated with the following entropy functional Accordingly, ξ = (gρ 1 , log ρ 2 ) T , F = diag(ρ) The problem is again associated with the logarithm entropy.
k N x L 1 error order L 2 error order L ∞ error order 1 Again ρ 1 and ρ 2 correspond to the film thickness and surfactant concentration respectively.
(5.2a) (5.2b) Here z = b(x, y) gives the impermeable interface between the seawater and the bedrock.
The saltwater sits on the bedrock between z = b(x, y) and z = b(x, y) + ρ 2 (x, y, t). From z = b(x, y) + ρ 2 (x, y, t) to z = b(x, y) + ρ 2 (x, y, t) + ρ 1 (x, y, t) is the freshwater, which is immersible with the saltwater. The parameter µ ∈ (0, 1) is the mass density ratio between the freshwater and the saltwater. (5.2) is associated with the energy functional It can be show that there exists a non-negative solution to (5.2) satisfying the energy decay We assume the zero-flux boundary condition and use a 20 × 20 square mesh on [0, 1] × [0, 1].
We compute to t = 12 with τ = 0.002(h x ) 2 . The parameter is set as µ = 0.9. The initial converge to a steady state as t becomes large.

Concluding remarks
In this paper, we extend the DG method in [41] to solve cross-diffusion systems with a gradient flow structure. The difficulty is that the non-negativity of solutions is usually an essential part for the entropy stability. One needs to design numerical schemes to preserve      the entropy structure, as well as to ensure the positivity of solutions. We adopt the Gauss-Lobatto quadrature rule in the DG method, so that the resulting semi-discrete scheme are subject to an entropy inequality consistent with that of the continuum system. Furthermore, for a class of problems, with a suitable choice of numerical fluxes, the numerical method is compatible with the positivity-preserving procedure established in [43] and [42]. The extension to two-dimensional problems on Cartesian meshes is also discussed. In our numerical tests, we observe that the methods achieve high-order accuracy. Though the convergence rate is reduced with odd-order piecewise polynomials, the optimal rate can be retrieved by imposing larger Lax-Friedrichs constant in the numerical fluxes. Numerical simulations to problems with positivity-preserving issues have also been performed.
A Fully discrete entropy inequality Theorem A.1. Consider the Euler forward time discretization of the one-dimensional scheme (2.2) with central and Lax-Friedrichs fluxes (2.3) on a uniform mesh. Suppose z · F z ≥ β F |z| 2 , z · Dξz ≥ β Dξ |z| 2 , |F | ≤ β F and |Dξ| ≤ β Dξ uniformly in ρ for some fixed constants β F , β Dξ , β F and β Dξ . If τ ≤ min i ( Here |F | and |Dξ| are l 2 matrix norms of F and Dξ respectively, c 0 is a constant for norm equivalence and c inv is the constant in the inverse estimate.
Proof. We omit all subscripts h and superscripts n in this proof.
Here we have applied the mean value theorem and ζ is on the line segment between ρ n+1 and ρ. Let η = ρ n+1 −ρ τ . Then Using the scheme (2.2), we obtain where ζ i+ 1 2 lies between ρ − i+ 1 2 and ρ + i+ 1 2 . Our main task is to estimate Here c 1 and c 2 are constants to be determined. Note that ∼ I i | · | 2 dx defines a norm on P k (I i ). Using norm equivalence and inverse estimates, we have where c inv is a constant not less than 1 in the inverse estimates. Taking c 1 = c 2 0 c inv h 2 and c 2 = c 0 c inv h , then we have Therefore, we conclude that The proof is completed by substituting the time step restriction into the inequality.