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.


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$.

  • 1
    $\begingroup$ If GMRES can't drive the residual to zero, then your matrix must have a nullspace and the right hand side is not in the range space. $\endgroup$ Jul 31, 2022 at 3:10
  • 1
    $\begingroup$ Are you just increasing the number of allowed iterations or are you also increasing the number of steps before the GMRES method restarts (assuming you are using a restart GMRES)? $\endgroup$ Jul 31, 2022 at 21:39
  • 1
    $\begingroup$ If your system is small enough, you can try using something like SuperLU to try and rule out the linear solver as the reason for convergence failures. $\endgroup$ Jul 31, 2022 at 21:42
  • $\begingroup$ I have increased both the number of GMRES iterations and the number of restart iterations with little success. I have also switched to a LU solver and have verified that the Jacobian matrix is not singular (using the cheaply computed determinant from the LU decomposition) at any iteration of the Newton method. Doesn't this imply that the kernel of the matrix is simply the trivial vector space containing only the zero vector? It does seem that the residual vector is outside the range of the Jacobian. In HDG, we can tune the stabilization matrix. I wonder if it is the culprit. $\endgroup$
    – Wil
    Aug 1, 2022 at 13:31
  • $\begingroup$ What data you specifying at the boundaries, and how are you imposing the BC? I get the impression that you're specifying the entire state at all boundaries, rather than just ingoing characteristics (or some approximation thereof) so it sounds to me like the issue has to do with the problem being over-constrained. $\endgroup$ Aug 2, 2022 at 11:10


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.