## 2 Methods  
### Spatial Discretization  
The differential form of a conservation law system in a three-dimensional space could be written as  
\[\partial_{t}\underline{U}+\partial_{x}\underline{F}^{x}+\partial_{y} \underline{F}^{y}+\partial_{z}\underline{F}^{z}=\underline{Q} \tag{1}\]  
in which  
* \(\partial_{t},\partial_{x},\partial_{y},\partial_{z}\) represent \(\frac{\partial}{\partial t},\frac{\partial}{\partial x},\frac{\partial}{ \partial y},\frac{\partial}{\partial z}\), respectively;
* \(\underline{U}\) is a \(K\times 1\) matrix of unknowns, each row of which is a scalar function depending on position \(\vec{x}\) and time \(t\);
* \(\underline{F}^{\mu}\) (where \(\mu=x,y,z\)) is also a \(K\times 1\) matrix, which is the dot-product of the flux \(\underline{\vec{F}}\) (whose value depends on \(\underline{U}\)) and \(\vec{e}_{\mu}\) (which is an unit vector along the positive direction of the \(\mu\)-axis);
* \(\underline{Q}\) is a \(K\times 1\) matrix of \(0\)'s. This differential equation could be turned into an integral equation, by multiplying both sides with an arbitrary function \(V(\vec{x})\) and integrating the product on an arbitrary control volume \(\Omega\): \[\int_{\Omega}\Big{(}\partial_{t}\underline{U}+\
abla\cdot\underline{\vec{F}} \Big{)}V=\underline{Q}.\]  
To weaken the smoothness requirements on \(\underline{\vec{F}}\), we apply integral by parts and [PERSON]'s divergence theorem to it, which will lead to the weak form of Equation (1):  
\[\int_{\Omega}\Big{(}\partial_{t}\underline{U}-\underline{\vec{F}}\cdot\
abla V \Big{)}+\int_{\partial\Omega}\Big{(}\vec{v}\cdot\underline{\vec{F}}\Big{)}V =\underline{Q}, \tag{2}\]  
where \(\vec{v}\) is the outer normal unit vector of the control surface \(\partial\Omega\) (which is the boundary of \(\Omega\)). To introduce spatial discretization for Equation (2), we choose the linear space spanned by polynomials up to the \(p\)-th degree over \(\Omega\), denoted as \(\mathcal{V}^{p}(\Omega)\), as the approximation space. Let \(\underline{\phi}(\vec{x})=\big{[}\phi_{1}(\vec{x})\cdots\phi_{L}(\vec{x}) \big{]}\) be a basis of \(\mathcal{V}^{p}(\Omega)\), in this article, we choose  
\[\phi_{1}(\vec{x})=1,\qquad \begin{bmatrix}\phi_{2}(\vec{x})\\ \phi_{3}(\vec{x})\\ \phi_{4}(\vec{x})\end{bmatrix} = \begin{bmatrix}x-x_{0}\\ y-y_{0}\\ z-z_{0}\end{bmatrix} ,\qquad \begin{bmatrix}\phi_{5}(\vec{x})\\ \phi_{6}(\vec{x})\\ \phi_{7}(\vec{x})\\ \phi_{8}(\vec{x})\\ \phi_{9}(\vec{x})\\ \phi_{10}(\vec{x})\end{bmatrix} = \begin{bmatrix}(x-x_{0})(x-x_{0})\\ (x-x_{0})(y-y_{0})\\ (x-x_{0})(z-z_{0})\\ (y-y_{0})(y-y_{0})\\ (y-y_{0})(z-z_{0})\\ (z-z_{0})(z-z_{0})\end{bmatrix} ,\qquad\ldots,\]  
in which \((x_{0},y_{0},z_{0})\) is the geometric center of \(\Omega\). Then \(\underline{U}\) and \(V\) could be approximated as  
\[\underline{U}(\vec{x},t)\approx\underline{U}^{h}(\vec{x},t)=\sum_{l=1}^{L} \underline{\underline{U}}_{l}(t)\,\phi_{l}(\vec{x}),\qquad V(\vec{x})\approx V ^{h}(\vec{x})=\sum_{l=1}^{L}\hat{V}_{l}\,\phi_{l}(\vec{x}),\]  
where each \(\underline{\underline{U}}_{l}(t)\) is a \(K\times 1\) matrix of temporal functions, and each \(\hat{V}_{l}\) is a constant number (which is arbitrary since \(V(\vec{x})\) is arbitrary). Substitute them into the weak form (Equation (1)), we obtain \[\sum_{l}\hat{V}_{l}\left[\sum_{k}\left(\int_{\Omega}\phi_{l}\phi_{k}\right)\frac{ \mathrm{d}\underline{\hat{U}}_{k}}{\mathrm{d}t}+\int_{\Omega}(\
abla\phi_{l}) \cdot\bar{E}\left(\underline{U}^{h}\right)+\oint_{\partial\Omega}\phi_{l}\, \underline{F}^{\
u}\left(\underline{U}^{h}_{l},\underline{U}^{h}_{O}\right) \right]=\underline{Q} \tag{3}\]  
where \(\vec{v}\cdot\underline{F}=:\underline{F}^{\
u}\) is the normal flux on \(\partial\Omega\), whose value could be solved from \(\underline{U}^{h}_{l}\) (the approximated inner-side state) and \(\underline{U}^{h}_{O}\) (the approximated outer-side state). We implement this procedure as an independent module called the Riemann solver of the conservation law (Equation (1)), see [31] for details.