## 2 Numerical Methods
### Flux Reconstruction Method  
In the FR approach, first the computational domain \(\Omega\) is divided into \(N\) non-overlapping elements \(\Omega_{h},h=1,\ldots,N\). Taking non-uniformity into consideration, a coordinate transformation is applied to map the physical coordinates \(\mathbf{x}=(x,y,z)\) to the computational coordinates \(\xi\)=\((\xi,\eta,\zeta)\) in a standard hexahedral computational element \(\Omega_{h}=[-1,1)^{3}\). The transformation Jacobian matrix \(\mathbf{J}\) is defined as \(J_{ij}=\partial x_{i}/\partial\xi_{j}\). Then for the following governing equations:\[\frac{\partial U}{\partial t}+\
abla\cdot\mathbf{F}=0\,, \tag{6}\]  
the transformed governing equations take the following form:  
\[\frac{\partial\hat{U}}{\partial t}+\hat{\
abla}\cdot\mathbf{F}=0\,, \tag{7}\]  
where  
\[\hat{U} =|\mathbf{J}|U\,, \tag{8}\] \[\hat{\mathbf{F}} =|\mathbf{J}|\mathbf{J}^{-1}\mathbf{F}\,,\] \[\hat{\
abla} =(\partial/\partial\xi,\partial/\partial\eta,\partial/\partial \zeta)\,.\]  
In the standard computational element, a set of points (solution points, SPs) are defined, on which the values of the flow field are stored and updated by solving the equations above. In this paper, the SPs are located at the Legendre-Gauss points. By adopting \(N_{p}\) solution points, for example, in \(\xi\) direction, the 1-D Lagrange basis polynomial of \(N_{p}-1\) degree is  
\[l_{i}(\xi)=\prod_{j=1,j\
eq i}^{N_{p}}\left(\frac{\xi-\xi_{j}}{\xi_{i}-\xi_{j} }\right)\,,\quad i=1,\dots,N_{p}\,. \tag{9}\]  
Once given the values at the SPs, the solution \(\hat{U}\) in the standard computational element \(\Omega_{h}\) is approximated by \(\hat{U}_{h}\), which is reconstructed by a tensor product of the three polynomials in \(\xi\), \(\eta\) and \(\zeta\) direction, respectively:  
\[\hat{U}_{h}(\xi,\eta,\zeta)=\sum_{k=1}^{N_{p}}\sum_{j=1}^{N_{p}}\sum_{i=1}^{N_ {p}}\hat{U}_{h|i,j,k}\cdot l_{i}(\xi)\cdot l_{j}(\eta)\cdot l_{k}(\xi)\,, \tag{10}\]  
where \(\hat{U}_{h|i,j,k}\) is the current value of the approximation function of solution at SP located by \((i,j,k)\) in the computational element \(\Omega_{h}\). Next, the polynomials of fluxes are constructed for time advancement of the solution in element \(\Omega_{h}\). First, the discontinuous polynomials of fluxes in each element are constructed:  
\[\hat{\mathbf{F}}_{h}^{D} =\{\hat{f}_{h}^{D},\hat{g}_{h}^{D},\hat{h}_{h}^{D}\}\,, \tag{11}\] \[\hat{f}_{h}^{D}(\xi,\eta,\zeta) =\sum_{k=1}^{N_{p}}\sum_{j=1}^{N_{p}}\sum_{i=1}^{N_{p}}\hat{f}_{h |i,j,k}\cdot l_{i}(\xi)\cdot l_{j}(\eta)\cdot l_{k}(\xi)\,,\] \[\hat{g}_{h}^{D}(\xi,\eta,\zeta) =\sum_{k=1}^{N_{p}}\sum_{j=1}^{N_{p}}\sum_{i=1}^{N_{p}}\hat{g}_{h |i,j,k}\cdot l_{i}(\xi)\cdot l_{j}(\eta)\cdot l_{k}(\xi)\,,\] \[\frac{\partial\hat{U}}{\partial t}+\hat{\
abla}\cdot\mathbf{F}=0\,, \tag{7}\]0  
where superscript D denotes discontinuous, \(\hat{f}_{h|i,j,k}\), \(\hat{g}_{h|i,j,k}\), \(\hat{h}_{h|i,j,k}\) are the current value of the approximation fuction of flux at SP located by \((i,j,k)\) in the computational element \(\Omega_{h}\), along \(\xi\), \(\eta\) and \(\zeta\) direction, respectively. Note that the polynomials \(\hat{f}_{h}^{D}\), \(\hat{g}_{h}^{D}\), \(\hat{h}_{h}^{D}\) are piecewise and might not be continuous across the interfaces between the elements. In the FR approach, the discontinuous polynomials of fluxes are then made continuous by introducing correction functions. The corrected fluxes are then defined as:\[