Discontinuous Galerkin Method Tutorial: Derivation, Fluxes, and Python Implementation

The discontinuous Galerkin method, usually shortened to DG, is a numerical technique for solving partial differential equations. It combines several useful properties of finite volume and finite element methods. Like the finite volume method, DG can preserve quantities locally through carefully defined fluxes. Like the finite element method, it supports flexible meshes and polynomial approximations inside each element.

The defining feature of DG is that the numerical solution does not need to remain continuous across element boundaries. Each element has its own polynomial representation. Neighboring elements communicate through numerical fluxes evaluated at their shared interfaces.

This structure makes DG useful for wave propagation, transport equations, compressible flow, shallow-water systems, electromagnetics, and problems that contain shocks or material interfaces. It is also attractive for parallel computing because many calculations can be performed independently inside each element.

How DG Differs From Other Methods

A finite volume method usually stores one average value in each cell and calculates the flow of a conserved quantity through cell faces. It provides strong local conservation, but many basic finite volume schemes use relatively low-order spatial approximations.

A continuous finite element method represents the solution with polynomial basis functions that connect continuously across element boundaries. This offers geometric flexibility and high-order approximation, but discontinuities cannot be represented directly without additional treatment.

DG uses polynomial basis functions inside every element while allowing separate values on either side of an interface. Numerical fluxes determine how those values interact. The method therefore provides local conservation, high-order approximation, and direct support for discontinuous solutions.

The Continuous Advection Equation

The one-dimensional linear advection equation provides a clear introduction to DG:

∂u/∂t + ∂f(u)/∂x = 0

For constant transport velocity a, the physical flux is:

f(u) = a u

The equation then becomes:

∂u/∂t + a ∂u/∂x = 0

This model describes a profile moving through the domain without changing shape. When a is positive, information travels from left to right. When it is negative, information travels from right to left.

Dividing the Domain Into Elements

Suppose the domain extends from zero to L. DG divides it into non-overlapping elements:

Ω = K₁ ∪ K₂ ∪ ... ∪ Kₙ

Within each element, the numerical solution is represented by a polynomial:

uₕ(x, t) = Σ Uᵢ(t) φᵢ(x)

The functions φᵢ are local basis functions, and the coefficients Uᵢ are the degrees of freedom that change with time.

A polynomial of degree zero stores one constant value in each element. A linear approximation uses two local degrees of freedom in one dimension. A quadratic approximation uses three.

Unlike continuous finite elements, DG does not force the polynomial from one element to match the polynomial in the next. At a shared interface, the solution may therefore have a left value and a right value.

Deriving the Element-Wise Weak Form

The DG formulation begins by multiplying the governing equation by a test function v and integrating over one element:

∫K (∂uₕ/∂t) v dx + ∫K (∂f(uₕ)/∂x) v dx = 0

Integration by parts moves the spatial derivative from the physical flux to the test function:

∫K (∂uₕ/∂t) v dx
- ∫K f(uₕ) ∂v/∂x dx
+ f(uₕ) v |∂K = 0

The boundary term is essential. It describes the flow through the left and right sides of the element.

Because the DG solution may be discontinuous at an interface, the physical flux is not uniquely defined there. One element provides a value from the left, while its neighbor provides another value from the right. DG replaces the ambiguous physical flux with a numerical flux:

f̂(u⁻, u⁺)

The final weak form for one element is:

∫K (∂uₕ/∂t) v dx
- ∫K f(uₕ) ∂v/∂x dx
+ f̂R vR
- f̂L vL = 0

This equation contains an element-local volume contribution and two interface contributions. The numerical fluxes are the only terms that directly connect neighboring elements.

Why Numerical Fluxes Matter

A numerical flux must be consistent. When the two interface values are equal, it should reproduce the physical flux:

f̂(u, u) = f(u)

It must also provide appropriate stability. For transport problems, the flux should respect the direction in which information moves.

The flux choice controls how strongly neighboring elements interact and how much numerical dissipation is introduced. An unsuitable flux can create oscillations, excessive smoothing, or incorrect wave speeds.

Upwind Flux

For linear advection with positive velocity, information comes from the left side of an interface. The upwind flux is therefore:

f̂(u⁻, u⁺) = a u⁻    when a > 0

For negative velocity, information comes from the right:

f̂(u⁻, u⁺) = a u⁺    when a < 0

A compact expression is:

f̂ = a⁺u⁻ + a⁻u⁺

a⁺ = max(a, 0)
a⁻ = min(a, 0)

The upwind flux is simple, stable, and widely used. It introduces some numerical dissipation, but the amount generally decreases as the mesh is refined or the polynomial order increases.

Local Lax-Friedrichs Flux

The local Lax-Friedrichs flux, also called the Rusanov flux, is frequently used for nonlinear conservation laws:

f̂(u⁻, u⁺)
= 0.5[f(u⁻) + f(u⁺)]
- 0.5 α(u⁺ - u⁻)

The parameter α is an estimate of the largest characteristic speed at the interface. For constant linear advection, it is usually |a|.

The first part averages the physical fluxes. The second part adds dissipation that helps control unstable interface modes. This flux is robust and easy to implement, although it may be more diffusive than a specialized approximate Riemann solver.

Riemann Solvers

Nonlinear hyperbolic systems may contain several waves traveling at different speeds. Examples include the Euler equations for compressible flow and the shallow-water equations.

A Riemann solver examines the left and right states at an interface and estimates the waves produced by their interaction. Exact Riemann solvers can be expensive, so practical DG codes commonly use approximate methods such as Roe, HLL, HLLC, or Rusanov fluxes.

The correct choice depends on the equation, desired accuracy, robustness requirements, and ability to preserve important physical properties.

The Semi-Discrete Matrix System

After selecting basis and test functions, the weak form can be expressed as a system of ordinary differential equations:

M dU/dt = R(U)

The mass matrix contains integrals of basis-function products:

Mᵢⱼ = ∫K φᵢ φⱼ dx

The residual R(U) includes the volume derivative and numerical flux contributions.

Because DG basis functions belong to individual elements, the global mass matrix has a block-diagonal structure. Each element contributes a small independent block. These blocks can be inverted separately.

The mass matrix is not automatically diagonal for every basis and integration rule. It becomes diagonal or approximately diagonal in common nodal DG formulations that use matching interpolation and quadrature points, such as Gauss-Lobatto collocation. Other formulations use small dense element matrices.

Choosing Basis Functions and Quadrature

DG methods commonly use modal or nodal polynomial bases. A modal basis represents the solution through polynomial modes, often based on Legendre polynomials. A nodal basis stores solution values at interpolation nodes inside the element.

Nodal formulations are convenient because interface values can be obtained directly when nodes are placed on element boundaries. Gauss-Lobatto points include both endpoints of the reference element, while Gauss points remain inside it.

Numerical quadrature evaluates the integrals in the weak form. The quadrature rule must be sufficiently accurate for the polynomial degree and any nonlinear terms. Insufficient quadrature can introduce aliasing error and instability.

Polynomial Order and Accuracy

For sufficiently smooth solutions, a well-designed DG method can achieve an error proportional to approximately h^(p+1). Increasing the polynomial degree can therefore improve accuracy without adding more elements.

Higher order is not always better. Shocks and sharp discontinuities can create oscillations near the jump. Limiters, artificial viscosity, filtering, or shock-capturing techniques may be needed.

A Linear DG Solver in NumPy

The following educational example implements a linear p = 1 DG method for periodic one-dimensional advection. Each element contains two degrees of freedom located at its endpoints.

The implementation uses the element weak form, an upwind flux, and a third-order strong-stability-preserving Runge-Kutta method.

import numpy as np

# Domain and model parameters
length = 1.0
number_of_elements = 80
velocity = 1.0
final_time = 0.5

element_width = length / number_of_elements
jacobian = element_width / 2.0

# Reference-element mass matrix for linear basis functions
mass_reference = np.array([
    [2.0 / 3.0, 1.0 / 3.0],
    [1.0 / 3.0, 2.0 / 3.0]
])

mass_matrix = jacobian * mass_reference
inverse_mass = np.linalg.inv(mass_matrix)

# S[i, j] = integral(phi_j * derivative(phi_i)) on [-1, 1]
volume_matrix = np.array([
    [-0.5, -0.5],
    [ 0.5,  0.5]
])

left_vector = np.array([1.0, 0.0])
right_vector = np.array([0.0, 1.0])

# Physical coordinates of local DG nodes
left_edges = np.arange(number_of_elements) * element_width
right_edges = left_edges + element_width
coordinates = np.column_stack((left_edges, right_edges))

# Smooth periodic initial condition
solution = (
    0.5
    + 0.5 * np.sin(2.0 * np.pi * coordinates / length)
)

def upwind_flux(left_state, right_state, speed):
    positive_speed = max(speed, 0.0)
    negative_speed = min(speed, 0.0)

    return (
        positive_speed * left_state
        + negative_speed * right_state
    )

def spatial_residual(values):
    residual = np.zeros_like(values)

    for element in range(number_of_elements):
        left_neighbor = (element - 1) % number_of_elements
        right_neighbor = (element + 1) % number_of_elements

        # States on the left interface
        left_inside = values[element, 0]
        left_outside = values[left_neighbor, 1]

        # States on the right interface
        right_inside = values[element, 1]
        right_outside = values[right_neighbor, 0]

        flux_left = upwind_flux(
            left_outside,
            left_inside,
            velocity
        )

        flux_right = upwind_flux(
            right_inside,
            right_outside,
            velocity
        )

        local_rhs = (
            velocity * volume_matrix @ values[element]
            + flux_left * left_vector
            - flux_right * right_vector
        )

        residual[element] = inverse_mass @ local_rhs

    return residual

# Conservative time-step estimate
time_step = 0.1 * element_width / abs(velocity)
current_time = 0.0

while current_time < final_time:
    dt = min(time_step, final_time - current_time)

    # SSP-RK3 stage 1
    stage_one = solution + dt * spatial_residual(solution)

    # SSP-RK3 stage 2
    stage_two = (
        0.75 * solution
        + 0.25 * (
            stage_one
            + dt * spatial_residual(stage_one)
        )
    )

    # SSP-RK3 stage 3
    solution = (
        (1.0 / 3.0) * solution
        + (2.0 / 3.0) * (
            stage_two
            + dt * spatial_residual(stage_two)
        )
    )

    current_time += dt

print("Simulation completed")
print("Final time:", current_time)
print("Minimum value:", solution.min())
print("Maximum value:", solution.max())

This example is intentionally limited to a smooth linear problem. Production DG solvers also need robust boundary treatments, multidimensional mappings, accurate quadrature, nonlinear flux evaluation, limiters, and more advanced time integration.

Time Integration and the CFL Condition

The spatial DG discretization creates a system of ordinary differential equations. An explicit method such as Runge-Kutta can then advance that system in time.

The maximum stable time step depends on wave speed, element size, polynomial degree, flux, and time integrator. A common scaling is:

Δt ∝ h / [(2p + 1)|a|]

The exact stability constant is method-dependent. It should be established through analysis, documentation, or numerical testing rather than treated as universal.

Increasing the polynomial order normally reduces the largest stable explicit time step. High-order DG may therefore require more steps even when it needs fewer elements.

DG for Advection-Dominated Problems

Standard continuous Galerkin methods can develop oscillations when advection dominates diffusion. Stabilized continuous methods such as SUPG modify the test functions to add control along streamlines.

DG handles transport through interface fluxes. Upwind or approximate Riemann fluxes introduce direction-aware stabilization and preserve local conservation.

SUPG remains effective for many smooth convection-dominated problems. DG is attractive when local conservation, discontinuous solutions, complex meshes, or element-wise adaptivity are central requirements.

DG vs. FVM and Continuous FEM

Advantages of the DG Method

• Local conservation is built into the interface flux balance.
• Polynomial degree can vary between elements.
• Discontinuities do not violate the approximation space.
• Element calculations are well suited to parallel hardware.
• Complex and unstructured meshes can be supported.
• Mesh refinement and polynomial enrichment can be combined.
• Different numerical fluxes can be selected for different equations.

Limitations of DG

DG normally uses more degrees of freedom than a continuous finite element method of the same polynomial degree because neighboring elements do not share interface values.

The method also requires careful flux design. Hyperbolic, elliptic, and mixed equations need different interface treatments. Diffusion operators require formulations such as interior penalty, local DG, or related techniques.

Explicit high-order DG methods can have restrictive time-step limits. Nonlinear problems may require limiters, entropy-stable fluxes, positivity preservation, or artificial viscosity.

These features make DG powerful but more difficult to implement correctly than basic finite difference or finite volume schemes.

Common Implementation Mistakes

One common mistake is using the physical flux directly at a discontinuous interface without defining how the two neighboring values should be combined.

Other frequent errors include:

• Using the wrong sign for the left or right boundary contribution
• Applying inconsistent normal-vector orientations
• Assuming the mass matrix is always diagonal
• Using quadrature that is too weak for nonlinear terms
• Ignoring the polynomial-degree dependence of the CFL limit
• Applying periodic or physical boundary conditions incorrectly
• Using high-order polynomials near shocks without a limiter
• Testing only visual output instead of convergence and conservation

How to Validate a DG Solver

Begin with a problem that has a known analytical solution. Periodic linear advection is useful because the exact profile is simply shifted by a t.

Run the solver with several mesh resolutions and measure an error norm. For a smooth solution, the observed convergence rate should approach the theoretical order of the method.

Check conservation by integrating the solution across the domain. For periodic advection, the total mass should remain nearly constant.

Test both positive and negative velocities to verify the upwind direction. Constant initial data should remain constant. Boundary states should also be tested separately before moving to nonlinear equations.

When to Choose DG

DG is a strong choice when the problem contains waves, shocks, material interfaces, or transport-dominated behavior. It is also useful when local conservation is essential or when the simulation benefits from high-order accuracy on an unstructured mesh.

A basic finite volume method may remain simpler for low-order conservation problems. Continuous finite elements may be more efficient for smooth elliptic or structural problems where discontinuities are not expected.

The numerical method should follow the mathematical structure of the PDE rather than current popularity or software availability.

Conclusion

The discontinuous Galerkin method represents the solution with independent polynomials inside each element. Its weak formulation produces natural boundary terms that are replaced by numerical fluxes. These fluxes control how information moves between elements and allow the complete method to preserve quantities locally.

DG combines high-order polynomial approximation, geometric flexibility, discontinuity handling, and strong element locality. These advantages make it valuable for wave propagation, hyperbolic conservation laws, convection-dominated systems, and large parallel simulations.

The method also introduces additional complexity. Fluxes, quadrature, time-step restrictions, limiters, and boundary conditions must be selected carefully. A reliable implementation should be developed gradually, beginning with a simple linear equation and verified through convergence and conservation tests.

For researchers who already understand finite volume or finite element methods, DG provides a natural next step toward high-order, locally conservative PDE solvers in scientific Python.