∇∇ ΔΔ αα γγ μμ ℕN ℝR ∂∂ νν ΦΦ ττ ωω ∞∞ ϕϕ title

# Flux-corrected transport stabilization of an evolutionary cross-diffusion
cancer invasion model

Shahin Heydari 111e-mail<EMAIL_ADDRESS>, corresponding author
Charles University, Faculty of Mathematics and Physics, Sokolovská 83, 18675
Praha 8, Czech Republic  Petr Knobloch 222e-mail<EMAIL_ADDRESS>Charles University, Faculty of Mathematics and Physics, Sokolovská 83, 18675
Praha 8, Czech Republic  Thomas Wick 333e-mail<EMAIL_ADDRESS>hannover.de Leibniz University Hannover, Institute of Applied Mathematics,
Welfengarten 1, 30167 Hannover, Germany

###### Abstract

In the present work, we investigate a model of the invasion of healthy tissue
by cancer cells which is described by a system of nonlinear PDEs consisting of
a cross-diffusion-reaction equation and two additional nonlinear ordinary
differential equations. We show that when the convective part of the system,
the chemotactic term, is dominant, then straightforward numerical methods for
the studied system may be unstable. We present an implicit finite element
method using conforming $P_{1}$ or $Q_{1}$ finite elements to discretize the
model in space and the $\theta$-method for discretization in time. The
discrete problem is stabilized using a nonlinear flux-corrected transport
approach. It is proved that both the nonlinear scheme and the linearized
problems used in fixed-point iterations are solvable and positivity
preserving. Several numerical experiments are presented in 2D using the
deal.II library to demonstrate the performance of the proposed method.

Key words: cancer invasion, cross-diffusion equation, FEM-FCT stabilization,
positivity preservation, existence of solutions
AMS Classification (2020): 65M22, 65M60, 92C17, 35Q92,

abstract=true

##  1 Introduction

Keller and Segel [29, 28] proposed the first mathematical model for
description of chemotactical processes. Chemotaxis refers to the motion in the
direction to (or away from) the position of higher concentration based on the
gradient of chemical substances and its chemotacitivity character which
controls the speed of this motion. Their model has been widely extended and
followed to develop more sophisticated and complex chemotaxis models and
played a vitally important role in many areas of science, in particular in
medical and biological applications, for example, bacteria and cell
aggregation [1, 39, 51], tumor angiogenisis and invasion [3, 11, 12, 13],
biological pattern formation [2, 51], and immune cell migration [54]. From the
analytical point of view, mathematical analysis for chemotaxis systems of
equations is a challenge and causes many questions especially in the context
of the existence and uniqueness of solutions. In the last three decades, many
researchers have been actively involved and answered some of these questions
[40, 16, 23, 49, 22, 10]. From the numerical point of view, so far a great
deal of research on chemotaxis models has been done in various areas,
including the finite difference method [11, 12, 31], discontinuous Galerkin
method [18, 35], finite element method [43, 56, 57], finite volume method
[20], operator-splitting methods [42], or fractional step algorithms [52].
However, many analytical and numerical aspects are still untouched and call
for further investigation.

The chemotaxis problems are usually strongly coupled nonlinear systems of
equations whose solutions represent concentrations or densities and need to be
non-negative in order to satisfy the physics behind the system. Hence, it is
difficult to construct an efficient and accurate numerical method that does
not produce solutions with negative values. Another interesting aspect is
singular, spiky and oscillatory behavior of the solutions. In particular, when
the chemotaxis term dominates the diffusion and reaction terms, in other words
large chemosensitivity is present which corresponds to large Reynolds numbers,
it may give rise to nonphysical oscillations in the solution. To overcome this
problem, stabilization methods can be applied. Up to now, many scientists used
flux-corrected transport (FCT) algorithms, i.e., nonlinear high-resolution
schemes introduced by Boris and Book [8, 7, 9], later developed based on
linear finite element discretizations by Kuzmin, Löhner et al. [36, 33, 34,
32], and further extended to linear and nonlinear space-time FEM-FCT in [19].
In [46], an implicit flux-corrected transport scheme was developed and applied
to three benchmark examples of the general Keller–Segal model in two spatial
dimensions. It was shown that the proposed method is positivity preserving and
sufficiently accurate, even in the cases where solutions blow up in the center
or at the boundary of the domain. The investigations of the blow-up behavior
of the solutions were further extended to three spatial dimensions in [47]. In
[45, 44], an FEM-FCT scheme was coupled with a level-set method to obtain
positivity preserving solutions on a stationary surface and evolving-in-time
surfaces. It was shown that the proposed method is able to produce accurate
numerical solutions, which makes it possible to couple the partial
differential equations defined on a specific domain with the PDEs that are
defined on the surface of this domain. This scheme was further used with
operator-splitting techniques to solve chemotaxis models in 3D. The operator-
splitting method splitted a 3D problem into a sequence of 1D subproblems and
the FEM-FCT algorithm was used to solve each 1D subproblem separately [25]. In
[48], the authors used an efficient adaptive moving mesh finite element
approach based on the parabolic Monge–Ampère method for determining the
coordinate transformation for the adaptive mesh combined with an FCT scheme
which guarantees the non-negativity of the solutions. As a result, the
computational cost was significantly reduced. All aforementioned techniques
were also applied to the same benchmark examples. A different case was studied
in [24], where the authors used the pressure-correction scheme and flux-
corrected transport algorithm to propose an efficient linear positivity-
preserving method for the solution of chemotaxis–Stokes equations.

In this work, we focus on a cancer-invasion model developed in [41], modeling
the motion of cancer cells, degradation of extracellular matrix, and certain
enzymes (e.g., protease). The extracellular matrix is degraded upon contact
with protease which is produced where cancer cells and extracellular matrix
meet and decay over the time. In [21], we extended the proposed model by a
diffusion term, gave a rigorous proof for the existence of the global
classical solution and presented numerical results for a Galerkin finite
element discretization. In the present paper, a diffusion term is not
considered, which makes the problem more challenging. In [30], one of the
authors of the present paper applied a positivity preserving non-standard
finite difference method to solve the nonlinear system in 1D, see also [14]
for related approaches. Here, we consider the finite element method and apply
the FCT technique to guarantee the positivity preservation. First, however, we
consider the more diffusive nonlinear low-order method. An additional
nonlinearity is then introduced by the flux correction. We prove that both
nonlinear problems are solvable and positivity preserving. To the best of our
knowledge, the current work is a first attempt to gain an insight into the
applicability of the FCT technique to the numerical solution of a chemotaxis
system without self-diffusion and to provide a rigorous analysis of the
solvability and positivity preservation. Note that the existence and
uniqueness for the FEM-FCT method applied to linear evolutionary convection-
diffusion equations has been addressed only recently in [27, 26]. We also
present a fixed-point algorithm for the iterative solution of the FCT
discretization and prove that it is well posed and provides a non-negative
solution at each step. Consequently, the non-negativity of the approximate
solution is guaranteed independently of the choice of a stopping criterion.
The properties of the proposed FCT scheme are illustrated by various numerical
simulations carried out using our newly designed algorithm in the deal.II
library [4, 5].

The outline of this paper is as follows. In Section 2, we formulate the
mathematical model which is discretized by the Galerkin method in Section 3.
Then, the FCT stabilization is introduced in Section 4, where also the
solvability and positivity preservation is proved. The fixed-point algorithm
is proposed and investigated in Section 5. In Section 6, we report several
numerical simulations in two spatial dimensions carried out for various
regimes. Finally, our results are summarized in Section 7.

##  2 Mathematical model

In this section, we discuss the following nondimensionalized continuous model
of a malignant cancer invasion proposed by Perumpanani et al. in [41, 37]. The
model contains three unknown variables, namely the cancer cell density
$u=u(x,t)$, connective tissue $c=c(x,t)$, and protease $p=p(x,t)$, and it
consists of the equations

$\displaystyle\dfrac{\partial u}{\partial
t}=\mu\,u(1-u)-\chi\,\nabla\cdot(u\nabla c)\qquad$ $\displaystyle\text{in
$\Omega\times(0,T]$}\,,$ (2.1) $\displaystyle\dfrac{\partial c}{\partial
t}=-pc$ $\displaystyle\text{in $\Omega\times(0,T]$}\,,$ (2.2)
$\displaystyle\dfrac{\partial p}{\partial t}=\epsilon^{-1}(uc-p)$
$\displaystyle\text{in $\Omega\times(0,T]$}\,,$ (2.3)

where $\Omega$ is a bounded polyhedral domain in $ℝ^{d},d\in\\{1,2,3\\}$,
$[0,T]$ is a time interval, and $\mu$, $\chi$, $\epsilon$ are positive
constants. Here, $\mu$ and $\chi$ denote the proliferation and haptotaxis rate
of cancer cells, respectively, and the parameter $\epsilon$ is supposed to be
small since the units of connective tissues and invasive cells are much larger
than the protease. In the process of invasion, the connective tissue is
affected by the invasive flux of $u\nabla c$ into its compartment. Since the
connective tissue does not contain any empty space large enough for passing of
passive cancer cells, it degrades by protease which is produced by invasive
cancer cells upon contact with connective tissue. It can be shown that if the
initial conditions of the above model are non-negative, then the computed
solutions stay non-negative at all times, for more details see [37, 38, 14]
and the references therein.

The system (2.1)–(2.3) is subjected to the homogeneous Neumann boundary
condition

$u\,\frac{\partial c}{\partial
n}=0\qquad\text{on}\,\,\,\partial\Omega\times[0,T]\,,$ (2.4)

where $n$ is the unit outward normal vector on $\partial\Omega$. The above
equations are endowed with the initial conditions

$u(x,0)=u^{0}(x)\,,\quad c(x,0)=c^{0}(x)\,,\quad p(x,0)=p^{0}(x)\,,\qquad
x\in\Omega\,,$ (2.5)

where $u^{0},\,c^{0},\,p^{0}:\Omega\to[0,1]$ are given functions.

In [21], we considered a modified version of (2.1)–(2.3) containing an extra
diffusion term in (2.1). Precisely, instead of the equation (2.1), we
considered

$\dfrac{\partial u}{\partial t}=\mu\,u(1-u)-\chi\,\nabla\cdot(u\nabla
c)+\alpha^{-1}\Delta u\qquad\quad\text{in $\Omega\times(0,T]$}$ (2.6)

with a positive constant $\alpha$. This required to replace the boundary
condition (2.4) by

$\alpha^{-1}\,\frac{\partial u}{\partial n}=\chi\,u\,\frac{\partial
c}{\partial n}\qquad\text{on}\,\,\,\partial\Omega\times[0,T]\,.$

Thus, the problem considered in this paper corresponds to the limit case
$\alpha\to\infty$ of the problem from [21]. In that paper, we proved the
existence of global classical solutions for two- and three-dimensional bounded
domains $\Omega$ with smooth boundaries and we proved that these solutions are
non-negative. Moreover, we showed that by fixing the proliferation rate $\mu$
and varying the haptotaxis $\chi$ one can make either the diffusion or the
transport of the cancer cells dominant. The domination of the convection term
can produce spurious oscillations and a blow-up in the solution of the system
as it is the case to be considered in here.

##  3 A Galerkin discretization

The solution of the problem (2.1)–(2.5) satisfies

$\displaystyle\left(\dfrac{\partial u}{\partial
t},v\right)=\mu\,\big{(}u(1-u),v\big{)}+\chi\,\big{(}u\nabla c,\nabla
v\big{)}\qquad$ $\displaystyle\text{in $(0,T]$ and for $v\in
H^{1}(\Omega)$}\,,$ (3.1) $\displaystyle
c(x,t)=c^{0}(x)\,{\mathrm{e}}^{-\int_{0}^{t}\,p(x,s)\,{\mathrm{d}}s}$
$\displaystyle\forall\,\,(x,t)\in\Omega\times[0,T]\,,$ (3.2) $\displaystyle
p(x,t)={\mathrm{e}}^{-t/\epsilon}\left[p^{0}(x)+\frac{1}{\epsilon}\,\int_{0}^{t}\,u(x,s)\,c(x,s)\,{\mathrm{e}}^{s/\epsilon}\,{\mathrm{d}}s\right]\quad$
$\displaystyle\forall\,\,(x,t)\in\Omega\times[0,T]\,,$ (3.3)

where $(\cdot,\cdot)$ denotes the inner product in $L^{2}(\Omega)$ or
$L^{2}(\Omega)^{d}$. To define an approximate solution of (2.1)–(2.5), we
first introduce a triangulation ${\mathcal{T}}_{h}$ of $\Omega$ consisting of
simplicial (for $d=1,2,3$), quadrilateral (for $d=2$) or hexahedral (for
$d=3$) shape-regular cells possessing the usual compatibility properties (see,
e.g., [15]). For any cell $K\in{\mathcal{T}}_{h}$, we denote by $h_{K}$ the
diameter of $K$ and assume that $h_{K}\leq h$. We denote by $V_{h}\subset
H^{1}(\Omega)$ the usual conforming $P_{1}$ or $Q_{1}$ finite element space
constructed using the triangulation ${\mathcal{T}}_{h}$. Let
$\phi_{1},\dots,\phi_{M}$ be the standard basis functions of $V_{h}$
associated with the vertices $x_{1},\dots,x_{M}$ of ${\mathcal{T}}_{h}$. Thus,
the basis functions are non-negative and satisfy $\phi_{i}(x_{j})=\delta_{ij}$
for $i,j=1,\dots,M$, where $\delta_{ij}$ is the Kronecker symbol. Any function
$v_{h}\in V_{h}$ can be identified with a coefficient vector
$\mathbf{v}=(v_{j})_{j=1}^{M}$ with respect to these basis functions.
Precisely, introducing the bijective operator $\pi_{h}:{\mathbb{R}}^{M}\to
V_{h}$ by

$\pi_{h}\mathbf{v}=\sum_{j=1}^{M}\,v_{j}\,\phi_{j}\,,$

one has $v_{h}=\pi_{h}\mathbf{v}$. The assumed shape regularity of
${\mathcal{T}}_{h}$ implies that

$\|\nabla\phi_{i}\|_{L^{2}(K)}\leq\kappa\,h_{K}^{d/2-1}\qquad\forall\,\,K\in{\mathcal{T}}_{h}\,,\,\,i=1,\dots,M\,,$
(3.4)

where $\kappa$ is a fixed constant independent of $i$, $K$, and $h$. Next, the
time interval $[0,T]$ is decomposed by $0=t_{0}<t_{1}<\dots<t_{N}=T$ and we
set $\tau_{n}=t_{n}-t_{n-1}$, $n=1,\dots,N$. At each time level $t_{n}$, the
solution of (2.1)–(2.5) will be approximated by functions
$u_{h}^{n},c_{h}^{n},p_{h}^{n}\in V_{h}$. These functions can be identified
with coefficient vectors $\mathbf{u}^{n}=(u_{j}^{n})_{j=1}^{M}$,
$\mathbf{c}^{n}=(c_{j}^{n})_{j=1}^{M}$,
$\mathbf{p}^{n}=(p_{j}^{n})_{j=1}^{M}$, respectively, satisfying
$u_{h}^{n}=\pi_{h}\mathbf{u}^{n}$, $c_{h}^{n}=\pi_{h}\mathbf{c}^{n}$,
$p_{h}^{n}=\pi_{h}\mathbf{p}^{n}$. Note that $u_{h}^{n}(x_{i})=u_{i}^{n}$,
$c_{h}^{n}(x_{i})=c_{i}^{n}$, and $p_{h}^{n}(x_{i})=p_{i}^{n}$ for
$i=1,\dots,M$. We set

$u^{0}_{i}=u^{0}(x_{i})\,,\qquad c^{0}_{i}=c^{0}(x_{i})\,,\qquad
p^{0}_{i}=p^{0}(x_{i})\,,\qquad\quad i=1,\dots,M\,.$ (3.5)

Using linear interpolation with respect to time between the time levels gives
functions $u_{h,\tau}$, $c_{h,\tau}$, $p_{h,\tau}$ defined on
$\overline{\Omega}\times[0,T]$. For example, $u_{h,\tau}$ satisfies

$u_{h,\tau}(x,t)=\frac{1}{\tau_{n+1}}\left[u_{h}^{n+1}(x)(t-t_{n})+u_{h}^{n}(x)(t_{n+1}-t)\right]\qquad\forall\,\,x\in\overline{\Omega}\,,\,\,t\in[t_{n},t_{n+1}]\,,\,\,n=0,\dots,N-1\,,$

or, equivalently,

$u_{h,\tau}(x_{i},t)=\frac{1}{\tau_{n+1}}\left[u_{i}^{n+1}(t-t_{n})+u_{i}^{n}(t_{n+1}-t)\right]\qquad\forall\,\,i=1,\dots,M\,,\,\,t\in[t_{n},t_{n+1}]\,,\,\,n=0,\dots,N-1\,.$

Replacing the space $H^{1}(\Omega)$ in (3.1) by $V_{h}$ and applying the
$\theta$-method for discretization in time (with $\theta\in[0,1]$), one
obtains

$\displaystyle\left(\frac{u_{h}^{n+1}-u_{h}^{n}}{\tau_{n+1}},v_{h}\right)=\theta\,\mu\,\big{(}u_{h}^{n+1}(1-u_{h}^{n+1}),v_{h}\big{)}+\theta\,\chi\,\big{(}u_{h}^{n+1}\nabla
c_{h}^{n+1},\nabla v_{h}\big{)}$ $\displaystyle\hskip
28.45274pt+(1-\theta)\,\mu\,\big{(}u_{h}^{n}(1-u_{h}^{n}),v_{h}\big{)}+(1-\theta)\,\chi\,\big{(}u_{h}^{n}\nabla
c_{h}^{n},\nabla v_{h}\big{)}\qquad\forall\,\,v_{h}\in
V_{h}\,,\,\,n=0,\dots,N-1\,.$ (3.6)

Defining the matrices $\mathbb{M}=(m_{ij})_{i,j=1}^{M}$ and
${\mathbb{A}}^{n}=(a^{n}_{ij})_{i,j=1}^{M}$ with

$m_{ij}=(\phi_{j},\phi_{i})\,,\qquad
a^{n}_{ij}=-\mu\,\big{(}\phi_{j}(1-u_{h}^{n}),\phi_{i}\big{)}-\chi\,\big{(}\phi_{j}\nabla
c_{h}^{n},\nabla\phi_{i}\big{)}\,,$

the discrete variational problem (3.6) can be written in the matrix form

$({\mathbb{M}}+\theta\,\tau_{n+1}\,{\mathbb{A}}^{n+1})\,\mathbf{u}^{n+1}=({\mathbb{M}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{A}}^{n})\,\mathbf{u}^{n}\,,\qquad
n=0,\dots,N-1\,.$ (3.7)

The relations (3.2) and (3.3) suggest to define the coefficients of
$c_{h}^{n}$ and $p_{h}^{n}$ by

$\displaystyle
c_{i}^{n}=c^{0}(x_{i})\,{\mathrm{e}}^{-\int_{0}^{t_{n}}\,p_{h,\tau}(x_{i},s)\,{\mathrm{d}}s}\,,$
$\displaystyle i=1,\dots,M\,,\,\,n=0,\dots,N\,,$ (3.8) $\displaystyle
p_{i}^{n}={\mathrm{e}}^{-t_{n}/\epsilon}\left[p^{0}(x_{i})+\frac{1}{\epsilon}\,\int_{0}^{t_{n}}\,u_{h,\tau}(x_{i},s)\,c_{h,\tau}(x_{i},s)\,{\mathrm{e}}^{s/\epsilon}\,{\mathrm{d}}s\right],\qquad$
$\displaystyle i=1,\dots,M\,,\,\,n=0,\dots,N\,.$ (3.9)

Then, for $i=1,\dots,M$ and $n=0,\dots,N-1$, one has

$\displaystyle
c_{i}^{n+1}=c_{i}^{n}\,{\mathrm{e}}^{-\int_{t_{n}}^{t_{n+1}}\,p_{h,\tau}(x_{i},s)\,{\mathrm{d}}s}\,,$
(3.10) $\displaystyle
p_{i}^{n+1}={\mathrm{e}}^{-\tau_{n+1}/\epsilon}\,p^{n}_{i}+\frac{1}{\epsilon}\,{\mathrm{e}}^{-t_{n+1}/\epsilon}\int_{t_{n}}^{t_{n+1}}\,u_{h,\tau}(x_{i},s)\,c_{h,\tau}(x_{i},s)\,{\mathrm{e}}^{s/\epsilon}\,{\mathrm{d}}s\,.$
(3.11)

A direct computation gives

$\displaystyle c_{i}^{n+1}$ $\displaystyle=$
$\displaystyle\,c_{i}^{n}\,{\mathrm{e}}^{-\tau_{n+1}\,(p^{n+1}_{i}+p^{n}_{i})/2}\,$
, (3.12) $\displaystyle p_{i}^{n+1}$ $\displaystyle=$
$\displaystyle\,\,{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\,p^{n}_{i}+\frac{1}{\tau_{n+1}^{2}}\Bigg{\\{}$
$\displaystyle\Big{(}u^{n+1}_{i}\,(\epsilon-\tau_{n+1})-u^{n}_{i}\,\epsilon\Big{)}\Big{(}c^{n+1}_{i}\,(\epsilon-\tau_{n+1})-c^{n}_{i}\,\epsilon\Big{)}$
$\displaystyle-\,$
$\displaystyle\Big{(}u^{n+1}_{i}\,\epsilon-u^{n}_{i}\,(\epsilon+\tau_{n+1})\Big{)}\Big{(}c^{n+1}_{i}\,\epsilon-c^{n}_{i}\,(\epsilon+\tau_{n+1})\Big{)}\,{\mathrm{e}}^{-\tau_{n+1}/\epsilon}$
$\displaystyle+\,$
$\displaystyle(u^{n+1}_{i}-u^{n}_{i})(c^{n+1}_{i}-c^{n}_{i})\,\epsilon^{2}\left(1-{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\right)\Bigg{\\}},$
(3.13)

for $i=1,\dots,M$ and $n=0,\dots,N-1$. Note that the effects described by the
model (2.1)–(2.5), such as chemotaxis, strongly rely on the nonlinear coupling
terms. Therefore, all nonlinearities are treated implicitly in the discrete
problem (3.7)–(3.9).

To compute a solution of the nonlinear problem (3.7)–(3.9) at time $t_{n+1}$
(assuming that the solution vectors $\mathbf{u}^{n}$, $\mathbf{c}^{n}$, and
$\mathbf{p}^{n}$ at the previous time instant $t_{n}$ are known), we apply
simple fixed-point iterations leading to sequences
$\mathbf{u}^{n+1}_{k}=(u_{j,k}^{n+1})_{j=1}^{M}$,
$\mathbf{c}^{n+1}_{k}=(c_{j,k}^{n+1})_{j=1}^{M}$, and
$\mathbf{p}^{n+1}_{k}=(p_{j,k}^{n+1})_{j=1}^{M}$. We set
$\mathbf{u}^{n+1}_{0}=\mathbf{u}^{n}$, $\mathbf{c}^{n+1}_{0}=\mathbf{c}^{n}$,
$\mathbf{p}^{n+1}_{0}=\mathbf{p}^{n}$ and then, for $k>0$ and $i=1,\dots,M$,
we define

$\displaystyle c_{i,k}^{n+1}$ $\displaystyle=$
$\displaystyle\,\,c_{i}^{n}\,{\mathrm{e}}^{-\tau_{n+1}\,(p^{n+1}_{i,k-1}+p^{n}_{i})/2}$
$\displaystyle\,,$ (3.14) $\displaystyle p_{i,k}^{n+1}$ $\displaystyle=$
$\displaystyle{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\,p^{n}_{i}+\frac{1}{\tau_{n+1}^{2}}\Bigg{\\{}$
$\displaystyle\Big{(}u^{n+1}_{i,k-1}\,(\epsilon-\tau_{n+1})-u^{n}_{i}\,\epsilon\Big{)}\Big{(}c^{n+1}_{i,k}\,(\epsilon-\tau_{n+1})-c^{n}_{i}\,\epsilon\Big{)}$
$\displaystyle-\,$
$\displaystyle\Big{(}u^{n+1}_{i,k-1}\,\epsilon-u^{n}_{i}\,(\epsilon+\tau_{n+1})\Big{)}\Big{(}c^{n+1}_{i,k}\,\epsilon-c^{n}_{i}\,(\epsilon+\tau_{n+1})\Big{)}\,{\mathrm{e}}^{-\tau_{n+1}/\epsilon}$
$\displaystyle+\,$
$\displaystyle(u^{n+1}_{i,k-1}-u^{n}_{i})(c^{n+1}_{i,k}-c^{n}_{i})\,\epsilon^{2}\left(1-{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\right)\Bigg{\\}}.$
(3.15)

The iterate $\mathbf{u}^{n+1}_{k}$ is computed by solving the linear system

$({\mathbb{M}}+\theta\,\tau_{n+1}\,{\mathbb{A}}^{n+1}_{k-1})\,\mathbf{u}^{n+1}_{k}=({\mathbb{M}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{A}}^{n})\,\mathbf{u}^{n}\,,$
(3.16)

where the matrix ${\mathbb{A}}^{n+1}_{k-1}$ is defined by

${\mathbb{A}}^{n+1}_{k-1}=\Big{(}-\mu\,\big{(}\phi_{j}(1-u_{h,k-1}^{n+1}),\phi_{i}\big{)}-\chi\,\big{(}\phi_{j}\nabla
c_{h,k}^{n+1},\nabla\phi_{i}\big{)}\Big{)}_{i,j=1}^{M}$ (3.17)

and $u_{h,k-1}^{n+1}=\pi_{h}\mathbf{u}^{n+1}_{k-1}$ and
$c_{h,k}^{n+1}=\pi_{h}\mathbf{c}^{n+1}_{k}$ are the finite element functions
corresponding to the coefficient vectors $\mathbf{u}^{n+1}_{k-1}$ and
$\mathbf{c}^{n+1}_{k}$, respectively.

The linear system (3.16) has the form

${\mathbb{B}}\,\mathbf{u}^{n+1}={\mathbb{K}}\,\mathbf{u}^{n}$ (3.18)

and it is desirable that this system is positivity preserving, i.e., that
$\mathbf{u}^{n+1}\geq 0$ if $\mathbf{u}^{n}\geq 0$. A necessary and sufficient
condition for this property is ${\mathbb{B}}^{-1}\,{\mathbb{K}}\geq 0$ but
this condition is difficult to verify. Sufficient conditions are formulated in
the following lemma. Note that, throughout the paper, an inequality of the
type $\mathbf{u}^{n}\geq 0$ means that the inequality holds for each component
of the vector $\mathbf{u}^{n}$. Similarly, the fact that all entries of a
matrix $\mathbb{K}$ are non-negative is expressed by ${\mathbb{K}}\geq 0$.

######  Lemma 3.1.

Let the matrices ${\mathbb{B}}=(b_{ij})_{i,j=1}^{M}$ and
${\mathbb{K}}=(k_{ij})_{i,j=1}^{M}$ satisfy

$b_{ii}\geq 0\,,\quad k_{ii}\geq 0\,,\quad b_{ij}\leqslant 0\,,\quad
k_{ij}\geqslant 0,\qquad\forall\,\,i,j=1,\dots,M\,,\,\,i\neq j\,,$

and let $\mathbb{B}$ be a strictly diagonally dominant or an irreducibly
diagonally dominant matrix. Then $\mathbb{B}$ is an M-matrix and the scheme
(3.18) is positivity preserving.

  * Proof.

According to [53, Theorem 3.27], $\mathbb{B}$ is an M-matrix. Thus,
${\mathbb{B}}^{-1}\geq 0$ and hence also ${\mathbb{B}}^{-1}\,{\mathbb{K}}\geq
0$, which implies the result. ∎

In general, the linear system (3.16) originating from a standard Galerkin
discretization does not satisfy the above constraints because the mass matrix
is non-negative and the stiffness matrix may contain positive off-diagonal
entries. Our numerical results in Section 6 show that indeed the concentration
$\mathbf{u}$ may become negative in some parts of the computational domain
$\Omega$.

##  4 FCT stabilization

As we will see in Section 6, the magnitude of the solutions gradients can be
extremely large in some regions. The solution of the Galerkin discretization
from the previous section may become negative especially in these regions. As
a remedy, in the following we will modify the Galerkin discretization to
guarantee a positivity preservation property. As shown by Kuzmin [33, 34, 32],
this property can be readily enforced at the discrete level using a
conservative manipulation of the mass and stiffness matrices. The former will
be approximated by its diagonal counterpart $\mathbb{M}_{\mathrm{L}}$
constructed using row-sum mass lumping, whereas the latter will be modified by
adding an artificial diffusion matrix. To limit the amount of the artificial
diffusion, the FEM-FCT approach will be applied following [34].

Since the methods considered in this section guarantee that the approximate
solutions are non-negative, it is possible to replace the matrix
${\mathbb{A}}^{n}$ from the previous section by
$\tilde{\mathbb{A}}^{n}=(\tilde{a}^{n}_{ij})_{i,j=1}^{M}$ with

$\tilde{a}^{n}_{ij}=-\mu\,\big{(}\phi_{j}(1-|u_{h}^{n}|),\phi_{i}\big{)}-\chi\,\big{(}\phi_{j}\nabla
c_{h}^{n},\nabla\phi_{i}\big{)}\,.$

The matrix $\tilde{\mathbb{A}}^{n}$ is more suitable for theoretical
considerations than the matrix ${\mathbb{A}}^{n}$. However, a non-negative
approximate solution $u^{n}_{h}$, $c^{n}_{h}$, $p^{n}_{h}$ satisfying a
discrete problem based on the matrix $\tilde{\mathbb{A}}^{n}$ will satisfy
also the corresponding discrete problem with the original matrix
${\mathbb{A}}^{n}$.

Using the matrix $\tilde{\mathbb{A}}^{n}$, we introduce a symmetric artificial
diffusion matrix ${\mathbb{D}}^{n}=(d^{n}_{ij})_{i,j=1}^{M}$ defined by

$d^{n}_{ij}=-\max\\{\tilde{a}^{n}_{ij},0,\tilde{a}^{n}_{ji}\\}\quad\mbox{for}\
i\neq j\,,\qquad\quad d^{n}_{ii}=-\sum_{j=1,j\neq i}^{M}d^{n}_{ij}\,,$

and we set ${\mathbb{L}}^{n}=\tilde{\mathbb{A}}^{n}+{\mathbb{D}}^{n}$. Note
that ${\mathbb{L}}^{n}=(l^{n}_{ij})_{i,j=1}^{M}$ is a Z-matrix (i.e., it has
non-positive off-diagonal entries). Furthermore, we introduce the lumped mass
matrix $\mathbb{M}_{\mathrm{L}}=\mbox{diag}(m_{1},\dots,m_{M})$ with

$m_{i}=\sum_{j=1}^{M}\,m_{ij}\,,\qquad i=1,\dots,M\,.$

Now, the simplest way to enforce the positivity preservation is to consider
the so-called low-order method corresponding to the so-called high-order
method (3.7) which is defined by

$(\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1})\,\mathbf{u}^{n+1}=(\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n})\,\mathbf{u}^{n}\,,\qquad
n=0,\dots,N-1\,.$ (4.1)

Note that the matrix ${\mathbb{L}}^{n+1}$ depends on $\mathbf{u}^{n+1}$ and
$\mathbf{c}^{n+1}$ so that the low-order problem is again nonlinear. In
contrast to the Galerkin discretization (3.7), it is now possible to assure
the positivity preservation for sufficiently small time steps.

######  Lemma 4.1.

Let the time step $\tau_{n+1}$ satisfy the conditions

$(1-\theta)\,\tau_{n+1}\,l^{n}_{ii}\leq
m_{i}\,,\qquad\theta\,\tau_{n+1}\,\Big{(}\mu\,m_{i}+\chi\,\big{(}\nabla
c_{h}^{n+1},\nabla\phi_{i}\big{)}\Big{)}<m_{i}\,,\qquad\quad i=1,\dots,M\,.$
(4.2)

Then the matrix
$\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n}$ has non-
negative entries and
$\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1}$ is an
M-matrix.

  * Proof.

The first condition in (4.2) implies that
$\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n}$ has non-
negative diagonal entries. The off-diagonal entries of this matrix are non-
negative as well, since $\mathbb{M}_{\mathrm{L}}$ is diagonal and
${\mathbb{L}}^{n}$ is a Z-matrix.

Denoting
$\mathbb{B}=\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1}$,
one has for any $i\in\\{1,\dots,M\\}$

$\sum_{j=1}^{M}\,b_{ij}=m_{i}+\theta\,\tau_{n+1}\,\sum_{j=1}^{M}\,\tilde{a}^{n+1}_{ij}=m_{i}-\theta\,\tau_{n+1}\,\Big{(}\mu\,\big{(}1-|u_{h}^{n+1}|,\phi_{i}\big{)}+\chi\,\big{(}\nabla{c}_{h}^{n+1},\nabla\phi_{i}\big{)}\Big{)}\,,$

where we used the fact that $\sum_{j=1}^{M}\,\phi_{j}=1$. Since
$(1,\phi_{i})=m_{i}$, it follows from the second condition in (4.2) that
$\sum_{j=1}^{M}\,b_{ij}>0$. Thus, $b_{ii}>\sum_{j\neq i}\,|b_{ij}|$, i.e.,
$\mathbb{B}$ is strictly diagonally dominant and hence non-singular. Moreover,
$\mathbb{B}$ is a matrix of non-negative type and hence it is an M-matrix
(see, e.g., [6, Corollary 3.13]). ∎

######  Corollary 4.2.

Let the time step $\tau_{n+1}$ satisfy the conditions (4.2). Then the low-
order scheme (4.1) is positivity preserving, i.e.,

$\mathbf{u}^{n}\geqslant 0\quad\Rightarrow\quad\mathbf{u}^{n+1}\geqslant 0\,.$
(4.3)

  * Proof.

According to Lemma 4.1, the matrix
$\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1}$ is non-
singular,
$(\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1})^{-1}\geq 0$,
and $\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n}\geq 0$,
which immediately implies (4.3). ∎

######  Remark 4.3.

The second condition in (4.2) involves ${c}^{n+1}_{h}$ which implicitly
depends on $\tau_{n+1}$ through $\mathbf{p}^{n+1}$ and hence also through
$\mathbf{u}^{n+1}$. Therefore, it is desirable to replace this condition by a
condition independent of $c^{n+1}_{h}$. This is possible since we will show
that the values of ${c}^{n+1}_{h}$ are in the interval $[0,1]$. Then,
employing (3.4), one gets

$\big{(}\nabla{c}_{h}^{n+1},\nabla\phi_{i}\big{)}=\sum_{j=1}^{M}\,{c}^{n+1}_{j}\,(\nabla\phi_{j},\nabla\phi_{i})\leq\sum_{K\ni
x_{i}}\,\sum_{j=1}^{M}\,\|\nabla\phi_{j}\|_{L^{2}(K)}\,\|\nabla\phi_{i}\|_{L^{2}(K)}\leq
n_{\mathrm{v}}\,\kappa^{2}\sum_{K\ni x_{i}}h_{K}^{d-2}\,,$

where $n_{\mathrm{v}}$ is the number of vertices of a cell in
${\mathcal{T}}_{h}$ ($n_{\mathrm{v}}=d+1$ for simplices, $n_{\mathrm{v}}=4$
for quadrilaterals, and $n_{\mathrm{v}}=8$ for hexahedra). Thus, if the time
step $\tau_{n+1}$ satisfies

$\theta\,\tau_{n+1}\left(\mu\,m_{i}+\chi\,n_{\mathrm{v}}\,\kappa^{2}\sum_{K\ni
x_{i}}h_{K}^{d-2}\right)<m_{i}\,,\qquad i=1,\dots,M\,,$ (4.4)

and ${c}_{h}^{n+1}\in[0,1]$, then the second condition in (4.2) holds. Note
that (4.4) may be significantly more restrictive than (4.2).

To prove that the low-order discretization consisting of the equations (4.1),
(3.12), and (3.13) has a solution, we shall use the following consequence of
Brouwer’s fixed-point theorem.

######  Lemma 4.4.

Let $X$ be a finite-dimensional Hilbert space with inner product
$(\cdot,\cdot)_{X}$ and norm $\|\cdot\|_{X}$. Let $P:X\to X$ be a continuous
mapping and $K>0$ a real number such that $(Px,x)_{X}>0$ for any $x\in X$ with
$\|x\|_{X}=K$. Then there exists $x\in X$ such that $\|x\|_{X}<K$ and $Px=0$.

  * Proof.

See [50, p. 164, Lemma 1.4].∎

######  Theorem 4.5.

Consider any $n\in\\{0,\dots,N-1\\}$ and let
$\mathbf{u}^{n},\mathbf{c}^{n},\mathbf{p}^{n}\in{\mathbb{R}}^{M}$ satisfy
$\mathbf{u}^{n}\geq 0$, $1\geq\mathbf{c}^{n}\geq 0$, $\mathbf{p}^{n}\geq 0$.
Let the time step $\tau_{n+1}$ satisfy the conditions

$(1-\theta)\,\tau_{n+1}\,l^{n}_{ii}\leq
m_{i}\,,\qquad\theta\,\tau_{n+1}\left(\mu\,m_{i}+\chi\,n_{\mathrm{v}}\,\kappa^{2}\sum_{K\ni
x_{i}}h_{K}^{d-2}\right)<m_{i}\,,\qquad i=1,\dots,M\,.$ (4.5)

Then there exist vectors
$\mathbf{u}^{n+1},\mathbf{c}^{n+1},\mathbf{p}^{n+1}\in{\mathbb{R}}^{M}$
satisfying (4.1), (3.12), (3.13) and $\mathbf{u}^{n+1}\geq 0$,
$1\geq\mathbf{c}^{n+1}\geq 0$, $\mathbf{p}^{n+1}\geq 0$.

  * Proof.

To get rid of the exponential dependence on $\mathbf{p}^{n+1}$ when estimating
the nonlinear terms in (4.1), we replace (3.12) by

$c_{i}^{n+1}=c_{i}^{n}\,{\mathrm{e}}^{-\tau_{n+1}\,(|p^{n+1}_{i}|+p^{n}_{i})/2}\,,\qquad
i=1,\dots,M\,.$ (4.6)

At the end of the proof, we will show that $\mathbf{p}^{n+1}\geq 0$ so that
the original relation (3.12) will be recovered.

For $u,p\in{\mathbb{R}}$ and $i=1,\dots,M$, we introduce the notation

$C_{i}(p)=c_{i}^{n}\,{\mathrm{e}}^{-\tau_{n+1}\,(|p|+p^{n}_{i})/2}$ (4.7)

and

$\displaystyle P_{i}(u,p)=\,\,\frac{1}{\tau_{n+1}^{2}}$
$\displaystyle\,u\,C_{i}(p)\,\Bigg{\\{}(\epsilon-\tau_{n+1})^{2}+\epsilon^{2}\left(1-2\,{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\right)\Bigg{\\}}$
$\displaystyle+\frac{\epsilon}{\tau_{n+1}^{2}}$
$\displaystyle\,u\,c^{n}_{i}\,\Bigg{\\{}\tau_{n+1}\left(1+{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\right)-2\,\epsilon\left(1-{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\right)\Bigg{\\}}$
$\displaystyle+\frac{\epsilon}{\tau_{n+1}^{2}}$
$\displaystyle\,u^{n}_{i}\,C_{i}(p)\,\Bigg{\\{}\tau_{n+1}\left(1+{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\right)-2\,\epsilon\left(1-{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\right)\Bigg{\\}}$
$\displaystyle+\frac{1}{\tau_{n+1}^{2}}$
$\displaystyle\,u^{n}_{i}\,c^{n}_{i}\,\Bigg{\\{}\epsilon^{2}\left(2-{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\right)-(\epsilon+\tau_{n+1})^{2}\,{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\Bigg{\\}}+{\mathrm{e}}^{-\tau_{n+1}/\epsilon}\,p^{n}_{i}\,.$
(4.8)

Then, the validity of (4.6) and (3.13) is equivalent to

$c^{n+1}_{i}=C_{i}(p^{n+1}_{i})\,,\quad
p^{n+1}_{i}=P_{i}(u^{n+1}_{i},p^{n+1}_{i})\,,\qquad i=1,\dots,M\,.$

Note that

$|P_{i}(u,p)|\leq(|u|+u^{n}_{i})\,c^{n}_{i}\left(\frac{2\,\epsilon+\tau_{n+1}}{\tau_{n+1}}\right)^{2}+p^{n}_{i}\qquad\forall\,\,u,p\in{\mathbb{R}}\,,\,\,i=1,\dots,M\,.$
(4.9)

Furthermore, for $\mathbf{u},\mathbf{p}\in{\mathbb{R}}^{M}$ and
$i,j=1,\dots,M$, we denote

$\displaystyle
A_{ij}(\mathbf{u},\mathbf{p})=-\mu\,\big{(}\phi_{j}(1-|\pi_{h}\mathbf{u}|),\phi_{i}\big{)}-\chi\,\big{(}\phi_{j}\nabla(\pi_{h}\mathbf{C}(\mathbf{p})),\nabla\phi_{i}\big{)}\,,$
(4.10) $\displaystyle
D_{ij}(\mathbf{u},\mathbf{p})=-\max\\{A_{ij}(\mathbf{u},\mathbf{p}),0,A_{ji}(\mathbf{u},\mathbf{p})\\}\quad\mbox{for}\
i\neq j\,,\qquad\quad D_{ii}(\mathbf{u},\mathbf{p})=-\sum_{j=1,j\neq
i}^{M}D_{ij}(\mathbf{u},\mathbf{p})\,,$ (4.11) $\displaystyle
S_{i}(\mathbf{u},\mathbf{p})=m_{i}\,u_{i}+\theta\,\tau_{n+1}\,\sum_{j=1}^{M}\,(A_{ij}(\mathbf{u},\mathbf{p})+D_{ij}(\mathbf{u},\mathbf{p}))\,u_{j}-[(\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n})\,\mathbf{u}^{n}]_{i}\,,$
(4.12)

where $\mathbf{C}(\mathbf{p})=(C_{i}(p_{i}))_{i=1}^{M}$. Then (4.1) with
$c^{n+1}_{h}$ defined by (4.6) is equivalent to

$S_{i}(\mathbf{u}^{n+1},\mathbf{p}^{n+1})=0\,,\qquad i=1,\dots,M\,.$

Therefore, defining the operator $P:{\mathbb{R}}^{2M}\to{\mathbb{R}}^{2M}$ by

$P\,\mathbf{U}=(S_{1}(\mathbf{u},\mathbf{p}),\dots,S_{M}(\mathbf{u},\mathbf{p}),p_{1}-P_{1}(u_{1},p_{1}),\dots,p_{M}-P_{M}(u_{M},p_{M}))\qquad\forall\,\,\mathbf{U}=(\mathbf{u},\mathbf{p})\in{\mathbb{R}}^{2M}\,,$
(4.13)

the vectors $\mathbf{u}^{n+1}$, $\mathbf{c}^{n+1}$, $\mathbf{p}^{n+1}$ are a
solution of (4.1), (4.6), (3.13) if and only if
$\mathbf{U}=(\mathbf{u}^{n+1},\mathbf{p}^{n+1})$ satisfies $P\,\mathbf{U}=0$
and $\mathbf{c}^{n+1}=\mathbf{C}(\mathbf{p}^{n+1})$.

To show that the equation $P\,\mathbf{U}=0$ has a solution, we will verify the
assumptions of Lemma 4.4. Since it is obvious that the operator $P$ is
continuous, it suffices to investigate the product
$(P\,\mathbf{U},\mathbf{U})$, where $(\cdot,\cdot)$ is the Euclidean inner
product in ${\mathbb{R}}^{2M}$. We will denote the corresponding norm by
$\|\cdot\|$. The Euclidean norm in ${\mathbb{R}}^{M}$ will be denoted by
$\|\cdot\|_{M}$. Since the matrix
$(D_{ij}(\mathbf{u},\mathbf{p}))_{i,j=1}^{M}$ is symmetric and has zero row
sums and non-positive off-diagonal entries, one obtains

$\sum_{i,j=1}^{M}\,u_{i}\,D_{ij}(\mathbf{u},\mathbf{p})\,u_{j}=-\frac{1}{2}\,\sum_{i,j=1}^{M}\,D_{ij}(\mathbf{u},\mathbf{p})(u_{i}-u_{j})^{2}\geq
0\qquad\forall\,\,\mathbf{u},\mathbf{p}\in{\mathbb{R}}^{M}\,.$

Furthermore, since $0\leq\mathbf{C}(\mathbf{p})\leq 1$, the expressions
$\big{(}\phi_{j}\nabla(\pi_{h}\mathbf{C}(\mathbf{p})),\nabla\phi_{i}\big{)}$
can be bounded independently of $\mathbf{p}$. Therefore, using the equivalence
of norms on finite-dimensional spaces, one obtains

$\theta\,\tau_{n+1}\,\sum_{i,j=1}^{M}\,u_{i}\,(A_{ij}(\mathbf{u},\mathbf{p})+D_{ij}(\mathbf{u},\mathbf{p}))\,u_{j}\geq\theta\,\tau_{n+1}\,\mu\,\|\pi_{h}\mathbf{u}\|_{L^{3}(\Omega)}^{3}-C_{1}\,\|\mathbf{u}\|_{M}^{2}\geq
C_{2}\,\|\mathbf{u}\|_{M}^{3}-C_{1}\,\|\mathbf{u}\|_{M}^{2}\,,$

where $C_{1}$ and $C_{2}$ are positive constants independent of $\mathbf{u}$
and $\mathbf{p}$. Thus,

$\sum_{i=1}^{M}\,u_{i}\,S_{i}(\mathbf{u},\mathbf{p})\geq
C_{2}\,\|\mathbf{u}\|_{M}^{3}-C_{1}\,\|\mathbf{u}\|_{M}^{2}-C_{3}\,\|\mathbf{u}\|_{M}\,,$
(4.14)

where
$C_{3}=\|(\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n})\,\mathbf{u}^{n}\|_{M}$.
Finally, using (4.9), it follows that

$\sum_{i=1}^{M}\,p_{i}\,(p_{i}-P_{i}(u_{i},p_{i}))\geq\|\mathbf{p}\|_{M}^{2}-C_{4}\,\|\mathbf{p}\|_{M}\,\|\mathbf{u}\|_{M}-C_{5}\,\|\mathbf{p}\|_{M}$

with positive constants $C_{4}$ and $C_{5}$ independent of $\mathbf{u}$ and
$\mathbf{p}$. Applying the Young inequality, the previous two inequalities
imply that there exist positive constants $C_{6}$ and $C_{7}$ such that

$(P\,\mathbf{U},\mathbf{U})\geq\frac{1}{2}\,\|\mathbf{U}\|^{2}+C_{2}\,\|\mathbf{u}\|_{M}^{3}-C_{6}\,\|\mathbf{u}\|_{M}^{2}-C_{7}\geq\frac{1}{2}\,\|\mathbf{U}\|^{2}-\frac{C_{6}^{3}}{C_{2}^{2}}-C_{7}\qquad\forall\,\,\mathbf{U}=(\mathbf{u},\mathbf{p})\in{\mathbb{R}}^{2M}\,.$

Thus, for any $K>\sqrt{2\,C_{6}^{3}/C_{2}^{2}+2\,C_{7}}$, one has
$(P\,\mathbf{U},\mathbf{U})>0$ for any $\mathbf{U}\in{\mathbb{R}}^{2M}$ with
$\|\mathbf{U}\|=K$. Therefore, according to Lemma 4.4, there exists a solution
$\mathbf{U}$ of the equation $P\,\mathbf{U}=0$ and hence also a solution
$\mathbf{u}^{n+1}$, $\mathbf{c}^{n+1}$, $\mathbf{p}^{n+1}$ of (4.1), (4.6),
and (3.13).

It immediately follows from (4.6) that $0\leq\mathbf{c}^{n+1}\leq 1$. Thus,
according to Corollary 4.2 and Remark 4.3, the solution satisfies
$\mathbf{u}^{n+1}\geq 0$. Since (3.13) is equivalent to (3.11), one also has
$\mathbf{p}^{n+1}\geq 0$ and hence (3.12) is satisfied as well. ∎

Although the solution of (4.1), (3.12), (3.13) does not possess negative
values under the time step restrictions (4.2), it is usually very inaccurate
since too much artificial diffusion is introduced by the modifications leading
to the low-order method (4.1), cf. Section 6.3. Therefore, in the FEM-FCT
methodology, a correction term $\overline{\mathbf{f}}^{n+1}$ is added in such
a way that the method becomes less diffusive while negative values are still
excluded. This leads to an extension of (4.1) in the form

$(\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1})\,\mathbf{u}^{n+1}=(\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n})\,\mathbf{u}^{n}+\overline{\mathbf{f}}^{n+1}\,.$

The high-order method (3.7) (with ${\mathbb{A}}^{n}$ replaced by
$\tilde{\mathbb{A}}^{n}$) is recovered if

$\overline{\mathbf{f}}^{n+1}=(\mathbb{M}_{\mathrm{L}}-{\mathbb{M}})(\mathbf{u}^{n+1}-\mathbf{u}^{n})+\theta\,\tau_{n+1}\,{\mathbb{D}}^{n+1}\,\mathbf{u}^{n+1}+(1-\theta)\,\tau_{n+1}\,{\mathbb{D}}^{n}\,\mathbf{u}^{n}\,.$
(4.15)

Since ${\mathbb{D}}^{n}$ has zero row sums, one can write

$({\mathbb{D}}^{n}\,\mathbf{u}^{n})_{i}=\sum_{j=1}^{M}\,d^{n}_{ij}\,(u^{n}_{j}-u^{n}_{i})\,,\qquad
i=1,\dots,M\,.$

For the terms with the matrices ${\mathbb{D}}^{n+1}$ and
$\mathbb{M}_{\mathrm{L}}-{\mathbb{M}}$ (which also have zero row sums), one
can proceed analogously and hence (4.15) holds if an only if

$\overline{\mathbf{f}}^{n+1}=\Bigg{(}\sum_{j=1}^{M}\,f^{n+1}_{ij}\Bigg{)}_{i=1}^{M}\,,$

where the algebraic fluxes $f^{n+1}_{ij}$ are given by

$f^{n+1}_{ij}=-m_{ij}\,(u_{j}^{n+1}-u_{i}^{n+1})+m_{ij}\,(u^{n}_{j}-u^{n}_{i})+\theta\,\tau_{n+1}\,d^{n+1}_{ij}\,(u_{j}^{n+1}-u_{i}^{n+1})+(1-\theta)\,\tau_{n+1}\,d^{n}_{ij}\,(u^{n}_{j}-u^{n}_{i})\,.$
(4.16)

Because $\mathbb{M}$, ${\mathbb{D}}^{n+1}$, and ${\mathbb{D}}^{n}$ are
symmetric matrices, one has $f^{n+1}_{ij}=-f^{n+1}_{ji}$. Note also that the
fluxes depend on (unknown) values of the approximate solution at time level
$t_{n+1}$.

Now, the idea of the FCT approach is to limit the fluxes $f^{n+1}_{ij}$ by
solution dependent correction factors $\alpha^{n+1}_{ij}\in[0,1]$ called
limiters so that the non-negativity of the approximate solution can be
guaranteed but less artificial diffusion is introduced than in case of the
low-order method. This leads to the discrete problem

$(\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1})\,\mathbf{u}^{n+1}=(\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n})\,\mathbf{u}^{n}+\Bigg{(}\sum_{j=1}^{M}\,\alpha^{n+1}_{ij}\,f^{n+1}_{ij}\Bigg{)}_{i=1}^{M}\,.$
(4.17)

The original Galerkin discretization is recovered for $\alpha_{ij}=1$ while
the largest amount of artificial diffusion is introduced for $\alpha_{ij}=0$.
The latter setting is appropriate in the neighborhood of steep fronts and
large gradients. The artificial diffusion can be removed in regions where the
solution is smooth and where non-positive off-diagonal entries of the
stiffness matrix do not pose any threat to non-negativity. The corrected
fluxes depend on the approximate solution in a nonlinear way but since the
problem in here is already nonlinear, we can treat both nonlinearities
simultaneously.

It is convenient to write the nonlinear problem (4.17) in the form

$\displaystyle\mathbb{M}_{\mathrm{L}}\,\overline{\mathbf{u}}$ $\displaystyle=$
$\displaystyle(\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n})\,\mathbf{u}^{n}\,,$
(4.18) $\displaystyle\mathbb{M}_{\mathrm{L}}\,\tilde{\mathbf{u}}$
$\displaystyle=$
$\displaystyle\mathbb{M}_{\mathrm{L}}\,\overline{\mathbf{u}}+\Bigg{(}\sum_{j=1}^{M}\,\alpha^{n+1}_{ij}\,f^{n+1}_{ij}\Bigg{)}_{i=1}^{M}\,,$
(4.19)
$\displaystyle(\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1})\,\mathbf{u}^{n+1}$
$\displaystyle=$ $\displaystyle\mathbb{M}_{\mathrm{L}}\,\tilde{\mathbf{u}}\,.$
(4.20)

According to Lemma 4.1, the steps (4.18) and (4.20) are positivity preserving
under the conditions (4.2). To guarantee the positivity preservation of the
second step, the limiters $\alpha^{n+1}_{ij}$ have to be defined
appropriately. We will apply the Zalesak algorithm [55] which will be
described next.

The solution of the nonlinear problem (4.18)–(4.20) is computed by fixed-point
iterations where the algebraic fluxes are calculated using the previous
iterate. Since the properties of the Zalesak algorithm do not depend on the
form of these fluxes, we will denote them simply by $f_{ij}$. Then, the aim is
to find limiters $\alpha_{ij}\in[0,1]$ such that the solution
$\tilde{\mathbf{u}}$ of

$\mathbb{M}_{\mathrm{L}}\,\tilde{\mathbf{u}}=\mathbb{M}_{\mathrm{L}}\,\overline{\mathbf{u}}+\Bigg{(}\sum_{j=1}^{M}\,\alpha_{ij}\,f_{ij}\Bigg{)}_{i=1}^{M}$

satisfies

$\overline{u}_{i}^{\mathrm{min}}\leq\tilde{u}_{i}\leq\overline{u}_{i}^{\mathrm{max}}\,,\qquad
i=1,\dots,M\,,$ (4.21)

where

$\overline{u}_{i}^{\mathrm{min}}=\min_{j\in{\mathcal{N}}_{i}\cup\\{i\\}}\overline{u}_{j}\,,\quad\overline{u}_{i}^{\mathrm{max}}=\max_{j\in{\mathcal{N}}_{i}\cup\\{i\\}}\overline{u}_{j}\,,\qquad
i=1,\dots,M\,,$

and ${\mathcal{N}}_{i}$ is the index set of neighbour vertices to the vertex
$x_{i}$ (note that two vertices of the triangulation ${\mathcal{T}}_{h}$ are
called neighbouring if they are contained in the same mesh cell). To preserve
conservativity, it is important that the limiters $\alpha_{ij}$ form a
symmetric matrix. The limiting process begins with cancelling all fluxes that
are diffusive in nature and tend to flatten the solution profiles, cf. [34].
The required modification is

$f_{ij}:=0\qquad\mbox{if}\,\,\,\,f_{ij}\,(\overline{u}_{j}-\overline{u}_{i})>0\,.$
(4.22)

The remaining fluxes are truly antidiffusive and the computation of
$\alpha_{ij}$ involves the following steps:

  * 1.

Compute the sum of positive/negative antidiffusive fluxes into node $i$

$P_{i}^{+}=\sum_{j\in{\mathcal{N}}_{i}}\,\max\\{0,{f}_{ij}\\}\,,\qquad
P_{i}^{-}=\sum_{j\in{\mathcal{N}}_{i}}\,\min\\{0,{f}_{ij}\\}\,.$ (4.23)

  * 2.

Compute the distance to a local extremum of the auxiliary solution
$\overline{\mathbf{u}}$

$Q_{i}^{+}=m_{i}\,(\overline{u}_{i}^{\mathrm{max}}-\overline{u}_{i})\,,\qquad
Q_{i}^{-}=m_{i}\,(\overline{u}_{i}^{\mathrm{min}}-\overline{u}_{i})\,.$ (4.24)

  * 3.

Compute the nodal correction factors for the net increment to node $i$

$R_{i}^{+}=\min\left\\{1,\dfrac{Q_{i}^{+}}{P_{i}^{+}}\right\\},\qquad
R_{i}^{-}=\min\left\\{1,\dfrac{Q_{i}^{-}}{P_{i}^{-}}\right\\}.$ (4.25)

If a denominator is zero, set the respective value of $R_{i}^{+}$ or
$R_{i}^{-}$ equal to $1$.

  * 4.

Check the sign of the antidiffusive flux and define the correction factor by

$\alpha_{ij}=\left\\{\begin{array}[]{ll}\min\\{R_{i}^{+},R_{j}^{-}\\}&\mbox{if}\quad\,f_{ij}>0\,,\\\\[2.84526pt]
1&\mbox{if}\quad\,f_{ij}=0\,,\\\\[2.84526pt]
\min\\{R_{i}^{-},R_{j}^{+}\\}&\mbox{if}\quad\,f_{ij}<0\,.\end{array}\right.$
(4.26)

It can be easily verified (see, e.g., [6]) that this algorithm leads to the
property (4.21).

Now we are in a position to prove the solvability and positivity preservation
for the above FCT discretization.

######  Theorem 4.6.

Consider any $n\in\\{0,\dots,N-1\\}$ and let
$\mathbf{u}^{n},\mathbf{c}^{n},\mathbf{p}^{n}\in{\mathbb{R}}^{M}$ satisfy
$\mathbf{u}^{n}\geq 0$, $1\geq\mathbf{c}^{n}\geq 0$, $\mathbf{p}^{n}\geq 0$.
Let the time step $\tau_{n+1}$ satisfy the conditions (4.5). Then there exist
vectors
$\mathbf{u}^{n+1},\mathbf{c}^{n+1},\mathbf{p}^{n+1}\in{\mathbb{R}}^{M}$
satisfying (4.17), (3.12), (3.13) where the fluxes $f^{n+1}_{ij}$ are given by
(4.16) and (4.22) and the limiters $\alpha^{n+1}_{ij}$ are computed using the
Zalesak algorithm (4.23)–(4.26) from the fluxes $f^{n+1}_{ij}$. Moreover,
these vectors satisfy $\mathbf{u}^{n+1}\geq 0$, $1\geq\mathbf{c}^{n+1}\geq 0$,
and $\mathbf{p}^{n+1}\geq 0$.

  * Proof.

The proof follows the lines of that of Theorem 4.5. Thus, we again start with
replacing (3.12) by (4.6). We again define $C_{i}$, $P_{i}$, $A_{ij}$, and
$D_{ij}$ by (4.7), (4.8), (4.10) and (4.11), respectively, whereas $S_{i}$ are
now defined by

$\displaystyle
S_{i}(\mathbf{u},\mathbf{p})=m_{i}\,u_{i}+\theta\,\tau_{n+1}\,\sum_{j=1}^{M}\,(A_{ij}(\mathbf{u},\mathbf{p})+D_{ij}(\mathbf{u},\mathbf{p}))\,u_{j}-\sum_{j=1}^{M}\,\alpha_{ij}(\mathbf{u},\mathbf{p})\,\tilde{f}_{ij}(\mathbf{u},\mathbf{p})$
$\displaystyle\hskip
241.84842pt-[(\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n})\,\mathbf{u}^{n}]_{i}\,,$

where $\alpha_{ij}(\mathbf{u},\mathbf{p})$ are defined by the Zalesak
algorithm (4.23)–(4.26) for the algebraic fluxes
$\tilde{f}_{ij}(\mathbf{u},\mathbf{p})$ defined by

$\tilde{f}_{ij}(\mathbf{u},\mathbf{p})=\left\\{\begin{array}[]{ll}f_{ij}(\mathbf{u},\mathbf{p})&\mbox{if}\quad\,f_{ij}(\mathbf{u},\mathbf{p})(\overline{u}_{j}-\overline{u}_{i})\leq
0\,,\\\\[2.84526pt]
0&\mbox{if}\quad\,f_{ij}(\mathbf{u},\mathbf{p})(\overline{u}_{j}-\overline{u}_{i})>0\,,\end{array}\right.$

with $\overline{\mathbf{u}}$ from (4.18) and

$f_{ij}(\mathbf{u},\mathbf{p})=\big{(}-m_{ij}+\theta\,\tau_{n+1}\,D_{ij}(\mathbf{u},\mathbf{p})\big{)}(u_{j}-u_{i})+\big{(}m_{ij}+(1-\theta)\,\tau_{n+1}\,d^{n}_{ij}\big{)}(u^{n}_{j}-u^{n}_{i})\,.$

Then, defining the operator $P:{\mathbb{R}}^{2M}\to{\mathbb{R}}^{2M}$ by
(4.13), the vectors $\mathbf{u}^{n+1}$, $\mathbf{c}^{n+1}$, $\mathbf{p}^{n+1}$
are a solution of (4.17), (4.6), (3.13) if and only if
$\mathbf{U}=(\mathbf{u}^{n+1},\mathbf{p}^{n+1})$ satisfies $P\,\mathbf{U}=0$
and $\mathbf{c}^{n+1}=\mathbf{C}(\mathbf{p}^{n+1})$.

The solvability of the equation $P\,\mathbf{U}=0$ will be again proved using
Lemma 4.4. To show the continuity of the operator $P$ at any point
$\widetilde{\mathbf{U}}\equiv(\widetilde{\mathbf{u}},\widetilde{\mathbf{p}})\in{\mathbb{R}}^{2M}$,
it suffices to consider the terms
$\alpha_{ij}(\mathbf{u},\mathbf{p})\,\tilde{f}_{ij}(\mathbf{u},\mathbf{p})$
since the remaining terms in the definition of $P$ are clearly continuous.
Moreover, $f_{ij}$ and hence also $\tilde{f}_{ij}$ are continuous. Thus, if
$\tilde{f}_{ij}(\widetilde{\mathbf{U}})\neq 0$, then the denominators in the
formulas defining $\alpha_{ij}(\mathbf{U})$ with
$\mathbf{U}=(\mathbf{u},\mathbf{p})$ do not vanish in a neighborhood of
$\widetilde{\mathbf{U}}$ and hence $\alpha_{ij}$ is continuous at
$\widetilde{\mathbf{U}}$. Consequently, also $\alpha_{ij}\,\tilde{f}_{ij}$ is
continuous at $\widetilde{\mathbf{U}}$. If
$\tilde{f}_{ij}(\widetilde{\mathbf{U}})=0$, then

$|(\alpha_{ij}\,\tilde{f}_{ij})(\mathbf{U})-(\alpha_{ij}\,\tilde{f}_{ij})(\widetilde{\mathbf{U}})|=|(\alpha_{ij}\,\tilde{f}_{ij})(\mathbf{U})|\leq|\tilde{f}_{ij}(\mathbf{U})|=|\tilde{f}_{ij}(\mathbf{U})-\tilde{f}_{ij}(\widetilde{\mathbf{U}})|\,,$

which shows that $\alpha_{ij}\,\tilde{f}_{ij}$ is again continuous at
$\widetilde{\mathbf{U}}$.

To estimate $(P\,\mathbf{U},\mathbf{U})$ from below, let us denote

$\widetilde{\alpha}_{ij}(\mathbf{u},\mathbf{p})=\left\\{\begin{array}[]{ll}\alpha_{ij}(\mathbf{u},\mathbf{p})&\mbox{if}\quad\,f_{ij}(\mathbf{u},\mathbf{p})(\overline{u}_{j}-\overline{u}_{i})\leq
0\,,\\\\[2.84526pt]
0&\mbox{if}\quad\,f_{ij}(\mathbf{u},\mathbf{p})(\overline{u}_{j}-\overline{u}_{i})>0\,.\end{array}\right.$

Then $\widetilde{\alpha}_{ij}$ again form a symmetric matrix and
$\alpha_{ij}\,\tilde{f}_{ij}=\widetilde{\alpha}_{ij}\,f_{ij}$. Therefore,
$S_{i}(\mathbf{u},\mathbf{p})$ can be written in the form

$\displaystyle S_{i}(\mathbf{u},\mathbf{p})=$
$\displaystyle\sum_{j=1}^{M}\,m_{ij}\,u_{j}+\theta\,\tau_{n+1}\,\sum_{j=1}^{M}\,A_{ij}(\mathbf{u},\mathbf{p})\,u_{j}$
$\displaystyle+$
$\displaystyle\sum_{j=1}^{M}\,\big{(}1-\widetilde{\alpha}_{ij}(\mathbf{u},\mathbf{p})\big{)}\big{(}-m_{ij}+\theta\,\tau_{n+1}\,D_{ij}(\mathbf{u},\mathbf{p})\big{)}(u_{j}-u_{i})$
$\displaystyle+$
$\displaystyle\sum_{j=1}^{M}\,\big{(}1-\widetilde{\alpha}_{ij}(\mathbf{u},\mathbf{p})\big{)}\big{(}m_{ij}+(1-\theta)\,\tau_{n+1}\,d^{n}_{ij}\big{)}(u^{n}_{j}-u^{n}_{i})$
$\displaystyle-$
$\displaystyle[({\mathbb{M}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{A}}^{n})\,\mathbf{u}^{n}]_{i}\,.$

Denoting
$B_{ij}=\big{(}1-\widetilde{\alpha}_{ij}(\mathbf{u},\mathbf{p})\big{)}\big{(}-m_{ij}+\theta\,\tau_{n+1}\,D_{ij}(\mathbf{u},\mathbf{p})\big{)}$,
one has

$\sum_{i,j=1}^{M}\,u_{i}\,\big{(}1-\widetilde{\alpha}_{ij}(\mathbf{u},\mathbf{p})\big{)}\big{(}-m_{ij}+\theta\,\tau_{n+1}\,D_{ij}(\mathbf{u},\mathbf{p})\big{)}(u_{j}-u_{i})=-\frac{1}{2}\,\sum_{i,j=1}^{M}\,B_{ij}\,(u_{i}-u_{j})^{2}\geq
0\,,$

since the matrix $(B_{ij})_{i,j=1}^{M}$ is symmetric and has non-positive off-
diagonal entries. Therefore, one again obtains (4.14) where the constants
$C_{1}$, $C_{2}$ are the same as in the proof of Theorem 4.5 and

$C_{3}=\|\mathbf{g}\|_{M}+\|({\mathbb{M}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{A}}^{n})\,\mathbf{u}^{n}\|_{M}\,,$

where

$g_{i}=\sum_{j=1}^{M}\,|m_{ij}+(1-\theta)\,\tau_{n+1}\,d^{n}_{ij}|\,|u^{n}_{j}-u^{n}_{i}|\,,\qquad
i=1,\dots,M\,.$

Thus, in the same way as in the proof of Theorem 4.5, one concludes that there
exists a solution $\mathbf{U}$ of the equation $P\,\mathbf{U}=0$ and hence
also a solution $\mathbf{u}^{n+1}$, $\mathbf{c}^{n+1}$, $\mathbf{p}^{n+1}$ of
(4.17), (4.6), and (3.13).

To prove the positivity preservation, we write (4.17) in the form
(4.18)–(4.20). Since
$\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n}\geq 0$
according to Lemma 4.1, one has $\overline{\mathbf{u}}\geq 0$. Applying
(4.21), one gets $\tilde{\mathbf{u}}\geq 0$. Since $0\leq\mathbf{c}^{n+1}\leq
1$ due to (4.6), it follows from Lemma 4.1 and Remark 4.3 that the matrix
$\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1}$ is an
M-matrix. Consequently, $\mathbf{u}^{n+1}\geq 0$ in view of (4.20). Since
(3.13) is equivalent to (3.11), one also has $\mathbf{p}^{n+1}\geq 0$ and
hence (3.12) is satisfied as well. ∎

##  5 Iterative solution of the FCT discretization

To compute a solution of the nonlinear problem (4.17), (3.12), (3.13) at time
$t_{n+1}$, we will proceed similarly as for the Galerkin discretization in
Section 3. Thus, given approximations $\mathbf{u}^{n+1}_{k-1}$,
$\mathbf{c}^{n+1}_{k-1}$, $\mathbf{p}^{n+1}_{k-1}$ (with some $k>0$) of
$\mathbf{u}^{n+1}$, $\mathbf{c}^{n+1}$, $\mathbf{p}^{n+1}$, respectively, we
compute $\mathbf{c}^{n+1}_{k}$, $\mathbf{p}^{n+1}_{k}$ using (3.14), (3.15).
The iterate $\mathbf{u}^{n+1}_{k}$ is computed by solving the linear system

$(\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1}_{k-1})\,\mathbf{u}^{n+1}_{k}=(\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n})\,\mathbf{u}^{n}+\Bigg{(}\sum_{j=1}^{M}\,\alpha^{n+1}_{ij,k-1}\,f^{n+1}_{ij,k-1}\Bigg{)}_{i=1}^{M}\,,$
(5.1)

where
${\mathbb{L}}^{n+1}_{k-1}={\mathbb{A}}^{n+1}_{k-1}+{\mathbb{D}}^{n+1}_{k-1}$
with the matrix ${\mathbb{A}}^{n+1}_{k-1}$ defined in (3.17) and the
artificial diffusion matrix ${\mathbb{D}}^{n+1}_{k-1}$ defined by

$d^{n+1}_{ij,k-1}=-\max\\{a^{n+1}_{ij,k-1},0,a^{n+1}_{ji,k-1}\\}\quad\mbox{for}\
i\neq j\,,\qquad\quad d^{n+1}_{ii,k-1}=-\sum_{j=1,j\neq
i}^{M}d^{n+1}_{ij,k-1}\,.$ (5.2)

The algebraic fluxes $f^{n+1}_{ij,k-1}$ are given by

$f^{n+1}_{ij,k-1}=\big{(}-m_{ij}+\theta\,\tau_{n+1}\,d^{n+1}_{ij,k-1}\big{)}(u_{j,k-1}^{n+1}-u_{i,k-1}^{n+1})+\big{(}m_{ij}+(1-\theta)\,\tau_{n+1}\,d^{n}_{ij}\big{)}(u^{n}_{j}-u^{n}_{i})$
(5.3)

and we again consider the prelimiting step

$f^{n+1}_{ij,k-1}:=0\qquad\mbox{if}\,\,\,\,f^{n+1}_{ij,k-1}\,(\overline{u}_{j}-\overline{u}_{i})>0\,,$
(5.4)

with $\overline{\mathbf{u}}$ from (4.18). The limiters $\alpha^{n+1}_{ij,k-1}$
are computed from the fluxes $f^{n+1}_{ij,k-1}$ using the Zalesak algorithm
(4.23)–(4.26). The following result shows that, under suitable time step
restrictions, the above-defined iterates are uniquely determined and preserve
non-negativity. This is important since, in practice, the fixed-point
iterations are usually terminated when a stopping criterion is met, i.e.,
typically before reaching the solution of the nonlinear problem (4.17),
(3.12), (3.13).

######  Theorem 5.1.

Consider any $n\in\\{0,\dots,N-1\\}$ and $k\in{\mathbb{N}}$ and let
$\mathbf{u}^{n},\mathbf{c}^{n},\mathbf{p}^{n}\in{\mathbb{R}}^{M}$ and
$\mathbf{u}^{n+1}_{k-1},\mathbf{p}^{n+1}_{k-1}\in{\mathbb{R}}^{M}$ be
arbitrary vectors satisfying $\mathbf{u}^{n}\geq 0$, $1\geq\mathbf{c}^{n}\geq
0$, $\mathbf{p}^{n}\geq 0$, and $\mathbf{u}^{n+1}_{k-1}\geq 0$,
$\mathbf{p}^{n+1}_{k-1}\geq 0$. Let $\mathbf{c}^{n+1}_{k}$,
$\mathbf{p}^{n+1}_{k}$ be given by (3.14), (3.15). Let the time step
$\tau_{n+1}$ satisfy the conditions

$(1-\theta)\,\tau_{n+1}\,l^{n}_{ii}\leq
m_{i}\,,\qquad\theta\,\tau_{n+1}\Big{(}\mu\,\big{(}1-\pi_{h}\mathbf{u}^{n+1}_{k-1},\phi_{i}\big{)}+\chi\,\big{(}\nabla(\pi_{h}\mathbf{c}^{n+1}_{k}),\nabla\phi_{i}\big{)}\Big{)}<m_{i}\,,\qquad
i=1,\dots,M\,.$ (5.5)

Then the linear system (5.1) has a unique solution $\mathbf{u}^{n+1}_{k}$ and
one has $\mathbf{u}^{n+1}_{k}\geq 0$, $1\geq\mathbf{c}^{n+1}_{k}\geq 0$, and
$\mathbf{p}^{n+1}_{k}\geq 0$.

  * Proof.

The formula (3.14) immediately implies that $1\geq\mathbf{c}^{n+1}_{k}\geq 0$.
Since (3.15) can be written in the form (3.11) with $u_{h,\tau}$ and
$c_{h,\tau}$ defined using $\mathbf{u}^{n+1}_{k-1}$ and
$\mathbf{c}^{n+1}_{k}$, respectively, at time $t_{n+1}$, one has
$\mathbf{p}^{n+1}_{k}\geq 0$. Since
$\mathbb{M}_{\mathrm{L}}-(1-\theta)\,\tau_{n+1}\,{\mathbb{L}}^{n}\geq 0$
according to Lemma 4.1, the solution of (4.18) satisfies
$\overline{\mathbf{u}}\geq 0$. Then (4.21) implies $\tilde{\mathbf{u}}\geq 0$
for the solution of

$\mathbb{M}_{\mathrm{L}}\,\tilde{\mathbf{u}}=\mathbb{M}_{\mathrm{L}}\,\overline{\mathbf{u}}+\Bigg{(}\sum_{j=1}^{M}\,\alpha^{n+1}_{ij,k-1}\,f^{n+1}_{ij,k-1}\Bigg{)}_{i=1}^{M}\,.$

Finally, we use the fact that $\mathbf{u}^{n+1}_{k}$ satisfies

$(\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1}_{k-1})\,\mathbf{u}^{n+1}_{k}=\mathbb{M}_{\mathrm{L}}\,\tilde{\mathbf{u}}\,.$
(5.6)

It follows from the proof of Lemma 4.1 that, under the second condition in
(5.5), the matrix
$(\mathbb{M}_{\mathrm{L}}+\theta\,\tau_{n+1}\,{\mathbb{L}}^{n+1}_{k-1})$ is an
M-matrix and hence $\mathbf{u}^{n+1}_{k}$ is uniquely determined and satisfies
$\mathbf{u}^{n+1}_{k}\geq 0$. ∎

######  Remark 5.2.

From the physical point of view, the quantities $u$, $c$, and $p$ should be
not only non-negative but also bounded by $1$ from above. We have proved that
this is the case for the approximations of $c$. Moreover, if this would be
true also for the approximations of $u$, the integral form (3.9) would provide
this property also for the approximations of $p$. Unfortunately, a proof of
the upper bound for the approximations of $u$ is not available and numerical
results suggest that this bound can be violated. Note that a standard proof of
upper bounds for FCT discretizations relies on the decomposition
(4.18)–(4.20). Then, in particular, one would need that the solution of (4.18)
satisfies $\overline{\mathbf{u}}\leq 1$ if ${\mathbf{u}}^{n}\leq 1$. Choosing
${\mathbf{u}}^{n}={\mathbf{1}}$ (a vector with all components equal to $1$),
this requirement implies that
$(1-\theta)\,{\mathbb{A}}^{n}\,{\mathbf{1}}=(1-\theta)\,{\mathbb{L}}^{n}\,{\mathbf{1}}\geq
0$, i.e., the row sums of the matrix $(1-\theta)\,{\mathbb{A}}^{n}$ have to be
non-negative. Similarly, to derive an upper bound from (5.6), one would need
that $\theta\,{\mathbb{A}}^{n+1}_{k-1}\,{\mathbf{1}}\geq 0$. It is clear that
the validity of these row sum conditions cannot be expected.

######  Remark 5.3.

The second condition on $\tau_{n+1}$ in (5.5) depends on
${\mathbf{c}}^{n+1}_{k}$ which itself depends on $\tau_{n+1}$. Consequently,
in general, one has to proceed iteratively to find $\tau_{n+1}$ which
satisfies (5.5). To avoid this and also the dependence of $\tau_{n+1}$ on the
fixed-point iteration index $k$, it is possible to replace (5.5) by (4.5), cf.
Remark 4.3.

We summarize the procedure for obtaining a high-resolution positivity
preserving scheme for solving (2.1)–(2.5) in Algorithm 5.1.

Algorithm 5.1 Iterative scheme for computing an approximation of the solution
to the nonlinear FCT discretization.

1: Choose a tolerance $\mathrm{Tol}>0$ and a damping factor $\beta\in(0,1]$.

2: Compute the initial values $\mathbf{c}^{0}$, $\mathbf{p}^{0}$, and
$\mathbf{u}^{0}$ by (3.5).

3: Compute the mass matrix $\mathbb{M}$ and the lumped mass matrix
$\mathbb{M}_{\mathrm{L}}$.

4: for $n=0,1,\ldots,N-1$ do

5: Compute the stiffness matrix ${\mathbb{A}}^{n}$ and the artificial
diffusion matrix ${\mathbb{D}}^{n}$ and set
${\mathbb{L}}^{n}={\mathbb{A}}^{n}+{\mathbb{D}}^{n}$.

6: Choose $\tau_{n+1}$ satisfying (4.5).

7: Compute the intermediate solution $\overline{\mathbf{u}}$ from (4.18).

8: Set $\mathbf{c}^{n+1}_{0}=\mathbf{c}^{n}$,
$\mathbf{p}^{n+1}_{0}=\mathbf{p}^{n}$, and
$\mathbf{u}^{n+1}_{0}=\mathbf{u}^{n}$.

9: for $k=1,2,\ldots$ do

10: Compute $\mathbf{c}^{n+1}_{k}$ from (3.14) using
$\mathbf{c}^{n},\mathbf{p}^{n}$ and $\mathbf{p}_{k-1}^{n+1}$.

11: Compute $\mathbf{p}^{n+1}_{k}$ from (3.15) using
$\mathbf{p}^{n},\mathbf{c}^{n},\mathbf{u}^{n}$, $\mathbf{c}^{n+1}_{k}$, and
$\mathbf{u}_{k-1}^{n+1}$.

12: Compute the stiffness ${\mathbb{A}}^{n+1}_{k-1}$ from (3.17) using
$\mathbf{c}^{n+1}_{k}$ and $\mathbf{u}_{k-1}^{n+1}$.

13: Compute the artificial diffusion matrix ${\mathbb{D}}^{n+1}_{k-1}$ from
(5.2) and set
${\mathbb{L}}^{n+1}_{k-1}={\mathbb{A}}^{n+1}_{k-1}+{\mathbb{D}}^{n+1}_{k-1}$.

14: Compute the algebraic fluxes $f^{n+1}_{ij,k-1}$ from (5.3) and (5.4).

15: Compute the limiters $\alpha^{n+1}_{ij,k-1}$ by the Zalesak algorithm
(4.23)–(4.26) using the fluxes $f^{n+1}_{ij,k-1}$ and the intermediate
solution $\overline{\mathbf{u}}$.

16: Compute $\mathbf{u}^{n+1}_{k}$ by solving the linear system (5.1).

17: if
$\max\left\\{\|\mathbf{c}_{k}^{n+1}-\mathbf{c}_{k-1}^{n+1}\|_{M},\|\mathbf{p}_{k}^{n+1}-\mathbf{p}_{k-1}^{n+1}\|_{M},\|\mathbf{u}_{k}^{n+1}-\mathbf{u}_{k-1}^{n+1}\|_{M}\right\\}<\mathrm{Tol}$
then

18: Go to line 23.

19: else

20: Set
$\mathbf{c}_{k}^{n+1}:=\beta\,\mathbf{c}_{k}^{n+1}+(1-\beta)\,\mathbf{c}_{k-1}^{n+1}\,,\quad\mathbf{p}_{k}^{n+1}:=\beta\,\mathbf{p}_{k}^{n+1}+(1-\beta)\,\mathbf{p}_{k-1}^{n+1}\,,\quad\mathbf{u}_{k}^{n+1}:=\beta\,\mathbf{u}_{k}^{n+1}+(1-\beta)\,\mathbf{u}_{k-1}^{n+1}$.

21: end if

22: end for

23: Set
$\mathbf{c}^{n+1}=\mathbf{c}_{k}^{n+1}\,,\quad\mathbf{p}^{n+1}=\mathbf{p}_{k}^{n+1}\,,\quad\mathbf{u}^{n+1}=\mathbf{u}_{k}^{n+1}$.

24: end for

##  6 Numerical results

In the following, we present several numerical experiments to verify the
positivity preserving properties of the proposed scheme for the model
(2.1)–(2.5).

The computations are performed on a square domain $\Omega=(0,20)^{2}$ which is
decomposed into quadrilateral mesh cells obtained by uniform refinements.
Precisely, after $r$ refinements, the triangulation ${\mathcal{T}}_{h}$
consists of $2^{2\,r}$ equal squares. If not otherwise stated, we consider
five refinements, i.e., ${\mathcal{T}}_{h}$ consists of $32\times 32$ mesh
cells. As explained above, conforming bilinear finite elements are used for
approximating all unknown variables. The final time is $T=50$ and the
parameter $\epsilon=0.2$ is used. The values of the remaining parameters of
the model will be specified for the particular computations. The initial
conditions are defined by

$u^{0}(x)={\mathrm{e}}^{-|x|^{2}}\,,\qquad
c^{0}(x)=1-\frac{1}{2}\,{\mathrm{e}}^{-|x|^{2}}\,,\qquad
p^{0}{(x)}=\frac{1}{2}\,{\mathrm{e}}^{-|x|^{2}}\,.$

If not otherwise stated, we apply the A-stable Crank-Nicolson method
corresponding to $\theta=0.5$ for the time discretization. In one case, we
will also discuss the application of the unconditionally stable backward Euler
method corresponding to $\theta=1$. Algorithm 5.1 is used with the tolerance
$\mathrm{Tol}=10^{-8}$ and the damping factor $\beta=0.5$. The linear system
(5.1) is solved using the sparse direct solver UMFPACK [17]. Our newly
developed algorithms are implemented in the open-source finite element library
deal.II [4, 5].

###  6.1 Comparison between the standard Galerkin FEM and the FEM-FCT scheme
in presence of diffusion

To begin with, in the first example we consider the modified model subjected
to an extra diffusion term in the equation (2.1) with diffusion coefficient
$\alpha^{-1}$, as considered in [21], i.e., the equation (2.1) is replaced by
(2.6). We consider $\alpha=10$, $\chi=1$, and $\mu=1$. As can be seen from
Figs. 1 and 3, the FEM-FCT scheme introduces slightly more artificial
diffusion than the standard Galerkin FEM. One can observe that the cancer
cells invade the extracellular matrix and occupy the whole domain completely
at the final time. Next, we decrease the amount of the diffusion by setting
$\alpha=1000$ and keep the proliferation and haptotaxis rate as before. As can
be seen from Figs. 5 and 6, the standard Galerkin FEM shows some oscillations
in the front layer and the numerical simulation breaks down when the solution
reaches the boundary of the computational domain, whereas applying the FEM-FCT
removes the oscillations and keeps the solution positive at all times. The
corresponding snapshots of the cancer cell density, extracellular matrix, and
protease are plotted along the line $y=x$ in Figs. 2, 4, 6, and 8.

(a) t = 0

(b) t = 10

(c) t = 20

(d) t = 30

Figure 1:  Cancer cell invasion $u$ at different time instants $t=0,10,20,30$,
obtained with the standard Galerkin FEM for $\alpha=10$, $\mu=1$ and $\chi=1$.

(a) t = 0

(b) t = 10

(c) t = 20

(d) t = 30

Figure 2:  Cancer cell invasion $u$, connective tissue $c$, and protease $p$
at different time instants $t=0,10,20,30$, obtained with the standard Galerkin
FEM for $\alpha=10$, $\mu=1$ and $\chi=1$.

(a) t = 0

(b) t = 10

(c) t = 20

(d) t = 30

Figure 3:  Cancer cell invasion $u$ at different time instants $t=0,10,20,30$,
obtained with the FEM-FCT scheme for $\alpha=10$, $\mu=1$ and $\chi=1$.

(a) t = 0

(b) t = 10

(c) t = 20

(d) t = 30

Figure 4:  Cancer cell invasion $u$, connective tissue $c$, and protease $p$
at different time instants $t=0,10,20,30$, obtained with the FEM-FCT scheme
for $\alpha=10$, $\mu=1$ and $\chi=1$.

(a) t = 0

(b) t = 10

(c) t = 20

(d) t = 30

Figure 5: Cancer cell invasion $u$ at different time instants $t=0,10,20,30$,
obtained with the standard Galerkin FEM for $\alpha=1000$, $\mu=1$ and
$\chi=1$.

(a) t = 0

(b) t = 10

(c) t = 20

(d) t = 30

Figure 6: Cancer cell invasion $u$, connective tissue $c$, and protease $p$ at
different time instants $t=0,10,20,30$, obtained with the standard Galerkin
FEM for $\alpha=1000$, $\mu=1$ and $\chi=1$.

(a) t = 0

(b) t = 10

(c) t = 20

(d) t = 30

Figure 7: Cancer cell invasion $u$ at different time instants $t=0,10,20,30$,
obtained with the FEM-FCT scheme for $\alpha=1000$, $\mu=1$ and $\chi=1$.

(a) t = 0

(b) t = 10

(c) t = 20

(d) t = 30

Figure 8: Cancer cell invasion $u$, connective tissue $c$, and protease $p$ at
different time instants $t=0,10,20,30$, obtained with the FEM-FCT scheme for
$\alpha=1000$, $\mu=1$ and $\chi=1$.

###  6.2 The FEM-FCT scheme in absence of diffusion for $\chi=1$, $\mu=1$

In this section, we consider the case without the diffusion term, i.e.,
utilizing (2.1) now, and again set $\chi=\mu=1$. This case was studied in [14,
30], where the authors applied a nonstandard finite difference (NSFD) scheme
using Mickens rules. The proposed methods were successful in comparison to
standard finite difference methods at obtaining positive solutions, however,
some wiggles still remained in the vicinity of the front layer. On the other
hand, deriving an efficient NSFD scheme heavily depends on the type of the
system and the discretization of different terms. Therefore, in this work, we
applied the FEM-FCT methodology to remove the oscillations in the front layer
while keeping the solutions positive at all times, see Figs. 9 and 10. Next,
we check numerically whether the approximate solutions converge. To this end,
we computed the integrals of the solutions at the final time $t=50$ for
different numbers of global refinements, see Table 1. The results correspond
to the situation where the tumor is completely malignant and invades the whole
extracellular matrix. In Table 2, we study the values of the solutions at the
point $(20,20)$, the differences between two consecutive iterative solutions,
and the numbers of fixed-point iterations for different time steps. In
particular, we observe that the proposed scheme is convergent with respect to
the time step size. The convergence of the cancer cell invasion $u$ with
respect to the time step and the mesh width at two different time instants is
also studied in Fig. 11 by means of solution graphs along the line $y=x$.

(a) t = 0

(b) t = 10

(c) t = 20

(d) t = 30

Figure 9: Cancer cell invasion $u$ at different time instants $t=0,10,20,30$,
obtained with the FEM-FCT scheme for $\mu=1$ and $\chi=1$.

(a) t = 0

(b) t = 10

(c) t = 20

(d) t = 30

Figure 10: Decay of the extracellular matrix $c$ at different time instants $t=0,10,20,30$, computed using the FEM-FCT scheme for $\mu=1$ and $\chi=1$.  Table 1: Convergence of the mean values with respect to global mesh refinement at the last time instant $t=50$. # of refinements | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---
# DOF | 81 | 289 | 1089 | 4225 | 16641
$\int_{\Omega}c_{h}(x)\,\mathrm{d}x$ | 0.02362193 | 0.02670467 | 0.03284535 | 0.03042843 | 0.03680451
$\int_{\Omega}p_{h}(x)\,\mathrm{d}x$ | 0.02373726 | 0.02685417 | 0.03308441 | 0.03062835 | 0.03712137
$\int_{\Omega}u_{h}(x)\,\mathrm{d}x$ | 0.99999999 | 0.99999998 | 0.99999976 | 0.99999943 | 0.99999889
Table 2: Convergence of the solutions at the point $(20,20)$ with respect to the time step $\tau$ and convergence of the fixed-point iterations, both at the last time instant $t=50$.  $\tau$ | $c^{n+1}_{k}$ | $\|c^{n+1}_{k}-c^{n+1}_{k-1}\|$ | $p^{n+1}_{k}$ | $\|p^{n+1}_{k}-p^{n+1}_{k-1}\|$ | $u^{n+1}_{k}$ | $\|u^{n+1}_{k}-u^{n+1}_{k-1}\|$ | # Iterations
---|---|---|---|---|---|---|---
1.0 | 0.0333129 | 4.6407721e-09 | 0.0335600 | 5.0321162e-09 | 1.0004557 | 6.9181459e-11 | 21
0.1 | 0.0388396 | 5.5772037e-09 | 0.0391805 | 5.7753161e-09 | 1.0007926 | 4.4710309e-11 | 18
0.01 | 0.0387508 | 9.0822211e-09 | 0.0390900 | 9.2197598e-09 | 1.0007004 | 1.4400819e-10 | 14
0.001 | 0.0387514 | 7.2775629e-09 | 0.0390907 | 7.3738174e-09 | 1.0007904 | 1.2164834e-10 | 11
0.0001 | 0.0387516 | 5.8227737e-09 | 0.0390908 | 5.8989104e-09 | 1.0007904 | 9.7728734e-11 | 8

(a) t = 10

(b) t = 20

(c) t = 10

(d) t = 20

Figure 11: Cancer cell invasion $u$ computed using the FEM-FCT scheme at time
instants $t=10$ and $t=20$ for different time steps (first two pictures) and
for different numbers of global refinements (last two pictures).

###  6.3 Effect of haptotactic domination

In this section, we investigate the effect of directional movement of cancer
cells inside the domain. This is a very important property in cancer modeling
which can lead to metastasis. In metastasis, the cancer cells are moving to
the other parts of the body and start proliferate, forming a new tumor in the
new part, and invade the surrounding tissues. In this case, it is very
difficult to detect the location of cancerous cells and this is one of the
predominant causes of most deaths due to cancer. In the following, we only
study a very simple case of haptotactic dominating mechanism of the cancer
cell motion. In addition to the absence of the diffusion effect in the system,
there is only a small amount of the proliferation rate: we set $\mu=0.0001$
and $\chi=1$ in the computations. As a result of the haptotactic migration
domination, a small cluster of cancer cells builds up at the beginning and
this initial amount is expected to move along the direction of the gradient of
the extracellular matrix. As Fig. 12 indicates, the numerical simulation by
the standard Galerkin FEM breaks down in a very short amount of time after the
time instant $t=15$. Next, we apply the FEM-FCT scheme and the low-order
method, see Figs. 13 and 14, respectively. We observe that, in both cases, the
stabilization prevents the blow-up in the system and leads to non-negative
solutions. However, some oscillations still remain in the interior layer.
These oscillations could be suppressed by adaptive mesh refinement, which is
however out of the scope of this paper. As expected, the low-order method
provides a more diffusive solution than the FEM-FCT scheme. It is interesting
that, combining the FEM-FCT scheme with the backward Euler method
($\theta=1$), oscillation-free solutions are obtained, see Fig. 15.

(a) t = 0

(b) t = 5

(c) t = 15

Figure 12: The effect of the haptotactic rate on the cancer cell invasion $u$
at different time instants $t=0,5,15$, computed using the standard Galerkin
FEM for $\mu=0.0001$ and $\chi=1$.

(a) t = 10

(b) t = 20

(c) t = 30

(d) t = 40

Figure 13: The effect of the haptotactic rate on the cancer cell invasion $u$
at different time instants $t=10,20,30,40$, computed using the FEM-FCT scheme
for $\mu=0.0001$ and $\chi=1$.

(a) t = 10

(b) t = 20

(c) t = 30

(d) t = 40

Figure 14: The effect of the haptotactic rate on the cancer cell invasion $u$
at different time instants $t=10,20,30,40$, computed using the low-order
method for $\mu=0.0001$ and $\chi=1$.

(a) t = 10

(b) t = 20

(c) t = 30

(d) t = 40

Figure 15: The effect of the haptotactic rate on the cancer cell invasion $u$
at different time instants $t=10,20,30,40$, computed using the FEM-FCT scheme
with $\theta=1$ for $\mu=0.0001$ and $\chi=1$.

##  7 Conclusions

In this paper, we proposed a fully discrete nonlinear high-resolution
positivity preserving FEM-FCT scheme for chemotaxis equations without self-
diffusion term describing a model of cancer invasion. We proved the
solvability and positivity preservation of both the nonlinear discrete problem
and the linear problems appearing in fixed-point iterations. A series of
numerical experiments are shown to verify the robustness of the proposed
method. Derivation of error estimates is left to future work.

## Acknowledgments

This work was initiated during a research stay of the first author at the
Institute of Applied Mathematics at the Leibniz University Hanover from
November 2021 to April 2022 for which hospitality is still gratefully
acknowledged. The work of Shahin Heydari was further supported through the
grant No. 396921 of the Charles University Grant Agency. The work of Petr
Knobloch was supported through the grant No. 22-01591S of the Czech Science
Foundation.

## References

  * [1] Masashi Aida, Tohru Tsujikawa, Messoud Efendiev, Atsushi Yagi, and Masayasu Mimura.  Lower estimate of the attractor dimension for a chemotaxis growth system.  J. London Math. Soc. (2), 74(2):453–474, 2006.
  * [2] Masashi Aida and Atsushi Yagi.  Target pattern solutions for chemotaxis-growth system.  Sci. Math. Jpn., 59(3):577–590, 2004.
  * [3] Alexander R.A. Anderson, Mark A.J. Chaplain, E. Luke Newman, Robert J.C. Steele, and Alastair M. Thompson.  Mathematical modelling of tumour invasion and metastasis.  Computational and mathematical methods in medicine, 2(2):129–154, 2000.
  * [4] Daniel Arndt, Wolfgang Bangerth, Denis Davydov, Timo Heister, Luca Heltai, Martin Kronbichler, Matthias Maier, Jean-Paul Pelteret, Bruno Turcksin, and David Wells.  The DEAL.II finite element library: Design, features, and insights.  Comput. Math. Appl., 81:407–422, 2021.
  * [5] Daniel Arndt, Wolfgang Bangerth, Marco Feder, Marc Fehling, Rene Gassmöller, Timo Heister, Luca Heltai, Martin Kronbichler, Matthias Maier, Peter Munch, Jean-Paul Pelteret, Simon Sticko, Bruno Turcksin, and David Wells.  The deal.II library, Version 9.4.  J. Numer. Math., 30(3):231–246, 2022.
  * [6] Gabriel R. Barrenechea, Volker John, and Petr Knobloch.  Finite element methods respecting the discrete maximum principle for convection-diffusion equations.  SIAM Rev., accepted for publication, 2023.
  * [7] D.L. Book, J.P. Boris, and K. Hain.  Flux-corrected transport II: Generalizations of the method.  J. Comput. Phys., 18(3):248–283, 1975.
  * [8] Jay P. Boris and David L. Book.  Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works.  J. Comput. Phys., 11(1):38–69, 1973.
  * [9] J.P. Boris and D.L. Book.  Flux-corrected transport. III. minimal-error FCT algorithms.  J. Comput. Phys., 20(4):397–431, 1976.
  * [10] Vincent Calvez, Lucilla Corrias, and Mohamed Abderrahman Ebde.  Blow-up, concentration phenomenon and global existence for the Keller–Segel model in high dimension.  Comm. Partial Differential Equations, 37(4):561–584, 2012.
  * [11] M. A. J. Chaplain and G. Lolas.  Mathematical modelling of cancer cell invasion of tissue: the role of the urokinase plasminogen activation system.  Math. Models Methods Appl. Sci., 15(11):1685–1734, 2005.
  * [12] M. A. J. Chaplain and G. Lolas.  Mathematical modelling of cancer invasion of tissue: dynamic heterogeneity.  Netw. Heterog. Media, 1(3):399–439, 2006.
  * [13] Mark A.J. Chaplain and Andrew M. Stuart.  A model mechanism for the chemotactic response of endothelial cells to tumour angiogenesis factor.  Mathematical Medicine and Biology: A Journal of the IMA, 10(3):149–168, 1993.
  * [14] Michael Chapwanya, Jean M.-S. Lubuma, and Ronald E. Mickens.  Positivity-preserving nonstandard finite difference schemes for cross-diffusion equations in biosciences.  Comput. Math. Appl., 68(9):1071–1082, 2014.
  * [15] P. G. Ciarlet.  The finite element method for elliptic problems.  North-Holland, Amsterdam, 1978.
  * [16] L. Corrias, B. Perthame, and H. Zaag.  Global solutions of some chemotaxis and angiogenesis systems in high space dimensions.  Milan J. Math., 72:1–28, 2004.
  * [17] Timothy A. Davis.  Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method.  ACM Trans. Math. Software, 30(2):196–199, 2004.
  * [18] Yekaterina Epshteyn and Alexander Kurganov.  New interior penalty discontinuous Galerkin methods for the Keller–Segel chemotaxis model.  SIAM J. Numer. Anal., 47(1):386–408, 2008/09.
  * [19] Dianlei Feng, Insa Neuweiler, Udo Nackenhorst, and Thomas Wick.  A time-space flux-corrected transport finite element formulation for solving multi-dimensional advection-diffusion-reaction equations.  J. Comput. Phys., 396:31–53, 2019.
  * [20] Francis Filbet.  A finite volume scheme for the Patlak–Keller–Segel chemotaxis model.  Numer. Math., 104(4):457–488, 2006.
  * [21] Mario Fuest, Shahin Heydari, Petr Knobloch, Johannes Lankeit, and Thomas Wick.  Global existence of classical solutions and numerical simulations of a cancer invasion model.  ESAIM Math. Model. Numer. Anal., 57(4):1893–1919, 2023.
  * [22] Dirk Horstmann and Marcello Lucia.  Uniqueness and symmetry of equilibria in a chemotaxis model.  J. Reine Angew. Math., 654:83–124, 2011.
  * [23] Dirk Horstmann and Michael Winkler.  Boundedness vs. blow-up in a chemotaxis system.  J. Differential Equations, 215(1):52–107, 2005.
  * [24] Xueling Huang, Xinlong Feng, Xufeng Xiao, and Kun Wang.  Fully decoupled, linear and positivity-preserving scheme for the chemotaxis–Stokes equations.  Comput. Methods Appl. Mech. Engrg., 383:Paper No. 113909, 19, 2021\.
  * [25] Xueling Huang, Xufeng Xiao, Jianping Zhao, and Xinlong Feng.  An efficient operator-splitting FEM-FCT algorithm for 3D chemotaxis models.  Engineering with Computers, 36(4):1393–1404, 2020.
  * [26] Volker John and Petr Knobloch.  Existence of solutions of a finite element flux-corrected-transport scheme.  Appl. Math. Lett., 115:Paper No. 106932, 6, 2021.
  * [27] Volker John, Petr Knobloch, and Paul Korsmeier.  On the solvability of the nonlinear problems in an algebraically stabilized finite element method for evolutionary transport-dominated equations.  Math. Comp., 90(328):595–611, 2021.
  * [28] Evelyn F. Keller and Lee A. Segel.  Initiation of slime mold aggregation viewed as an instability.  Journal of theoretical biology, 26(3):399–415, 1970.
  * [29] Evelyn F. Keller and Lee A. Segel.  Model for chemotaxis.  Journal of theoretical biology, 30(2):225–234, 1971.
  * [30] M. Mehdizadeh Khalsaraei, Sh. Heydari, and L. Davari Algoo.  Positivity preserving nonstandard finite difference schemes applied to cancer growth model.  J. Cancer Treat. Res., 4(4):27–33, 2016.
  * [31] Mikhail K. Kolev, Miglena N. Koleva, and Lubin G. Vulkov.  An unconditional positivity-preserving difference scheme for models of cancer migration and invasion.  Mathematics, 10(1):131, 2022.
  * [32] D. Kuzmin and S. Turek.  Flux correction tools for finite elements.  J. Comput. Phys., 175(2):525–558, 2002.
  * [33] Dmitri Kuzmin.  Explicit and implicit FEM-FCT algorithms with flux linearization.  J. Comput. Phys., 228(7):2517–2534, 2009.
  * [34] Dmitri Kuzmin.  Algebraic flux correction I. Scalar conservation laws.  In Dmitri Kuzmin, Rainald Löhner, and Stefan Turek, editors, Flux-corrected transport. Principles, algorithms, and applications, pages 145–192. Springer, Dordrecht, second edition, 2012.
  * [35] Xingjie Helen Li, Chi-Wang Shu, and Yang Yang.  Local discontinuous Galerkin method for the Keller–Segel chemotaxis model.  J. Sci. Comput., 73(2-3):943–967, 2017.
  * [36] Rainald Löhner, Ken Morgan, Jaime Peraire, and Mehdi Vahdati.  Finite element flux-corrected transport (FEM–FCT) for the Euler and Navier–Stokes equations.  Int. J. Numer. Methods Fluids, 7(10):1093–1109, 1987.
  * [37] B. P. Marchant, J. Norbury, and A. J. Perumpanani.  Travelling shock waves arising in a model of malignant invasion.  SIAM J. Appl. Math., 60(2):463–476, 2000.
  * [38] B. P. Marchant, J. Norbury, and J. A. Sherratt.  Travelling wave solutions to a haptotaxis-dominated model of malignant invasion.  Nonlinearity, 14(6):1653–1671, 2001.
  * [39] Masayasu Mimura and Tohru Tsujikawa.  Aggregating pattern dynamics in a chemotaxis model including growth.  Physica A: Statistical Mechanics and its Applications, 230(3-4):499–543, 1996.
  * [40] Vidyanand Nanjundiah.  Chemotaxis, signal relaying and aggregation morphology.  Journal of Theoretical Biology, 42(1):63–105, 1973.
  * [41] Abbey J. Perumpanani, Jonathan A. Sherratt, John Norbury, and Helen M. Byrne.  A two parameter family of travelling waves with a singular barrier arising from the modelling of extracellular matrix mediated cellular invasion.  Phys. D, 126(3-4):145–159, 1999.
  * [42] David L. Ropp and John N. Shadid.  Stability of operator splitting methods for systems with indefinite operators: advection-diffusion-reaction systems.  J. Comput. Phys., 228(9):3508–3516, 2009.
  * [43] Norikazu Saito.  Conservative upwind finite-element method for a simplified Keller–Segel system modelling chemotaxis.  IMA J. Numer. Anal., 27(2):332–365, 2007.
  * [44] Andriy Sokolov, Ramzan Ali, and Stefan Turek.  An AFC-stabilized implicit finite element method for partial differential equations on evolving-in-time surfaces.  J. Comput. Appl. Math., 289:101–115, 2015.
  * [45] Andriy Sokolov, Robert Strehl, and Stefan Turek.  Numerical simulation of chemotaxis models on stationary surfaces.  Discrete Contin. Dyn. Syst. Ser. B, 18(10):2689–2704, 2013.
  * [46] R. Strehl, A. Sokolov, D. Kuzmin, and S. Turek.  A flux-corrected finite element method for chemotaxis problems.  Comput. Methods Appl. Math., 10(2):219–232, 2010.
  * [47] Robert Strehl, Andriy Sokolov, Dmitri Kuzmin, Dirk Horstmann, and Stefan Turek.  A positivity-preserving finite element method for chemotaxis problems in 3D.  J. Comput. Appl. Math., 239:290–303, 2013.
  * [48] M. Sulman and T. Nguyen.  A positivity preserving moving mesh finite element method for the Keller–Segel chemotaxis model.  J. Sci. Comput., 80(1):649–666, 2019.
  * [49] Youshan Tao and Mingjun Wang.  A combined chemotaxis-haptotaxis system: the role of logistic source.  SIAM J. Math. Anal., 41(4):1533–1558, 2009.
  * [50] Roger Temam.  Navier-Stokes equations. Theory and numerical analysis.  North-Holland, Amsterdam, 1977.
  * [51] R. Tyson, S.R. Lubkin, and James D. Murray.  A minimal mechanism for bacterial pattern formation.  Proceedings of the Royal Society of London. Series B: Biological Sciences, 266(1416):299–304, 1999.
  * [52] Rebecca Tyson, L. G. Stern, and Randall J. LeVeque.  Fractional step methods applied to a chemotaxis model.  J. Math. Biol., 41(5):455–475, 2000.
  * [53] Richard S. Varga.  Matrix iterative analysis.  Springer-Verlag, Berlin, 2000.
  * [54] Dianqing Wu.  Signaling mechanisms for regulation of chemotaxis.  Cell research, 15(1):52–56, 2005.
  * [55] Steven T. Zalesak.  Fully multidimensional flux-corrected transport algorithms for fluids.  J. Comput. Phys., 31(3):335–362, 1979.
  * [56] Jiansong Zhang, Jiang Zhu, and Rongpei Zhang.  Characteristic splitting mixed finite element analysis of Keller–Segel chemotaxis models.  Appl. Math. Comput., 278:33–44, 2016.
  * [57] Shubo Zhao, Xufeng Xiao, Jianping Zhao, and Xinlong Feng.  A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces.  Comput. Math. Appl., 79(11):3189–3205, 2020.