# Bound-preserving discontinuous Galerkin methods with modified Patankar time
integrations for chemical reacting flows111The first and third authors were
supported by the NSF grant DMS-1818467 and Simons Foundation 961585.

Fangyao Zhu222Department of Mathematical Sciences, Michigan Technological
University, Houghton, MI 49931. E-mail<EMAIL_ADDRESS>Juntao
Huang333Department of Mathematics and Statistics, Texas Tech University,
Lubbock, TX 79409. Email<EMAIL_ADDRESS>Yang Yang444Department of
Mathematical Sciences, Michigan Technological University, Houghton, MI 49931.
E-mail<EMAIL_ADDRESS>
###### Abstract

In this paper, we develop bound-preserving discontinuous Galerkin (DG) methods
for chemical reactive flows. There are several difficulties in constructing
suitable numerical schemes. First of all, the density and internal energy are
positive, and the mass fraction of each species is between 0 and 1. Secondly,
due to the rapid reaction rate, the system may contain stiff sources, and the
strong-stability-preserving explicit Runge-Kutta method may result in limited
time step sizes. To obtain physically relevant numerical approximations, we
apply the bound-preserving technique to the DG methods. Though traditional
positivity-preserving techniques can successfully yield positive density,
internal energy and mass fractions, it may not enforce the upper bound 1 of
the mass fractions. To solve this problem, we need to (1) make sure the
numerical fluxes in the equations of the mass fractions are consistent with
that in the equation of the density; (2) choose conservative time integrations
such that the summation of the mass fractions is preserved. With the above two
conditions, the positive mass fractions have summation 1, then they are all
between 0 and 1. For time discretization, we apply the modified Runge-
Kutta/multi-step Patankar methods, which are explicit for the flux while
implicit for the source. Such methods can handle stiff sources with relatively
large time steps, preserve the positivity of the target variables, and keep
the summation of the mass fractions to be 1. Finally, it is not
straightforward to combine the bound-preserving DG methods and the Patankar
time integrations. The positivity-preserving technique for DG method requires
positive numerical approximations at the cell interfaces, while Patankar
methods can keep the positivity of the pre-selected point-values of the target
variables. To match the degree of freedom, we use $Q^{k}$ polynomials on
rectangular meshes for problems in two space dimensions. To evolve in time, we
first read the polynomials at the Gaussian points. Then suitable slope
limiters can be applied to enforce the positivity of the solutions at those
points, which can be preserved by the Patankar methods, leading to positive
updated numerical cell averages. In addition, we use another slope limiter to
get positive solutions used for the bound-preserving technique for the flux.
Numerical examples are given to demonstrate the good performance of the
proposed schemes.

Key Words: compressible Euler equations, chemical reacting flows, bound-
preserving, discontinuous Galerkin method, modified Patankar method

##  1 Introduction

In this paper, we develop numerical methods for the chemical reacting flows
[30]

$\displaystyle\rho_{t}+m_{x}+n_{y}=0$ (1.1a) $\displaystyle
m_{t}+(mu+p)_{x}+(mv)_{y}=0$ (1.1b) $\displaystyle
n_{t}+(nu)_{x}+(nv+p)_{y}=0$ (1.1c) $\displaystyle
E_{t}+((E+p)u)_{x}+((E+p)v)_{y}=0$ (1.1d)
$\displaystyle(c_{1})_{t}+(mz_{1})_{x}+(nz_{1})_{y}=s_{1}$ (1.1e)
$\displaystyle...$
$\displaystyle(c_{M-1})_{t}+(mz_{M-1})_{x}+(nz_{M-1})_{y}=s_{M-1}.$ (1.1f)

Here $\rho$ is the density, $(u,v)$ is the velocity field, $m=\rho u$ and
$n=\rho v$ are the momentum in $x$ and $y$ directions, $E$ is the total energy
which is the summation of the kinetic and internal energy of the fluid, $p$ is
the pressure, $z_{i}$ is the mass fraction of the $i$-th species with
$i=1,2,...,M$ and $M$ being the number of species. For $1\leq i\leq M$, we
have $c_{i}=\rho z_{i}$ and $\sum_{i=1}^{M}z_{i}=1$. Therefore, the total
density satisfies $\rho=\sum_{i=1}^{M}c_{i}$ and $0\leq z_{i}\leq 1$. The
pressure is

$p=RT\sum_{i=1}^{M}\frac{\rho_{i}}{M_{i}}$

where $R$ is the universal gas constant, $T=\rho/p$ is the temperature, and
$M_{i}$ is the molar mass of the $i$-th species. The total energy $E$ is
expressed as

$E=\sum_{i=1}^{M}\rho_{i}e_{in,i}(T)+\sum_{i=1}^{M}\rho_{i}q_{i}+\frac{1}{2}\rho(u^{2}+v^{2})$

where $q_{i}$ is the enthalpy for the $i$th species and $e_{in,i}(T)=C_{i}T$
is the internal energy of the $i$th species with $C_{i}=3R/2M_{i}$ and
$5R/2M_{i}$ for the monoatomic and diatomic species, respectively.

We can write (1.1) in a compact form :

$\displaystyle\bm{U}_{t}+\bm{F}\left(\bm{U}\right)_{x}+\bm{G}\left(\bm{U}\right)_{y}=\bm{S}\left(\bm{U}\right)$
(1.2)

where $\bm{S}\left(\bm{U}\right)$ denotes the source describing the chemical
reactions of the form

$\nu_{1,j}^{{}^{\prime}}X_{1}+\nu_{2,j}^{{}^{\prime}}X_{2}+...+\nu_{M,j}^{{}^{\prime}}X_{M}\leavevmode\nobreak\
{\rightleftharpoons}\leavevmode\nobreak\
\nu_{1,j}^{{}^{\prime\prime}}X_{1}+\nu_{2,j}^{{}^{\prime\prime}}X_{2}+...+\nu_{M,j}^{{}^{\prime\prime}}X_{M},\quad
j=1,...,J,$

where $J$ is the total number of reactions, $\nu_{i,j}^{{}^{\prime}}$ and
$\nu_{i,j}^{{}^{\prime\prime}}$ are the stoichiometric coefficients of the
reactants and productions of the $i$th species in the $j$th reaction. For non-
equilibrium chemistry, the rate of production in (1.1) can be written as

$s_{i}=M_{i}\sum_{j=1}^{J}(\nu_{i,j}^{{}^{\prime\prime}}-\nu_{i,j}^{{}^{\prime}})(k_{f,j}\prod_{s=1}^{M}\left(\frac{c_{s}}{M_{s}}\right)^{\nu_{s,j}^{{}^{\prime}}}-k_{b,j}\prod_{s=1}^{M}\left(\frac{c_{s}}{M_{s}}\right)^{\nu_{s,j}^{{}^{\prime\prime}}}),\quad
i=1,...,M.$

where $k_{f,j}$ and $k_{b,j}$ are functions of temperature indicating the
forward and backward reaction rates.

The solution to chemical reacting flows may contain shocks. Moreover,
numerical schemes may result in non-physical numerical approximations, e.g.
the density and pressure are negative, and the mass fractions are out of the
interval $[0,1]$. The non-physical numerical approximations may further lead
to ill-posed problems and eventually blow-up of the numerical simulations.
Hence, constructing a bound-preserving scheme is essential. In this paper, we
apply discontinuous Galerkin (DG) method that is high order accurate and
flexible on geometry. The method was first introduced in 1973 [25] for the
neutron transport equation, a time independent hyperbolic equation. Later,
Cockburn and Shu extended the DG method for solving time dependent problems
such as nonlinear convection problems and Euler equations. The framework was
given in a series of papers [3, 4, 5, 6, 7], where the DG method was coupled
with Runge-Kutta (RK) time integration along with TVB nonlinear limiters to
achieve non-oscillatory properties for strong shocks. However, the TVB limiter
is not sufficient to maintain the positivity of the numerical approximations.
In [29], high order DG methods for two-dimensional gaseous detonations were
constructed to preserve the positivity of density, pressure, and all the mass
fractions. The idea is to apply first order Euler forward time discretization
and find a sufficient condition for the cell averages of the DG numerical
approximations to be positive. Then a slope limiter is applied to construct
new physically relevant numerical approximations, keeping the original cell
averages. The time discretizations can be extended to high order strong-
stability-presering (SSP) RK or multi-step (MS) methods [11, 26, 27] since
they are convex combinations of the Euler forward method. For our problem, we
need to preserve not only the lower bound 0 but also the upper bound 1 for the
mass fractions. The first work preserving the physical bounds of the mass
fractions was given in [12, 2], where the compressible miscible displacements
in porous media were discussed. Later, the idea was adapted to multi-species
and multi-reaction detonations to construct high order DG schemes in [9, 10].
The basic strategy is to apply the positivity-preserving techniques to each
mass fraction $z_{i}$ and enforce $\sum_{i=1}^{M}z_{i}=1$ to obtain physically
relevant numerical approximations by using conservative time integrations. The
extension to finite difference methods was also given in [8].

Another difficulty is how to deal with the stiff source terms due to the rapid
reaction rate in chemical reactive flows. Direct application of the explicit
time integration may result in limited time steps, and one may consider
Implicit-Explicit (IMEX) methods which treats the source term implicitly while
the flux term explicitly. By doing so, the positivity-preserving technique can
also be applied to the flux terms. Unfortunately, it is not easy to construct
implicit solvers for the source terms if the reactions are complicated and the
system turns out to be fully nonlinear. Therefore, one may try to linearize
the source term, keeping the conservation of mass. One method that can
preserve the physical bounds and treat the source term implicitly is the
modified Patankar (MP) time integration, which was developed from the Patankar
trick that first introduced in [24].

To demonstrate the idea of the MP time integration, we consider a class of
system of ODEs which describes the chemical reactions by ignoring the
convection terms in (1.2):

$\displaystyle\frac{dc_{i}}{dt}=M_{i}\sum_{r=1}^{R}(\nu_{i,r}^{{}^{\prime\prime}}-\nu_{i,r}^{{}^{\prime}})(k_{f,r}\prod_{s=1}^{M}\left(\frac{c_{s}}{M_{s}}\right)^{\nu_{s,r}^{{}^{\prime}}}-k_{b,r}\prod_{s=1}^{M}\left(\frac{c_{s}}{M_{s}}\right)^{\nu_{s,r}^{{}^{\prime\prime}}}),\quad
i=1,...,M.$ (1.3)

Equation (1.3) can be seen as production-destruction equations which have the
form

$\displaystyle\frac{dc_{i}}{dt}=P_{i}(c)-D_{i}(c),\quad i=1,2,...,M,$ (1.4)

with

$P_{i}(c)=\sum_{j=1}^{M}p_{ij}(c),\quad D_{i}(c)=\sum_{j=1}^{M}d_{ij}(c)$

and

$p_{ij}(c)=d_{ji}(c)\geq 0$

where $c=(c_{1},c_{2},...,c_{M})^{T}$ denotes the concentrations, the
production function $p_{ij}(c)$ describes the transformation rate of $j$-th
component to $i$-th component; whereas $d_{ij}(c)$, the destruction function,
denotes the transformation rate of $i$th component to $j$-th component. The
solution of (1.4) is conservative which means that $\sum_{i=1}^{N}c_{i}(t)$
remains unchanged with respect to time. In addition, the solution should be
positive if the initial condition is positive and $d_{ij}(c)=0$ for $c_{i}=0$.

The Patankar trick [24] was to multiply the destruction term by the ratio of
the numerical approximations between two time levels, leading to a linear
implicit scheme. Unfortunately, such a method fails the mass conservation. In
contrast, the MP scheme in [1] guarantees both positivity and conservation
properties with any time step size, and the method reads as:

$\displaystyle c_{i}^{n+1}=c_{i}^{n}+\Delta
t\left(\sum_{j}p_{ij}(c^{n})\frac{c_{j}^{n+1}}{c_{j}^{n}}-\sum_{j}d_{ij}(c^{n})\frac{c_{i}^{n+1}}{c_{i}^{n}}\right).$
(1.5)

However, the scheme is first order accurate only. In [19, 20], MP RK schemes
of second and third order were introduced. In [23], the MP trick is adapted to
deferred correction (DeC) schemes and developed MPDeC schemes of arbitrary
order of accuracy. More recently, the stability analysis of the MP schemes are
given in [16, 17, 13, 18]. However, it is not easy to couple the convection
terms into the schemes in [19], [20] and [23]. Therefore, instead of using the
classical form of RK methods, Huang and Shu applied RK schemes of Shu-Osher
form [26] and constructed another class of MP RK schemes with the SSP
structure. These conservative and unconditionally positivity-preserving MP RK
methods were given in [14, 15] along with the necessary and sufficient
conditions derived to obtain the desired order of accuracy. In these two
works, the time integration was combined with the positivity-preserving finite
difference weighted essentially non-oscillatory (WENO) scheme.

In this paper, we first construct the MS MP time integration and couple the
RK/MS MP time integration with DG spatial discretization. The combination is
not straightforward due to the inconsistency of the collocation points
required in the positivity-preserving technique and the MP time integrations.
Due to the SSP structure of the MP time integration to be used in this paper,
we use (1.5) to demonstrate the full algorithm, and the fully-discrete scheme
with convection term can be written as

$c_{i}^{n+1}=c_{i}^{n}+\Delta tf(c^{n})+\Delta
t\left(\sum_{j}p_{ij}(c^{n})\frac{c_{j}^{n+1}}{c_{j}^{n}}-\sum_{j}d_{ij}(c^{n})\frac{c_{i}^{n+1}}{c_{i}^{n}}\right),$

where $f(c)$ is the convection term after DG spatial discretization. Given the
physically relevant numerical approximations at time level $n$, we first treat
the convection term only, and define $\tilde{c}_{i}=c_{i}^{n}+\Delta
tf(c^{n})$. The positivity-preserving technique yields positive numerical cell
average of $\tilde{c}$. Then we apply a slope limiter to obtain positive
$\tilde{c}$ at Gaussian points in each cell as this is required by the MP time
integration. Next, we update the numerical approximation at each Gaussian
point by using the MP time integration, and the resulted numerical
approximations are also positive, leading to positive numerical cell averages.
Finally, another slope limiter is called to have positive numerical
approximations at the cell interfaces which will be used in the next time
level.

The paper is organized as follows. In Section 2, we first demonstrate the MP
RK method and then develop the second and third order MS MP scheme and show
the necessary and sufficient conditions to obtain the desired order of
accuracy. In Section 3, we show the positivity-preserving property of the DG
method as well as the bound-preserving technique for the mass fractions. In
Section 4, we give numerical examples to demonstrate the performance of the
scheme. Some concluding remarks will be given in Section 5.

##  2 The ODE solver

In this section, we demonstrate the time integrations to be used in this
paper. We first review the MP RK method given in [14, 15] and then construct
the MP MS methods.

###  2.1 The modified Patankar Runge-Kutta methods

In this subsection, we briefly present the MP RK methods proposed in [14, 15].
The general production-destruction system can be written as

$\displaystyle\frac{dc_{i}}{dt}=\sum_{j=1}^{N}p_{ij}(c)-\sum_{j=1}^{N}d_{ij}(c),\quad
i=1,2,\dots,N,$ (2.1)

where $c_{i}$ is the concentration of the $i$-th component ${\bf
c}=(c_{1},c_{2},\cdots,c_{M})^{T}$. The second order MP RK scheme [14] is

$\displaystyle c_{i}^{(0)}$ $\displaystyle=$ $\displaystyle c_{i}^{n},$ (2.2a)
$\displaystyle c_{i}^{(1)}$ $\displaystyle=$
$\displaystyle\alpha_{10}c_{i}^{(0)}+\Delta
t\beta_{10}\left(\sum_{j}p_{ij}(c^{(0)})\frac{c_{j}^{(1)}}{c_{j}^{(0)}}-\sum_{j}d_{ij}(c^{(0)})\frac{c_{i}^{(1)}}{c_{i}^{(0)}}\right),$
(2.2b) $\displaystyle c_{i}^{n+1}$ $\displaystyle=$
$\displaystyle\alpha_{20}c_{i}^{(0)}+\alpha_{21}c_{i}^{(1)}$ (2.2c)
$\displaystyle+\Delta
t\left(\sum_{j}\left(\beta_{20}p_{ij}(c^{(0)})+\beta_{21}p_{ij}(c^{(1)})\right)\frac{c_{j}^{n+1}}{(c_{j}^{(1)})^{s}(c_{j}^{(0)})^{1-s}}\right.$
$\displaystyle\left.-\sum_{j}\left(\beta_{20}d_{ij}(c^{(0)})+\beta_{21}d_{ij}(c^{(1)})\right)\frac{c_{i}^{n+1}}{(c_{i}^{(1)})^{s}(c_{i}^{(0)})^{1-s}}\right)$

with the parameters satisfying

$\alpha_{10}=1,\quad\alpha_{20}=1-\alpha,\quad\alpha_{21}=\alpha,$
$\beta_{10}=\beta,\quad\beta_{20}=1-\frac{1}{2\beta}-\alpha\beta,\quad\beta_{21}=\frac{1}{2\beta},$
$s=\frac{1-\alpha\beta+\alpha\beta^{2}}{\beta(1-\alpha\beta)}$

where

$0\leq\alpha\leq 1,\quad\beta>0,\quad\alpha\beta+\frac{1}{2\beta}\leq 1.$

The third order MP RK scheme [15] is

$\displaystyle c_{i}^{(0)}$ $\displaystyle=$ $\displaystyle c_{i}^{n},$ (2.3a)
$\displaystyle c_{i}^{(1)}$ $\displaystyle=$
$\displaystyle\alpha_{10}c_{i}^{(0)}+\Delta
t\beta_{10}\left(\sum_{j}p_{i,j}(c^{(0)})\frac{c_{j}^{(1)}}{c_{j}^{(0)}}-\sum_{j}d_{i,j}(c^{(0)})\frac{c_{i}^{(1)}}{c_{i}^{(0)}}\right),$
(2.3b) $\displaystyle\rho_{i}$ $\displaystyle=$ $\displaystyle
n_{1}c_{i}^{(1)}+n_{2}c_{i}^{n}\left(\frac{c_{i}^{(1)}}{c_{i}^{n}}\right)^{2},$
$\displaystyle c_{i}^{(2)}$ $\displaystyle=$
$\displaystyle\alpha_{20}c_{i}^{(0)}+\alpha_{21}c_{i}^{(1)}$ (2.3c)
$\displaystyle+\Delta
t\sum_{j}\left(\left(\beta_{20}p_{ij}^{(0)}+\beta_{21}p_{ij}^{(1)}\right)\frac{c_{j}^{(2)}}{\rho_{j}}-\left(\beta_{20}d_{ij}^{(0)}+\beta_{21}d_{ij}^{(1)}\right)\frac{c_{i}^{(2)}}{\rho_{i}}\right),$
$\displaystyle\mu_{i}$ $\displaystyle=$ $\displaystyle
c_{i}^{n}\left(\frac{c_{i}^{(1)}}{c_{i}^{n}}\right)^{s},$ $\displaystyle
a_{i}$ $\displaystyle=$
$\displaystyle\eta_{1}c_{i}^{n}+\eta_{2}c_{i}^{(1)}+\Delta
t\sum_{j}\left(\left(\eta_{3}p_{ij}^{n}+\eta_{4}p_{ij}^{(1)}\right)\frac{a_{j}}{\mu_{j}}-\left(\eta_{3}d_{ij}^{n}+\eta_{4}d_{ij}^{(1)}\right)\frac{a_{i}}{\mu_{i}}\right),$
(2.3d) $\displaystyle\sigma_{i}$ $\displaystyle=$ $\displaystyle
a_{i}+zc_{i}^{n}\frac{c_{i}^{(2)}}{\rho_{i}},$ $\displaystyle c_{i}^{n+1}$
$\displaystyle=$
$\displaystyle\alpha_{30}c_{i}^{(0)}+\alpha_{31}c_{i}^{(1)}+\alpha_{32}c_{i}^{(2)}$
(2.3e) $\displaystyle+\Delta
t\sum_{j}\left(\left(\beta_{30}p_{ij}^{(0)}+\beta{31}p_{ij}^{(1)}+\beta_{32}p_{ij}^{(2)}\right)\frac{c_{j}^{n+1}}{\sigma_{j}}-\left(\beta_{30}d_{ij}^{(0)}+\beta_{31}d_{ij}^{(1)}+\beta{32}d_{ij}^{(2)}\right)\frac{c_{i}^{n+1}}{\sigma_{i}}\right)$

where the coefficients $\alpha_{ij}$ and $\beta_{ij}$, and the parameters
$n_{1}$, $n_{2}$, $\eta_{1}$, $\eta_{2}$, $\eta_{3}$, $\eta_{4}$, $s$ and $z$
are given in [15]. This scheme is third-order accurate and unconditionally
positivity-preserving.

###  2.2 The modified Patankar multi-step methods

We focus on the production-destruction equations (2.1) and derive MP MS scheme
in this subsection. In [11], the explicit SSP MS methods for the nonlinear ODE

$\displaystyle\frac{du}{dt}=L(u)$ (2.4)

are given as follows: for the second-order,

$\displaystyle u^{n+1}=\frac{1}{4}u^{n-2}+(\frac{3}{4}u^{n}+\frac{3}{2}\Delta
tL(u^{n})),$ (2.5)

and for the third-order,

$\displaystyle u^{n+1}=(\frac{11}{27}u^{n-3}+\frac{4}{9}\Delta
tL(u^{n-3}))+(\frac{16}{27}u^{n}+\frac{16}{9}\Delta tL(u^{n})).$ (2.6)

We can construct second and third order MP MS method according to it.

####  2.2.1 Second-order scheme

The second-order explicit SSP multistep method for (2.1) is

$\displaystyle
c_{i}^{n+1}=\frac{1}{4}c_{i}^{n-2}+\frac{3}{4}c_{i}^{n}+\frac{3}{2}\Delta
t(\sum_{j}p_{ij}^{n}-\sum_{j}d_{ij}^{n}).$ (2.7)

Following [14], we make some modification to preserve conservation and
positivity unconditionally. The MP MS method is rewritten as

$\displaystyle
c_{i}^{n+1}=\frac{1}{4}c_{i}^{n-2}+\frac{3}{4}c_{i}^{n}+\frac{3}{2}\Delta
t(\sum_{j}p_{ij}^{n}\frac{c_{j}^{n+1}}{\sigma_{j}}-\sum_{j}d_{ij}^{n}\frac{c_{i}^{n+1}}{\sigma_{i}}),$
(2.8)

where $\sigma_{i}\geq 0$ is an undetermined function of $c_{i}^{n-2}$,
$c_{i}^{n-1}$ and $c_{i}^{n}$.

If we assume

$\displaystyle\sigma_{i}=c_{i}^{n+1}+O(\Delta t^{2})$ (2.9)

then by (2.8) we have

$\displaystyle c_{i}^{n+1}$
$\displaystyle=\frac{1}{4}c_{i}^{n-2}+\frac{3}{4}c_{i}^{n}+\frac{3}{2}\Delta
t(\sum_{j}p_{ij}^{n}\frac{c_{j}^{n+1}}{c_{j}^{n+1}+O(\Delta
t^{2})}-\sum_{j}d_{ij}^{n}\frac{c_{i}^{n+1}}{c_{i}^{n+1}+O(\Delta t^{2})})$
$\displaystyle=\frac{1}{4}c_{i}^{n-2}+\frac{3}{4}c_{i}^{n}+\frac{3}{2}\Delta
t(\sum_{j}p_{ij}^{n}(1+O(\Delta t^{2}))-\sum_{j}d_{ij}^{n}(1+O(\Delta
t^{2})))$
$\displaystyle=\frac{1}{4}c_{i}^{n-2}+\frac{3}{4}c_{i}^{n}+\frac{3}{2}\Delta
t(\sum_{j}p_{ij}^{n}-\sum_{j}d_{ij}^{n})+O(\Delta t^{3})$

which shows that (2.8) is second-order accurate by using the fact that the
explicit MS method (2.7) is second-order accurate. Therefore, (2.9) is a
sufficient condition that (2.8) is second-order accurate.

In the second-order scheme (2.8), if we take

$\displaystyle\sigma_{i}=(c_{i}^{n})^{s}(c_{i}^{n-1})^{r}(c_{i}^{n-2})^{1-r-s}$
(2.10)

with $s$ and $r$ to be determined later, then the Taylor expansions for
$c_{i}^{n-k}$ at $t^{n+1}$ is

$\displaystyle c_{i}^{n-k}=c_{i}^{n+1}-(k+1)\Delta
t(P_{i}^{n+1}-D_{i}^{n+1})+O(\Delta t^{2})$ (2.11)

for $k=0,1,2$. For convenience of notation, we denote

$\displaystyle Q_{1}:=P_{i}^{n+1}-D_{i}^{n+1}$ (2.12)

for further analysis and then (2.11) becomes

$\displaystyle c_{i}^{n-k}=c_{i}^{n+1}-(k+1)\Delta tQ_{1}+O(\Delta t^{2}).$
(2.13)

Therefore, the Taylor expansions for $\sigma_{i}$ at $t^{n+1}$ is

$\displaystyle\sigma_{i}$
$\displaystyle=(c_{i}^{n})^{s}(c_{i}^{n-1})^{r}(c_{i}^{n-2})^{1-r-s}$
$\displaystyle=(c_{i}^{n+1}-\Delta tQ_{1})^{s}(c_{i}^{n+1}-2\Delta
tQ_{1})^{r}(c_{i}^{n+1}-3\Delta tQ_{1})^{1-r-s}+O(\Delta t^{2})$
$\displaystyle=c_{i}^{n+1}+(r+2s-3)\Delta tQ_{1}+O(\Delta t^{2}),$

and now we have

$\displaystyle r+2s-3=0$ (2.14)

which is required for (2.10) to be satisfied.

Before stating the main theorem, we give the following lemma:

######  Lemma 2.1.

The MP MS scheme (2.8) is conservative. If
$\sum_{i=1}^{M}c_{i}^{n-2}=\sum_{i=1}^{M}c_{i}^{n}=1$, then
$\sum_{i=1}^{M}c_{i}^{n+1}=1$ for $i=1,2,...,M$.

###### Proof.

We add up (2.8) over $i$ and use the fact that the modified Patankar scheme is
conservative, i.e.
$\sum_{i,j}(p_{ij}^{n}\frac{c_{j}^{n+1}}{\sigma_{j}}-d_{ij}^{n}\frac{c_{i}^{n+1}}{\sigma_{i}})=0$
to obtain

$\sum_{i}c_{i}^{n+1}=\frac{1}{4}\sum_{i}c_{i}^{n-2}+\frac{3}{4}\sum_{i}c_{i}^{n}.$

∎

######  Lemma 2.2.

The MP MS scheme (2.8) is unconditionally positivty-preserving. That is, for
all $\Delta t>0$ and $c_{i}^{n-2}$,$c_{i}^{n-1}$ and $c_{i}^{n}>0$, we have
$c_{i}^{n+1}>0$ for $i=1,2,...,M$.

The proof is similar to Lemma 2.7 in [19], so we skip it. This lemma directly
results in the following theorem.

######  Theorem 2.3.

The MP MS scheme (2.8) is second-order accurate with

$\sigma_{i}=(c_{i}^{n})^{s}(c_{i}^{n-1})^{r}(c_{i}^{n-2})^{1-r-s}$

and

$r+2s-3=0.$

Moreover, it is conservative in the sense that

$\sum_{i=1}^{M}c_{i}^{n-2}=\sum_{i=1}^{M}c_{i}^{n-1}=\sum_{i=1}^{M}c_{i}^{n}=\sum_{i=1}^{M}c_{i}^{n+1}$

and unconditionally positivity-preserving: if $c_{i}^{n}\geq 0$ for
$i=1,2,...M$, then $c_{i}^{n+1}\geq 0$ for $i=1,2,...,M$.

####  2.2.2 Third-order scheme

Following the same approach in the previous subsection, we have the modified
third order scheme

$\displaystyle
c_{i}^{n+1}=\frac{11}{27}c_{i}^{n-3}+\frac{16}{27}c_{i}^{n}+\Delta
t(\sum_{j}(\frac{4}{9}p_{ij}^{n-3}+\frac{16}{9}p_{ij}^{n})\frac{c_{j}^{n+1}}{\sigma_{j}}-(\frac{4}{9}d_{ij}^{n-3}+\frac{16}{9}d_{ij}^{n})\frac{c_{i}^{n+1}}{\sigma_{i}}),$
(2.15)

where $\sigma_{i}\geq 0$ are undetermined functions of $c_{i}^{n-k}$,
$k=0,1,2,3$. If we assume that

$\displaystyle\sigma_{i}=c_{i}^{n+1}+O(\Delta t^{3}),$ (2.16)

then (2.15) is third-order accurate.

We will derive explicit expressions of $\sigma_{i}$ in (2.15) for the
sufficient condition (2.16) to be satisfied. Following the analysis for the
second-order scheme (2.7), we try to make

$\displaystyle\sigma_{i}=(c_{i}^{n})^{s}(c_{i}^{n-1})^{r}(c_{i}^{n-2})^{p}(c_{i}^{n-3})^{1-r-s-p}$
(2.17)

with $s$, $r$ and $p$ to be determined. Taylor expansions for $c_{i}^{n-k}$ at
$t^{n+1}$ with $k=0,1,2,3$ give

$\displaystyle c_{i}^{n-k}=c_{i}^{n+1}-(k+1)\Delta
t(P_{i}^{n+1}-D_{i}^{n+1})+\frac{(k+1)^{2}}{2}\Delta
t^{2}\frac{\partial(P_{i}^{n+1}-D_{i}^{n+1})}{\partial
c}(P^{n+1}-D^{n+1})+O(\Delta t^{3}).$ (2.18)

For convenience of notation, we denote

$\displaystyle Q_{1}:=P_{i}^{n+1}-D_{i}^{n+1},\quad
Q_{2}:=\frac{\partial(P_{i}^{n+1}-D_{i}^{n+1})}{\partial c}(P^{n+1}-D^{n+1}),$
(2.19)

and then

$\displaystyle c_{i}^{n-k}=c_{i}^{n+1}-(k+1)\Delta
tQ_{1}+\frac{(k+1)^{2}}{2}\Delta t^{2}Q_{2}+O(\Delta t^{3}).$ (2.20)

Taylor expansion for $\sigma_{i}$ at $t^{n+1}$ yields

$\displaystyle\sigma_{i}$
$\displaystyle=(c_{i}^{n})^{s}(c_{i}^{n-1})^{r}(c_{i}^{n-2})^{p}(c_{i}^{n-3})^{1-r-s-p}$
$\displaystyle=(c_{i}^{n+1}-\Delta tQ_{1}+\frac{1}{2}\Delta
t^{2}Q_{2})^{s}(c_{i}^{n+1}-2\Delta tQ_{1}+2\Delta t^{2}Q_{2})^{r}$
$\displaystyle\qquad(c_{i}^{n+1}-3\Delta tQ_{1}+\frac{9}{2}\Delta
t^{2}Q_{2})^{p}(c_{i}^{n+1}-4\Delta tQ_{1}+8\Delta
t^{2}Q_{2})^{1-r-s-p}+O(\Delta t^{3})$
$\displaystyle=c_{i}^{n+1}+(p+2r+3s-4)\Delta
tQ_{1}-(7p+12r+15s-16)\frac{\Delta t^{2}}{2}Q_{2}$
$\displaystyle\qquad+((p+2r+3s)^{2}-4r-9s-p)\frac{\Delta
t^{2}Q_{1}^{2}}{2c_{i}^{n+1}}+O(\Delta t^{3})$

Letting the coefficients to be zero, we can solve out

$\displaystyle r=-3s+6,\quad p=3s-8.$ (2.21)

######  Theorem 2.4.

The MP MS scheme (2.15) is third-order accurate with

$\sigma_{i}=(c_{i}^{n})^{s}(c_{i}^{n-1})^{r}(c_{i}^{n-2})^{p}(c_{i}^{n-3})^{1-r-s-p}$

and

$r=-3s+6,\quad p=3s-8.$

Moreover, it is conservative which means that

$\sum_{i=1}^{N}c_{i}^{n-3}=\sum_{i=1}^{N}c_{i}^{n-2}=\sum_{i=1}^{N}c_{i}^{n-1}=\sum_{i=1}^{N}c_{i}^{n}=\sum_{i=1}^{N}c_{i}^{n+1}$

and unconditionally positivity-preserving: if $c_{i}^{n}\geq 0$ for
$i=1,2,...M$, then $c_{i}^{n+1}\geq 0$ for $i=1,2,...,M$.

##  3 Bound-preserving DG scheme

In this section, we discuss the positivity-preserving DG scheme. We first
review the DG scheme and then discuss its positivity-preserving property. Then
we demonstrate the bound preserving technique for mass fractions. Finally, we
apply the ODE solver to our problem.

###  3.1 The DG scheme

In this subsection, we demonstrate the DG scheme for (1.2). We define the
finite element space $V_{h}^{k}$ as

$V_{h}^{k}=\\{v:v|_{K}\in Q^{k}(K),\>\>\forall K\in\Omega_{h}\\},$

where $Q^{k}(K)$ is the set of tensor product polynomials of degree at most
$k$ in cell $K$. Then the DG scheme is to find the numerical solution
${\bm{U_{h}}}\in[V_{h}^{k}]^{M+3}$ such that for all $v\in[V_{h}^{k}]^{M+3}$
we have

$\displaystyle\frac{d}{dt}\int_{K}{\bm{U_{h}}}vd\bm{x}+\sum_{e\in\partial
K}\int_{e}\bm{H}(U_{h}^{int},U_{n}^{ext},\nu_{e})vd\Gamma-\int_{K}{\bm{F}(\bm{U_{h}})}v_{x}d\bm{x}-\int_{K}{\bm{G}(\bm{U_{h}})}v_{y}d\bm{x}=\int_{K}{\bm{S}(\bm{U_{h}})}vd\bm{x},$
(3.1)

where $\nu_{e}$ is the outward normal of the edge $e$ on element $K$. We use
Lax-Friedrichs flux in this paper

$\displaystyle\bm{H}(U_{1},U_{2},\bm{\nu}_{e})=\frac{1}{2}[f(U_{1})\bm{\nu}_{e}+f(U_{2})\bm{\nu}_{e}-\alpha(U_{1}-U_{2})],$
(3.2)

with $\alpha=\||\langle u,v\rangle|+c\|_{\infty}$, where $c=\sqrt{\frac{\gamma
p}{\rho}}$ is the sound speed. Note that $f=({\bm{F},\bm{G}})$ and
$\bm{H}=(h_{\rho},h_{m},h_{n},h_{E},h_{1},\cdots,h_{M-1})^{T}$. The following
definition [9, 10] is used in the bound-preserving technique.

######  Definition 3.1.

We say the elements in numerical flux $\bm{H}$ in (3.2) are consistent if
$h_{\rho}=h_{i}$ if we take $z_{i}=1$ for all $1\leq i\leq M-1$.

As discussed in [9, 10], the Lax-Friedrichs flux is consistent. To present the
positivity-preserving technique, we consider Euler forward time
discretization. We take the test function to be 1, then the equation satisfied
by the cell averages is

$\displaystyle\bar{{\bm{U}}}_{K}^{n+1}=\bar{U}_{K}^{n}-\frac{\Delta
t}{|K|}\sum_{e\in\partial
K}\int_{e}{\bm{H}(\bm{U_{h}}^{int},\bm{U_{h}}^{ext}},\bm{\nu}_{e})d\Gamma+\frac{\Delta
t}{|K|}\int_{K}{\bm{S}(\bm{U_{h}}})d\bm{x},$ (3.3)

where $\bar{U}_{K}^{n}=\frac{1}{|K|}\int_{K}Ud\bm{x}$ is the cell average of
$U$ in cell $K$ at time level $n$ and $\Delta t$ is the time step size. The
integrals are approximated by proper quadrature rules which will be discussed
in the next section.

###  3.2 Bound-Preserving technique

In this section, we develop the bound-preserving technique for the convection
term. First, ignoring the source term from (3.3), we have

$\displaystyle\bar{\bm{U}}_{K}^{n+1}=\bar{\bm{U}}_{K}^{n}-\frac{\Delta
t}{|K|}\sum_{e\in\partial
K}\int_{e}\bm{H}(\bm{U_{h}}^{int},\bm{U_{h}}^{ext},\bm{\nu}_{e})d\Gamma.$
(3.4)

We are seeking numerical approximations chosen from an admissible set $G$
defined as

$G=\left\\{\bm{U}=\left(\begin{array}[]{cc}\rho\\\ m\\\ n\\\ E\\\ c_{1}\\\
...\\\ c_{M-1}\\\ \end{array}\right),\quad\rho>0,\quad p>0,\quad z_{i}>0,\quad
i=1,...,M,\quad\sum_{i=1}^{M}z_{i}=1\right\\}.$

As demonstrated in [9], $G$ is a convex set.

In this paper, we consider rectangular meshes. The spatial domain
$\Omega=[a,b]\times[c,d]$ is partitioned into $N\times N$ cells. The
partitions are assumed to be uniform throughout this paper for simplicity.
However, this assumption is not essential. We denote the cells in the $x$ and
$y$ directions by $I_{i}=[x_{i-1/2},x_{i+1/2}],$ and
$J_{j}=[y_{j-1/2},y_{j+1/2}],$ respectively. for each cell $K\in\Omega$ and
let
$K_{ij}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$
be the $(i,j)$-th cell. We then define a set of quadrature points $S_{K}$ on
cell $K$. We use $L$ points Gauss quadrature with $L\geq k+1$ for the
integrals in (3.4). The Gauss quadrature points on
$[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$ and
$[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$ are denoted as

$p_{i}^{x}=\\{x_{i}^{\beta}:\beta=1,...,L\\}\quad and\quad
p_{j}^{y}=\\{y_{j}^{\beta}:\beta=1,...,L\\}.$

In addition, we denote $\hat{L}$ Gauss-Lobatto points with $2\hat{L}-1\geq k$
on $[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$ and
$[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$ as

$\hat{p}_{i}^{x}=\\{\hat{x}_{i}^{\alpha}:\alpha=0,...,\hat{L}\\}\quad
and\quad\hat{p}_{j}^{y}=\\{\hat{y}_{j}^{\alpha}:\alpha=0,...,\hat{L}\\}.$

Define

$\displaystyle
S_{K_{i,j}}=(p_{i}^{x}\otimes\hat{p}_{j}^{y})\cup(\hat{p}_{i}^{x}\otimes
p_{j}^{y})\cup(p_{i}^{x}\otimes p_{j}^{y}).$ (3.5)

After defining the quadrature points, we can now state the theorem given in
[9, 10] below.

######  Theorem 3.2.

If the numerical approximation $\bm{U}_{K}(x,y)\in G\quad\forall(x,y)\in
S_{K_{ij}}$, where $S_{K_{ij}}$ is defined in (3.5) for rectangular meshes,
then the DG scheme (3.4) is positivity-preserving, namely,
$\bar{\bm{U}}_{K}^{n+1}\in G$ under the time step restriction

$\displaystyle\alpha(\lambda_{1}+\lambda_{2})\leq\hat{\omega}_{1}$ (3.6)

where $\lambda_{1}=\frac{\Delta{t}}{\Delta x}$ and
$\lambda_{2}=\frac{\Delta{t}}{\Delta y}.$

###  3.3 Bound preserving technique for MP method

In this subsection, we focus on the bound preserving technique for MP time
integrations. We apply DG scheme as in section 3.1 to the spatial
discretization except source term. We use $k$-point Gaussian quadrature for
the integral in each direction in the convection terms and take Lagrangian
basis so that the point values at the quadrature points are known. Therefore,
we denote $c_{k,i,l}=c_{k,i,l}(t)$ as the point value at the $l$-th Gaussian
quadrature point of the $i$-th species at the $k$-th cell after the spatial
discretization, and define $c_{k,l}=(c_{k,1,l},\cdots,c_{k,M,l})^{T}$. To
include the convection term and construct the semi-discrete scheme, we take
the test function to be the $l$-th Lagrangian basis and formulate the system
of ODEs base on the destruction and production equation (2.1) in the following
form:

$\displaystyle\frac{dc_{k,i,l}}{dt}=F_{k,i,l}(c)+P_{k,i,l}(c)-D_{k,i,l}(c),\quad
k=1,\cdots,N^{2},\quad i=1,\cdots,M.$ (3.7)

$F_{k,i,l}=F_{k,i,l}(c)$ denotes the contributions of the convection terms
after spatial discretizations in the PDEs. The production and destruction
terms are $P_{k,i,l}=P_{k,i,l}(c)=\sum_{j=1}^{M}p_{i,j}(c_{k,l})$ and
$D_{k,i,l}=D_{k,i,l}(c)=\sum_{j=1}^{M}d_{i,j}(c_{k,l})$ which satisfy

$p_{i,j}(c_{k,l})=d_{j,i}(c_{k,l})\geq 0,\quad\forall i,j,k,l\quad and\quad
c\geq 0$

For simplicity of presentation, we drop the subscript $l$ in the rest of this
section. We make the following assumption on (3.7) :
Assumption 2.1. The Euler forward method for the convection term satisfies the
positivity-preserving property: if $c_{k,i}^{n}\geq 0$ for all $k,i$, then
there exists $\Delta t_{0}>0$ such that

$c_{k,i}^{n}+\Delta tF_{k,i}(c^{n})\geq 0$

for all $k,i$ and $\Delta t\leq\Delta t_{0}$.

The bound-preserving technique can be extended to SSP RK/MS methods which are
convex combinations of forward Euler. For simplicity, we consider the second-
order MPMS scheme (2.8) only and incorporate the convection term $F_{k,i}$.
The fully discrete scheme is given as

$\displaystyle
c_{k,i}^{n+1}=\frac{1}{4}c_{k,i}^{n-2}+\frac{3}{4}c_{k,i}^{n}+\frac{3}{2}\Delta
tF_{k,i}(c^{n})+\frac{3}{2}\Delta
t\left(\sum_{j}p_{k,i,j}^{n}\frac{c_{k,j}^{n+1}}{\sigma_{k,j}}-\sum_{j}d_{k,i,j}^{n}\frac{c_{k,i}^{n+1}}{\sigma_{k,i}}\right)$
(3.8)

where

$\sigma_{i}=(c_{i}^{n})^{s}(c_{i}^{n-1})^{r}(c_{i}^{n-2})^{1-r-s}.$

Clearly, the scheme is positivity-preserving if the time step satisfies

$\Delta t\leq\frac{1}{2}\Delta t_{0}.$

Similarly for the third-order MP MS scheme (2.8), incorporating the convection
term $F_{k,i}$ we have

$\displaystyle
c_{k,i}^{n+1}=\frac{11}{27}c_{k,i}^{n-3}+\frac{16}{27}c_{k,i}^{n}+\Delta
t(\frac{4}{9}F_{k,i}(c^{n-3})+\frac{16}{9}F_{k,i}(c^{n}))$
$\displaystyle+\Delta
t(\sum_{j}(\frac{4}{9}p_{k,i,j}^{n-3}+\frac{16}{9}p_{k,i,j}^{n})\frac{c_{k,j}^{n+1}}{\sigma_{j}}-(\frac{4}{9}d_{k,i,j}^{n-3}+\frac{16}{9}d_{k,i,j}^{n})\frac{c_{k,i}^{n+1}}{\sigma_{i}})$
(3.9)

where

$\sigma_{i}=(c_{i}^{n})^{s}(c_{i}^{n-1})^{r}(c_{i}^{n-2})^{p}(c_{i}^{n-3})^{1-r-s-p}.$

The scheme is positivity-preserving if the time step satisfies

$\Delta t\leq\frac{1}{3}\Delta t_{0}.$

###  3.4 The limiter

In this subsection, We discuss the limiter to be used to deal with negative
numerical approximations with positive numerical cell averages.

To enforce the positivity, a suitable limiter [32] can be added to keep
density, pressure and mass fraction being within the physical bounds.
Following [9, 10], the main idea of the positivity-preserving limiter is to
modify the DG polynomial $w_{K}(x)$ into polynomial

$\tilde{w}_{K}=\theta_{K}(w_{K}-\bar{w}_{K})+\bar{w}_{k}$

where $\bar{w}_{K}\in G$ is the cell average and $\theta_{K}\in[0,1]$. This
limiter does not degenerate the accuracy of the DG polynomial $w_{K}(x)$ [33].
The algorithm of the limiter on each fixed element $K$ is given below:

  1. 1.

For a small number $\epsilon=10^{-13}$, if $\bar{\rho}>\epsilon$, we proceed
to the next step. Otherwise, we simply take $U=\bar{U}$

  2. 2.

Enforce the positivity of the density $\rho$: Compute the minimum value as

$\rho_{min}=\min_{(x,y)\in S_{K}}\rho(x,y),$

where $S_{K}$ is defined in (3.5). If $\rho_{min}<0$, then take

$\hat{\rho}=\bar{\rho}+\theta_{K}(\rho-\bar{\rho}),$
$\hat{c}_{i}=\bar{c}_{i}+\theta_{K}(c_{i}-\bar{c}_{i}),i=1,...,M-1$

with

$\theta_{K}=\frac{\bar{\rho}-\epsilon}{\bar{\rho}-\rho_{min}}.$

Here $\hat{r}_{M}=\bar{r}_{M}+\theta(r_{M}-\bar{r}_{M})$ to keep
$\sum_{i}^{M}\hat{r}_{i}=\hat{\rho}$.

  3. 3.

Enforce positivity of the mass fractions: For $1\leq i\leq M$, define
$\hat{S}_{i}=\\{(x,y)\in S_{K}:\hat{c}_{i}(x,y)\leq 0\\}$. Take

$\displaystyle\tilde{c}_{i}=\hat{c}_{i}+\theta(\frac{\bar{c}_{i}}{\bar{\rho}}\hat{\rho}-\hat{c}_{i}),\
1\leq i\leq M-1,$ $\displaystyle\>\>\>\theta=\max_{1\leq i\leq
M}\max_{(x,y)\in\hat{S}_{i}}\\{\frac{-\hat{c}_{i}(x,y)\bar{\rho}}{\bar{c}_{i}\hat{\rho}(x,y)-\hat{c}_{i}(x,y)\bar{\rho}},0\\}.$
(3.10)

  4. 4.

Modify the pressure: Denote
$\tilde{{\bm{U}}}=(\hat{\rho},m,n,E,\tilde{c}_{1},...,\tilde{c}_{M-1})^{T}$.
For each $\bm{x}\in S$, if $\tilde{\bm{U}}\in G$, then take
$\theta_{\bm{x}}=1$. Otherwise, take

$\theta_{\bm{x}}=\frac{p(\bar{\bm{U}})}{p(\bar{\bm{U}})-p(\tilde{\bm{U}}(x))}.$

Then, we use

$\bm{U}^{new}=\bar{\bm{U}}+\theta(\tilde{\bm{U}}-\bar{\bm{U}}),\quad\theta=\min_{\bm{x}\in
S_{K}}\theta_{\bm{x}},$

as the new DG approximation.

###  3.5 Full algorithm and the main theorem

We have discussed the DG scheme and MP MS time integration. Putting them
together, we have a conservative and positivity-preserving scheme for solving
Euler equations. Now we demonstrate the full algorithm of our method in the
flow chart. We use the MPMS2 method on rectangular meshes as an example:

  1. 1.

Consider the convection term only

${\bm{U}}_{k,i}^{(1)}=\frac{1}{4}\bm{U}_{k,i}^{n-2}+\frac{3}{4}\bm{U}_{k,i}^{n}+\frac{3}{2}\Delta
t\bm{F}_{k,i}(\bm{U}^{n}),$

where $\bm{U}=(\rho z_{1},\cdots,\rho z_{M},m,E)^{T}.$

  2. 2.

Apply limiter in Section 3.4 to polynomials at the Gaussian points with
$S_{k_{1}}=p_{i}^{x}\otimes p_{j}^{y}$.

  3. 3.

Use Patankar to deal with source

$c_{k,i}^{n+1}=c_{k,i}^{(1)}+\frac{3}{2}\Delta
t\left(\sum_{j}p_{k,i,j}^{n}\frac{c_{k,j}^{n+1}}{\sigma_{k,j}}-\sum_{j}d_{k,i,j}^{n}\frac{c_{k,i}^{n+1}}{\sigma_{k,i}}\right),$

where $c=(\rho z_{1},\cdots,\rho z_{M})^{T}.$

  4. 4.

Apply limiter in Section 3.4 to polynomials at points in
$S_{k_{2}}=(\hat{p}_{i}^{x}\otimes
p_{j}^{y})\cup(p_{i}^{x}\otimes\hat{p}_{j}^{y})$.

  5. 5.

Let $\bm{U}^{n-2}=\bm{U}^{n-1}$, $\bm{U}^{n-1}=\bm{U}^{n}$, and
$\bm{U}^{n}=\bm{U}^{n+1}$, then go back to step 1 and repeat the process 1-4
until the final time is reached.

Remark 3.1. The reason we apply the limiter twice in the above flow chart is
because each limiter aims for different set of points. In Step 2, the limiter
is applied on the quadrature points. In this way, we ensure that the point
values on quadrature points are kept positive in each cell $K$. Fortunately,
Patankar in Step 3 will not destroy the positivity of the cell average and the
quadrature point values, but it may not yield positive numerical
approximations on the cell interfaces. So we apply the limiter again in Step 4
to those points. The numerical approximations are guaranteed to be in the
admissible set.
The main theorem is stated below.

######  Theorem 3.3.

Consider the DG scheme (3.1) coupled with the second-order MP MS method (3.8)
or the third-order one (3.9).The schemes are bound preserving: if
$\bm{U}_{k,i}^{n}\in G$ for all $1\leq k\leq N^{2}$,$1\leq i\leq M$, then
$\bm{U}_{k,i}^{n+1}\in G$ under the condition $\Delta
t\leq\min_{i}\frac{\alpha_{i}}{\beta_{i}}\Delta\tilde{t}$ with
$\Delta\tilde{t}$ satisfies the condition in (3.6).

###### Proof.

Let us first prove the positivity preserving property for the DG scheme with
the first order explicit Euler forward time integration. Without source term,
the equation satisfied by cell average is given in (3.3). The positivity-
preserving technique yields positive numerical cell average. Then applying
slope limiter, we obtain the positive numerical approximation at Gaussian
points. The momentum and energy equations are trivial since there is no source
for them. For the equations of mass fractions, we approximate the source terms
using MP trick

$\displaystyle\sum_{j}p_{ij}^{n}\frac{c_{j}^{n+1}}{\sigma_{j}}-\sum_{j}d_{ij}^{n}\frac{c_{i}^{n+1}}{\sigma_{i}}.$
(3.11)

The goal is to show that $\bar{c}_{k}^{n+1}\geq 0$. A sufficient condition is
to obtain positive numerical approximations at the Gaussian quadrature points
and this can be achieved by the MP time integration. We can refer to [1] for
the proof of positivity of the technique. Extend to second and third order
multi-step time discretization will keep the positivity since they are convex
combinations of Euler forward scheme.

To start the discussion about the upper bound of the mass fraction $z_{i}\leq
1$, we first subtract (1.1e)-(1.1f) from (1.1a) to obtain a new equation

$\displaystyle(c_{M})_{t}+(mz_{M})_{x}+(nz_{M})_{y}=s_{M}$ (3.12)

Then we apply the positivity-preserving technique to each $c_{k,i}^{n}$ for
$i=1,\dots,M$. Therefore, $c_{i}>0(z_{i}>0)$ are obtained at all time levels.
Moreover, we use consistent fluxes and conservative time integration to
preserve the total mass $\sum_{i=1}^{M}c_{i}=\rho$. As a result, the numerical
approximation of $z_{i}$ is bounded in interval $[0,1]$. ∎

##  4 Numerical examples

In this section, we provide numerical experiments to show the performance of
the numerical scheme.

###  4.1 Test of the ODE solver

In this part, we test the accuracy of MP MS methods in solving linear and
nonlinear ODEs

######  Example 1.

This example is to test the accuracy on a linear problem:

$\frac{dc_{1}}{dt}=c_{2}-ac_{1},$ $\frac{dc_{2}}{dt}=ac_{1}-c_{2},$

with constant $a>0$, and initial value $c_{1}(0)=c_{1}^{0}$,
$c_{2}(0)=c_{2}^{0}$. The exact solutions are

$c_{1}(t)=(1+b\exp(-(a+1)t))c_{1}^{\infty},\quad
c_{2}(t)=c_{1}^{0}+c_{2}^{0}-c_{1}(t),$

with the parameters $c_{1}^{\infty}$ and $b$ determined by

$c_{1}^{\infty}=\frac{c_{1}^{0}+c_{2}^{0}}{a+1},\quad
b=\frac{c_{1}^{0}}{c_{1}^{\infty}}-1.$

In the numerical experiment, we take $c_{1}^{0}=4.5$, $c_{2}^{0}=3.2$, $a=2.7$
and the final time $t=1$.

| MPMS2 | MPMS3
---|---|---
$\Delta t$ | Error | Order | Error | Order
1/20 | 1.88e-02 | – | 1.03e-03 | –
1/40 | 4.56e-03 | 2.04 | 1.27e-04 | 3.02
1/80 | 1.09e-03 | 2.07 | 1.58e-05 | 3.00
1/160 | 2.75e-04 | 1.99 | 2.04e-06 | 2.96
1/320 | 6.97e-05 | 1.98 | 2.61e-07 | 2.97
Table 1: Example 1: Accuracy test of MPMS2 and MPMS3 for linear ODEs

The errors between the numerical and exact solutions at the final time are
listed in Table 1. From the table, we can observe optimal convergence rates
for second and third order MP MS methods.

######  Example 2.

This example is to test the accuracy for solving a nonlinear problem given in
[14]:

$\displaystyle\frac{dc_{1}}{dt}$ $\displaystyle=$ $\displaystyle
F_{1}(c)-\frac{c_{1}c_{2}}{c_{1}+1},$ $\displaystyle\frac{dc_{2}}{dt}$
$\displaystyle=$ $\displaystyle F_{2}(c)+\frac{c_{1}c_{2}}{c_{1}+1}-ac_{2},$
$\displaystyle\frac{dc_{3}}{dt}$ $\displaystyle=$ $\displaystyle
F_{3}(c)+ac_{2},$

where $(F_{1}(c),F_{2}(c),F_{3}(c))$ denotes "convection terms".

To express this system of ODEs in the form of production-destruction
equations, we set

$p_{21}=d_{12}(c)=\frac{c_{1}c_{2}}{c_{1}+1},\quad
p_{32}(c)=d_{23}(c)=ac_{2},$

and $p_{ij}=d_{ij}=1$ for other sets of $i$, $j$. The initial conditions are
set as $c_{1}^{0}=9.98$, $c_{2}^{0}=0.01$ and $c_{3}^{0}=0.01$. The convection
terms are

$(F_{1}(c),F_{2}(c),F_{3}(c))=\left(c_{1}c_{2}c_{3},\frac{c_{3}}{c_{2}},c_{2}c_{2}c_{3}^{2}\right).$

The final time $t=1$ and the parameter $a=1$.

The errors between the numerical and exact solutions at the final time are
listed in Table 2, from which we can also observe the optimal convergence
rates.

| MPMS2 | MPMS3
---|---|---
$\Delta t$ | Error | Order | Error | Order
1/20 | 2.13e-03 | – | 3.67e-04 | –
1/40 | 4.43e-04 | 2.26 | 3.68e-05 | 3.32
1/80 | 1.14e-04 | 1.96 | 4.115e-06 | 3.16
1/160 | 2.72e-05 | 2.06 | 4.87e-07 | 3.07
1/320 | 6.76e-06 | 2.01 | 5.924e-08 | 3.04
Table 2: Example 2: Accuracy test of MPMS2 and MPMS3 for nonlinear ODEs

Remark 4.1. The parameters $s$, $r$, and $p$ were chosen based on the
relationships in (2.14) and (2.21) for the numerical schemes to be at their
desired accuracy. We tested for different values of $s$ from $-1$ to $5$ hence
different values of $r$ and $p$ with fixed time step in example 4.1 and
example 4.2. We want to know the influence of different set of parameters on
the error. The results are the four plots given in Figure 1. At this point, we
could only observe and conclude that there seems to be a set of parameters
which gives the smallest error. This is an interesting topic and we leave it
in the future work.

(a) Example 4.1 using MPMS2

(b) Example 4.1 using MPMS3

(c) Example 4.2 using MPMS2

(d) Example 4.2 using MPMS3

Figure 1: Parameter $s$ and error plots

###  4.2 Euler Equations with Three species Reactions and General Equation of
State

######  Example 3.

We solve the three species model of the one-dimensional Euler system with a
more general equation of state in [30]

$\displaystyle U$ $\displaystyle=$
$\displaystyle(\rho_{1},\rho_{2},\rho_{3},\rho u,E)^{T},$ $\displaystyle F(U)$
$\displaystyle=$ $\displaystyle(\rho_{1}u,\rho_{2}u,\rho_{3}u,\rho
u^{2}+P,(E+p)u)^{T},$ $\displaystyle S(U)$ $\displaystyle=$
$\displaystyle(2M_{1}\omega,-M_{2}\omega,0,0,0)^{T}.$

The rate of the chemical reaction is given by

$\omega=\left(k_{f}(T)\frac{\rho_{2}}{M_{2}}-k_{b}(T)\left(\frac{\rho_{1}}{M_{1}}\right)^{2}\right)\sum_{s=1}^{3}\frac{\rho_{s}}{M_{s}},\quad
k_{f}=CT^{-2}e^{-E/T},$
$k_{b}=k_{f}/exp(b_{1}+b_{2}logz+b_{3}z+b_{4}z^{2}+b_{5}z^{3}),\quad
z=10000/T.$

The parameters are $M_{1}=0.016$, $M_{2}=0.032$, $M_{3}=0.028$,
$q_{1}=1.558\times 10^{7}$, $R=8.31447215$, $C_{0}=2.9\times 10^{17}m^{3}$,
$E_{0}=59750K$, and $b_{1}=2.855$, $b_{2}=0.988$, $b_{3}=-6.181$,
$b_{4}=-0.023$, $b_{5}=-0.001$. For this model, we split the source terms
$2M_{1}\omega$ into two parts:

$2M_{1}\omega=\omega_{+}-\omega_{-},$

with

$\omega_{+}=2M_{1}k_{f}(T)\frac{\rho_{2}}{M_{2}}\sum_{s=1}^{3}\frac{\rho_{s}}{M_{s}}\geq
0\quad
and\quad\omega_{-}=2M_{1}k_{b}(T)\left(\frac{\rho_{1}}{M_{1}}\right)^{2}\sum_{s=1}^{3}\frac{\rho_{s}}{M_{s}}\geq
0.$

The production and destruction terms are

$p_{2,1}=d_{1,2}=\omega_{-}\geq 0,\quad p_{1,2}=d_{2,1}=\omega_{+}\geq 0.$

The eigenvalues of the Jacobian are $(u,u,u,u-c,u+c)$ where
$c=\sqrt{\frac{\gamma p}{\rho}}$ with
$\gamma=1+\frac{p}{T\sum_{s=1}^{3}\rho_{s}e_{s}^{{}^{\prime}}(T)}$. The
initial conditions are: the densities $\rho_{1}$,$\rho_{2}$ and $\rho_{3}$ are
$5.251896311257204\times 10^{-5}$, $3.748071704863518\times
10^{-5}$,$2.962489471973072\times 10^{-4}$ on the left, and
$8.341661837019181\times 10^{-8}$, $9.455418692098664\times 10^{-11}$,
$2.748909430004963\times 10^{-7}$ on the right. The velocities are zero. The
pressures are 1000 on the left and 1 on the right. The final time is
$t=0.0001$. The densities, velocity and pressure are presented in Figure2 for
MPRK2, Figure 3 for MPMS2 and Figure 4 for MPMS3 methods. We can observe some
oscillations in the numerical approximations. This is because we did not apply
oscillation suppressors such as WENO algorithm in the numerical scheme. Though
the oscillations exist, the bound-preserving technique is enough to stabilize
our numerical scheme.

(a) $\rho_{1}$ (b) $\rho_{2}$

(c) $\rho_{3}$ (d) $u$

(e) $p$

Figure 2: Example 4.3: three species reaction problem at t=0.0001 using MPRK2

(a) $\rho_{1}$ (b) $\rho_{2}$

(c) $\rho_{3}$ (d) $u$

(e) $p$

Figure 3: Example 4.3: three species reaction problem at t=0.0001 using MPMS2

(a) $\rho_{1}$ (b) $\rho_{2}$

(c) $\rho_{3}$ (d) $u$

(e) $p$

Figure 4: Example 4.3: three species reaction problem at t=0.0001 using MPMS3

###  4.3 Reactive Euler Equations

Consider the reactive Euler equations in 2D with

$\displaystyle\bm{U}$ $\displaystyle=$ $\displaystyle(\rho,m,n,E,\rho Y)^{T},$
$\displaystyle\bm{F}(\bm{U})$ $\displaystyle=$ $\displaystyle(m,\rho
u^{2}+p,\rho uv,(E+p)u,\rho uY)^{T},$ $\displaystyle\bm{G}(\bm{U})$
$\displaystyle=$ $\displaystyle(n,\rho uv,\rho v^{2}+p,(E+p)v,\rho vY)^{T},$
$\displaystyle\bm{S}(\bm{U})$ $\displaystyle=$
$\displaystyle(0,0,0,0,\omega),$ (4.1)

where

$m=\rho u,\quad n=\rho v,\quad
E=\frac{1}{2}\rho(u^{2}+v^{2})+\frac{p}{\gamma-1}+\rho qY,$

with $Y$ denoting the reactant mass fraction and $\omega$ being the source
term which is expressed as

$\omega=-\tilde{K}\rho Ye^{-\tilde{T}/T}.$

To fit into our framework, we rewrite (4.1) into

$\displaystyle\bm{U}$ $\displaystyle=$ $\displaystyle(\rho Y,\rho
Z,m,n,E)^{T},$ $\displaystyle\bm{F}(\bm{U})$ $\displaystyle=$
$\displaystyle(\rho uY,\rho uZ,\rho u^{2}+p,\rho uv,(E+p)u)^{T},$
$\displaystyle\bm{G}(\bm{U})$ $\displaystyle=$ $\displaystyle(\rho vY,\rho
vZ,\rho uv,\rho v^{2}+p,(E+p)v)^{T},$ $\displaystyle\bm{S}(\bm{U})$
$\displaystyle=$ $\displaystyle(\omega,-\omega,0,0,0),$

where $Z$ denotes the unreacted mass fraction.

######  Example 4.

In this example, the numerical convergence of our scheme is tested. The domain
is set to be $[0,2]\times[0,2]$. The initial condition is given as follows: if
$x^{2}+y^{2}\leq 0.36$, then $(\rho,u,v,p,Y)=(1,0,0,80,0)$; otherwise,
$(\rho,u,v,p,Y)=(1,0,0,10^{-9},1)$. The boundary conditions on the bottom and
left are reflective. The terminal time is t = 0.2. We use uniform rectangular
meshes with mesh sizes $\Delta x=\Delta y=1/120$. The numerical results with
MPMS2 and MPRK2 are shown in Figure 5 and Figure 6, respectively, which are
comparable to the results in [29] where the converged solutions are observed.
Also, we can observe some spurious oscillations in the numerical
approximations, and this is mainly due to the lack of mechanisms to suppress
oscillations.

######  Example 5.

We test the detonation diffraction in this example. The same parameters and
initial conditions in [29] are applied. The initial condition are as follows:
if $x<0.5$, then $(\rho,u,v,E,Y)=(11,6.18,0,970,1)$; otherwise,
$(\rho,u,v,E,Y)=(1,0,0,5,5,1)$. The boundary conditions are reflective except
that at $x=0$, where we take $(\rho,u,v,E,Y)=(11,6.18,0,970,1)$. The terminal
time is $t=0.6$. The parameters are $\gamma=1.2$, $q=50$, $\tilde{T}=50$ and
$\tilde{K}=2566.4$. The numerical schemes may produce negative density and or
pressure for this example which can lead to blow-up of the numerical
simulations. DG coupled with the time integration MPRK2, and MPMS2 are tested
and the numerical results with $\Delta x=\Delta y=1/48$ are shown in Figures
7-8. Our numerical results agree with previous works.

(a) Contour of Density (b) Contour of Pressure

(c) Cut along $y=0$ (d) Cut along $y=0$

Figure 5: Example4.4 numerical convergence study MPMS2

(a) Contour of Density (b) Contour of Pressure

(c) Cut along $y=0$ (d) Cut along $y=0$

Figure 6: Example4.4 numerical convergence study MPRK2

(a) Colored Contour of Density (b) Contour Line of Density

(c) Colored Contour of Pressure (d) Contour Line of Pressure

Figure 7: Example4.5 detonation diffraction problem using MPRK2

(a) Colored Contour of Density (b) Contour Line of Density

(c) Colored Contour of Pressure (d) Contour Line of Pressure

Figure 8: Example 4.5 detonation diffraction problem using MPMS2

##  5 Concluding remarks

In this paper, we constructed MP MS method up to third order accuracy for
production-destruction equations. Coupled with bound-preserving DG methods,
the scheme we obtained have both conservative and bound-preserving properties.
Numerical experiments have shown the accuracy and effectiveness of the
proposed schemes. In the numerical experiments, we take appropriately small
time step. To the best of our knowledge, there is no modified Patankar
multistep method in the literature. Therefore, our method is a supplement to
the existing method. The multistep methods generally require less function
evaluations in each step than the Runge-Kutta method. Therefore, for ODEs or
PDEs without characteristic speed changing wildly, the multistep method works
better than the Runge-Kutta method. However, if the velocity changes wildly,
the Runge-Kutta method would perform better.

## Conflict of interest

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

## References

  * [1] H. Burchard, E. Deleersnijder, A. Meister, A high-order conservative Patankar-type discretisation for stiff systems of production–destruction equations. Applied Numerical Mathematics, v47 (2003), pp. 1–30
  * [2] N. Chuenjarern, Z. Xu and Y. Yang, High-order bound-preserving discontinuous Galerkin methods for compressible miscible displacements in porous media on triangular meshes, Journal of Computational Physics, v378 (2019), pp.110-128.
  * [3] B. Cockburn and C.-W. Shu, The Runge-Kutta local projection P1-discontinuous- Galerkin finite element method for scalar conservation laws, Mathematical Modelling and Numerical Analysis (M2AN), v25 (1991), pp. 337-361.
  * [4] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework, Mathematics of Computation, v52 (1989), pp. 411-435.
  * [5] B. Cockburn, S.-Y. Lin and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one dimensional systems, Journal of Computational Physics, v84 (1989), pp. 90-113.
  * [6] B. Cockburn, S. Hou and C.-W. Shu, The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case, Mathematics of Computation, v54 (1990), pp. 545-581.
  * [7] B. Cockburn and C.-W. Shu, The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems, Journal of Computational Physics, v141 (1998), pp. 199-224.
  * [8] J. Du and Y. Yang, High-order bound-preserving finite difference methods for multispecies and multireaction detonations, Communications on Applied Mathematics and Computation, accepted.
  * [9] J. Du, C. Wang, C. Qian, Y. Yang, _High-order bound-preserving discontinuous Galerkin methods for stiff multispecies detonation_ , SIAM Journal on Scientific Computing, v41 (2019), pp. B250–B273.
  * [10] J. Du and Y. Yang, Third-order conservative sign-preserving and steady-state-preserving time integrations and applications in stiff multispecies and multireaction detonations, Journal of Computational Physics, v395 (2019), pp. 489-510.
  * [11] S. Gottlieb, C.-W. Shu, and E. Tadmor. _Strong stability-preserving high-order time discretization methods_ , SIAM review, 43 (2001), pp. 89–112.
  * [12] H. Guo and Y. Yang, Bound-preserving discontinuous Galerkin method for compressible miscible displacement problem in porous media, SIAM Journal on Scientific Computing, v39 (2017), pp. A1969– A1990.
  * [13] J. Huang, T. Izgin, S. Kopecz, A. Meister and C.-W. Shu, On the stability of strong-stability-preserving modified Patankar Runge-Kutta schemes. arXiv preprint arXiv:2205.01488 (2022).
  * [14] J. Huang, C.-W. Shu, _Positivity-preserving time discretizations for production-destruction equations with applications to non-equilibrium flows_. Journal of Scientific Computing, v78 (2019), pp. 1181-1839.
  * [15] J. Huang, W. Zhao and C.-W. Shu, A third-order unconditionally positivity-preserving scheme for production-destruction equations with applications to non-equilibrium flows , Journal of Scientific Computing, v79 (2019), pp. 1015-1056.
  * [16] T. Izgin, S. Kopecz and A. Meister, On Lyapunov stability of positive and conservative time integrators and application to second order modified Patankar–Runge–Kutta schemes, ESAIM: Mathematical Modelling and Numerical Analysis, v56(2022), 1053-1080.
  * [17] T. Izgin, S. Kopecz and A. Meister, On the Stability of Unconditionally Positive and Linear Invariants Preserving Time Integration Schemes. arXiv preprint arXiv:2202.11649 (2022).
  * [18] T. Izgin, and P. Öffner. On the Stability of Modified Patankar Methods. arXiv preprint arXiv:2206.07371 (2022).
  * [19] S. Kopecz, A. Meister, _On order conditions for modified Patankar–Runge–Kutta schemes_ , Applied Numerical Mathematics, v123 (2018), pp. 159-179.
  * [20] S. Kopecz, A. Meister, Unconditionally positive and conservative third order modified Patankar–Runge– Kutta discretizations of production–destruction systems. BIT Numer. Math. v58 (2018), pp. 691–728.
  * [21] Y. Lv and M. Ihme, Discontinuous Galerkin method for multicomponent chemically reacting flows and combustion, Journal of Computational Physics, v270 (2014), pp. 105–137.
  * [22] Y. Lv and M. Ihme, High-order discontinuous Galerkin method for applications to multicomponent and chemically reacting flows, Acta Mechanica Sinica, v33 (2017), pp. 486–499.
  * [23] P. Öffner and D. Torlo, Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes. Applied Numerical Mathematics, v153 (2020), pp. 15-34.
  * [24] S. Patankar, Numerical Heat Transfer and Fluid Flow.CRCPress,London(1980)
  * [25] W. H. Reed and T. R. Hill, Triangular mesh methods for the Neutron transport equation, Los Alamos Scientific Laboratory Report LA-UR-73-479, Los Alamos, NM, 1973.
  * [26] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics, v77(2), pp. 439–471.
  * [27] C.-W. Shu, Total-variation-diminishing time discretizations, SIAM Journal on Statistical and Scientific Computing, v9 (1988), pp. 1073–1084.
  * [28] G. Strang, On the construction and comparison of difference schemes, SIAM Journal on Numerical Analysis, v5 (1968), pp. 506–517.
  * [29] C. Wang, X. Zhang, C.-W. Shu, and J. Ning, _Robust high order discontinuous Galerkin schemes for two-dimensional gaseous detonations_ , Journal of Computational physics v231 (2012), pp. 653–665
  * [30] W. Wang, C.-W. Shu, H.C. Yee and B. Sjogreen, High order well-balanced schemes and applications to non-equilibrium flow, Journal of Computational Physics, v228 (2009), pp. 6682-6702.
  * [31] X. Zhang, C.-W. Shu, _Positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations with source terms_ , Journal of Computational physics, v230 (2011), 1238–1248.
  * [32] X. Zhang and C.-W. Shu, _On maximum-principle-satisfying high order schemes for scalar conservation laws_ , Journal of Computational Physics, v229 (2010), 3091-3120.
  * [33] X. Zhang, _On positivity-preserving high order discontinuous Galerkin schemes for compressible Navier-Stokes equations_ , Journal of Computational Physics, v328 (2017), 301-343.