## 2 High-Order Accurate Numerical Method
### Spatial Discretization  
In this work, the discontinuous [PERSON] finite element method is considered to be implemented on each element with polynomial functions of a maximum degree equal to \(p\), which results in the method being \((p+1)\)th order accurate. Let \(d\) be the number of space dimensions, and the number of degrees of freedom per equation per element is \(N_{p}=\prod_{i=1}^{d}(p+i)/i\). Assume that the computational domain \(\Omega\subset R^{d}\) is subdivided into a collection of non-overlapping finite elements \(\{\Omega_{i}\}\) and define the following space for both solution and test function:  
\[\mathbf{V}_{h}^{p}=\Big{\{}\mathbf{v}_{h}\in[L_{2}(\Omega)]^{m}:\left.v_{h}\right|_{ \Omega_{i}}\in\left[V_{p}^{m}\right]\ \ \forall\Omega_{i}\in\Omega\Big{\}}, \tag{11}\]  
where \(V_{p}^{m}\) is the space of vector-valued polynomials of a degree, at most \(p\), and \(m\) is the dimension of the conservative state vector (note: \(m=4\) for present 2D case and \(m=5\) for 3D). Multiplying Equation (9) by a test function \(\mathbf{W}_{h}\in\mathbf{V}_{h}^{p}\) and integrating by parts on each element \(\varOmega_{i}\), we obtain the weak formulation of the DG method:  
\[\int_{\varOmega_{i}}W_{h}^{T}\varGamma\frac{\partial q_{h}}{\partial t}d\varOmega +\int_{\partial\varOmega_{i}}\mathbf{W}_{h}^{T}\mathbf{F}(\mathbf{q}_{h})\cdot\mathbf{n}d \sigma-\int_{\varOmega_{i}}\
abla\mathbf{W}_{h}^{T}\cdot\mathbf{F}(\mathbf{q}_{h})d\varOmega =0\quad\forall\mathbf{W}_{h}\in\mathbf{V}_{h}^{p}, \tag{12}\]  
where \(\mathbf{q}_{h}\) is the approximate solution at element \(\varOmega_{i}\) and can be represented by basis functions \(\phi=\left\{\phi_{j}\right\}\in\mathbf{V}_{h}^{p}\) as follows:  
\[\mathbf{q}_{h}=\sum_{j=1}^{N_{p}}\widetilde{\mathbf{q}}_{i,j}(t)\phi_{j}(\mathbf{x}), \tag{13}\]  
where the unknowns \(\widetilde{\mathbf{q}}_{i,j}(t)\) are the solution coefficients of the element to be solved. According to Equation (12), we have the following system of \(N_{p}\) equations for each \(\varOmega_{i}\)  
\[\int_{\varOmega_{i}}\phi_{j}\mathbf{\Gamma}\frac{\partial q_{h}}{\partial t}d \varOmega+\int_{\partial\varOmega_{i}}\phi_{j}\mathbf{F}(\mathbf{q}_{h})\cdot\mathbf{n}d \sigma-\int_{\varOmega_{i}}\
abla\phi_{j}\cdot\mathbf{F}(\mathbf{q}_{h})d\varOmega=0 \quad 1\leq j\leq N_{p}. \tag{14}\]  
Thus, the unknowns \(\widetilde{\mathbf{q}}_{i,j}(t)\) can be determined through solving system (14). Keep in mind that the flux terms are not uniquely defined at element interface due to discontinuous function approximation. In order to treat this issue properly, the flux term \(\mathbf{F}(\mathbf{q}_{h})\cdot\mathbf{n}\), appearing in the second term of Equation (14), is usually required to be replaced by a numerical flux \(\mathbf{H}\big{(}\mathbf{q}_{h}^{+},\mathbf{q}_{h}^{-},\mathbf{n}\big{)}\), where \(\mathbf{q}_{h}^{+}\) is the neighboring element interface state, \(\mathbf{q}_{h}^{-}\) the internal interface state, and \(\mathbf{n}\) the direction normal to the interface. It is reported that a Roe-type flux function [24] satisfies the flux consistency conditions [25] associated with ensuring the formal accuracy of the scheme, and hence is adopted in the present paper. As mentioned above, the wave speeds of the governing equations were changed by preconditioning, and the related dissipative part of the numerical fluxes \(\mathbf{H}\big{(}\mathbf{q}_{h}^{+},\mathbf{q}_{h}^{-},\mathbf{n}\big{)}\) are supposed to be modified accordingly. Following the work of Reference [21], the resultant numerical flux can be written as  
\[\mathbf{H}\big{(}\mathbf{q}_{h}^{+},\mathbf{q}_{h}^{-},\mathbf{n}\big{)}=\frac{1}{2}\big{(} \mathbf{F}\big{(}\mathbf{q}_{h}^{+}\big{)}+\mathbf{F}\big{(}\mathbf{q}_{h}^{-}\big{)}\big{)} \cdot\mathbf{n}-\frac{1}{2}\sigma\mathbf{\Gamma}|\mathbf{A}_{\mathbf{\Gamma}}|\Delta\mathbf{q} \cdot\mathbf{n}. \tag{15}\]  
It can be noted that the effect of preconditioning is considered in the dissipative part of the last term of the Equation (15).