Convergence of interior penalty DG methods

Question

I’m currently having some issues with my routine for the linear advection-diffusion problem. The model problem is as follows: $$ \nabla\cdot(\mathbf{s} u) - \nabla\cdot(\kappa\nabla u) = f, \;\;\;\text{in}\;\Omega \subseteq \mathbb{R}^{2} $$ $$ u = g_{D} \;\;\;\text{on}\;\Gamma_{D} $$

The boundary is purely Dirichlet in this case. The formulation is: $$ a(u,v) = \sum_{K \in \mathcal{M}_{h}}\int_{K}\kappa\nabla u\cdot\nabla v - \sum_{e \in \mathcal{E}_{h}}\int_{e}\left\{\kappa\nabla u\cdot\mathbf{n}\right\}[[v]] + \sum_{e \in \mathcal{E}_{h}}\left\{\kappa\nabla v\cdot\mathbf{n}\right\}[[u]] + \sum_{e \in \mathcal{E}_{h}}\frac{\sigma}{|e|}\int_{e}[[u]][[v]] $$ $$ b(u,v) = -\sum_{K \in \mathcal{M}_{h}}\int_{K}\mathbf{s}u\cdot\nabla v + \sum_{e \in \mathcal{E}_{\text{int}}}\int_{e}\left\{\mathbf{s}\cdot\mathbf{n}\right\}u^{up}[[v]] + \sum_{e \in \Gamma_{\text{out}}}\int_{e}(\mathbf{s}\cdot\mathbf{n})uv $$ $$ \ell(v) = \int_{\Omega}fv + \sum_{e \in \Gamma_{D}}\int_{e}\left( \kappa\nabla v\cdot\mathbf{n} + \frac{\sigma}{|e|}v\right)g_{D} - \sum_{e \in \Gamma_{\text{in}}}\int_{e}(\mathbf{s}\cdot\mathbf{n})g_{D}v $$ So, we have $$ a(u,v) + b(u,v) = \ell(v) $$

The inflow and outflow boundaries are defined as $\Gamma_{\text{in}} = \left\{e \in \partial\Omega| \;\mathbf{s}\cdot\mathbf{n} < 0\right\}$ and $\Gamma_{\text{in}} = \partial\Omega\backslash\Gamma_{\text{in}}$. And $\mathcal{M}_{h}$ is the set of all mesh elements (triangular in my case), $\mathcal{E}_{h}$ the set of all edges of the mesh elements, and $\mathcal{E}_{\text{int}}$ just the set of interior edges.

For the linear diffusion problem, my routine has no issues at all. Now if I use the formulation above for the exact solution $u(x,y) = x(x-1)y(y-1)$, for some values of $\kappa$ and $\mathbf{s}$ yield suboptimal convergence when I’m quite positive that they shouldn’t.

The source I’m primarily using is “Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations” by B. Riviere. For this formulation, the convergence rate (for penalization of $\sigma \ge 0$) in the $L^{2}$ norm should be $p$ if $p$ is even and $p+1$ if $p$ is odd, where $p$ is the degree of the basis functions used. In the energy norm, we should have convergence rate of $p$ in all cases.

Now here are the examples I am testing and the results. The exact solution is taken to be $u(x,y) = x(x-1)y(y-1)$ and the source term is chosen accordingly. I take $\Omega = [0,1]^{2}$.

1) $\kappa = 10^{-10}$, $\mathbf{s} = [1,1]^{T}$. $L^{2}$ convergence is $p$ for all $p$ (suboptimal). Energy norm convergence is $p-1$ for all $p$ (suboptimal) (with the exception of $p=1$, which has convergence rate $p$). This is for $\sigma = 1$.

2) Same example as above, but this time with $\sigma = 0$. I obtain the proper convergence rates that I stated in the paragraph above.

3) Same example as in 1), with $\sigma = 1$, $\kappa=10^{-10}$ and $\mathbf{s} = [1000,1000]^{T}$. I obtain the proper convergence rates here too.

4) In general, for $\mathbf{s} = [1,1]^{T}$, if I take $\kappa=10^{-k}$ for some positive integer $k > 3$ (approximately speaking), I attain the suboptimal convergence rate.

5) If I multiply the entire model equation by $10^{10}$ on both sides and use the resulting formulation (so $\kappa = 1$ and $\mathbf{s} = [10^{10}, 10^{10}]^{T}$ and the source term $f$ is modified accordingly), I obtain proper convergence rates.

Any ideas what is going on here? Is this to be expected or a bug in the code? I’ve tested on a number of a examples and can’t really find other examples where the rate is suboptimal. It’s very perplexing to me that in examples 1) and 2), if I increase the advection it actually gives me the expected convergence.

Can any DG experts chime in here?

Answer 1

Since you haven't had an answer yet, I'll reformulate my comment.

Saying that a method is $p$-th order accurate implies that a polynomial manufactured solution of lesser order can be captured exactly. For example, a 2nd order method will represent a linear solution up to machine precision. This has often helped me in finding implementation issues. It's easier and faster than checking convergence rates, although you need to do that afterwards for non-polynomial solutions (sines and such).

Now consider the case of constant basis functions, which is similar but not equal to finite volume methods. In that case, note that the only remaining diffusion term is the penalty term! So it's highly unlikely that a constant $\sigma$ can give satisfactory results for such a wide range of $\kappa$. Search the literature for suggestions about $\sigma(\kappa)$.