Space-time adaptive ADER-DG finite element method with LST-DG predictor and a posteriori sub-cell WENO finite-volume limiting for simulation of non-stationary compressible multicomponent reactive flows

In this work, a classic mathematical problem formulation has been used to describe non-stationary gas-dynamic flow of multicomponent reactive media. The general properties of the gas-dynamic flow were described using a system of non-stationary one-dimensional Euler equations. The multicomponent composition of the reacting medium and the reactions occurring in it were described by expanding the original system of Euler equations by a system of non-stationary one-dimensional convection-reaction equations. The transfer velocities of the components in the convection-reaction equations were chosen to be the same and equal to the gas-dynamic velocity of the medium (one-speed approximation of the transfer). The density of the components was chosen by introducing into consideration the mass concentrations of the components, while the density of the medium was also determined from the Euler equations. In this case, the source term associated with the energy yield from the reactions occurring in the medium was included in the system of Euler equations on the right side of the energy equation. Thus, the system of non-stationary one-dimensional Euler equations and the system of non-stationary one-dimensional convection-reaction equations formed a quasilinear system of hyperbolic equations for describing non-stationary compressible reacting flows, which takes the following form:

where $\mathbf{U}$ is the vector of conserved values, $\mathbf{F}$ is the flux terms and $\mathbf{S}$ is the source term:

where $\rho$ is the mass density; $u$ is the velocity; $p$ is the pressure; $\varepsilon$ is the total energy density including the thermal $e$ and the kinetic contributions $\varepsilon=e+\frac{1}{2}\rho v^{2}$; $\mathbf{c}^{T}=[c_{1},\ldots,c_{R}]$ is a vector of mass concentrations $c_{k}$ of the component of the reacting medium, which determines the mass fraction of the $k$-th component in the mixture: the density of the $k$-th component can be obtained by the formula $\rho_{k}=\rho c_{k}$, and thus the relation $\sum_{k}c_{k}=1$; $R$ is the amount of components in the mixture; $\mathrm{S}_{e}$ is the source term that determine energy yield associated with reactions occurring in the medium; $\mathbf{S}_{r}$ is the source terms that determines the rates of reactions. The multicomponent and reaction properties of the flow are considered in the form of non-stationary one-dimensional convection-reaction equations

where, as can be seen from the system of equations (1), the transfer velocity of the mixture components and the total mass density are determined from the Euler equations. The reaction rates included in the term $\mathbf{S}_{r}$ were chosen in the form of the law of mass action, where the reaction rate constants could have an arbitrary dependence on the temperature $T$ of the medium. The equation of state of the multicomponent mixture was chosen in the form of the perfect gas equation. The caloric equation of state was specified in the form $p=(\gamma-1)e$, where $\gamma$ was calculated from the ratio of the weighted average, by mass concentrations $c_{k}$, specific heat capacities of components. Empirical values of the constants, $R$ and $\gamma$, which are used in the local approximation of the wide-range equations of state by the equation of state of a perfect gas Zeldovich_Raizer_2002 ; Lunev_2017 ; Anderson_1989 , can be used instead of the classical values.

Thus, in the present work, a computational scheme was used for the simplest mathematical model of unsteady compressible reacting flows, however, it is quite in demand in solving applied problems Oran_Boris_2005 ; Lunev_2017 ; Nagnibeda_2009 ; Anderson_1989 .

The numerical scheme used in this work is based on the space-time adaptive ADER finite-element DG method, which has a high accuracy in resolution the smooth components of the solution. Nonphysical anomalies of the numerical solution that arise in the areas of discontinuities and strong gradients of the solution, associated with the linearity of the basic numerical scheme and Godunov’s theorem, are corrected a posteriori by a sufficiently stable high-precision ADER-WENO finite-volume scheme. An excellent detailed description of this computational scheme is given in the basic works ader_dg_ideal_flows ; ader_dg_dev_1 ; ader_dg_dev_2 ; ader_weno_lstdg ; ader_dg_diss_flows ; ader_dg_ale ; ader_dg_grmhd ; ader_dg_gr_prd ; ader_dg_PNPM ; PNPM_DG ; ader_dg_eff_impl ; fron_phys ; exahype ; ader_dg_hpc_impl_1 ; ader_dg_hpc_impl_2 ; ader_dg_hpc_impl_3 ; ader_dg_hpc_impl_4 of the developers of this method. Details of the internal structure of the ADER-DG are presented in the works ader_dg_dev_1 ; ader_dg_dev_2 ; PNPM_DG . Details of the internal structure of the ADER-DG and ADER-WENO finite volume methods, in which the LST-DG prediction method is used, are presented in the work ader_weno_lstdg . Peculiarities of mathematical formulation and efficiency possibility of implementation are discussed in ader_dg_eff_impl ; fron_phys and exahype ; ader_dg_hpc_impl_1 ; ader_dg_hpc_impl_2 ; ader_dg_hpc_impl_3 ; ader_dg_hpc_impl_4 . The modern state of development of the space-time adaptive ADER finite element DG method with a posteriori correction technique of solutions on sub-cells by the finite-volume limiter using AMR for use on unstructured meshes and using the ALE approach is presented in the works ader_dg_mod_1 ; ader_dg_mod_2 . For this reason, the description of the computational method in this paragraph will be given only briefly enough to understand the general structure of the method used in this particular case.

The space-time adaptive ADER-DG finite element method with LST-DG predictor and a posteriori sub-cell WENO finite-volume limiting involves a sequence of steps ader_dg_ideal_flows ; ader_dg_dev_1 ; ader_dg_dev_2 ; ader_weno_lstdg ; ader_dg_diss_flows ; ader_dg_ale ; ader_dg_grmhd ; ader_dg_gr_prd ; ader_dg_PNPM ; PNPM_DG ; ader_dg_eff_impl ; fron_phys ; exahype ; ader_dg_hpc_impl_1 ; ader_dg_hpc_impl_2 ; ader_dg_hpc_impl_3 ; ader_dg_hpc_impl_4 :

• •

  a LST-DG predictor, using which a local discrete space-time solution in the small is obtained;

• •

  a pure ADER discontinuous Galerkin $\mathbb{P}_{N}\mathbb{P}_{N}$ scheme, using which a preliminary high accuracy solution is obtained;

• •

  a determination of the admissibility of the obtained high accuracy solution and identification of “troubled” cells;

• •

  a recalculation of the solution in “troubled” cells by a stable ADER-WENO finite-volume limiter.

The one-dimensional computational domain $\Omega=[a,b]$ was represented by the mesh $\Omega=\cup_{i}\Omega_{i}$, where $\Omega_{i}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$. The piecewise polynomials DG-representation $\mathbf{U}_{h}(x,t^{n})$ are based on the Legendre interpolation polynomials of degree $N$, specified on the space cell $\Omega_{i}$ (mapped into the space reference element $\omega=[0;1]$):

where $\hat{\mathbf{u}}_{p}^{n}$ are the DG-representation coefficients; $\varphi_{p}=\varphi_{p}(\xi)$ are the basis functions; and index $p\in[0,N]$. The interpolation nodes of the Legendre polynomials $\varphi_{p}(\xi)$ were chosen at the nodal points of the Gauss-Legendre quadrature formula.

The solution on a new $(n+1)$-th time step is obtained using the fully discrete one-step ADER-DG scheme ader_dg_ideal_flows ; ader_dg_dev_1 ; ader_dg_dev_2 ; ader_weno_lstdg ; ader_dg_diss_flows ; ader_dg_ale ; ader_dg_grmhd ; ader_dg_gr_prd ; ader_dg_PNPM ; PNPM_DG ; ader_dg_eff_impl ; fron_phys ; exahype ; ader_dg_hpc_impl_1 ; ader_dg_hpc_impl_2 ; ader_dg_hpc_impl_3 ; ader_dg_hpc_impl_4 , that was chosen in the classical form of projection onto the elements $\varphi_{h}(\xi)$ of the set of basis functions $\left\{\varphi_{h}(\xi)\right\}$:

where $\mathbf{q}_{h}=\mathbf{q}_{h}(x,t)$ is the discrete space-time solution obtained by using the LST-DG-predictor; $\mathbf{\mathfrak{G}}=\mathbf{\mathfrak{G}}(\mathbf{q}_{h}^{-},\mathbf{q}_{h}^{+})$ is the Riemann solver; $\mathbf{q}_{h}^{L}$ and $\mathbf{q}_{h}^{R}$ is the discrete space-time solution in left and right cells; left and right states $(\mathbf{q}_{h}^{-},\mathbf{q}_{h}^{+})$ of the Riemann problem $\mathbf{\mathfrak{G}}$ are chosen as a solution at cell $\Omega_{i}$ interfaces. In this work, the Rusanov flux Rusanov_solver (as noted in the work ader_dg_ideal_flows , the Rusanov flux sometimes referred to as the local Lax-Friedrichs flux) and HLLE Riemann solver HLLE_solver_1 ; HLLE_solver_2 were used as a Riemann solver Toro_solvers_2009 ; HLLE Riemann solver is characterized by high stability when used in high rarefaction ranges HLLE_solver_2 ; Toro_solvers_2009 . In this paper, the author would also like to note a very interesting shock-stable modification of the HLLC Riemann solver with reduced numerical dissipation HLLC_solver_shock_stable_2020 , which may arouse interest in using it with the ADER-DG schemes.

The computational scheme was based on space-time adaptive ADER finite-element discontinuous Galerkin method, which is linear in the sense of Godunov’s theorem. Therefore, when using a DG-representation with a degree $N\geqslant 1$, with a formal order of accuracy not less than $2$, nonphysical oscillations arise in the numerical solution. Nonphysical anomalies of the solution were found in the solution obtained at a new time step, using the physical admissibility criteria (PAC) and the relaxed discrete maximum principle (DMP) in the sense of polynomials (numerical admissibility criteria; NAC). The physical admissibility criteria in this work were chosen in a form similar to the work ader_dg_ideal_flows : density, pressure (more precisely, internal energy) must be positive; the concentrations of all components must be non-negative; the normalization condition $\sum_{k}c_{k}=1$ for the mass concentrations of the medium components must be satisfied (this automatically includes a check for the condition $c_{k}\leqslant 1$); the values calculated on the new time step must not be nan, which can occur due to a serious loss of accuracy. The numerical admissibility criteria, in the form of DMP, without significant changes are taken from works ader_dg_ideal_flows ; ader_dg_diss_flows :

where $V_{i}=\Omega_{i-1}\,\cup\Omega_{i}\,\cup\Omega_{i+1}$, $\mathbf{\delta}$ is chosen to be a solution-dependent tolerance given by

where $\delta_{0}=10^{-4}$ and $\epsilon=10^{-3}$, in accordance with the recommendations of the works ader_dg_ideal_flows ; ader_dg_diss_flows . The inequality (18) is checked separately for each value of the vector of conserved values. The $\min$ and $\max$ functions return vectors of conserved values. Violation of inequality for at least one value leads to violation of the criterion.

Simultaneous verification by the PAC and the NAC for a finite-element cell $\Omega_{i}$ determines that the numerical solution obtained by the space-time adaptive ADER finite-element DG method is admissible, while the cell is marked with indicator $\beta=0$; otherwise, the cell is marked as “troubled”, with indicator $\beta=1$ ader_dg_ideal_flows . In such “troubled” cells, a sub-grid with $N_{s}$ spatial sub-cells $\Omega_{i,k}$ was generated, where the piecewise-constant alternative data representation $\mathbf{v}_{h}=\hat{\mathit{P}}\mathbf{U}_{h}$ was constructed using the reversible conservative $L_{2}$-interpolation procedure with operator $\hat{\mathit{P}}$:

The verification of the PAC and the NAC, in algorithmic and software implementations, was carried out not on the DG-representation $\mathbf{U}_{h}$ of the solution, but immediately on the piecewise-constant alternative data representation $\mathbf{v}_{h}$; for which a small tolerance $\mathbf{\delta}$ was included in the inequalities of the DMP (see details in ader_dg_ideal_flows ; ader_dg_diss_flows ).

The solution in “troubled” cells was recalculated using the stable high-precision ADER-WENO-FV scheme, that was chosen in the classical form ader_dg_ideal_flows ; ader_dg_dev_1 ; ader_dg_dev_2 ; ader_weno_lstdg ; ader_dg_diss_flows ; ader_dg_ale ; ader_dg_grmhd ; ader_dg_gr_prd ; ader_dg_PNPM ; PNPM_DG ; ader_dg_eff_impl ; fron_phys ; exahype ; ader_dg_hpc_impl_1 ; ader_dg_hpc_impl_2 ; ader_dg_hpc_impl_3 ; ader_dg_hpc_impl_4 :

where the introduced notation corresponds to the notation for the formula (17); $\tau^{n}=t^{n+1}-t^{n}$ is the time step. The main difference is that the discrete space-time solution $\mathbf{q}_{h}(x,t)$ is obtained by the LST-DG-predictor based on the WENO-reconstruction $\mathbf{w}_{h}(x,t^{n})$ as an initial condition, instead of DG-representation $\mathbf{U}^{n}_{h}$. The resulting finite-volume solution, in the form of the piecewise-constant alternative data representation $\mathbf{v}^{n+1}_{h}$, is converted into a DG-representation $\mathbf{U}^{n+1}_{h}$ using the inverse transformation specified by the suitable high order accurate reconstruction operator $\hat{\mathit{R}}$, which determines the relation $\mathbf{U}^{n+1}_{h}=\hat{\mathit{R}}\mathbf{v}^{n+1}_{h}$. The suitable high order accurate reconstruction operator $\hat{\mathit{R}}$ is implemented using the least squares method implemented using a pseudo-inverse matrix.

It is important to note that the operators $\hat{\mathit{P}}$ and $\hat{\mathit{R}}$ satisfy the reversibility condition $\hat{\mathit{R}}\circ\hat{\mathit{P}}=1$, however, this condition is satisfied only in one direction: direct transformation is reversible $\mathbf{U}_{h}\rightarrow\mathbf{v}_{h}\rightarrow\mathbf{U}_{h}$, but the reverse transformation is irreversible, in the general case, $\mathbf{v}_{h}\rightarrow\mathbf{U}_{h}\rightarrow\mathbf{v}^{*}_{h}\neq\mathbf{v}_{h}$; so there is the relation $\hat{\mathit{P}}\circ\hat{\mathit{R}}\neq 1$. This is due to the fact that the amount of sub-cells $N_{s}$ of sub-grid in cell is greater than the amount of coefficients $N+1$ of the DG-representation, and the system of equations for obtaining the DG-representation $\mathbf{U}_{h}$ from the piecewise-constant representation $\mathbf{v}_{h}$ is overdetermined (thus, it is solved by the least squares method, using a pseudo-inverse matrix), and the system of equations for obtaining piecewise-constant representation $\mathbf{v}_{h}$ from a DG-representation $\mathbf{U}_{h}$ is underdetermined. Therefore, in algorithmic and software implementations, each cell is associated with the status (such as is_dg_the_primary_representation) that determines which of the DG-representations of the solution at the time step $t^{n}$ is primary — whether the DG-representation $\mathbf{U}^{n}_{h}$ is obtained as a result of using the space-time adaptive ADER finite-element DG scheme and is its solution, or the DG-representation $\mathbf{U}^{n}_{h}$ is the result of a transformation from a piecewise-constant representation $\mathbf{v}_{h}$, using an operator $\hat{\mathit{R}}$. This is an important aspect of the implementation, since the next time step $t^{n+1}$ of ADER-DG scheme may result in non-admissible DG-representation $\mathbf{U}^{n+1}_{h}$ (“troubled” cell) of the solution, so it will be necessary to recalculate the solution by the ADER-WENO-FV scheme, which will require a solution $\mathbf{v}^{n}_{h}$ at this time step $t^{n}$. Depending on which solution $\mathbf{v}^{n}_{h}$ is primary at this time step $t^{n}$, it will be used as a piecewise-constant representation for the ADER-WENO-FV scheme. Thus, if the cell remains “troubled” for several time steps, all these time steps will be performed only by the ADER-WENO-FV scheme, and there will be no problem with the irreversibility of the transformation operators $\hat{\mathit{P}}\circ\hat{\mathit{R}}\neq 1$.

The WENO-reconstruction in high-precision ADER-WENO-FV scheme was carried out according to the piecewise-constant alternative data representation, for each sub-cell, based on the expansion in basis functions (similar to the representation (16), however, in this case, mapped into the space reference element $\omega=[0;1]$ is performed separately for each sub-cell $\Omega_{i,k}$; so each “troubled” cell $\Omega_{i}$ contains $N_{s}$ separate representations of the WENO-reconstructions):

where $\hat{\mathbf{w}}_{p}^{n}$ are the WENO-reconstruction coefficients. A detailed description of the reconstruction procedure is presented in the work ader_weno_lstdg . The piecewise constant representation used to perform the reconstruction is taken from the piecewise constant representation of the given problem cell and two neighboring cells — “ghost” subgrids are built in neighboring cells. Technically, in algorithmic and software implementations, the piecewise-constant representation of two neighboring cells is already computed, whether they are “troubled” or not; the piecewise-constant representation calculation is still performed for the verification of the PAC and the NAC of the cells; at the same time, NAC still requires the calculation of the piecewise-constant representation in neighboring cells. The resulting WENO-reconstruction $\mathbf{w}_{h}(x,t^{n})$ is used as an initial condition in the LST-DG-predictor to obtain a discrete space-time solution $\mathbf{q}_{h}(x,t)$ to be used in the high-precision ADER-WENO-FV scheme (21) for each sub-cell $\Omega_{i,k}$.

ADER paradigm and the solution of a Generalized Riemann Problem (GRP) in this computational scheme are based on the use of a local discrete space-time solution $\mathbf{q}_{h}(x,t)$. The space-time adaptive ADER-DG scheme (17) and ADER-WENO-FV scheme (21) use the discrete space-time solution $\mathbf{q}_{h}(x,t)$ in the space-time cell $\Omega_{i}\times[t^{n},t^{n+1}]$ for ADER-DG scheme and in the space-time cell $\Omega_{i,k}\times[t^{n},t^{n+1}]$ for ADER-WENO-FV scheme (mapped into the space-time reference element $[0;1]\times[0;1]$). The discrete space-time solution $\mathbf{q}_{h}(x,t)$ is chosen in the form ader_dg_ideal_flows ; ader_dg_dev_1 ; ader_dg_dev_2 ; ader_weno_lstdg ; ader_dg_diss_flows ; ader_dg_ale ; ader_dg_grmhd ; ader_dg_gr_prd ; ader_dg_PNPM ; PNPM_DG ; ader_dg_eff_impl ; fron_phys ; exahype ; ader_dg_hpc_impl_1 ; ader_dg_hpc_impl_2 ; ader_dg_hpc_impl_3 ; ader_dg_hpc_impl_4 :

where basis functions $\Phi_{\mathfrak{p}}(\xi,\tau)$ are chosen as tensor-products of Lagrange interpolation polynomials $\varphi_{n}(\xi)$ with degree $N$, which were also used in the DG-representation (16) and WENO-reconstruction (22):

where tensor index $\mathfrak{p}=(p_{0},p_{1})$, and $0\leqslant p_{0},\ p_{1}\leqslant N$. In this case, the following method of solving the GRP is used: first evolve the data locally in the small inside each element and then interact the evolved data at the element interfaces via a classical Riemann solver PNPM_DG ; ADER_GRP_LST_DG_1 ; ADER_GRP_LST_DG_2 ; ADER_GRP_LST_DG_3 ; ADER_GRP_LST_DG_4 ; ADER_GRP_LST_DG_5 (references in work ader_dg_ideal_flows ).

The LST-DG predictor for obtaining the discrete space-time solution $\mathbf{q}_{h}(\xi(x),\tau(t))$ in the space-time control volume $\Omega\times[t^{n},t^{n+1}]$ is chosen ader_dg_ideal_flows ; ader_dg_dev_1 ; ader_dg_dev_2 ; ader_weno_lstdg ; ader_dg_diss_flows ; ader_dg_ale ; ader_dg_grmhd ; ader_dg_gr_prd ; ader_dg_PNPM ; PNPM_DG ; ader_dg_eff_impl ; fron_phys ; exahype ; ader_dg_hpc_impl_1 ; ader_dg_hpc_impl_2 ; ader_dg_hpc_impl_3 ; ader_dg_hpc_impl_4 as

where summation is assumed over the tensor index $\mathfrak{q}=(q_{0},q_{1})$ and scalar index $k$ (Einstein notation); $\hat{\mathbf{u}}_{k}$ is the representation coefficients of the initial conditions (at this time step $t^{n}$): for the case of the ADER-DG scheme $\hat{\mathbf{u}}_{k}$ is the DG-representation coefficients — $\hat{\mathbf{u}}_{k}=\hat{\mathbf{U}}^{n}_{k}$; for the case of the ADER-WENO-FV scheme $\hat{\mathbf{u}}_{k}$ is the WENO-reconstruction coefficients — $\hat{\mathbf{u}}_{k}=\hat{\mathbf{w}}^{n}_{k}$. Functions $\tilde{\mathbf{F}}$ and $\tilde{\mathbf{S}}$ are rescaled flux $\mathbf{F}$ and source $\mathbf{S}$ terms due to mapping $\Omega\times[t^{n},t^{n+1}]\leftrightarrow[0;1]\times[0;1]$; it is important to note here that in the cases of the ADER-DG scheme and the ADER-WENO-FV scheme, the spatial step of the cells differs — the cell step $h$ and the sub-cell step $h_{s}$ differ by $N_{s}$ times; so the rescaling factors will be different. The LST-DG predictor is used to obtain the discrete space-time solution $\mathbf{q}_{h}(\xi(x),\tau(t))$ in each cell $\Omega_{i}$ in the case of a ADER-DG scheme, as well as in each subcell $\Omega_{i,k}$ of each “troubled” cell $\Omega_{i}$ (and also, in certain cases, in nearby sub-cells to the “troubled” cells, to calculate the left flux terms in the left sub-cell and right flux terms in the right sub-cell, if neighboring cells are not “troubled” and WENO-reconstruction is not required in them). LST-DG predictor (25) takes a compact notation Zanotti_lectures_2016 ; Jackson_2017

where $\mathbb{K}_{\mathfrak{p}\mathfrak{q}}^{\tau}$, $\mathbb{K}_{\mathfrak{p}\mathfrak{q}}^{\xi}$, $\mathbb{F}_{\mathfrak{p}\mathfrak{h}}$ and $\mathbb{M}_{\mathfrak{p}\mathfrak{q}}$ are the matrices of integrals of $\varphi_{p}(\xi)$ and $d\varphi_{p}(\xi)/d\xi$, which have a convenient index structure due to the orthogonality of the set of basis functions, and can be pre-computed in the code. The process of obtaining a solution $\hat{\mathbf{q}}_{\mathfrak{p}}$ to system of nonlinear algebraic equations (26) with stiff terms can be organized using one of the iterative scheme ader_dg_ideal_flows ; ader_dg_dev_1 ; ader_dg_dev_2

where matrix products on the right side have certain properties: the matrix $[(\mathbb{K}^{\tau})^{-1}]_{\mathfrak{p}\mathfrak{q}}\mathbb{K}_{\mathfrak{q}\mathfrak{r}}^{\xi}$ has zero eigenvalues Zanotti_lectures_2016 ; Jackson_2017 , therefore, in the case of a problem without source terms, a simple iterative procedure is converged to a fixed point (Banach fixed point theorem); the matrix $[(\mathbb{K}^{\tau})^{-1}]_{\mathfrak{p}\mathfrak{q}}\mathbb{F}_{\mathfrak{q},k}$ also has an interesting property:

where $\delta_{p_{1},k}$ is the Kronecker delta symbol and $\mathbb{I}_{p_{0}}=1$ is just a one; which is a consequence of the correspondence principle: in the case of null fluxes terms $\mathbf{F}\equiv 0$ and null sources terms $\mathbf{S}\equiv 0$, the solution $\mathbf{q}(x,t)$ should not depend on time $t$ and for the representation coefficients the solution $\hat{\mathbf{q}}_{p_{0},k}=\hat{\mathbf{u}}_{p}$ can be obtained. In the present work, the system of equations was solved by an iterative method, the left side of the system of equations was solved by an internal iterative procedure. In the case of a weakly stiff terms of the problem, the internal iterative procedure was a simple fixed point method. In the case of a strong stiff terms of the problem, the internal iterative procedure was implemented by Newton’s method.

Despite the fact that the methods used in this paper are locally implicit, the resulting ADER-DG and ADER-WENO-FV schemes are explicit (local implicitness of the LST-DG predictor step). Therefore, the Courant-Friedrichs-Lewy (CFL) condition is imposed on the time discretization step $\tau^{n}=t^{n+1}-t^{n}$, which takes the form ader_dg_ideal_flows ; ADER_DG_time_step_1 ; ADER_DG_time_step_2 :

where $|\lambda_{max}|$ is the maximum signal velocity, $h$ is the spatial characteristic mesh size and $h_{s}=h/N_{s}$ is the spatial characteristic sub-grid size. If $N_{s}=2N+1$ then the discretization time steps $\tau^{n}_{DG}=\tau^{n}_{FV}$ ader_dg_ideal_flows , therefore, in this work, the condition $N_{s}=2N+1$ was chosen.

The software implementation was made using the C++ programming language. Multithreading execution was provided using the OpenMP standard. Parallelization by threads was performed primitively — loops were parallelized using OpenMP pragmas, with the allocation of critical code sections and private and shared variables. Multiprocessing execution and inter-process communication was provided within the framework of the MPI standard. Parallelization by processes was also performed primitively — the full spatial mesh was divided into ranges, in each of which calculations were performed by separate processes; at each step, data was transferred from the boundary cells of the selected ranges. Load balancing of processes and threads in the parallelized version of the program was not performed, which in some cases led to a slight idle time of computing resources — the iterative processes of obtaining discrete space-time solutions by the LST-DG predictor for various initial states often converged for a different number of iterations, but this difference was not significant in terms of computation time.