## 3 Numerical Framework
### Discontinuous Galerkin Discretization  
We are interested in the space discretization of the VDI flow model (3) using a high-order discontinuous Galerkin (DG) method. Let \(\mathcal{K}_{h}=\{\mathrm{K}\}\) denote the collection of non-overlapping mesh elements \(\mathrm{K}\) of the domain \(\Omega\in\mathbb{R}^{d}\), such that  
\[\overline{\Omega}_{h}=\bigcup_{\mathrm{K}\in\mathcal{K}_{h}}\overline{\mathrm{ K}}.\]  
Let \(\mathcal{F}_{h}=\{\mathrm{F}\}\) denote the set of faces of the computational mesh \(\Omega_{h}\). \(\mathcal{F}_{h}^{b}\) collects the faces at the boundary of \(\Omega_{h}\) and \(\mathcal{F}_{h}^{i}=\mathcal{F}_{h}\setminus\mathcal{F}_{h}^{b}\) the internal faces, i.e., for any \(\mathrm{F}\in\mathcal{F}_{h}^{i}\), there exist two elements \(\mathrm{K}^{+},\mathrm{K}^{-}\in\mathcal{K}_{h}\), such that \(\mathrm{F}\in\partial\mathrm{K}^{+}\cap\partial\mathrm{K}^{-}\). Then, let us define the discrete polynomial space  
\[\mathbb{P}_{d}^{k}(\mathcal{K}_{h})=\Big{\{}v_{h}\in L^{2}(\Omega)\,|\,v_{h| \mathrm{K}}\in\mathbb{P}_{d}^{k}(\mathrm{K}),\,\forall\mathrm{K}\in\mathcal{K }_{h}\Big{\}},\]  
where \(k\geq 0\) is an integer and \(\mathbb{P}_{d}^{k}(\mathrm{K})\) denotes the restriction to \(\mathrm{K}\) of the polynomial functions of \(d\) variables and total degree \(\leq k\). For this polynomial space, a set of orthonormal and hierarchical basis functions is computed according to [29]. By multiplying Equation (3) by an arbitrary smooth test function \(\mathbf{v}=\{v_{1},\ldots,v_{m}\}\), and integrating by parts, we obtain the weak formulation  
\[\int_{\Omega}\mathbf{v}\cdot\left(\mathbf{D}\mathbf{P}(\mathbf{w})\frac{\partial \mathbf{w}}{\partial t}\right) \mathrm{d}\mathbf{x}-\int_{\Omega}\
abla\mathbf{v}:\left[\mathbf{F}^{c}(\mathbf{w} )+\mathbf{F}^{v}(\mathbf{w},\
abla\mathbf{w})\right]\mathrm{d}\mathbf{x}\\ +\oint_{\partial\Omega}\mathbf{v}\otimes\mathbf{n}:\left[\mathbf{F}^{c}(\mathbf{ w})+\mathbf{F}^{v}(\mathbf{w},\
abla\mathbf{w})\right]\mathrm{d}\mathbf{x}+\int_{\Omega}\mathbf{v} \cdot\mathbf{s}(\mathbf{w})\,\mathrm{d}\mathbf{x}=0, \tag{4}\]  
where \(\mathbf{n}\) is the unit vector normal to the boundary. In order to discretize Equation (4) by means of discontinuous Galerkin methods, the solution \(\mathbf{w}\) and the test function \(\mathbf{v}\) are replaced with a finite element approximation \(\mathbf{w}_{h}\) and a discrete test function \(\mathbf{v}_{h}\), respectively, both of which belong to \([\mathbb{P}_{d}^{k}(\mathcal{K}_{h})]^{m}\). Accordingly, within each mesh element \(\mathrm{K}\in\mathcal{K}_{h}\) each component \(w_{hj}\), \(j=1,\ldots,m\), of \(\mathbf{w}_{h}\) can be expressed as  
\[w_{hj}=\phi_{l}\mathcal{W}_{jl},\]  
where \(\phi_{l}\) are the orthonormal and hierarchical basis functions and \(W_{jl}\) are the components of the global vector \(\mathbf{W}\) of the unknown degrees of freedom (dof), being \(l=1,\ldots,N_{dof}^{k}\). As a generic function, \(\phi_{h}\in\mathbb{P}_{d}^{k}(\mathcal{K}_{h})\) is double-valued over an internal mesh face and the jump trace operator \(\llbracket\phi_{h}\rrbracket=\phi_{h}|_{\mathrm{K}^{+}}-\phi_{h}|_{\mathrm{K}^ {-}}\) is introduced. Whenever applied to a vector, this operator acts component-wise. is the Riemann problem we ought to find a solution to. In (7), \(x_{n0}\in\mathrm{F}\) denotes the position of the discontinuity at initial time \(t_{0}\), and the subscripts \(L\) (left) and \(R\) (right) denote the states of the two mesh elements \(\mathrm{K}^{+}\) and \(\mathrm{K}^{-}\) sharing \(\mathrm{F}\). According to the Godunov scheme approach, the solution at the locus \(\mathcal{S}=(x_{n}-x_{n0})/t=0\) (vertical axis of Figure 1) is directly employed to evaluate the convective numerical fluxes. Since the system (6) does not fit into the framework of first-order hyperbolic conservation laws, the hyperbolic wave pattern solving the Riemann problem cannot be found.