I am solving the steady state compressible Euler equations in conservation form in 2D in a rectangular domain with a discontinuous Galerkin (DG) code I have developed. As a test, the boundary conditions are Dirichlet conditions set to the same far-field values. The Ma number is much less than 0.3.
I expect the solution to simply be the far-field values at every point in the domain. Indeed, for constant (q = 0) polynomials defined on the domain, I obtain this solution across both coarse and fine meshes. However, when the polynomial order is increased (q >= 1), the preconditioned Krylov solver (GMRES) can not satisfy the Newton system which describes the DG implementation of the Euler equations. Instead, the solution for all conservation variables is highly oscillatory. The solution on the boundary is correct, however.
Finally, my question is, assuming I have implemented the Newton system correctly (which I have painstakingly checked many, many times), what might be cause of this behavior? Allowing the Krylov solver more iterations decreases the residual of the solution to the Newton system slightly, but doesn't provide indication that the solver would converge given enough iterations. Considering conservation equations for other flux formulations, e.g. scalar advection in a constant velocity field, will converge under some p-refinement.
I suppose I should mention that I have implemented a hybridized DG scheme, where the Riemann solver is implicitly defined by a stabilization term on the numerical flux and I haven't seen flux limiters mentioned. The literature reports that this stabilization should have less effect on the final solution as the polynomial order increases. I have found little change in the behavior of the solution for various stabilization values beyond a small threshold value.
Edit:
In HDG, we introduce additional variables defined on the interfaces of the elements. We also define a numerical flux function defined on the interfaces written in terms of the trace and element variables which is given by
\begin{equation} f_h \left( u, \hat{u}, n \right) = F(\hat{u}) \cdot n + \sigma \left( u, \hat{u}, n\right) \left( u - \hat{u} \right) \end{equation}
where $u$ and $\hat{u}$ are the conservation variables defined on the element and trace, respectively, $F(q)$ is the Euler flux, and $\sigma$ is the stabilization matrix which implicitly defines the Riemann solver. I found the most success setting $\sigma = \tau I$ for large values of $\tau$, though the literature presents other formulations which include characteristics of $\frac{\partial F(q)}{\partial q} \cdot n$. However, these have precluded convergence in my case. The boundary conditions are imposed by solving
\begin{equation} \int_{\partial \Omega} b_h \left(u, \hat{u}, n \right) \cdot \mu = 0 \end{equation}
on the domain boundaries where $\mu$ is a test function in the same space to which $\hat{u}$ belongs and $b_h \left(u, \hat{u}, n \right) = f_h \left( u, \hat{u}, n \right) + f_h \left( u_b, \hat{u}, -n \right)$ with $u_b$ denoting the boundary state. After non-dimensionalization, $u_b$ is set to $\left[ \rho: 1, \rho u: 0.1, \rho v: 0, E: 1.78 \right]^T$.
