{'fraction_non_alphanumeric': 0.05251032729583731, 'fraction_numerical': 0.018390530664124563, 'mean_word_length': 4.038423053832299, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An unconditionally energy stable and positive upwind DG scheme for the Keller-Segel model

The well-suited discretization of the Keller-Segel equations for chemotaxis has become a very challenging problem due to the convective nature inherent to them. This paper aims to introduce a new upwind, mass-conservative, positive and energy-dissipative discontinuous Galerkin scheme for the Keller-Segel model. This approach is based on the gradient-flow structure of the equations. In addition, we show some numerical experiments in accordance with the aforementioned properties of the discretization. The numerical results obtained emphasize the really good behaviour of the approximation in the case of chemotactic collapse, where very steep gradients appear.


Introduction
Since the introduction in the 70's of the biological Keller-Segel model for chemotaxis phenomena [27,28], it and many related variants have attracted a great deal of interest in the mathematical community. Chemotaxis, a biological process through which organisms (e.g. cells) migrate in response to a chemical stimulus, is modelled by means of nonlinear systems of partial differential equations (PDE). The classical one can be written as follows: find two real valued functions, u = u(x, t) and v = v(x, t), defined in Ω × [0, T ] such that: in Ω × (0, T ), (1a) ∂ n u := ∇u · n = 0, ∂ n v = 0, on ∂Ω × (0, T ), (1c) Herein Ω is a bounded and smooth domain of R d , with d ∈ N the spatial dimension, and the parameters are k i > 0 for i ∈ {0, 1, 2, 3, 4}. The mathematical formulation (1) can be interpreted in biological terms as follows: u denotes a certain cell distribution (or population of organisms, in general) at the position x ∈ Ω and time t ∈ [0, T ], whereas v stands for the concentration of chemoattractant (i.e. a chemical signal towards which cells are induced to migrate). Both cells and chemoatractanct experiment some diffusion in the spatial domain. This auto-diffusion phenomena (experimented by cells and chemoattractant) are designed by the terms −k 0 ∆u and −k 2 ∆v, while the migration mechanism is modeled by the nonlinear crossdiffusion term −k 1 ∇ · (u∇v). This term is the major difficulty for theoretical analysis and also for numerical modelling of system (1). Further, the degradation and production of chemoattractant are associated with the terms −k 3 v and k 4 u, respectively. Note that the production of chemoattractant by the cells, to which cells are attracted, may eventually result in a chemotactic collapse, a phenomenon in which uncontrolled aggregation for u give rise to blowing up or exploding in finite time. This feature is well known and constitutes one of the outstanding characteristics of classical Keller-Segel model, and also one of its main challenges, specifically for numerical methods. Finally, the coefficient τ ∈ {0, 1} is considered to write at the same time the parabolic system when τ = 1, or the parabolic-elliptic for τ = 0.
Regarding the mathematical analysis for system (1): some results on sufficient conditions to ensure global existence and boundedness of solutions along time can be shown (see e.g. the review of Bellomo et al [5] and the references therein). They are based on mass conservation for u and on an energy dissipation law for this model (see section 3). For dimension d ≥ 2, these results require the initial density of cells, Ω u 0 , to be bigger than certain threshold. On the other hand, considerable research has been done in the direction of finding cases where chemotactic collapse arise. Among them, it is worth mentioning the first result in this direction, due to Herrero and Velázquez [23], where radially symmetric two-dimensional solutions which finite-time blow up are found. Other authors shed light on more general cases, for instance Horstmann and Wang [24] (non symmetric blow-up solutions) or Winkler [36] (higher dimensional case).
In the last decades, a lot of papers have been published dealing with these kind of theoretical issues both for the classical model (1) and for other models based on some extensions or generalizations. In general, they start from the Keller-Segel classical equations and modify them with the purpose of avoiding the non-physical blow up of solutions, producing solutions which are closer to the "real chemotaxis" phenomena observed in biology. Models include logistic, non-linear diffusion or production terms, chemo-repulsion effects or coupling with fluid equations [8,20,33,34,37]. See e.g. [3,5] for more examples. Thus, it is hoped that understanding the classical Keller-Segel equations may open new insights for dealing in depth with those other chemotaxis models.
On the other hand, taking into account the considerable efforts of the mathematical community in the theoretical analysis of Keller-Segel models, the number of papers dedicated to numerical analysis and simulation of chemotaxis equations is much lower, in relative terms. The main difficulty is the numerical approximation of the cross-diffusion term. Not only because of its non-linearity, but due to its convective nature, which makes particularly difficult to deal with using the finite element (FE) method (see, for instance, [6,17] for more details about this method). This difficulty is specifically significant in steep-gradient regions for v, which are precisely relevant in blow-up settings. Furthermore, preserving the physical properties of the continuous model (mass conservation, positivity and energy dissipation) in the discrete case adds an extra level of complication when it comes to designing a well-suited approximation.
Despite that, many interesting works have been published on numerical simulation of chemotaxis equations using different kind of approaches. For instance, the work by Saito [31] uses FE with upwind stabilization for the parabolic-elliptic Keller-Segel model (τ = 0 in (1)) showing mass conservation, positivity and error estimates. Also, Gutiérrez-Santacreu and Rodríguez-Galván demonstrate positivity, an energy law and a priori bounds for their FE scheme on acute meshes in [22]. In addition, other sophisticated techniques have been applied to this problem. This is the case of the finite volume (FV) method (we recommend [29] on this topic) which has become a very popular and successful approach as we can observe in papers like [10] by Chertock and Kurganov in which the authors devised positive preserving methods for the parabolic-parabolic formulation (τ = 1 in (1)), with demonstrated high accuracy and robustness, specifically on chemotactic collapse. It might also be pointed out the works of Saad and others, for instance in [26], where a volume finite element scheme for the capture of spatial patterns for a volume-filling chemotaxis is analyzed.
In this sense, it is also worth mentioning the very recent works [4,9,25,32] that show and analyze different techniques to approximate the solution of (1) while achieving mass conservation, positivity and energy-dissipation for a strictly positive initial cell condition. In the paper by Badía et al. [4] a discrete scheme using stabilized FE with a graph-Laplacian operator and a shock detector is proposed. This discretization satisfies, in general, both the mass conservation and positivity properties, and, in the case of acute meshes, it is also energy-dissipative. Otherwise, Huang and Shen develop in [25] a time-discrete approximation, admitting any spatial discretization, that preserves the positivity for the cell distribution u and is energy stable for a modified energy. This latest approach uses a suitable transformation of the solution for the positivity and the scalar auxiliary variable (SAV) technique for the energy stability. Alternatively, Shen and Xu introduced in [32] another general, positive (for the cell distribution u) and energy-dissipative, time-discrete scheme based on the gradient flow structure of the continuous model that admits different kind of spatial discretization such as FE, spectral methods or even finite-differences. The order of the latest approach and the blow-up phenomenon of the discrete solution, under a CFL condition, is studied in [9] by Chen et al.
Furthermore, discontinuous Galerkin (DG) methods (we refer the reader to [12,14] for a further insight) aroused the interest of researchers in recent years due to their flexibility for approximating, using standard meshes and computer libraries, different types of PDEs: elliptic, parabolic, hyperbolic. In the chemotaxis context, it is worth mentioning the paper of Y. Epshteyn and A. Kurganov [16], where the FV scheme given in [10] for the 2D Keller-Segel model (1) is extended to a DG scheme on cartesian meshes, obtaining good approximations even on blow-up regimes. Y. Epshteyn introduced two other related schemes in [15]. In all cases, different discontinuous Interior Penalty (IP) methods and upwinding techniques were considered for defining DG approximations. The schemes are even applied to the simulation of an haptotaxis model of tumor invasion into healthy tissue. Error estimates are shown but no energy property or maximum principle for the schemes is proven. In fact, spurious oscillations and negative values in the solution are reported.
More recent works make further progress in this direction, for instance in [38], where classical equations (1) are approximated by a positivity-preserving DG method with strong stability preserving (SSP) high order time discretizations. Error order estimates as well as positivity are shown in this work. Finally, in [21,30], the local discontinuous Galerkin method is applied, showing respectively positivity and energy dissipation.
In this work, we propose a new upwind DG scheme for the Keller-Segel model (1) that preserves the mass-conservation, positivity and energy stability properties of the continuous problem. As in [32], the proposed discretization takes advantage of the gradient flow structure of the model. First, Section 2 sets the notation that we are going to consider throughout the paper. In Section 3 we discuss the physical properties of the continuous model. Section 4 is the main part of the paper, in which we introduce the upwind DG scheme (8). In particular, in Subsection 4.1, we define the upwind approximation based on the ideas introduced in [1] along with some geometrical considerations for the mesh family T h , and we discuss the properties of the scheme in Subsection 4.2. Finally, in Section 5, we show several numerical tests in which we reproduce some blow-up results shown in the literature with one steady peak, [10], a peak moving towards the corner of the domain, [31], and several peaks, [2,7,18,35]. These numerical experiments endorse the good behaviour of the approximation obtained with the new scheme, which allows us to capture peaks reaching values up to the order of 10 7 . These kind of numerical results are rare in the literature due to the steep-gradients inherent to such sort of tests.

Notation
In this section we introduce the notation used throughout the paper. First, we consider a shaperegular triangular mesh T h = {K} K∈T h of size h over a bounded polygonal domain Ω ⊂ R d . Moreover, we note the set of edges or faces of T h by E h , which can be splitted into the interior edges E i h and the boundary edges h . Now, we fix the following orientation for the unit normal vector n e associated to and edge e ∈ E h of the mesh T h : • If e ∈ E i h is shared by the elements K and L, i.e. e = ∂K ∩ ∂L, then n e is exterior to K pointing to L (see Figure 1).
h , then n e points outwards of the domain Ω.
L K e n e Moreover, we denote the baricenter of the triangle K ∈ T h by C K . Now, we define the approximation spaces of discontinuous, P disc k (T h ), and continuous, P cont k (T h ), finite element functions over T h whose restriction to K ∈ T h are polynomials of degree k ≥ 0. In addition, the average { {·} } and the jump [[·]] of a scalar function v on an edge e ∈ E h are defined as follows: Regarding the time discretization, we take an equispaced partition Finally we set the following notation for the positive and negative parts of a scalar function v:
To the best knowledge of the authors, the existence of global solutions in time is still not clear in the literature.
and the initial conditions By taking the chemical potential of u, we can rewrite (2) as a gradient flux in which the flux direction is given by −∇µ which contains the effect of both the diffusion and the chemotaxis terms. This variational formulation consists of and the initial conditions Remark 3.2. By taking u = 1 in (2a) (or (4a)), any solution u conserves the mass, because (4), and adding the resulting expressions, one has that any solution (u, v) satisfies the following energy law where E : and H 1 (Ω) + = {u ∈ H 1 (Ω) : u ≥ 0}.

Fully discrete scheme
First, we regularize the chemical potential of u, defined in (3), by for some ε > 0. Note that ε is a regularization parameter since log(u + ε) is regular for all u ≥ 0. We will take ε = ε(h, ∆t) such that ε(h, ∆t) → 0 if (h, ∆t) → 0, being h > 0 the mesh size and ∆t > 0 the time step. Then, we propose the following decoupled fully discrete first order in time and upwind DG in space scheme for the model (1): Let v m ∈ P cont 1 (T h ) and u m ∈ P disc 0 (T h ) be given.
Step 1: for all v ∈ P cont 1 (T h ).
In order to preserve the positivity, we have introduced a truncation of the function u taking its positive part u ⊕ in the upwind part of (8b), which is consistent as the solution of the continuous model (1) satisfies u ≥ 0. Moreover, we take mass lumping in the terms ∂ t v m+1 , v h and k 3 v m+1 , v h in (8a).

Definition of upwind bilinear form a upw
h (·; ·, ·) Now, we are going to define the upwind bilinear form a upw h (·; ·, ·), introduced in the scheme (8). In order to achieve the energy stability with the scheme (8), we must consider the following hypothesis that will let us approximate the flux −∇µ accordingly: Hypothesis 1. The mesh T h of Ω is structured in the sense that the line between the baricenters of the triangles K and L is orthogonal to the interface e = K ∩ L ∈ E i h .
Then, we define the following upwind bilinear form to be applied to the flux −∇µ which may be discontinuous over E i h : where, for every e ∈ E i h with e = K ∩ L, with Π 0 being the projection on P disc 0 (T h ) and D e (T h ) the distance between the baricenters of the triangles K and L of the mesh T h that share e ∈ E i h , denoted by C K and C L , respectively. This way, we can rewrite (9) as Remark 4.2. Since the quadrature formula of the baricenter (or centroid) is exact for polynomials of order 1, if µ ∈ P disc 1 (T h ), we have that where C K is the baricenter of K ∈ T h . Hence, if µ ∈ P disc 1 (T h ), the expression (10) is the slope of the line between the points (C K , µ(C K )) and (C L , µ(C L )), which, under the Hypothesis 1, this line is parallel to the vector n e , with e = K ∩ L ∈ E i h . This expression is considered as an approximation of the discontinuous numerical normal flux ∇µ · n e for µ ∈ P disc 1 (T h ). In addition, observe that, since the baricenters are located 1/3 of the median from the side and 2/3 of the median from the vertex of the triangle, the expression (10) does not degenerate when h → 0, i.e., D e (T h ) > 0 for every e ∈ E i h and h > 0. A visual representation of the regular polygonal structure given by the lines between the adjacent baricenters is given in Figure 2.
Furthermore, in order to preserve the positivity of the variable v in our fully discrete scheme, we assume the following hypothesis.
Hypothesis 2. The mesh T h is acute, i.e., the angles of the triangles of T h are less than or equal to π/2. We give some examples, the meshes represented in Figure 3, which satisfy both Hypotheses 1 and 2. Proof. We will only prove the case a) Mesh 1 since the case b) Mesh 2 is anologous. Observe Mesh 1 in Figure 3. The baricenter of the triangles OAB, OBC and OCD are, respectively O+A+B . Hence, if we denote by e 1 the edge between OAB and OBC and by e 2 the edge between OBC and OCD, then Now, since l/|e 1 | = 1 and l/|e 2 | = 1/ √ 2, we can define D e (T h ) for any e ∈ E i h as  Proof. Just need to take u = 1 in (8b). Proof. Proving that if u m ≥ 0 and v m ≥ 0 (when τ > 0) then v m+1 ≥ 0 using that the mesh is acute is a classic result which can be found, for example, in [11,19]. Moreover, if u m ≥ 0 and v m ≥ 0 (when τ > 0) it follows from the equation (8b) that u m+1 ≥ 0 using the same arguments that are shown to prove the positivity result in [1] since the proof given in [1] is independent of the flux ∇µ m+1 ε . Proposition 4.6. There is a unique solution v m+1 of the linear equation (8a).

Properties of the scheme
Proof. Since we are dealing with a discrete linear problem, existence and unicity of the solution are equivalent. Hence, we just need to assume that there are two solutions of (8a), v 1 and v 2 , substract the expressions resulting of evaluating both solutions and test with v 1 − v 2 to prove unicity of the solution.
Proof. Consider the following well-known theorem: Theorem 4.8 (Leray-Schauder fixed point theorem). Let X be a Banach space and let T : X −→ X be a continuous and compact operator. If the set is bounded (with respect to α), then T has at least one fixed point.
Given u m ∈ P disc 0 (T h ) with u m ≥ 0 and the unique solution v m+1 of (8a), we define the map is the unique solution of the linear (and decoupled) problem: To check that T is well defined, it is straightforward to see that the solutions u of (13a) and µ ε of (13b) are unique, which involves their existence as P disc 0 (T h ) is a finite-dimensional space. Secondly, we will check that T is continuous. Let { u j } j∈N , { µ εj } j∈N ⊂ P disc 0 (T h ) be sequences such that lim j→∞ u j = u and lim j→∞ µ εj = µ ε . Taking into account that all norms are equivalent in P disc 0 (T h ) since it is a finite-dimensional space, the convergences u j → u and µ εj → µ ε are equivalent to the elementwise convergences ( u j ) K → u K and ( µ εj ) K → µ εK for every K ∈ T h (this may be seen, for instance, by using the norm · L ∞ (Ω) ). Taking limits when j → ∞ in (13) (with u . . = u j , µ ε . . = µ εj and (u, µ ε ) . . = T ( u j , µ εj )), using the notion of elementwise convergence and the fact that log(( u) ⊕ + ε) is continuous, we get that hence T is continuous. In addition, T is compact since P disc 0 (T h ) have finite dimension. Finally, let us prove that the set is bounded (independent of α). The case α = 0 is trivial so we will assume that α ∈ (0, 1]. If (u, µ ε ) ∈ B, then u ∈ P disc 0 (T h ) is the solution of 1 ∆t (u − αu m , u) = −α a upw h (µ ε ; u ⊕ , u), ∀u ∈ P disc 0 (T h ).
Now, testing (14) with u = 1, we get that and, as u m ≥ 0 and it can be proved that u ≥ 0 using the same arguments than in Theorem 4.5, we get that u L 1 (Ω) ≤ u m L 1 (Ω) .
Moreover, since u ≥ 0, µ ε ∈ P disc 0 (T h ) is the solution of the equation Hence, Thus, taking into account that u is bounded in P disc 0 (T h ), we conclude that µ ε is bounded in P disc 0 (T h ). Since P disc 0 (T h ) is a finite-dimensional space where all the norms are equivalent, we have proved that B is bounded.
where (8) and consider the equalities Then, adding the resulting expressions for (8a) and (8c) and substracting (8b), we obtain Now, using that δ t (u m+1 )F (u m+1 ) ≥ δ t (F (u m+1 )) for F (u m+1 ) = log(u m+1 + ε) (owing to the fact that F (u) is convex) we have that Hence, using Proposition 4.4, Thus, taking into account (19) and (20), we obtain the discrete energy law (4.10). In consequence, the scheme (8) is unconditionally energy stable since Proof. To show δ t E(u m+1 , v m+1 ) ≤ 0 it suffices to prove that a upw h (µ m+1 ε ; (u m+1 ) ⊕ , µ m+1 ε ) ≥ 0 since we know that the discrete energy satisfies (17). Now, take v = µ m+1 ε and use the definition (11) of the upwind bilinear form to get the following: Remark 4.12. Notice that the energy stability of the scheme (8) is obtained thanks to the approximation of the flux −∇µ ε made in the upwind bilinear form a upw h (·; ·, ·). In addition, the approximation −∇ 0 ne µ ε on the edges e ∈ E i h requires the assumption of the Hypothesis 1 for the mesh T h as discussed in the Remark 4.2.
In the test 5.1 we reproduce the first numerical experiment shown in the paper [10] by A. Chertock and A. Kurganov. In this paper, they use a scheme that preserves the positivity of both variables u and v, although they do not show any energy related result. Hence, in our case, we can improve the results shown in the aforementioned paper assuring that our scheme preserves the energy law of the continuous Keller-Segel model.
Then, in the test 5.2 we simulate the qualitative behaviour of the solution in the numerical experiment made by N. Saito in [31]. However, since the initial conditions used for the experiment in [31] are not specified we cannot reproduce the exact same test shown in this paper. In the case of our numerical test, the qualitative behaviour of the solution is similar to the one in [31] until the mesh is refined enough so that we capture the blow-up phenomenon.
Finally, in the test 5.3 we reproduce the qualitative behaviour of the results in [2,7,35] for different variations of chemotaxis equations and in [18] for the Keller-Segel equations, where blowup phenomena with multiple peaks are shown.
Remark 5.1. The scheme (8) preserves the positivity and conserves the mass of u m+1 which implies u m+1 ∈ L 1 (Ω). Therefore, we cannot expect an actual blow-up in the discrete case as it occurs in the continuous model. However, we observe in the numerical tests how the mass accumulates in some elements of T h leading to the formation of peaks.
In fact, we are able to capture peaks that reach values up to the order of 10 7 using the approximation shown in (8). These kind of numerical results are not usual in the literature due to the difficulties when approximating the steep gradients that this process involves.
The numerical results shown in this paper have been obtained using the Python library FEniCS. In order to improve the efficiency of the code, these have been run in parallel using several CPUs.
For the sake of a better visualization of the results, a P cont 1 -projection of u is represented in 3D using Paraview.

One bulge of cells
First, we reproduce the results shown in [10]. On this purpose, we consider the initial conditions which are plotted in Figure 4.
As stated in [10], the u and v components of the solution are expected to blow up in a finite time due to the initial conditions chosen. The result of the test using the scheme (8) with h ≈ 1.41 · 10 −3 and ∆t = 10 −6 is shown in Figure 5. In fact, we observe a blow-up phenomenon for a certain finite time in the range conjectured by A. Chertock and A. Kurganov in [10], t * ∈ (4.4 · 10 −5 , 10 −4 ), as our discrete approximation reaches values of order 10 6 in this time interval.
Moreover, the positivity is preserved for both u and v as stated in Theorem 4.5 and, unlike the scheme presented in [10], we are certain that the discrete energy decreases (both for E(·, ·) and E ε (·, ·)) using the scheme (8)

Three bulges of cells
Now, we show the results for a similar test to the one that appears in [31]. In this case, we take the parabolic-elliptic case (τ = 0) so that the characteristic speed of v is much faster than the characteristic speed of u. Moreover, we consider the initial condition which is plotted in Figure 8.
In Figure 9 we can observe the result of the test with h ≈ 2.83 · 10 −3 and ∆t = 10 −5 . In this case, the qualitative behavior of the solution is similar to the one shown in [31], with the peak of cells moving towards a corner of the domain. However, the qualitative behavior of the solution is different if we take h ≈ 7.07 · 10 −3 and ∆t = 10 −5 . Now, as represented in Figure 10, a blow-up phenomenon seems to occur in finite time and the peak of cells remains motionless far away from the corners of the domain.
In both cases, the positivity is preserved and the energy decreases in the discrete case. See Figures 11 and 13 (left) for the case h ≈ 2.83 · 10 −3 and Figures 12 and 13 (right) for the case h ≈ 7.07 · 10 −3 .

Blow-up phenomenon with multiple peaks
Finally, we show the results for a test in which it seems to happen a blow-up phenomenon with multiple peaks as it occurs, for instance, in the cases that appear in [2,7,35] for different variations of chemotaxis equations and in [18] for the Keller-Segel equations. On this purpose, we consider the initial conditions u0 = 1000(cos(2πx) cos(2πy) + 1), v0 = 500(sin(3πx) sin(3πy) + 1), which are plotted in Figure 14.
In Figure 15 we can observe the result of the test with h ≈ 3.54 · 10 −3 and ∆t = 10 −7 . Notice that we obtain 8 peaks of cells that reach very high values (up to values of order 10 7 ). Therefore, we can infer that we obtain a blow-up phenomenon at a certain finite time t * close to 10 −4 with multiple peaks.
Again, the positivity is preserved as shown in Figure 16 and the energy decreases in the discrete case as in Figure 17.