Bound-preserving finite element approximations of the Keller-Segel equations

March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations Santiago Badia School of Mathematics, Monash University, Clayton, Victoria 3800, Australia CIMNE, Centre Internacional de M`etodes, Num`erics a l’Enginyeria, Barcelona 08034, Spain santiago.b[email protected]du Jes´us Bonilla Los Alamos National Laboratory, Los Alamos, NM 87545, USA jbonil[email protected] Juan Vicente Guti´errez-Santacreu∗ Departamento de Matem´atica Aplicada I, E. T. S. I. Inform´atica, Universidad de Sevilla, Avda. Reina Mercedes, s/n. E-41012 Sevilla, Spain [email protected] Received 22 July 2022 Revised 19 September 2022 Accepted 9 January 2023 Published 13 March 2023 Communicated by F. Brezzi and N. Bellomo This paper aims to develop numerical approximations of the Keller–Segel equations that mimic at the discrete level the lower bounds and the energy law of the continuous problem. We solve these equations for two unknowns: the organism (or cell) density, which is a positive variable, and the chemoattractant density, which is a non-negative variable. We propose two algorithms, which combine a stabilized finite element method and a semi-implicit time integration. The stabilization consists of a nonlinear artificial diffusion that employs a graph-Laplacian operator and a shock detector that localizes local extrema. As a result, both algorithms turn out to be nonlinear and can generate cell and chemoattractant numerical densities fulfilling lower bounds. However, the first algorithm requires a suitable constraint between the space and time discrete parameters, whereas the second one does not. We design the latter to attain a discrete energy law on acute meshes. We report some numerical experiments to validate the theoretical results ∗Corresponding author 609 March 17, 2023 10:25 WSPC/103-M3AS 2350014 610 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu on blowup and nonblowup phenomena. In the blowup setting, we identify a locking phenomenon that relates the L∞(Ω)-norm to the L1(Ω)-norm limiting the growth of the singularity when supported on a macroelement. Keywords: Keller–Segel equations; nonlinear parabolic equations; stabilized finiteelement approximation; shock detector; lower bounds; energy law. AMS Subject Classification 2020: 35K20, 35K55, 65N30 1. Introduction 1.1. The Keller–Segel equations Dictyostelium is an ameba, a unicellular eukaryotic organism that lives in the soil, feeds on bacteria, and reproduces by bipartition. However, what truly sticks out about its life cycle is its behavior as a social ameba. When amebas lack food and therefore cannot divide, Dictyostelium opts for an alternate life cycle and joins with its congeners to lead to a stage of cell development and differentiation. Thus, Dictyostelium amebas are capable of movement, thanks to a molecular mechanism, in response to certain chemicals released by themselves. This aggregation phenomenon is known as chemotaxis. Keller and Segel17, 18 were the pioneers in deriving the first mathematical model to predict aggregation phenomena in populations of Dictyostelium discoideum or Escherichia coli. The model they proposed reads as follows. Let Ω ⊂Rd,d= 2 or 3, be a bounded domain, with nbeing its outward-directed unit normal vector to Ω, and let [0, T] be a time interval. Take Q= (0, T]×Ω and Σ = (0, T]×∂Ω. Then, the boundary-value problem for the Keller–Segel equations consists of finding u:¯ Q→ (0,∞), the organism (or cell) density, and v:¯ Q→[0,∞), the chemoattractant density, satisfying (∂tu−∆u=−∇ · (u∇v) in Q, ∂tv−∆v=u−vin Q, (1.1) subject to the initial conditions u(0) = u0and v(0) = v0in Ω,(1.2) and the no-flux boundary conditions ∇u·n= 0 and ∇v·n= 0 on Σ.(1.3) The dynamics of solutions to (1.1)–(1.3) is governed by an energy law, which can be formally deduced as follows. Suppose we have a classical solution (u, v) to system (1.1)–(1.3) on Q; multiply (1.1)1by log u−vand (1.1)2by ∂tuand integrate over Qand, finally, add the two resulting equations to get d dtE(u(t), v(t)) = −D(u(t), v(t)),∀t∈[0, T ],(1.4) March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations 611 where E(u, v). =1 2ZΩ |∇v(x)|2dx+1 2ZΩ v2(x)dx−ZΩ u(x)v(x)dx +ZΩ u(x)log u(x)dx, and D(u, v). =ZΩ |∆v(x)−v(x) + u(x)|2dx+ZΩ ∇u(x) pu(x)−pu(x)∇v(x) 2 dx. We note that the energy law (1.4) is only well-posed on the condition that u > 0. The importance of (1.4) is twofold. On the one hand, existence theory20 of twodimensional (2D) bounded solutions draws heavily on (1.4) and a Moser–Trudingerlike inequality, together with the condition RΩu0(x)dx∈(0,4π); if Ω ⊂R2is a ball and (u0, v0) are radially symmetric, one needs RΩu0(x)dx∈(0,8π). In dimension three, smallness conditions for (u0, v0) are required to prove bounded solutions.26 On the other hand, solutions, which may blowup either in finite or infinite time,16, 26 are also present in system (1.1)–(1.3). These unbounded solutions emerge from careful asymptotic analysis of (1.4), i.e. if Tmax ∈(0,∞] is a blowup time, then E(u(t), v(t)) → −∞ as t→Tmax.(1.5) Another essential fact13 is that there exist times when both the cell and chemoattractant densities should simultaneously blowup. In the particular case that Ω is a ball and (u0, v0) are radially symmetric, it is well known that there exist 2D solutions that blowup in finite time provided that RΩu0(x)dx>8πholds.12 A sufficient condition,27 which stems from (1.4) for initial data, generates three-dimensional solutions blowing up within finite time, but with RΩu0(x)dxof any size. To complete the bibliography of this paper on account of modeling and analytical results on the original Keller–Segel equations and their generalizations in cross-diffusion contexts, even with interaction in fluid environments, we refer the interested reader to Refs. 1, 3, 4, 14, and 15, where our numerical algorithms could have room for their possible applications. In the light of the above discussion, it is of great interest to construct numerical solutions that satisfy lower bounds — positivity for the organism density and nonnegativity for the chemoattractant density — and a discrete counterpart of (1.4). A direct numerical simulation is a simple approach to gaining insight into the issues of potential singularity formations for system (1.1)–(1.3). Nevertheless, this is a rudimentary approach because an artificial blowup might arise from purely numerical oscillations. After all, the numerical scheme does not guarantee lower bounds. On the contrary, adding numerical diffusion ensuring lower bounds might prevent the detection of a singularity. These are subtle issues for system (1.1)–(1.3) because March 17, 2023 10:25 WSPC/103-M3AS 2350014 612 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu specific solutions always have singular rapid-growth behavior and it is difficult to distinguish them from blowup formation. Several numerical discretizations of (1.1)–(1.3) have been proposed to satisfy lower bounds. Among them, we can find geometrical conditions11, 21 on meshes being weakly acute in two dimensions and acute in three dimensions, flux limiters,23–25 discontinuous Galerkin methods,19 and finite volume methods.7, 8 As for numerical solutions satisfying (1.4) at the discrete level, very few numerical algorithms are available. For example, a discontinuous Galerkin method was proposed,10 but without proving lower bounds. On the other hand, a standard finite element method11 has recently been proved to hold both lower bounds and a discrete energy law for acute meshes. This paper aims to construct finite element approximations satisfying lower bounds (as at the continuous level) without any requirement on the mesh in addition to quasi-uniformity. We shall provide two stabilized finite element methods2, 5 consisting of adding a diffusion artifact. The artificial diffusion combines a graph- Laplacian operator and a shock detector, which minimizes the amount of numerical diffusion introduced in the system. For the first algorithm, we can prove lower bounds for both unknowns. The second algorithm enjoys lower bounds and, on acute meshes, a discrete energy law. In doing so, some terms are lumped, and the algorithm uses a new discretization of the chemotaxis term. Nevertheless, lower bounds might be attained, even without lumping as for the first algorithm. Besides, the stabilization term added in the discrete chemoattractant equation to preserve lower bounds on general meshes impede obtaining a discrete energy inequality. On acute meshes, we can switch off this stabilization keeping the lower bounds and prove a discrete energy law. In any case, both algorithms are nonlinear. When using both algorithms in the context of approximating solutions that might blowup in finite time, we find a relation between the L1(Ω)- and L∞(Ω)- norms leading to a locking structure in the growth of Dirac-like potential singularities. This limitation appears when the singularity at hand is supported on a macroelement because the L∞(Ω)-norm is controlled by the L1(Ω)-norm, which remains globally bounded. 1.2. Notation Here Lp(Ω) and W1,p(Ω), p∈[1,∞], are the Lebesgue and Sobolev spaces endowed with the usual norm k · kLp(Ω) and k · kW1,p(Ω). In the particular case of p= 2, it is denoted as H1(Ω). The inner product of L2(Ω) is denoted by (·,·). Finally, by C we mean a constant which may vary with context but is always independent of the discrete parameters. 1.3. Layout Section 2 states the two numerical algorithms proposed in this work. We first give the features needed for constructing the finite element spaces and then announce March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations 613 our two algorithms, where we describe the stabilizing terms in detail. In particular, we motivate an expression of the shock detector for each algorithm. Moreover, the second algorithm incorporates a new discretization of the chemotaxis term. Then, in Sec. 3, we prove lower bounds for both algorithms and a discrete energy law for the second one. Finally, in Sec. 4, we end up with some numerical experiments to test our two algorithms, especially on lower bounds. 2. Construction of Numerical Schemes The section presents our two stabilized finite element methods jointly with the hypotheses and notation required for developing the stabilizing terms. 2.1. Hypotheses Let Ω be bounded domain in Rd, with d= 2 or 3. Its boundary is assumed to be polygonal or polyhedral. For Ω, consider {Kh}h>0to be a quasi-uniform family of conforming subdivisions of ¯ Ω into closed, convex subsets K, triangles or quadrilaterals (d= 2), and tetrahedra or hexahedra (d= 3), with hK:= diam(K) and h:= maxK∈KhhK, so that ¯ Ω = PK∈Kh. Moreover, let Nh={ai}I i=1 be the coordinates of the nodes of Kh. Associated to Khis the finite element space Xh={xh∈C0(¯ Ω) : xh|K∈P1or Q1∀K∈ Kh}, where P1is the set of linear polynomials on K, with Kbeing a triangle or tetrahedron, and Q1is the set of bilinear polynomials on K, with Kbeing a quadrilateral or hexahedron. Let {ϕai}I i=1 be the global shape functions for Xh, that is, ϕai(ai) = δij for i, j = 1, . . . , I, for which Ωai= supp ϕaidenotes the macroelement associated to each ϕaiand Nh(Ωai) is the set of the nodes belonging to Ωai. Furthermore, one defines the set of indices I(Ωai) = {j∈I:aj∈Ωai}. Let us introduce ih:C(¯ Ω) →Xhthe linear interpolation operator such that ihx(ai) = xh(ai) for i= 1, . . . , I. A discrete inner product is defined as (xh,¯xh)h=ZΩ ih(xh(x)¯xh(x))dx. The nodal values {xh(ai)}I i=1 are denoted {xi}I i=1. We now recall a well-known inverse estimate concerning Xh; see e.g. Lemma 4.5.3 in Ref. 6 for a proof. There exists Cinv >0, independent of h, such that kxhkW1,r (Ω) ≤Cinvh−1+dmin{1 r−1 p,0}kxhkLp(Ω),∀xh∈Xh.(2.1) Let us finally consider the following averaged interpolation operator defined as follows. Take, for each node ai∈ Nh, an element Kaisuch that ai∈Kai. Thus, we have Ihψ=X i∈I 1 |Kai|ZKai φ(x)dx!ϕai. March 17, 2023 10:25 WSPC/103-M3AS 2350014 614 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu We know9, 22 that there exists Csta >0, independent of h, such that kIhψkWs,p(Ω) ≤CstakψkWs,p(Ω),for s= 0,1 and 1 ≤p≤ ∞,(2.2) and Ihψ≥or >0,if ψ≥or >0.(2.3) 2.2. Finite element approximation It is assumed that (u0, v0)∈L1(Ω) ×H1(Ω) with u0>0 and v0≥0. Then, it is considered u0h=Ih(u0) and v0h=Ih(v0) satisfying (2.2) and (2.3); namely ku0hkL1(Ω) ≤Cstaku0kL1(Ω) and u0h>0,(2.4) and kv0hkH1(Ω) ≤Cstakv0kH1(Ω) and v0h≥0.(2.5) The discretization of problem (1.1) will be based on its variational formulation. As a time integration we use an implicit Euler time-stepping method lagging in time the value of vin (1.1)1in order to decouple the computation between the unknowns. Given N∈N, we let 0 = t0< t1<· · · < tN−1< tN=Tbe a uniform partitioning of [0, T] with time step k=T N. Moreover, we denote δtwn+1 h . =k−1(wn+1 h−wn h). 2.3. Algorithm 1 The first proposed algorithm reads as follows. Algorithm 1. Let u0 h=u0hand v0 h=v0h. Step (n+ 1): Known (un h, vn h)∈Xh×Xh, find (un+1 h, vn+1 h)∈Xh×Xh such that, for all xh∈Xh (δtun+1 h, xh)+(∇un+1 h,∇xh) −(un+1 h∇vn h,∇xh)+(Bu 1(un+1 h)un+1 h, xh)=0,(2.6) and (δtvn+1 h, xh)+(∇vn+1 h,∇xh) + (vn+1 h, xh)+(Bv 1(vn+1 h)vn+1 h, xh)=(un+1 h, xh).(2.7) The stabilizing terms above are defined as follows. For Ξ = uor v, we consider (BΞ 1(wh)uh, xh) = X i∈IX j∈I(Ωai) νΞ ij(wh, vh)ujxi`(i, j),(2.8) with νΞ ij(wh, vh) = max{αi(wh)fΞ ij, αj(wh)fΞ ji,0},for i6=j, March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations 615 and νΞ ii(wh) = X j∈I(Ωai)\{i} νΞ ij(wh), where fΞ ij is given by fu ij = (∇ϕaj,∇ϕai)−(ϕaj∇vh,∇ϕai) + k−1(ϕaj, ϕai), and fv ij = (∇ϕaj,∇ϕai)+(k−1+ 1)(ϕaj, ϕai), respectively. Additionally, `(i, j). = 2δij −1 is the graph-Laplacian operator, where δij is the Kronecker delta. The shock detector αi(wh) modulates the action of the stabilizing terms. In this work, we use the shock detector proposed in Ref. 2. Let ai∈ Nhand aj∈ Nh(Ωai)\{ai}. Define rij =aj−aias being the vector pointing from nodes aito aj, with ˆ rij . =rij |rij |being its normalized vector. Construct asym ij as the point at the intersection between the line that passes through aiand ajand ∂Ωai, that is, not aj(see Fig. 1). The set of all symmetric nodes with respect to node aiis represented with Nsym h(Ωai). Finally, define rsym ij . =asym ij −ai, and usym j . =uh(asym ij ). Then, one can define the gradient’s linear approximation of the jump and the mean at node aiin direction rij as J∇whKij . =wj−wi |rij|+wsym j−wi |rsym ij |,(2.9) {{|∇wh·ˆ rij|}}ij . =1 2 |wj−wi| |rij|+|wsym j−wi| |rsym ij |!.(2.10) Making use of these definitions, the shock detector at node ai∈ Nhfor whreads αi(wh). =       "|Pj∈I(Ωai)J∇whKij| Pj∈I(Ωai)2{{|∇wh·ˆ rij|}}ij #q ,if X j∈I(Ωai) {{|∇wh·ˆ rij|}}ij 6= 0, 0,otherwise, (2.11) for some q∈R+. Fig. 1. Representation of the symmetric node asym ij of ajconcerning aiin a macroelement of triangles (left) and right quadrilaterals (right). March 17, 2023 10:25 WSPC/103-M3AS 2350014 616 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu One key property of the shock detector is that it localizes local extremes. See Lemma 3.1 in Ref. 2 for a proof. Lemma 2.1. It follows that 0≤αi(wh)≤1and that αi(wh)=1for any extreme value at ai. 2.4. Algorithm 2 It is not obvious at all whether discrete solutions to Algorithm 1 defined by Eqs. (2.6) and (2.7) have the discrete counterpart of (1.4). Thus, our next task is designing a new algorithm that leads to discrete solutions fulfilling (1.4). The second algorithm is a variant of Algorithm 1 in which the chemotaxis and stabilizing terms are modified as follows. Let xh,˜xh,¯xh∈Xhand assume xh>0. Then, as Pi∈I∇ϕai= 0, one finds (xh∇˜xh,∇¯xh) = X k,j,i∈I xk˜xj¯xi(ϕak∇ϕaj,∇ϕai) =X k,j,i∈I xk(˜xj−˜xj)¯xi(ϕak∇ϕaj,∇ϕai) = 2 X k,j,i∈I xk(˜xj−˜xj)(¯xi−¯xj)(ϕak∇ϕaj,∇ϕai) =X k∈I i<j∈I xk(˜xj−˜xi)(¯xi−¯xj)(ϕak∇ϕaj,∇ϕai). From the simple observation that lim xj→xi xj−xi log xj−log xi =xi, is an approximation of the identity, we now approximate xkby xij where xij =     xj−xi log xj−log xi ,if xj6=xi, xi,if xj=xi. (2.12) Observe as well that ak,aj,aimust belong to Ωak∩Ωaj∩Ωaiso that (ϕak∇ϕaj,∇ϕai)6= 0; therefore, as h→0, one has ai,aj→ak, and hence xi:= x(ai)→xk:= x(ak). This way xij →xkas h→0. Thus (xh∇˜xh,∇¯xh)≈X i<j∈I xji(˜xj−˜xi)(¯xi−¯xj)(∇ϕaj,∇ϕai). Since we do not know a priori whether the discrete solution will be positive, we need to construct an auxiliary function that somehow extends the logarithmic March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations 617 function to a nonpositive value so that the coefficient (2.12) makes sense. Let ε > 0 and gε(s) =    slog s−s, if s > ε, s2−ε2 2ε+ (log ε−1)s, if s≤ε, and hence g0 ε(s) =    log s, if s > ε, s ε+ log ε−1,if s≤ε. Let us thus define (xh∇˜xh,∇¯xh)∗=X i<j∈I γji(xh)(˜xj−˜xi)(¯xi−¯xj)(∇ϕaj,∇ϕai),(2.13) with γji(xh) =      xj−xi g0 ε(xj)−g0 ε(xi),if xj6=xi, [xi]+,if xj=xi, (2.14) where the operator [x]+= max{0, x}stands for the positivity part. It should be noted that g0 εis bijective so g0 ε(xj) = g0 ε(xi) implies xj=xi. Algorithm 2 reads as follows. Algorithm 2. Let u0 h=u0hand v0 h=v0h. Step (n+ 1): Known (un h, vn h)∈Xh×Xh, find (un+1 h, vn+1 h)∈Xh×Xhsuch that, for all xh∈Xh (δtun+1 h, xh)h+ (∇un+1 h,∇xh)−(un+1 h∇vn h,∇xh)∗ + (Bu 2(un+1 h, vn h)un+1 h, xh) = 0,(2.15) and (δtvn+1 h, xh)h+ (∇vn+1 h,∇xh)+(vn+1 h, xh)h +γ(Bv 2(vn+1 h)vn+1 h, xh)=(un+1 h, xh)h,(2.16) with γ∈ {0,1}. Here the stabilizing term in (2.15) is given by (Bu 2(wh, vh)uh, xh) = X i<j∈I νu ji(wh, vh)(uj−ui)(xj−xi),(2.17) where νu ji(wh, vh) = max{¯αi(wh)fij,¯αj(wh)fji,0},for i6=j, March 17, 2023 10:25 WSPC/103-M3AS 2350014 624 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu satisfies (a) Lower bounds un h>0and vn h≥0. (b) L1(Ω)-bounds kun hkL1(Ω) =ku0hkL1(Ω), and kvn hkL1(Ω) ≤ kv0hkL1(Ω) +ku0hkL1(Ω). for all n= 0, . . . , N. Proof. The induction argument requires to verify that (3.1) holds for n= 0. Then, according to Lemma 3.1, the lower bounds (3.3) are satisfied for n= 0 from (2.4) and (2.5). In Lemma 3.2, we see that the L1(Ω)-bounds (3.8) and (3.9) hold for n= 0. Finally, Lemma 3.3 provides the L2(Ω)-bound (3.13) for n= 0. The general case follows from sequentially applying Lemmas 3.1–3.3, respectively. 3.2. Algorithm 2 We now address the question of lower bounds for Algorithm 2. Lemma 3.4. (Lower bounds) Let γ= 1. Assume that un h>0and vn h≥0. Then, it follows that the discrete solution (un+1 h, vn+1 h)to (2.15) and (2.16) satisfies un+1 h>0and vn+1 h≥0.(3.16) Proof. Suppose that there exists ai∈ Nhbeing the minimum such that un+1 i:= un+1 h(ai)≤0. Moreover, without loss of generality, we may assume that i < j for all j∈I(Ωai); otherwise, one should only take j∈I(Ωai) such that i < j. Let I∗(Ωai) = {j∈I(Ωai) : un+1 j=un+1 i}and let I∗ c(Ωai) be its complementary. Then, choose xh=ϕaiin (2.15) to get, using (un+1 h∇vn h,∇ϕi)∗=X j∈I∗ c(Ωai) vn j−vn i g0 ε(un+1 j)−g0 ε(un+1 i)(∇ϕaj,∇ϕai) ×(un+1 j−un+1 i)(ϕai(ai)−ϕai(aj)), that k−1(1, ϕai)un+1 i+X j∈I∗ c(Ωai) aji(un+1 h, vn h)(un+1 j−un+1 i)(ϕi(ai)−ϕi(aj)) +X j∈I∗ c(Ωai) νu ji(un+1 h, vn h)(un+1 j−un+1 i)(ϕai(aj)−ϕai(ai)) =k−1un i(1, ϕai), March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations 625 where we used X j∈I(Ωai) un+1 j(∇ϕaj,∇ϕai) =X j∈Ic(Ωai) (∇ϕaj,∇ϕai)(un+1 j−un+1 i)(ϕi(ai)−ϕi(aj)), and defined aji(un+1 h, vn h) = 1−vn j−vn i g0 ε(un+1 j)−g0 ε(un+1 i)!(∇ϕaj,∇ϕai). Therefore k−1(1, ϕai)un+1 i+X j∈I∗ c(Ωai)aji(un+1 h, vn h)−νu ji(un+1 h, vn h)(un+1 j−un+1 i) =k−1(1, ϕai)un i.(3.17) Since ¯αi(un+1 h) = 1 from Lemma 2.2, we have aji(un+1 h, vn h)−νu ji(un+1 h, vn h)≤0,for all j∈I∗ c(Ωai).(3.18) Applying (3.18) to (3.17) yields 0< un i≤un+1 i, which contradicts our assumption. An argument similar to that in the proof of Lemma 3.1 shows that vn+1 h≥0. Remark 3.1. Once we proved that the discrete solution un+1 his positive, the use of the truncating operator in (2.14) can be neglected. Remark 3.2. Algorithm 2 for γ= 0 produces non-negative discrete solutions to (2.16) provided that acuteness is assumed for Kh. The following lemma is the counterpart of Lemma 3.2. Lemma 3.5. (L1(Ω)-bounds) Under the conditions of Lemma 3.4, the discrete solution pair (un+1 h, vn+1 h)∈X2 hcomputed via (2.15) and (2.16) fulfills kun+1 hkL1(Ω) =ku0 hkL1(Ω),(3.19) and kvn+1 hkL1(Ω) ≤ kv0 hkL1(Ω) +ku0 hkL1(Ω).(3.20) Proof. Using the reasoning of Lemma 3.2 proves (3.19) and (3.20) from testing (2.15) and (2.16) by xh= 1 and noting (Bu 2(un+1 h, vn h)un+1 h,1) = 0, March 17, 2023 10:25 WSPC/103-M3AS 2350014 626 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu and (Bv 2(vn+1 h)vn+1 h,1) = 0. It is now shown that system (2.15) and (2.16) for γ= 0 possesses a discrete energy law associated with Eh(uh, vh). =1 2k∇vhk2+1 2kvhk2 h−(uh, vh)h+ (uh,log uh)h. Lemma 3.6. (A discrete energy law) Let γ= 0. Assume that This a weakly acute triangulation for d= 2 and acute triangulation for d= 3 and define I∗={(i, j)∈ I×I:un+1 i=un+1 j}with I∗ cbeing its complement. Then,the discrete solution (un+1 h, vn+1 h)∈X2 hcomputed via (2.15) and (2.16) satisfies Eh(un+1 h, vn+1 h)− Eh(un h, vn h) + ND(un+1 h, vn+1 h) + kkδtvn+1 hk2 h −kX i<j∈I∗ c (1 −¯α#(un+1 h))  g0 ε(un+1 j)−g0 ε(un+1 i) un+1 j−un+1 i!−1 2 (un+1 j−un+1 i) − g0 ε(un+1 j)−g0 ε(un+1 i) un+1 j−un+1 i!1 2 (vn j−vn i) 2 (∇ϕaj,∇ϕai) −kX i<j∈I∗ un+1 i(vn j−vn i)2(∇ϕaj,∇ϕai) = 0,(3.21) where ND(un+1 h, vn+1 h) =k∇(vn+1 h−vn h)k2 L2(Ω) +kvn+1 h−vn hk2 h+k(g00 ε(un+θ h),(δtun+1 h)2)h, and ¯α#∈ {αi, αj,0}is chosen such that νu ji(un+1 h, vn h) = ¯α#(un+1 h) 1−vn ¯ #−vn # g0 ε(un+1 ¯ #)−g0 ε(un+1 #)!(∇ϕa¯ #,∇ϕa#), with ¯ # = jwhen # = iand ¯ # = iwhen # = j. Proof. Select xh=ihg0 ε(un+1 h)−vn hin (2.15) and xh=δtvn+1 hin (2.16) to get (δtun+1 h, ihg0 ε(un+1 h)−vn h)h+ (∇un+1 h,∇(ihg0 ε(un+1 h)−vn h)) −(un+1 h∇vn h,∇(ihg0 ε(un+1 h)−vn h))∗ + (Bu 2(un+1 h, vn h)un+1 h, ihg0 ε(un+1 h)−vn h)=0,(3.22) and 1 2kk∇vn+1 hk2 L2(Ω) +1 2kkvn+1 hk2 h−1 2kk∇vn hk2 L2(Ω) −1 2kkvn hk2 h March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations 627 +1 2kk∇(vn+1 h−vn h)k2 L2(Ω) +1 2kkvn+1 h−vn hk2 h +kδtvn+1 hk2 h−(un+1 h, δtvn+1 h)h= 0. (3.23) We next pair some terms from (3.22) and (3.23) in order to handle them together. It is not hard to see that −(δtun+1 h, vn h)h−(un+1 h, δtvn+1 h)h=−1 k(un+1 h, vn+1 h)h+1 k(un h, vn h)h.(3.24) Now write (∇un+1 h,∇(ihg0 ε(un+1 h)−vn h)) =−X i<j∈I∗ c g0 ε(un+1 j)−g0 ε(un+1 i) un+1 j−un+1 i (un+1 j−un+1 i)2(∇ϕaj,∇ϕai) +X i<j∈I∗ c (un+1 j−un+1 i)(vn j−vn i)(∇ϕaj,∇ϕai), and (un+1 h∇vn h,∇(ihg0 ε(un+1 h)−vn h))∗ =−X i<j∈I∗ c (un+1 j−un+1 i)(vn j−vn i)(∇ϕaj,∇ϕai) +X i<j∈I∗ c un+1 j−un+1 i g0 ε(un+1 j)−g0 ε(un+1 i)(vn j−vn i)2(∇ϕaj,∇ϕai) +X i<j∈I∗ un+1 i(vn j−vn i)2(∇ϕaj,∇ϕai). Therefore (∇un+1 h,∇(ihg0 ε(un+1 h)−vn h)) −(un+1 h∇vn h,∇(ihg0 ε(un+1 h)−vn h))∗ =−X i<j∈I∗ c  un+1 j−un+1 i g0 ε(un+1 j)−g0 ε(un+1 i)!1 2 (un+1 j−un+1 i) − g0 ε(un+1 j)−g0 ε(un+1 i) un+1 j−un+1 i!1 2 (vn j−vn i) 2 (∇ϕaj,∇ϕai) −X i<j∈I∗ un+1 i(vn j−vn i)2(∇ϕaj,∇ϕai)>0, March 17, 2023 10:25 WSPC/103-M3AS 2350014 628 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu since (∇ϕaj,∇ϕai)≤0 due to the acuteness of Th. Moreover (B2(un+1 h, vn h)un+1 h, ihg0 ε(un+1 h)−vn h) =X i<j∈I∗ νu ji(un+1 h, vn h)(un+1 i−un+1 j)(g0 ε(un+1 j)−g0 ε(un+1 i)−(vn i−vn j)). Let us assume that νu ji(un+1 h, vn h) = ¯αi(un+1 h) 1−vn j−vn i g0 ε(un+1 j)−g0 ε(un+1 i)!(∇ϕaj,∇ϕai). Then 1−vn j−vn i g0 ε(un+1 j)−gε(un+1 i)!(un+1 j−un+1 i)(g0 ε(un+1 j)−g0 ε(un+1 i)−(vn i−vn j)) =g0 ε(un+1 j)−g0 ε(un+1 i) un+1 j−un+1 i (un+1 j−un+1 i)2−2(vn j−vn i)(un+1 j−un+1 i) +un+1 j−un+1 i g0 ε(un+1 j)−g0 ε(un+1 i)(vn j−vn i)2 = un+1 j−un+1 i g0 ε(un+1 j)−g0 ε(un+1 i)!1 2 (un+1 j−un+1 i) − g0 ε(un+1 j)−g0 ε(un+1 i) un+1 j−un+1 i!1 2 (vn j−vn i) 2 . As a result, we have (∇un+1 h,∇(ihg0 ε(un+1 h)−vn h)) −(un+1 h∇vn h,∇(ihg0 ε(un+1 h)−vn h))∗ + (Bu 2(un+1 h, vn h)un+1 h, ihg0 ε(un+1 h)−vn h) =−X i<j∈I∗ (1 −¯αi(un+1 h))  g0 ε(un+1 j)−g0 ε(un+1 i) un+1 j−un+1 i!−1 2 (un+1 j−un+1 i) − g0 ε(un+1 j)−g0 ε(un+1 i) un+1 j−un+1 i!1 2 (vn j−vn i) 2 (∇ϕaj,∇ϕai)>0,(3.25) since 0 ≤¯αi(un+1 h)≤1. Analogously, when νu ji(un+1 h, vn h) = ¯αj(un+1 h) 1−vn i−vn j g0 ε(un+1 i)−g0 ε(un+1 j)!(∇ϕai,∇ϕaj), holds. March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations 629 A Taylor polynomial of gεround un+1 hevaluated at un hyields gε(un h) = gε(un+1 h)−g0 ε(un+1 h)(un+1 h−un h) + g00 ε(un+θ h) 2(un+1 h−un h)2, where θ∈(0,1) such that un+θ h=θun+1 h+ (1 −θ)un h. Hence (δtun+1 h, g0 ε(un+1 h))h=1 k(gε(un+1 h),1)h−1 k(gε(un h),1)h+k 2(g00 ε(un+θ h),(δtun+1 h)2)h. The equality (3.21) follows by adding (3.22) and (3.23) and invoking (3.24) and (3.25). Remark 3.3. Lemma 3.6 can readily be (at least) extended to meshes with right quadrilaterals (d= 2) or hexahedra (d= 3) when bilinear finite elements are used, since (∇ϕai,∇ϕaj)≤0 for all i6=j∈I. We summarize the results of Algorithm 2 for γ= 0 in the following. Theorem 3.2. Let γ= 0. Assume that This a weakly acute triangulation for d= 2 and acute triangulation for d= 3. Then,it follows that the sequence of the discrete solution pairs {(un h, vn h)}N n=0 defined by Algorithm 2satisfies (a) Lower bounds un h>0and vn h≥0. (b) L1(Ω)-bounds kun hkL1(Ω) =ku0hkL1(Ω), and kvn hkL1(Ω) ≤ kv0hkL1(Ω) +ku0hkL1(Ω). (c) Energy law Eh(un h, vn h)− Eh(u0 h, v0 h) + ND(vn+1 h) + kkδtvn+1 hk2 L2(Ω) −k n−1 X m=0 X i<j∈I∗ c (1 −¯α#(um+1 h))  g0 ε(um+1 j)−g0 ε(um+1 i) um+1 j−um+1 i!1 2 ×(um+1 j−um+1 i)− um+1 j−um+1 i g0 ε(um+1 j)−g0 ε(um+1 i)!−1 2 (vm j−vm i) 2 ×(∇ϕaj,∇ϕai)−kX i<j∈I∗ un+1 i(vn j−vn i)2(∇ϕaj,∇ϕai)=0, for all n= 1, . . . , N. March 17, 2023 10:25 WSPC/103-M3AS 2350014 630 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu Remark 3.4. Observe that lim ε→0gε(s) = log s, for all s > 0. Therefore, if we take the limit as ε→0, we formally obtain, from (3.16), Algorithm 2 with the following definitions. The stabilizing term is now given by (Bu 2(un+1 h, vn h)un+1 h, xh) = X i<j∈I νu ji(un+1 h, vn h)(un+1 j−un+1 i)(xj−xi), where νu ji(un+1 h, vn h) = max{¯αi(un+1 h)fij,¯αj(un+1 h)fji,0},for i6=j, with fu ij =       1−vn j−vn i log(un+1 j)−log(un+1 i)!(∇ϕaj,∇ϕai),if un+1 j6=un+1 i, 0,if un+1 j=un+1 i, and the chemotaxis term is approximated as (xh∇˜xh,∇¯xh)∗=X i<j∈I γji(xh)(˜xj−˜xi)(¯xi−¯xj)(∇ϕaj,∇ϕai), where γji(xh) =    xj−xi log(xj)−log(xi),if xj6=xi, xi,if xj=xi. 4. Numerical Results Two new algorithms have been proposed and analyzed in the previous sections for approximating problem (1.1)–(1.3). Discrete solutions to Algorithms 1 and 2 for γ= 1 have lower bounds and preserve mass, whereas discrete solutions to Algorithm 2 for γ= 0 satisfy a discrete energy law as well provided that the mesh be weakly acute in dimension two or acute in dimension three. In this section, we present numerical results for two test problems using both algorithms to validate the theoretical results: The first one has a smooth solution and is used to investigate the effect of the numerical diffusion stemming from the stabilizing terms. The second example concerns blowup phenomena; thus, we assess if the two stabilized algorithms can deal with solutions that might develop singularities in finite time. As Algorithms 1 and 2 are nonlinear, Picard’s method with backtracking is used to carry out the iterations for solving the nonlinear system at each time step. The two stopping criteria used are to require the L∞(Ω)-norm of the residual to be less than 10−6or the increment to be less than 10−16. Moreover, we take q= 2 for Algorithms 1 and 2 in (2.11) and (2.18), respectively. Recall that the shock detector (2.11) for Algorithm 1 acts on both maxima and minima, whereas the March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations 631 shock detector (2.18) acts on minima only. In addition, for Algorithm 2, we set γ= 1 and ε= 10−6. 4.1. Smooth coalescence In this test, we consider the approximation of the Keller–Segel model (1.1)–(1.3) for Ω = (−π, π)2and the initial conditions u0= sin2(x)sin2(y) and v0= cos(x) + cos(y)+2. For both algorithms, we consider Xhto be a bilinear finite element space constructed over a uniformly structured 40 ×40 grid having mesh size h= 0.11107. The time step size is k= 0.02. Since RΩu0(x)dx=π2∈(0,4π), the expected dynamics is a smooth coalescence of the cell density because the chemoattractant density initially concentrates around the center of the domain. Figure 2 shows the evolution in time of the L1(Ω)-norm for the chemoattractant and cell densities. In particular, the mass conservation is numerically verified, as Lemmas 3.1 and 3.4 predict for the cell density. It also shows the competition between production and degradation of chemoattractant, where the latter prevails over the former, which implies a mass loss. As can be seen from Fig. 3 (first column), cells are initially distributed in four groups with the highest concentrations in the center of each quadrant. Instead, the highest concentration of chemoattractant is located at the origin as shown in Fig. 4. Then, the four groups of cells move toward the origin attracted by the chemoattractant while being diffused as displayed in Figs. 3 and 4 at times t= 1.3 and 2. Furthermore, the profile along the diagonal of the domain and the evolution of the energy functional Eh(uh, vh) can be found in Fig. 5 at times t= 0 and 2, and Fig. 6, respectively. It is evident that Eh(uh, vh) decays over time; therefore, the system evolves toward a steady-state solution. To somehow ascertain the numerical diffusion introduced by the stabilizing terms, we carry out the same numerical simulation via the algorithm proposed in Ref. 11, which comprises a standard finite method and a semi-implicit/implicit Fig. 2. Evolution of kuhkL1(Ω) and kvhkL1(Ω). March 17, 2023 10:25 WSPC/103-M3AS 2350014 632 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu Fig. 3. (Color online) Colormaps of uhat times t= 0,0.3 and 2. Euler time-stepping; that is, Algorithm 1 without the stabilizing terms. No significant differences are highlighted in Figs. 3 and 4 (bottom), implying that the shock-capturing works only in the presence of steep gradients. Figure 5 (bottom) shows the profiles at times t= 0 and 2 for the nonstabilized discrete solution. 4.2. Blowup phenomena In the following two numerical tests, we assess our two algorithms with initial data which generate solutions to the Keller–Segel problem (1.1)–(1.3) that blowup in finite time. In this fashion, the robustness and reliability of discrete solutions computed by Algorithms 1 and 2 are examined. The need of using stabilizing terms is a matter of the utmost importance. This was made manifest11 when a nonstabilized finite element method failed to prevent discrete solutions from being negative. In March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations 633 Fig. 4. (Color online) Colormaps of vhat times t= 0,0.3 and 2. this case, the usage of a positivity-preserving scheme is mandatory. Otherwise, the high gradients yield an oscillatory behavior of discrete solutions taking negative density values. 4.2.1. Blowup in the center We consider the Keller–Segel equation (1.1)–(1.3) on Ω = (0,1)2. As opposed to the previous test, both initial conditions u0and v0are now concentrated around about the point (0.5,0.5). In particular, we use initial conditions based on Example 5.2 in Ref. 19, i.e. u0= 840 exp(−84((x−0.5)2+ (y−0.5)2)), and v0= 420 exp(−42((x−0.5)2+ (y−0.5)2)). March 17, 2023 10:25 WSPC/103-M3AS 2350014 640 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu Fig. 20. Algorithm 2: Profiles of uhand vhalong the plane x= 0 at times t= 0, 0.15 and 0.2. Fig. 21. Algorithm 2: Evolution of kuhkL∞(Ω) and kvhkL∞(Ω), and kuhkL1(Ω) and kvhkL1(Ω). Fig. 22. Algorithm 2: Evolution of Eh(uh, vh). discrete solution at t= 2 just as the maximum is reached at the boundary whose support is one macroelement consisting of three triangles. In Figs. 18 and 21 (left), we depict the evolution of the L∞(Ω)-norm of the cell and chemoattractant densities whose maxima overgrow simultaneously on the boundary as seen in Figs. 17 and 21. It is observed that the blowup time is about 0.19, reaching a cell-density value of approximately 9.1·105for Algorithm 1 and 7.2·105for Algorithm 2, which arrives at the locking stage earlier. For the mesh used in this experiment, the locking cell density corresponds to maxn∈{0,...,N}kun hkL∞(Ω) ≤4 h2 min√3ku0kL1(Ω) ≈9.6·105. Here the L∞(Ω)-norm of the chemoattractant density undergoes a rapid growth at the singular point as well, which does not occur in the center-blowup March 17, 2023 10:25 WSPC/103-M3AS 2350014 Bound-preserving finite element approximations of the Keller–Segel equations 641 setting. The dynamics of the L1(Ω)-norm plotted in Figs. 18 and 21 (right) indicates that the chemoattractant mass increases; this contrasts with Figs. 10 and 13 (right), where the L1(Ω)-norm keeps constant practically. As for the energy functional Eh(uh, vh), a considerably pronounced drop in magnitude arises before the computational blowup time in Fig. 22. Acknowledgments Juan Vicente Guti´errez-Santacreu was partially supported by the Spanish Grant No. PGC2018-098308-B-I00 from Ministerio de Ciencias e Innovaci´on — Agencia Estatal de Investigaci´on with the participation of FEDER and by the Andalusian Grant No. P20 01120 from Junta de Andaluc´ıa (Consejer´ıa de Econom´ıa, Conocimiento, Empresas y Universidad). Los Alamos National Laboratory, an affirmative action/equal opportunity employer, is operated by Triad National Security, LLC for the National Nuclear Security Administration of U.S. Department of Energy under contract 89233218CNA000001. Los Alamos National Laboratory strongly supports academic freedom and a researcher’s right to publish; as an institution, however, the laboratory does not endorse the viewpoint of a publication or guarantee its technical correctness (LA-UR-22-27083). References 1. G. Arumugam and J. Tyagi, Keller–Segel chemotaxis models: A review, Acta Appl. Math. 171 (2021) 82. 2. S. Badia and J. Bonilla, Monotonicity-preserving finite element schemes based on differentiable nonlinear stabilization, Comput. Methods Appl. Mech. Engrg. 313 (2017) 133–158. 3. N. Bellomo, A. Bellouquid, Y. Tao and M. Winkler, Toward a mathematical theory of Keller–Segel models of pattern formation in biological tissues, Math. Models Methods Appl. Sci. 25 (2015) 1663–1763. 4. N. Bellomo, N. Outada, J. Soler, Y. Tao and M. Winkler, Chemotaxis and crossdiffusion models in complex environments: Models and analytic problems toward a multiscale vision, Math. Models Methods Appl. Sci. 32 (2022) 713–792. 5. J. Bonilla de Toro, Monotonicity-preserving finite element methods for hyperbolic problems, Doctoral thesis, Universitat Polit`ecnica de Catalunya (2019). 6. S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, Texts in Applied Mathematics, Vol. 15, 3rd edn. (Springer, 2008). 7. A. Chertock, Y. Epshteyn, H. Hu and A. Kurganov, High-order positivity-preserving hybrid finite-volume-finite-difference methods for chemotaxis systems, Adv. Comput. Math. 44 (2018) 327–350. 8. A. Chertock and A. Kurganov, A second-order positivity preserving central-upwind scheme for chemotaxis and haptotaxis models, Numer. Math. 111 (2008) 169–205. 9. V. Girault and J.-L. Lions, Two-grid finite-element schemes for the transient Navier– Stokes problem, M2AN Math. Model. Numer. Anal. 35 (2001) 945–980. 10. L. Guo, X. H. Li and Y. Yang, Energy dissipative local discontinuous Galerkin methods for Keller–Segel chemotaxis model, J. Sci. Comput. 78 (2019) 1387–1404. March 17, 2023 10:25 WSPC/103-M3AS 2350014 642 S. Badia, J. Bonilla & J. V. Guti´errez-Santacreu 11. J. V. Guti´errez-Santacreu and J. R. Rodr´ıguez-Galv´an, Analysis of a fully discrete approximation for the classical Keller–Segel model: Lower and a priori bounds, Comput. Math. Appl. 85 (2021) 69–81. 12. M. A. Herrero and J. J. L. Vel´azquez, A blow-up mechanism for a chemotaxis model, Ann. Scuola Norm. Sup. Pisa Cl. Sci. (4)24 (1997) 633–683. 13. D. Horstmann, The nonsymmetric case of the Keller–Segel model in chemotaxis: Some recent results, NoDEA Nonlinear Differential Equations Appl. 8(2001) 399–423. 14. D. Horstmann, From 1970 until present: The Keller–Segel model in chemotaxis and its consequences. I, Jahresber. Deutsch. Math.-Verein. 105 (2003) 103–165. 15. D. Horstmann, From 1970 until present: The Keller–Segel model in chemotaxis and its consequences. II, Jahresber. Deutsch. Math.-Verein. 106 (2004) 51–69. 16. D. Horstmann and G. Wang, Blow-up in a chemotaxis model without symmetry assumptions, European J. Appl. Math. 12 (2001) 159–177. 17. E. F. Keller and L. A. Segel, Initiation of slide mold aggregation viewed as an instability, J. Theor. Biol. 26 (1970) 399–415. 18. E. F. Keller and L. A. Segel, Model for chemotaxis, J. Theor. Biol. 30 (1971) 225–234. 19. X. H. Li, C.-W. Shu and Y. Yang, Local discontinuous Galerkin method for the Keller–Segel chemotaxis model, J. Sci. Comput. 73 (2017) 943–967. 20. T. Nagai, T. Senba and K. Yoshida, Application of the Trudinger-Moser inequality to a parabolic system of chemotaxis, Funkcial. Ekvac. 40 (1997) 411–433. 21. N. Saito, Error analysis of a conservative finite-element approximation for the Keller– Segel system of chemotaxis, Commun. Pure Appl. Anal. 11 (2012) 339–364. 22. L. R. Scott and S. Zhang, Finite element interpolation of non-smooth functions satisfying boundary conditions, Math. Comp. 54 (1990) 483–493. 23. R. Strehl, A. Sokolov, D. Kuzmin, D. Horstmann and S. Turek, A positivity-preserving finite element method for chemotaxis problems in 3D, J. Comput. Appl. Math. 239 (2013) 290–303. 24. R. Strehl, A. Sokolov, D. Kuzmin and S. Turek, A flux-corrected finite element method for chemotaxis problems, Comput. Methods Appl. Math. 10 (2010) 219–232. 25. M. Sulman and T. Nguyen, A positivity preserving moving mesh finite element method for the Keller–Segel chemotaxis model, J. Sci. Comput. 80 (2019) 649–666. 26. M. Winkler, Aggregation versus global diffusive behavior in the higher-dimensional Keller–Segel model, J. Differential Equations 248 (2010) 2889–2905. 27. M. Winkler, Finite-time blow-up in the higher-dimensional parabolic–parabolic Keller–Segel system, J. Math. Pures Appl. 100 (2013) 748–767.