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?

  • $\begingroup$ How do you solve your problem? $\endgroup$ Nov 26, 2013 at 13:23
  • 1
    $\begingroup$ your manufactured solution is quadratic, so you should capture the solution exactly for $p=2$ (zero discretization error) on any grid. Do you? If you don't, switch off the terms one by one to find the problem. If you still can't find the problem, use a linear manufactured solution, and eventually a constant. Once fixed, use a manufactured solution with sines to check convergence rates for any $p$. $\endgroup$
    – chris
    Nov 26, 2013 at 15:49
  • $\begingroup$ @chris Isn’t the manufactured solution quartic and requires $p=4$ to capture the solution exactly? For $p=4$, the solution is captured with error on the order of $10^{-14}$, which is close enough to machine precision I think. $\endgroup$ Nov 26, 2013 at 16:34
  • $\begingroup$ The solution is quadratic in each direction. $\endgroup$
    – Jesse Chan
    Nov 26, 2013 at 16:37
  • 1
    $\begingroup$ @JLC It would be unconventional (at best) to use a biquadratic space on a triangular mesh. Those spaces are generally used on quadrilateral meshes. $\endgroup$
    – Jed Brown
    Nov 29, 2013 at 15:59

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

  • $\begingroup$ your error is likely in the assembling routine. You should do a step by step debugging. Firstly by setting $\mathbf{s} = 0,$ you should be able to check from the $\| \cdot \|_{L_2}$ easily. If this however does succeed, then you can add the upwind terms and again a simple $L_2-$norm. It is also known from literature that the penalty parameter depends quadratically on the polynomial degree $(p)$ and also on the dimension $(d)$ of the problem, i.e. $\sigma = 2*(p+d)(p+1).$ The penalty cannot depend on the diffusion coefficient as this affects the boundedness and coercivity results. $\endgroup$
    – uli.xu
    Jan 4, 2017 at 2:40

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.