Question about dual- P1DCP_1^{DC} DG formulation in FreeFem++ documentation (nTonEdge logic & parameter al)

General Discussion
1

I am currently studying the section “Solution by Discontinuous-Galerkin FEM” in the FreeFem++ documentation, which presents a dual-P1DCP_1^{DC} formulation for solving a scalar advection problem.

The weak formulation presented is:

∫Ω(cn+1−cnδt+u⋅∇c)w+∫E(α∣n⋅u∣−12n⋅u)[c]w=∫EΓ−∣n⋅u∣cw∀w\int_{\Omega}\left(\frac{c^{n+1}-c^n}{\delta t} + u \cdot \nabla c \right) w + \int_E \left(\alpha |n \cdot u| - \frac{1}{2} n \cdot u \right) [c] w = \int_{E_{\Gamma}^{-}} |n \cdot u| c w \quad \forall w

In the associated code, the problem is defined as:

// Problem
problem Adual(cc, w)
    = int2d(Th)(
        (cc/dt + (v1*dx(cc) + v2*dy(cc)))*w
    )
    + intalledges(Th)(
        (1 - nTonEdge)*w*(al*abs(n) - n/2)*jump(cc)
    )
    - int2d(Th)(
        ccold*w/dt
    );

I have two questions regarding this implementation:

:one: About (1 - nTonEdge):

From the documentation and nTonEdge behavior, it seems that:

  • nTonEdge = 1 if the triangle has a boundary edge,

  • nTonEdge = 2 otherwise (i.e., internal edge with two adjacent elements).

In this case, the term (1 - nTonEdge) becomes:

  • 0 on boundary edges,

  • -1 on interior edges.

This seems counterintuitive to me, as it introduces a negative sign on interior edges. Could someone please explain why this negative sign is necessary, and how it fits into the weak formulation?

:two: About the parameter al:

The formulation involves a parameter al in the numerical flux:

(α∣n⋅u∣−12n⋅u)\left(\alpha |n \cdot u| - \frac{1}{2} n \cdot u\right)

Could you please clarify:

  • What is the typical choice for al?

  • What is the theoretical rationale behind it?

  • Should al = 1 always be used for full upwinding, or can it vary?

Any insights or references would be greatly appreciated.

2

Dear Tan,

You are right about the sign: there is a problem.

Indeed there is a mistake in the documentation: there should be a minus sign in front of the integral over E.

But the code is correct.

About the choice of alpha I don’t know.

3

why not /nTonEdge in + intalledges(Th)( (1 - nTonEdge)*w*(al*abs(n) - n/2)*jump(cc)? Do you know. I think the interior edges are computed twice, so should /nTonEdge. Am I right?

4

Dear Professor, why not /nTonEdge in + intalledges(Th)( (1 - nTonEdge)*w*(al*abs(n) - n/2)*jump(cc)? Do you know. I think the interior edges are computed twice, so should /nTonEdge. Am I right?

5

Dear Tan,

I agree that there is a mismatch between the formula and the code (factor 2). To fix the problem you should look at the mentioned reference :

ERN, A. and GUERMOND, J. L. Discontinuous Galerkin methods for Friedrichs’ symmetric systems. I. General theory. SIAM J. Numer. Anal., 2006

in particular section 5.2

6

I checked the reference, I find that the code is right. Indeed the correct formula is that there are two integrals for each internal edge: one for each of the two adjacent triangles.

For the choice of \alpha, according to the theory any positive value works, but there is a particularly interesting choice \alpha=1/2.