I am learning the theory behind DG-FEM methods using the Hesthaven/Warburton book and I am a bit confused about the role of the 'numerical flux.' I apologize if this is a basic question, but I have looked and not found a satisfactory answer to it.

Consider the linear scalar wave equation: $$\frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0$$ where the linear flux is given as $f(u) = au$.

As introduced in Hesthaven's book, for each element $k$, we end up with $N$ equations, one for each basis function, enforcing that the residual vanishes weakly:

$$R_h(x,t) = \frac{\partial u_h}{\partial t} + \frac{\partial au_h}{\partial x}$$

$$\int_{D^k} R_h(x,t) \psi_n(x) \, dx = 0 $$

Fine. So we go through integration by parts once to arrive at the 'weak form' (1) and integrate by parts twice to get the 'strong form' (2). I will adopt Hesthaven's sort-of-overkill but easily generalized surface integral form in 1D:

(1) $$ \int_{D^k} \left( \frac{\partial u_h^k}{\partial t} \psi_n-au_h^k\frac{d \psi_n}{d x} \right)\, dx = - \int_{\partial D^k} \hat{n}\cdot(au_h)^*\psi_n \,dx \qquad 1 \leq n\leq N$$

(2) $$ \int_{D^k} R_h \psi_n\, dx = \int_{\partial D^k} \hat{n}\cdot \left( au_h^k-(au_h)^* \right)\psi_n \,dx \qquad 1 \leq n\leq N $$

Why do we choose a numerical flux? Why don't we use the value of $au_h^k$ at the boundary in (1) instead of using a flux? Yes, it's true that the value of this quantity may be multiply defined across elements, but each equation is only over 1 element $D^k$, so why does this matter?

Further, the boundary term of the second integration by parts clearly yields a different quantity $au_h^k$ the second time in (2), which makes no sense to me. We are doing the same operation! Why wouldn't the two boundary terms just cancel, making (2) useless? How have we introduced new information?

Clearly I am missing something crucial to the method, and I would like to fix this. I have done some real and functional analysis, so if there is a more theory-based answer regarding the formulation, I would like to know!

  • 6
    $\begingroup$ One reason you choose a numerical flux so that you ensure conservation of $u$. If the flux at the boundary were not the same for each element that shares the boundary, the amount of $u$ flowing out of one element would be different than the amount flowing into the adjacent element. This is generally undesirable, since you're modeling a conservative transport equation. $\endgroup$ Commented May 18, 2016 at 23:14
  • 8
    $\begingroup$ Related to Tylers comment, but IMO even more important: the flux also introduces a coupling between the different subproblems. Otherwise there could be no propagation of information in a discrete sense. $\endgroup$ Commented May 19, 2016 at 19:20

3 Answers 3


The numerical flux is chosen to ensure that information in the problem travels in the direction of the characteristic curves of the equation (upwinding). As mentioned in the comments, the numerical flux is necessary in order to couple the subproblems defined on each element.

One way to get an intuition for the role of the numerical flux is to consider the following simple example.

Consider the scalar advection equation (where for simplicity $a=1$) $$ \frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} = 0 \qquad\text{on $\Omega$}, $$ where the domain is given by $\Omega = [0,1]$. Because this is a hyperbolic equation, and information is propagating from left to right, we need to enforce a boundary condition at $x = 0$ (but not at $x = 1$). For concreteness, suppose we enforce the Dirichlet condition $u(0,t) = g_D$ for some given $g_D$.

Suppose now we discretize this equation using the DG method, and we use two elements, $D_1 = [0,1/2]$ and $D_2 = [1/2,1]$. We could equally-well be discretizing the following set of two coupled PDEs, \begin{align*} \text{(PDE 1):}&& \quad v_t + v_x &= 0 \quad\text{on $D_1$},\\ \text{(PDE 2):}&& \quad w_t + w_x &= 0 \quad\text{on $D_2$}, \end{align*} where we will couple these equations to make them equivalent to the original equation.

To make the above equations well-posed, we need to enforce boundary conditions. As before, each equation is hyperbolic, and information is traveling from left to right. Therefore, we need to enforce a boundary condition for (PDE 1) on the left endpoint of $D_1$, and a boundary condition for (PDE 2) on the left endpoint of $D_2$.

The boundary condition on the left endpoint of $D_1$ must be chosen to be $v(0,t) = g_D$ in order to remain consistent with the original problem. We also look for smooth solutions, so the boundary condition on the left endpoint of $D_2$ must be chosen to enforce continuity. This condition reads $w(1/2,t) = v(1/2,t)$.

The DG method in this case chooses the numerical fluxes precisely to enforce the above boundary conditions. If we multiply by a test function $\psi$ and integrate by parts over each element $D_k$, we obtain boundary terms of the form \begin{align*} \int_{\partial D_1} \hat{n} \cdot v \psi \, dx &= \left[ v\psi \right]_0^{1/2}\\ \int_{\partial D_2} \hat{n} \cdot w \psi \, dx &= \left[ w\psi\right]_{1/2}^1 \\ \end{align*} In order to "weakly" enforce the boundary conditions, we replace $v$ and $w$ with the prescribed values at those points where boundary conditions are specified (i.e. the left endpoints of $D_1$ and $D_2$). This means we replace $v(0,t)$ by $g_D$ and $w(1/2,t)$ by $v(1/2,t)$ in the boundary integrals.

In other words, we define $u_h^* = g_D$ at $x = 0$, and $u_h^* = v(1/2,t)$ at $x=1/2$, and we recover exactly the standard upwind flux that is used in the DG method.

Looking at things this way, we can consider the numerical flux functions as weakly enforcing the boundary conditions on each element that are required to couple the equations in such a way that respects the characteristic structure of the equations.

For equations more complicated than constant-coefficient advection, information may not propagate always in the same direction, and so the numerical flux must be determined by solving (or approximating the solution to) a Riemann problem at the interface. This is discussed for linear problems in Section 2.4 of Hesthaven's book.


Very loosely speaking there are two things most discretization techniques need in order to converge to the actual solution of your PDE as you increase their approximation quality, regardless if you're using DG or not:

  1. Consistency (If a function $u$ satisfies the PDE, then it also satisfies your weak formulation)
  2. Stability (small changes in data result in small changes in answer)

The first steps of a DG derivation where you integrate by parts on each mesh element preserves (1) because you are starting with the PDE and only applying legal operations from there.

This does not give you (2) though. You can see this yourself by trying to assemble the matrix of a partially formulated DG weak form and looking at its eigenvalues - for time dependent problem we want them all on the left half plane, but without a proper numerical flux they will be everywhere. This leads to a solution that explodes exponentially in time even if the physical problem does not.

Thus you need to add terms to your formulation so that (2) is satisfied, but without hurting (1). This is tricky to do, but not impossible. You can replace function values with cell wall averages without hurting consistency, and you can always add the jump in at cell walls without hurting consistency (because for a solution $u$ with suitable smoothness property, the jump is just 0!)

The trick is to take combinations of jumps and averages and combine them in a way that your scheme is still consistent but also stable. After that a convergence theorem usually reveals itself.

This is the basics, but you can also often bring in additional physics into the numerical flux so that it doesn't simply satisfy these mathematical requirements but also plays nicely with conservation principles.


When you choose the test function equal to the trial function in the DG method you are creating an optimization problem. That is, you have a Galerkin rather than a Petrov-Galerkin method. You are seeking the time derivatives of the trial function amplitudes that will minimise the element residual in the L2 norm, and you make this miminimisation on the assumption of a given flux function at inflow.


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.