Derivation of weak form in undocumented DG Poisson demo

Question

In an attempt to become familiar with discontinuous Galerkin methods, I tried to derive the weak form for the simple Poisson equation, as used in the code here

# Define variational problem (ds(1) is Dirichlet, ds(2) is Neumann)
a = dot(grad(v), grad(u))*dx \
   - dot(avg(grad(v)), jump(u, n))*dS \
   - dot(jump(v, n), avg(grad(u)))*dS \
   + alpha/h_avg*dot(jump(v, n), jump(u, n))*dS \
   - dot(grad(v), u*n)*ds(1) \
   - dot(v*n, grad(u))*ds(1) \
   + (gamma/h)*v*u*ds(1)
L = v*f*dx - u0*dot(grad(v), n)*ds(1) + (gamma/h)*u0*v*ds(1) + g*v*ds(2)

I used a FEniCS lecture as a guide, which can be found here, and I understand the derivations made there. I will repeat the weak form from those slides below:

$$
\int_{\Omega} \nabla u \cdot \nabla v dx + \int_{\Gamma_i} \left( - \langle\nabla u\rangle[v] - \langle\nabla v\rangle[u] + \frac{\alpha}{h}[u][v] \right) dS + $$ $$\int_{\Gamma_e} \left( - \nabla u \cdot n v - \nabla v \cdot n u + \frac{\alpha}{h} u v \right) ds = \int_{\Omega} f v dx
$$

My question now is how to go from here to the form above? This is what I got so far:

1. Add the terms containing gamma on both sides to weakly enforce the Dirichlet boundary condition.
2. Split the external boundary in a Dirichlet and a Neumann part, and use dot(grad u, n) = g on the Neumann boundary

Doing so also results in term $\int_{\Gamma_N} \nabla v \cdot n u$ that is not accounted for in the code. I also miss a term with alpha on the external boundaries.

On the other hand, I cannot explain where the second term on the RHS comes from (u0*dot(grad(v), n)*ds(1)), neither why terms containing a testfunction v on the Dirichlet boundary don't disappear.

Any help on this would be very welcome.

Answer 1

The symmetry terms are only added on the Dirichlet boundary and the penalty term with alpha is a term with gamma here. This one is also not added at the Neumann part of the boundary. The penalty terms are essential for stability and enforcing the boundary conditions and everything with the normal derivative of the test function is added solely for symmetry and should be consistent, i.e. vanish for the strong solution.

Comments

Thank you for helping but I am still left with questions:

The symmetry terms are only added on the Dirichlet boundary

Why can you choose to only add them on the Dirichlet boundary? Following the derivation on the slides, you end up with a term over the whole external boundary (last term in slide 10), so if you add the symmetric term, wouldn't it be necessary over the whole boundary again?

and the penalty term with alpha is a term with gamma here.

There is a term with alpha in the code. It was my understanding that both terms with gamma (on the LHS and RHS) where added to weakly enforce the Dirichlet BC. If the terms with gamma in the code come from the term with alpha on the external boundary in the formula, then I don't see how?

I still don't understand where u0*dot(grad(v), n)*ds(1) comes from in the RHS, neither why it would vanish for the strong solution, since there is no u in here.

---

You can always add consistent terms, since in the derivation they are zero for exact solutions of your problem (u can be seen as continuous then, sufficient regularity assumed). They only matter in a discrete sense when you want to show stability of the finite-dimensional problem. Since u=u0 on the boundary, you can also add terms like (u-u0) times something, it doesn't matter in the derivation. Where it matters is in the discrete sense, since they penalize the jump at the boundaries, thus weakly enforce Dirichlet conditions for sufficiently large gamma (something like 10 times the polynomial degree squared is often sufficient). The terms with alpha are only present on interior facets in the example code (dS is for interior, ds for exterior facets, a small but important difference in the caps).

---

I think I understand, but then again I think I don't :-)
The way I see it one has to start with the original equation --> multiply by a testfunction and integrate over an element --> sum over all elements --> integrate by parts, resulting in the average and jump terms. It is only from here that you can add (or remove) consistent terms, no?

The two terms containing gamma meet the '(u - u0) times something' criterium, so those I understand.

The RHS term u0*dot(grad(v), n)*ds(1) on the other hand does not IMHO, because its opposite on the LHS dot(grad(v), u*n)*ds(1) was derived from the original PDE, and hence not added afterwards (second term on the second line of the formula above). So I still fail to see where that RHS term comes from, or why one can just add it here?

---

The term with grad(v) is also on the left with u instead of u0 (third but last on the LHS). It is just written differently for some reason. Put the u out of the brackets, subtract the RHS, and observe that it is again a term that contains (u-u0). You can add terms in strong or weak forms, it is always only a question of regularity and desireable properties you want to achieve for the discrete setting. Obtaining symmetry is the only reason why one would want the grad(v).n terms. The penalty is enough for stability, and you can actually also omit the symmetry term and obtain another variant or flip the sign and obtain unconditional stability in the DG norm (wrt gamma). However, the latter options sacrifice adjoint consistency which is crucial for optimal L2 error estimates.

---

As I mentioned in my previous comment, I realize that apart from the two terms containing gamma, there is another pair that can be rewritten into an '(u - u0) x ...' form. I also understand that one can add consistent terms like that at will.

But for this second pair, the term containing u was already present from the derivation. Hence only a term with u0 was added at some point, i.e. not consistent. I obviously am mistaken somewhere, but I can't put my finger on it.

---

Well, how do you think dot(grad(v), un)ds(1) is derived? A term dot(grad(u), vn)ds(1) pops out naturally using integration by parts, but how would you think one could end up with a term containing dot(grad(v),n) on the boundaries? See slide 7 of the PDF you linked to. The grad(v) terms are added on slide 12 for symmetry... They play no role in the derivation up to this point.

---

Ok, I believe I found my mistake. Thanks for your patience! I have a derivation now that reconstructs the weak form from the code, and that makes sense in my head.

I must say though, that if I'm correct, I find the slides very misleading (or maybe even wrong?).
Slide 11 mentions continuity of the solution [u] = 0 over all facets. It then adds a term $\int_{\Gamma} \langle \nabla v \rangle [u]$ over all facets, subsequently splitted in an integral over interior and exterior facets in slide 12. That is where I thought the boundary term containing dot(grad(v),n) came from, and hence why there is no matching term with u0 instead of u in the RHS.

In the derivation as I have it now, I assume that the continuity of the solution, [u] = 0, is something that only has meaning in the interior facets, so the correct term to add for symmetry is $\int_{\Gamma_{i}} \langle \nabla v \rangle [u] dS$ for the interior facets only, and then another $\int_{\Gamma_D} (u - u0) \nabla v \cdot n$ for the Dirichlet boundary, as you tried to tell me.

---

Glad you have it. Note that the derivation in the slides is for homogeneous Dirichlet BCs. This simplifies presentation, because you do not have to worry about the boundaries so much.

---

Oh, I forgot the fact that the example was for homogeneous Dirichlet BCs only. Then the slides are correct—although I think the attempt to keep things simple rather obscures the derivation here. Perhaps a matter of preference.

Thanks again for the help.