## 3 Stabilized Formulation  
Using any well-defined set of variables \(\mathbf{Y}\), it is possible to rewrite (2) in quasi-linear form as  
\[\mathbf{A}_{0}\mathbf{Y}_{,t}+\mathbf{A}_{i}\mathbf{Y}_{,i}=(\mathbf{K}_{ij}\mathbf{Y}_{,j})_{,i}+\mathbf{ S}, \tag{10}\]  
where \(\mathbf{A}_{0}=\mathbf{U}_{,\mathbf{Y}}\), \(\mathbf{A}_{i}=\mathbf{F}_{i,\mathbf{Y}}^{\text{adv}}\) is the \(i\)th Euler Jacobian matrix, and \(\mathbf{K}=[\mathbf{K}_{ij}]\) is the diffusivity matrix, where \(\mathbf{K}_{ij}\mathbf{Y}_{,j}=\mathbf{F}_{i}^{\text{diff}}\). In this paper, special attention is paid to the choice of pressure-primitive variables,  
\[\mathbf{Y}=\left\{ \begin{array}{l}p\\ \mathbf{u}\\ T,\end{array} \right\} \tag{11}\]  
where \(p\) is the pressure, \(u_{i}\) is the Cartesian velocity components, and \(T\) is the absolute temperature. The matrices and vectors for these variables can be found in [20]. This set of variables is endowed with the property that the incompressible limit is well behaved. The quasi-linear form can be written as  
\[\mathbf{\mathcal{L}Y}=S_{0}, \tag{12}\]  
with the differential operators defined below  
\[\mathbf{\mathcal{L}Y} = A_{0}\mathbf{Y}_{,t}+A_{i}\mathbf{Y}_{,i}-(\mathbf{K}_{ij}\mathbf{Y}_{,j})_{,i}- \mathbf{C}\mathbf{Y}\] \[\mathbf{\mathcal{L}}^{T}\mathbf{Y} = A_{0}^{T}\mathbf{Y}_{,i}+A_{i}^{T}\mathbf{Y}_{,i}-(\mathbf{K}_{ij}^{T}\mathbf{Y} _{,j})_{,j}-\mathbf{C}^{T}\mathbf{Y} \tag{13}\] \[-\mathbf{\mathcal{L}}^{*}\mathbf{Y} = A_{0}^{T}\mathbf{Y}_{,t}+A_{i}^{T}\mathbf{Y}_{,i}+(\mathbf{K}_{ij}^{T}\mathbf{Y} _{,i})_{,j}+\mathbf{C}^{T}\mathbf{Y}\] \[\mathbf{\mathcal{L}}_{\text{SUPG}}^{T}\mathbf{Y} = \mathbf{A}_{i}^{T}\mathbf{Y}_{,i}\]  
and the source term written as  
\[\mathbf{S} = \mathbf{C}\mathbf{Y}+\mathbf{S}_{0} \tag{14}\]  
Note that the difference between \(\mathbf{\mathcal{L}}^{T}\) and \(-\mathbf{\mathcal{L}}^{*}\) is the sign of the diffusive and source terms. In this paper, the time-discontinuous space-time Streamline Upwind Petrov-Galerkin (SUPG) and Subgrid-Scale (SGS) stabilized methods [26; 27] are considered. In the absence of source terms, and for linear shape functions, both methods coincide. Consider a space-time domain, where the time interval \(I=]0,T[\) is subdivided into \(N\) intervals \(I_{n}=]t_{n},t_{n+1}[\), \(n=0,1,\ldots,N-1\). We define for each time interval \(Q_{n}=\Omega\times I_{n}\) and \(P_{n}=\Gamma\times I_{n}\), where \(\Omega\) is the spatial domain and \(\Gamma\) its boundary. Finally, the \"slab\" \(Q_{n}\) is decomposed into the elements \(Q_{n}^{e},e=1,2,\ldots,(n_{\text{el}})_{n}\). Following [19; 20], the variational formulation is defined for the set of variables \(\mathbf{Y}\). Within each \(Q_{n}\), \(n=0,1,\ldots,N-1\), find \(\mathbf{Y}\in\mathcal{S}_{\mathbf{Y}}\) such that \(\forall\ \mathbf{W}\in\mathcal{V}_{\mathbf{Y}}\) :  
\[\int_{Q_{n}}\biggl{(}-\mathbf{W}_{,t}\cdot\mathbf{U}(\mathbf{Y})-\mathbf{W}_{,i} \cdot\mathbf{F}_{i}^{\text{adv}}(\mathbf{Y})+\mathbf{W}_{,i}\cdot\mathbf{K}_{ij}\mathbf{Y}_{,j}- \mathbf{W}\cdot\mathbf{S}\biggr{)}dQ\] \[\qquad\qquad+\int_{\Omega}\biggl{(}\mathbf{W}(t_{n+1}^{-})\cdot\mathbf{U }(\mathbf{Y}(t_{n+1}^{-}))-\mathbf{W}(t_{n}^{+})\cdot\mathbf{U}(\mathbf{Y}(t_{n}^{-}))\biggr{)}d\Omega\] \[\qquad\qquad\qquad+\sum_{e=1}^{(n_{\text{el}})_{n}}\int_{Q_{n}^{ e}}\biggl{(}\mathbb{L}\mathbf{W}\biggr{)}\cdot\mathbf{\tau}\biggl{(}\mathbf{\mathcal{L}}\mathbf{Y}- \mathbf{S}\biggr{)}dQ\] \[\qquad\qquad=\int_{P_{n}}\mathbf{W}\cdot\biggl{(}-\mathbf{F}_{i}^{\text{ adv}}(\mathbf{Y})+\mathbf{F}_{i}^{\text{diff}}(\mathbf{Y})\biggr{)}n_{i}dP \tag{15}\]  
The first and last integrals constitute the Galerkin terms expressed as functions of the variables \(\mathbf{Y}\), again written in conservative form to ensure that the weak solution is bestowed with the correct Rankine-Hugoniot relations.