1-D Linear Convection with DG1

Building upon our previous introduction, this example is the next natural step in the Discontinuous Galerkin (DG) framework. We move from piecewise-constant (DG0) approximations to linear polynomials (DG1), capturing the solution's gradient within each element.

Why move from DG0 to DG1?

In the DG0 method, we represented the solution with a single constant value per element. While simple and conservative, it acts exactly like a first-order finite volume scheme and suffers from heavy numerical diffusion. By moving to DG1, we enrich our local approximation space. Now, each element stores both a mean value and a slope. This allows the solution to vary linearly inside each element, drastically reducing numerical diffusion and paving the way toward high-order spectral accuracy.

Conservative form of linear convection

We continue to solve the 1-D linear convection equation in conservative form. Assuming a constant, positive wave speed $c > 0$, the equation is:

Elements and local coordinates

We partition our 1-D spatial domain into elements of uniform size $\Delta x$. To generalize the mathematics, we map each physical element $K_i$ onto a standard reference element with a local coordinate $\xi \in [-1, 1]$. Here, $\xi = -1$ corresponds to the left boundary, $\xi = 1$ to the right boundary, and $\xi = 0$to the element's center.

Linear polynomial approximation

Inside the $i$-th element, the numerical solution $u_h$ is approximated as a linear polynomial:

The two time-dependent coefficients dictate the shape: $a_i(t)$ is the cell average (mean mode), and $b_i(t)$ controls the gradient (slope mode).

Modal basis functions

This formulation uses a set of hierarchical modal basis functions (specifically, the first two Legendre polynomials):

• $\phi_0(\xi) = 1$
• $\phi_1(\xi) = \xi$

Because these basis functions are orthogonal on the interval $[-1, 1]$, the resulting mass matrix is diagonal. This means the evolution equations for the mean $\dot{a}_i$ and the slope $\dot{b}_i$ decouple beautifully, avoiding the need for an expensive matrix inversion.

Weak form and interface fluxes

To derive the semi-discrete scheme, we multiply the governing equation by our test functions ($\phi_0$ and $\phi_1$) and integrate over the element. Using integration by parts, we shift the spatial derivative from the solution onto the test functions, generating boundary terms evaluated at $\xi = \pm 1$. Because the polynomial approximation is discontinuous across element boundaries, we must introduce a single-valued numerical flux to link adjacent elements.

Upwind numerical flux

Information travels strictly from left to right ($c > 0$). Therefore, the physical flux at an interface is entirely determined by the upwind state—the value approaching the boundary from the left element. For the interface at the right boundary of element $i$ (where $\xi = 1$), the left state is simply the polynomial evaluated at that boundary:

SSP RK3 time integration

Unlike DG0, the DG1 solution varies spatially within the element. This internal gradient makes the scheme much more sensitive to temporal errors. A standard forward Euler explicit step is often unstable or heavily restricted. Instead, we use the 3-stage Strong Stability Preserving Runge-Kutta (SSP RK3) method. SSP RK3 is an explicit integrator built from convex combinations of forward Euler steps, ensuring that bounds and stability properties of the spatial discretization are maintained in time.

Understanding the CFL number for DG1

Stability requires limiting the time step size according to the Courant-Friedrichs-Lewy (CFL) condition:

Crucially, DG1 has a stricter practical stability limit than DG0 or first-order FVM. The higher-order polynomial content means signals can effectively cross an element faster than the mean flow $c$. For our standard DG1 formulation with SSP RK3, theoretical stability limits the CFL to $\le 1/3$ (though some specific problem setups allow slightly higher limits depending on exact definitions).

CFL sweep interpretation

In the companion Python code, we project a smooth Gaussian pulse onto the DG1 basis using numerical quadrature and simulate its transport under different CFL conditions:

• CFL = 0.4: The scheme is stable. Compared to DG0, the initial Gaussian pulse shape is preserved much better as it travels across the periodic domain, proving that the added slope degrees of freedom drastically reduce numerical diffusion.
• CFL = 0.8: The stability limit is severely violated. In the graphical output, the plotted values are aggressively clipped for readability, but inspecting the minimum, maximum, and L2 error values reveals massive exponential numerical growth ($\approx 10^{58}$).

Relationship to DG0, FVM, and FDM

To understand where we are in the landscape:

• FDM: Stores point values and approximates spatial derivatives using neighboring points.
• FVM: Stores cell averages and updates them strictly through fluxes entering and leaving the volume.
• DG0: Stores one constant coefficient (the mean) per element. Interface fluxes couple elements. It is algebraically equivalent to first-order upwind FVM.
• DG1: Stores two coefficients per element (mean and slope). It tracks internal gradients explicitly, bridging the gap between conservative finite volumes and high-order finite elements.

What comes next: higher-order DG and DGSEM

DG1 demonstrates the fundamental mechanics of the Discontinuous Galerkin method. The next logical evolution is to increase the polynomial degree $p$ arbitrarily. However, extending the modal basis approach to very high orders becomes computationally cumbersome. To achieve spectral-like accuracy efficiently, we will transition to the Discontinuous Galerkin Spectral Element Method (DGSEM), which employs nodal polynomials (like Lagrange basis functions), Legendre-Gauss-Lobatto (LGL) quadrature nodes, and collocation to build a blazing-fast, highly accurate solver.