Stability of Nonlinear Convection–Diffusion–Reaction Systems in Discontinuous Galerkin Methods

Abstract

In this work we provide an extension of the classical von Neumann stability analysis for high-order accurate discontinuous Galerkin methods applied to generalized nonlinear convection–reaction–diffusion systems. We provide a partial linearization under which a sufficient condition emerges that guarantees stability in this context. The stability behavior of these systems is then closely analyzed relative to Runge–Kutta Chebyshev (RKC) and strong stability preserving (RKSSP) temporal discretizations over a nonlinear system of reactive compressible gases arising in the study of atmospheric chemistry.

Appendix

1.1 Computation of the $G^{\star }$ Operator:

In this section we briefly outline the definition of the operators appearing in an RK scheme. Note that at the ith stage of an RK scheme we can write ${\varvec{u}}^{{(i)}} = G_i u^{(0)}$ where $G_{i}$ is an i th degree polynomial. Identifying the $G_i$ polynomials with the s-vectors of their coefficients, we have for the RKSSP scheme:

$$\begin{aligned} G_{0}= (1,0, \cdots ,0)^T, \quad G_{1} = (0,1,0, \cdots ,0)^T, \end{aligned}$$

and for $2 \le i \le s$,

$$\begin{aligned} G_{i,k} = {\left\{ \begin{array}{ll} \displaystyle \sum _{l=0}^{i-1} \alpha _{i-1,l}G_{l,0} -\frac{\alpha _{0,0}}{\beta _{0,0}}\beta _{i-1,l}G_{l,0} \quad &{}\text {for } k = 0,\\ \displaystyle \sum _{l=0}^{i-1} [\alpha _{i-1,l}G_{l,k} -\frac{a_{0,0}}{\beta _{0,0}}\beta _{i-1,l}G_{l,0} +\frac{1}{\beta _{1,1}}\beta _{i-1,l}G_{l,k-1}] \quad &{}\text {for } 2\le k \le i,\\ \displaystyle 0 \quad &{}\text {for } i+1\le k \le s.\\ \end{array}\right. } \end{aligned}$$

(6.1)

For the RKC scheme we have,

$$\begin{aligned} G_{i,k} = {\left\{ \begin{array}{ll} \displaystyle (1 - \mu _{i} - \nu _{i}) + \mu _{i}G_{i-1,0} + \nu _{i}G_{i-2,0} - \frac{1}{\tilde{\mu }_{1}}(\tilde{\mu }_{i}G_{i-1,0}+\tilde{\gamma }_{i}) \quad &{}\text {for } k = 0, \\ \displaystyle \mu _{i}G_{i-1,1} + \nu _{i}G_{i-2,1} + \frac{1}{\tilde{\mu }_{1}}(\tilde{\gamma }_{i} + \tilde{\mu }_{i}G_{i-1,0} -\tilde{\mu }_{i}G_{i-1,1}) \quad &{}\text {for } k = 1,\\ \displaystyle \mu _{i}G_{i-1,k} + \nu _{i}G_{i-2,k} + \frac{1}{\tilde{\mu }_{1}}(\tilde{\mu }_{i}G_{i-1,k-1} - \tilde{\mu }_{i}G_{i-1,k}) \quad &{}\text {for } 2\le k \le i,\\ \displaystyle 0 \quad &{}\text {for } i+1\le k \le s.\\ \end{array}\right. } \end{aligned}$$

(6.2)

Finally, we use $G^{\star }=G_s$ as the definition of the $G^\star$ operator.

1.2 Simple illustrative example: Burger’s equation

Consider the one-dimensional Burger’s equation

$$\begin{aligned} u_t + \frac{1}{2}(u^2)_{x} = 0, \end{aligned}$$

with periodic boundary data. The weak form of such a system satisfies (2.7), such that:

$$\begin{aligned} \frac{d}{dt}\int _{\Omega _{e_{i}}} \varphi _{h}u_{h} \,dx = F_{i} - \int _{\Omega _{e_{i}}} f(u_{h})\partial _{x}\varphi _{h}\,dx. \end{aligned}$$

(6.3)

Note here that $f = u^2/2$.

Now, to proceed a numerical flux is specified that is associated to $F_{i}$. For illustration, let us choose the local Lax-Friedrich’s flux

$$\begin{aligned} \hat{f}_{\mathrm {LLF}} = \frac{1}{2} (f|_{\Gamma _{i\ell }}+f|_{\Gamma _{\ell i}} - \alpha (u_{h}|_{\Gamma _{i\ell }}-u_{h}|_{\Gamma _{\ell i}}), \end{aligned}$$

such that,

$$\begin{aligned} F_{i} = \sum _{\ell \in \Xi (i)}\int _{\Gamma _{i\ell }} (\hat{f}_{\mathrm {LLF}})_{i\ell }\varphi _{h}|_{\Gamma _{i\ell }} n dS. \end{aligned}$$

Now that the flux has been specified, we proceed with the partial linearization from (3.3). The first step in doing so is to rewrite the representation with respect to the lth degree of freedom and replace f with the partial linearization in the volume integral, such that:

$$\begin{aligned} \frac{d}{dt} u_{l}^{i}(t) = \frac{2l+1}{h} \bigg \{ F_{il} - \int _{\Omega _{e_{i}}} (J_uf) u_{h} \partial _{x}\varphi _{l}(x)\,dx \bigg \}. \end{aligned}$$

(6.4)

Next the partial linearization is performed on the flux representations, such that:

$$\begin{aligned} F_{il} = \frac{1}{2} \sum _{\ell \in \Xi (i)}\int _{\Gamma _{i\ell }} \left( f|_{\Gamma _{i\ell }}+f|_{\Gamma _{\ell i}} - \alpha (u_{h}|_{\Gamma _{i\ell }}-u_{h}|_{\Gamma _{\ell i}})\right) n_{i\ell } \varphi _{l} dS, \end{aligned}$$

(6.5)

where the Jacobian in this case is trivially $(J_{u}f)|_{\Gamma _{i\ell }} = u_{h}|_{\Gamma _{i\ell }}$, and $\alpha = \max (u_{h}|_{\Gamma _{i\ell }},u_{h}|_{\Gamma _{\ell i}})$. This means we now have the system:

$$\begin{aligned} \frac{d}{dt} u_{l}^{i} = \frac{2l+1}{h} \bigg \{ \frac{1}{2} \sum _{\ell \in \Xi (i)}\int _{\Gamma _{i\ell }} ( [(J_{u}f)u_{h}]|_{\Gamma _{i\ell }}&+[(J_{u}f)u_{h}]|_{\Gamma _{\ell i}} - \alpha (u_{h}|_{\Gamma _{i\ell }}-u_{h}|_{\Gamma _{\ell i}})) \varphi _{l} n_{i\ell } dS \nonumber \\&- \int _{\Omega _{e_{i}}} (J_uf) u_{h} \partial _{x}\varphi _{l}(x)\,dx \bigg \}, \end{aligned}$$

(6.6)

so that using the same argument from Sect. 3.2, the flux can now be split relative to contributions determined by the stencil generated by cell “ownership,” as in (3.4).

To make this explicit, let us consider only the flux contribution,

$$\begin{aligned} \frac{2l+1}{2h} \sum _{\ell \in \Xi (i)}\int _{\Gamma _{i\ell }} ( [(J_{u}f)u_{h}]|_{\Gamma _{i\ell }} +[(J_{u}f)u_{h}]|_{\Gamma _{\ell i}} - \alpha (u_{h}|_{\Gamma _{i\ell }}-u_{h}|_{\Gamma _{\ell i}})) \varphi _{l} n_{i\ell } dS. \end{aligned}$$

This term is rewritten relative to cell ownership, such that:

$$\begin{aligned}&\frac{2l+1}{2h} \left( \sum _{k,m}^{n_{p}}F'_{km}u_{m}\varphi _{l}\right) \bigg |_{\Gamma _{i\ell }} + \frac{2l+1}{2h} \left( \sum _{k,m}^{n_{p}}F'_{km}u_{m}\varphi _{l}\right) \bigg |_{\Gamma _{\ell i}} \nonumber \\&= \left( \frac{2l+1}{2h} \sum _{\ell \in \Xi (i)}\int _{\Gamma _{i\ell }} (J_{u}f|_{\Gamma _{i\ell }}-\alpha )u_{h}|_{\Gamma _{i\ell }} \varphi _{l} n_{i\ell } dS \right) \nonumber \\&\quad + \left( \frac{2l+1}{2h} \sum _{\ell \in \Xi (i)}\int _{\Gamma _{\ell i}} (J_{u}f|_{\Gamma _{\ell i}}+\alpha )u_{h}|_{\Gamma _{\ell i}} \varphi _{l} n_{i\ell } dS \right) . \end{aligned}$$

(6.7)

This makes it clear now that the $F'_{km}$ denotes an evaluated form of $J_{u}f +\alpha$.

As stated in Sect 3.2, $F'_{jkm}$ is an $N\times (n_{p}+1)^{2}$ tensor, meaning for a single component flow such as Burger’s, we have that $F'_{km}$ is a $(n_{p}+1)\times (n_{p}+1)$ matrix, in this case defined by:

$$\begin{aligned} F'_{km} = \sum _{\ell \in \Xi (i)}\int _{\Gamma _{i\ell }} (J_{u}f|_{\Gamma _{i\ell }}+\alpha )_{k}\varphi _{m} n_{i\ell } dS. \end{aligned}$$

To see this we recognize that $(J_{u}f) = u_{h} = \sum _{l=0}^{n_{p}} u_{l}\varphi _{l}$ is a vector of length $n_{p}+1$, as is $\alpha$. Thus simply factoring out the basis function from the solution, and evaluating at the support points yields:

$$\begin{aligned}&\left( \frac{2l+1}{2h} \sum _{\ell \in \Xi (i)}\int _{\Gamma _{i\ell }} (J_{u}f|_{\Gamma _{i\ell }}-\alpha )u_{h}|_{\Gamma _{i\ell }} \varphi _{l} n_{i\ell } dS \right) \nonumber \\&\quad + \left( \frac{2l+1}{2h} \sum _{\ell \in \Xi (i)}\int _{\Gamma _{i\ell }} (J_{u}f|_{\Gamma _{\ell i}}+\alpha )u_{h}|_{\Gamma _{\ell i}} \varphi _{l} n_{i\ell } dS \right) \nonumber \\&\quad = \left( \frac{2l+1}{2h} \sum _{\ell \in \Xi (i)}\sum _{k,m}^{n_q}\int _{\Gamma _{i\ell }} (J_{u}f|_{\Gamma _{i\ell }}-\alpha )_{k}(\varphi _{m}u_{m})|_{\Gamma _{i\ell }} \varphi _{l} n_{i\ell } dS \right) \nonumber \\&\quad + \left( \frac{2l+1}{2h} \sum _{\ell \in \Xi (i)}\sum _{k,m}^{n_q}\int _{\Gamma _{\ell i}} (J_{u}f|_{\Gamma _{\ell i}}+\alpha )_{k} (\varphi _{m} u_{m})|_{\Gamma _{\ell i}} \varphi _{l} n_{i\ell } dS \right) . \end{aligned}$$

(6.8)

Notice that (6.6), after substituting the explicit form (6.8) in for the surface terms, comprise the matrices in (3.11) that end up forming the $G_{j}$ matrix from (3.12).

The factored contribution $F'_{km}$, being an $(n_p +1) \times (n_p + 1)$ matrix, clearly has dependencies on the degrees of freedom (e.g. the spatial nodes) of the cell that owns it. However, the assumption of the partial linearization (3.2) is that within the cell distance $\Delta x < h$, these spatial dependencies can be decoupled (by scaling arguments) from the dependencies in the Fourier expansion and subsequent shift. Numerical examples seem to reinforce the validity of this assumption.