## 2 High-Order Accurate Numerical Method
### Shock Capturing  
In order to simulate high-speed shocked flows, a kind of shock capturing technique [8] is adopted through adding a small amount of artificial viscosity to the preconditioned governing Equation (9) in order to deal with possible strong shocks, which reads:  
\[\mathbf{\Gamma}\frac{\partial q}{\partial t}+\
abla\cdot\mathbf{F}=\
abla\cdot(\mu \
abla\mathbf{q}), \tag{16}\]  
where \(\
abla\mathbf{q}\) is the gradient of the primitive variable vector and is computed by the BR2 scheme [27] in the present work, and \(\mu\) is the artificial viscosity coefficient. Keep in mind that the artificial viscosity introduced is intended to suppress non-physical oscillations in the vicinities of the strong shocks without affecting the solution accuracy in smooth regions. The smoothness of solution is measured by \(s_{i}=log_{10}S_{i}\) with \(S_{i}=\frac{\int_{\varOmega_{i}}(c_{h}-\hat{c}_{h})^{2}d\varOmega}{\int_{ \varOmega_{i}}c_{h}^{2}d\varOmega}\), in which \(c_{h}=\sum_{i=1}^{N_{p}}c_{i}(t)\phi_{i}(\mathbf{x})=\hat{c}_{h}+\sum_{i=(N_{p-1}+ 1)}^{N_{p}}c_{i}(t)\phi_{i}(\mathbf{x})\) is the selected numerical primitive variable, which could be pressure or density. In the present paper, pressure is used. Following the work of Reference [8], the resultant \(\mu_{i}\) associated with \(\Omega_{i}\) is computed by  
\[\mu_{i}= \begin{cases}0&s_{i}\leq s_{0}-\kappa\\ \frac{\mu_{0}}{2}\Big{(}1+\sin\frac{\pi(s_{i}-s_{0})}{2\kappa}\Big{)}&s_{0}- \kappa<s_{i}\leq s_{0}+\kappa\\ \mu_{0}&s_{i}>s_{0}+\kappa\end{cases} , \tag{17}\]  
where \(\mu_{0},s_{0},\) and \(\kappa\) are the control parameters that are usually chosen empirically by the users. It can be noted that the coefficients for each component in expansion (13) decay very quickly in smooth regions, while the strength of the discontinuity in shock regions will slow down the rate of such decay [8]. In other words, the elemental artificial viscosity coefficient \(\mu_{i}\) is close to zero in smooth regions and increases in shock regions. By adopting such an artificial viscosity coefficient, the shock capturing technique can restrain the spurious numerical oscillation in shock regions, causing almost no degrading of solution accuracy in the smooth regions. However, the above piecewise constant artificial viscosities are discontinuous at the element interfaces, which may still cause small numerical oscillations in the vicinities of the strong shocks [28]. The artificial viscosities, made to be continuous, are usually preferred. In the present work, such continuity is achieved through adopting a post-smoothing process, realized in two steps:  
Step 1. Transferring the elemental viscosity to the vertex of the element. The viscosity \(v_{k}\) at vertex \(k\) is determined by  
\[v_{k}=\sum_{i=1}^{n}w_{i}\mu_{i},\]  
where \(n\) is the number of the adjacent elements with weights \(w_{i}=\frac{|\Omega_{i}|}{\sum_{j=1}^{n}|\Omega_{j}|}\). Step 2. Computing the viscosity at each integral point of the element through a linear interpolation. It can be noted that a kind of smoothing is realized through a weighted average process in Step 1. In so doing, the semi-discrete weak form of the Equation (16) can be obtained from Equation (14) through adding the corresponding artificial viscosity terms, which reads:  
\[ \begin{split}&\int_{\Omega_{i}}\phi_{I}\Gamma\frac{\partial \mathbf{q}_{h}}{\partial t}d\Omega+\int_{\partial\Omega_{i}}\phi_{j}\mathbf{H}\big{(} \mathbf{q}_{h}^{+},\mathbf{q}_{h}^{-},\mathbf{n}\big{)}d\sigma-\int_{\Omega_{i}}\
abla \phi_{j}\cdot\mathbf{F}(\mathbf{q}_{h})d\Omega\\ &-\int_{\partial\Omega_{i}}\phi_{j}\
u\
abla\mathbf{q}_{h}\cdot nd \sigma+\int_{\Omega_{i}}\
abla\phi_{j}\cdot(\
u\
abla\mathbf{q}_{h})d\Omega=0 \quad\quad 1\leq j\leq N_{p}.\end{split}  \tag{18}\]  
Following the treatment of BR2 scheme [27], the additional flux term associated with artificial viscosity in Equation (18) is determined by  
\[\
u\
abla\mathbf{q}_{h}\cdot\mathbf{n}|_{\partial\Omega_{i}}=\left\{