Discontinuous Galerkin Discretization

GUERNICA discretizes configuration space using a discontinuous Galerkin (DG) method. For each discrete velocity \(\mathbf{v}_j\) from the DVM, DG produces an independent advection equation in configuration space, with coupling across velocities entering only through collisions.

This page serves as the reference for the DG weak form and the tensor-product / sum-factorized algebra used by GUERNICA’s matrix-free partial-assembly implementation.

PDE and notation

For each discrete velocity \(\mathbf{v}_j\), the transport equation is

\[ \frac{\partial f_j}{\partial t} + \nabla \cdot (\mathbf{v}_j f_j) = R_j, \]

where \(f_j(\mathbf{x},t) \equiv f_n(\mathbf{x},\mathbf{v}_j,t)\) and \(R_j\) contains collisions and sources evaluated at \(\mathbf{v}=\mathbf{v}_j\). For transport-only derivations, we temporarily drop the subscript \(j\) and treat \(\mathbf{v}\) as constant.

Let \(\Omega\) be partitioned into elements \(K\). On each element, approximate

\[ f(\mathbf{x},t)\big|_K = \sum_{a=1}^{N_p} f_a^K(t)\, \phi_a^K(\mathbf{x}), \]

with local basis \(\{\phi_a^K\}\). Test functions are taken from the same space.

Weak form with numerical flux

Multiply by a test function \(\phi_a\), integrate over an element \(K\), and integrate the divergence term by parts:

\[ \int_K \phi_a \frac{\partial f}{\partial t}\, d\mathbf{x} \;-\; \int_K \nabla \phi_a \cdot (\mathbf{v} f)\, d\mathbf{x} \;+\; \int_{\partial K} \phi_a \, \widehat{\mathbf{v} f}\cdot \mathbf{n}\, ds = \int_K \phi_a R\, d\mathbf{x}. \]

For constant \(\mathbf{v}\), the numerical flux is chosen as upwind:

\[ \widehat{\mathbf{v} f}\cdot \mathbf{n} = (\mathbf{v}\cdot \mathbf{n})\, f^{\mathrm{up}}, \qquad f^{\mathrm{up}} = \begin{cases} f^{-} & \mathbf{v}\cdot \mathbf{n} \ge 0,\\ f^{+} & \mathbf{v}\cdot \mathbf{n} < 0, \end{cases} \]

where \((-)\) and \((+)\) denote traces from the interior and exterior of \(K\), respectively. Boundary faces are handled by replacing the exterior trace \(f^{+}\) with a prescribed boundary state.

Semi-discrete matrix form (element-local)

Define the element-local coefficient vector \(\mathbf{f}^K(t)\in \mathbb{R}^{N_p}\) with entries \(f_b^K(t)\). The weak form yields the standard semi-discrete DG system

\[ \mathbf{M}^K \frac{d \mathbf{f}^K}{dt} = \mathbf{A}^K(\mathbf{v})\, \mathbf{f}^K + \mathbf{g}^K(\mathbf{v}; f^{-}, f^{+}) + \mathbf{r}^K, \]

where:

Mass matrix

\[ M^K_{ab} = \int_K \phi_a \phi_b \, d\mathbf{x}. \]

Volume advection operator

For constant \(\mathbf{v}\),

\[ A^K_{ab}(\mathbf{v}) = \int_K \nabla \phi_a \cdot (\mathbf{v}\, \phi_b)\, d\mathbf{x} = \sum_{q=1}^{d} v_q \int_K \frac{\partial \phi_a}{\partial x_q}\, \phi_b \, d\mathbf{x}. \]

(With the sign convention above, this term appears on the RHS with a \(+\) sign because the weak form has \(-\int_K \nabla\phi_a\cdot(\mathbf{v}f)\).)

Face flux contribution

The face term can be written as a linear operator applied to traces:

\[ g_a^K = - \int_{\partial K} \phi_a \, (\mathbf{v}\cdot \mathbf{n})\, f^{\mathrm{up}} \, ds. \]

This term couples neighboring elements through the upwind choice of trace values.

Source/collision projection

\[ r_a^K = \int_K \phi_a R \, d\mathbf{x}. \]

In GUERNICA, collisions are handled separately; the purpose of this page is the transport operator and its efficient evaluation.

Quadrature-based operator evaluation (matrix-free / partial assembly)

GUERNICA evaluates the DG transport operator without assembling global matrices by using quadrature rules on each element.

Let \(\{\xi_i\}_{i=1}^{n_q}\) and weights \(\{w_i\}\) define a quadrature rule on the reference element \(\hat{K}\). Using an isoparametric mapping \(\mathbf{x}=\mathbf{x}(\boldsymbol{\xi})\) with Jacobian \(J(\boldsymbol{\xi})\), volume integrals take the form

\[ \int_K (\cdot)\, d\mathbf{x} = \sum_{i=1}^{n_q} (\cdot)\big|_{\boldsymbol{\xi}_i}\, WJ_i, \qquad WJ_i = w_i \, \det(J(\boldsymbol{\xi}_i)). \]

Define the reference-element basis evaluated at quadrature points:

Interpolation matrix \(B \in \mathbb{R}^{n_q \times N_p}\),
\( B_{i b} = \hat{\phi}_b(\boldsymbol{\xi}_i), \)
Reference derivatives \(D^{(r)} \in \mathbb{R}^{n_q \times N_p}\) for each reference coordinate \(\xi_r\),
\( D^{(r)}_{i b} = \frac{\partial \hat{\phi}_b}{\partial \xi_r}(\boldsymbol{\xi}_i). \)

Physical-space derivatives use the chain rule:

\[ \nabla_{\mathbf{x}} \phi_b = J^{-T} \nabla_{\boldsymbol{\xi}} \hat{\phi}_b. \]

At quadrature point \(i\), define the geometric factors

\[ G^{(q,r)}_i = \frac{\partial \xi_r}{\partial x_q}\Big|_{\boldsymbol{\xi}_i}, \]

so that

\[ \frac{\partial \phi_b}{\partial x_q}\Big|_i = \sum_{r=1}^{d} G^{(q,r)}_i\, D^{(r)}_{i b}. \]

With these definitions, the action of the volume operator on a coefficient vector \(\mathbf{f}\) can be evaluated as:

Interpolate \(f\) to quadrature points:
\( f_i = \sum_b B_{i b} f_b \quad\Longleftrightarrow\quad \mathbf{f}_q = B\, \mathbf{f}. \)
Compute reference derivatives at quadrature points:
\( (\partial_{\xi_r} f)_i = \sum_b D^{(r)}_{i b} f_b \quad\Longleftrightarrow\quad \mathbf{d}^{(r)} = D^{(r)}\, \mathbf{f}. \)
Form the physical flux divergence contribution at quadrature points:
\( \big(\nabla \phi_a \cdot (\mathbf{v} f)\big) \;\Rightarrow\; \sum_{q=1}^{d} v_q \sum_{r=1}^{d} G^{(q,r)}_i\, d^{(r)}_i, \) weighted by \(WJ_i\).
Project back to coefficients (apply \(B^T\) or derivative transpose structure, depending on how the weak form is implemented).

This “operator application” viewpoint is what enables partial assembly and sum factorization.

Sum factorization on tensor-product elements (quads)

On quadrilateral elements (2D), GUERNICA uses tensor-product bases and quadrature rules. Let:

1D basis size \(p = N_{1D}\) (e.g., \(p = \text{order}+1\)),
1D quadrature size \(q = Q_{1D}\).

Then:

\(N_p = p^2\),
\(n_q = q^2\).

Tensor-product factorization of \(B\)

Let \(B^{(1D)} \in \mathbb{R}^{q \times p}\) be the 1D interpolation matrix. Then the 2D interpolation matrix is the Kronecker product

\[ B = B_y \otimes B_x, \qquad B_x = B^{(1D)},\;\; B_y = B^{(1D)}. \]

Naively forming \(B\) is unnecessary. Instead, interpolation is performed as successive 1D operations:

Given coefficients \(f_{a,b}\) arranged as a \(p\times p\) array:

Interpolate in \(x\):

\( \tilde{f}_{i,b} = \sum_{a=1}^{p} B_{x}(i,a)\, f_{a,b} \quad (i=1,\dots,q;\; b=1,\dots,p), \)

Interpolate in \(y\):

\( f_{i,j} = \sum_{b=1}^{p} B_{y}(j,b)\, \tilde{f}_{i,b} \quad (i=1,\dots,q;\; j=1,\dots,q). \)

This computes values at all \(q^2\) quadrature points in \(\mathcal{O}(p^2 q + p q^2)\) operations rather than \(\mathcal{O}(p^2 q^2)\).

Tensor-product factorization of derivatives

Let \(D^{(1D)} \in \mathbb{R}^{q \times p}\) be the 1D derivative matrix with entries \((D^{(1D)})_{i,a} = d\ell_a/d\xi(\xi_i)\) for the chosen 1D basis \(\{\ell_a\}\).

Then the reference derivatives can be written as:

\[ D^{(\xi)} = B_y \otimes D_x, \qquad D^{(\eta)} = D_y \otimes B_x, \]

where \(D_x = D^{(1D)}\) and \(D_y = D^{(1D)}\).

Operationally:

To compute \(\partial_{\xi} f\) at quadrature points:
apply \(D_x\) in \(x\), then \(B_y\) in \(y\)
To compute \(\partial_{\eta} f\) at quadrature points:
apply \(B_x\) in \(x\), then \(D_y\) in \(y\)

Both are evaluated via sequential 1D kernels with the same cost scaling as interpolation.

Incorporating geometry and advection speed

At each quadrature point \((i,j)\), physical-space derivatives are computed using geometric factors:

\[ \frac{\partial f}{\partial x} = G^{(x,\xi)}\, \frac{\partial f}{\partial \xi} + G^{(x,\eta)}\, \frac{\partial f}{\partial \eta}, \qquad \frac{\partial f}{\partial y} = G^{(y,\xi)}\, \frac{\partial f}{\partial \xi} + G^{(y,\eta)}\, \frac{\partial f}{\partial \eta}. \]

The advection contribution \(\mathbf{v}\cdot \nabla f\) is then a quadrature-pointwise combination:

\[ (\mathbf{v}\cdot \nabla f)_{i,j} = v_x \left( G^{(x,\xi)}_{i,j} f_{\xi,i,j} + G^{(x,\eta)}_{i,j} f_{\eta,i,j} \right) + v_y \left( G^{(y,\xi)}_{i,j} f_{\xi,i,j} + G^{(y,\eta)}_{i,j} f_{\eta,i,j} \right). \]

Finally, the result is projected back to coefficient space using the transpose sequence (e.g., applying \(B_x^T\), \(B_y^T\), \(D_x^T\), \(D_y^T\) as appropriate for the chosen weak-form implementation).

Face operators in matrix form (tensor-product factorization)

This section writes the DG face term in a matrix-like form and shows how it is evaluated using tensor-product (sum-factorized) face interpolation. The goal is not to derive DG again, but to provide a clean algebraic reference that maps directly onto the kernels in DG_Advection2D.

We continue to assume constant advection velocity \(\mathbf{v}\) on each application (true for each discrete velocity in the DVM).

Face weak form (recall)

For an element \(K\), the DG surface contribution is

\[ g_a^K(\mathbf{v}; f^-, f^+) = - \int_{\partial K} \phi_a \, (\mathbf{v}\cdot \mathbf{n})\, f^{\mathrm{up}} \, ds, \]

where \(f^{\mathrm{up}} = f^{-}\) if \(\mathbf{v}\cdot \mathbf{n} \ge 0\) and \(f^{\mathrm{up}} = f^{+}\) otherwise.

It is often convenient to write this using the standard flux-splitting notation:

\[ \mathbf{v}\cdot \mathbf{n} = (\mathbf{v}\cdot \mathbf{n})^{+} + (\mathbf{v}\cdot \mathbf{n})^{-}, \qquad a^{+} = \max(a,0),\;\; a^{-}=\min(a,0), \]

so that

\[ (\mathbf{v}\cdot \mathbf{n})\, f^{\mathrm{up}} = (\mathbf{v}\cdot \mathbf{n})^{+}\, f^{-} + (\mathbf{v}\cdot \mathbf{n})^{-}\, f^{+}. \]

Then the face term for element \(K\) can be written as

\[ g_a^K = - \int_{\partial K} \phi_a \left[ (\mathbf{v}\cdot \mathbf{n})^{+}\, f^{-} + (\mathbf{v}\cdot \mathbf{n})^{-}\, f^{+} \right] ds. \]

This makes it explicit that: - the interior element contributes through \((\mathbf{v}\cdot \mathbf{n})^{+}\), - the neighbor (or boundary state) contributes through \((\mathbf{v}\cdot \mathbf{n})^{-}\).

Face quadrature and face interpolation matrices

Consider a single face \(F \subset \partial K\). Let \(\{\zeta_\ell\}_{\ell=1}^{q_f}\) be a 1D quadrature rule on the reference face (for 2D quads, each face is 1D), with weights \(\{w^f_\ell\}\). The physical face integral becomes

\[ \int_{F} (\cdot)\, ds = \sum_{\ell=1}^{q_f} (\cdot)\big|_{\zeta_\ell}\, WJ^f_\ell, \qquad WJ^f_\ell = w^f_\ell \, J_F(\zeta_\ell), \]

where \(J_F\) is the face Jacobian (metric term).

Define the face interpolation matrix \(B_f \in \mathbb{R}^{q_f \times N_p}\) that maps element coefficients to face quadrature values:

\[ (B_f)_{\ell b} = \phi_b\big|_{F}(\zeta_\ell). \]

Given the element coefficient vector \(\mathbf{f}\in \mathbb{R}^{N_p}\), the trace values at face quadrature points are

\[ \mathbf{f}_F = B_f\, \mathbf{f}, \qquad (\mathbf{f}_F)_\ell = f(\zeta_\ell). \]

Similarly, for the neighbor element (or boundary state), we denote the values at the same physical face quadrature points by \(\mathbf{f}^+_F\).

Face term as a matrix application

At each face quadrature point \(\ell\), define the scalar normal speed

\[ \lambda_\ell = \mathbf{v}\cdot \mathbf{n}(\zeta_\ell). \]

For straight-sided elements with constant normals per face, \(\lambda_\ell\) is constant on the face; for curved geometry it may vary with \(\ell\).

Define diagonal matrices:

Quadrature weight/metric matrix:
\( W_f = \mathrm{diag}(WJ^f_\ell), \)
Upwind split normal-speed matrices:
\( \Lambda_f^{+} = \mathrm{diag}(\max(\lambda_\ell, 0)), \qquad \Lambda_f^{-} = \mathrm{diag}(\min(\lambda_\ell, 0)). \)

Then the face contribution to the element residual can be written compactly as:

\[ \mathbf{g}_F^K = - B_f^T \, W_f \left( \Lambda_f^{+}\, \mathbf{f}_F^{-} + \Lambda_f^{-}\, \mathbf{f}_F^{+} \right), \]

where \(\mathbf{f}_F^{-} = B_f \mathbf{f}^K\) is the interior trace and \(\mathbf{f}_F^{+}\) is the exterior trace (neighbor or boundary state) sampled on the same face quadrature points.

Expanding the interior term explicitly:

\[ \mathbf{g}_{F,\mathrm{int}}^K = - B_f^T \, W_f \Lambda_f^{+} B_f \, \mathbf{f}^K, \]

and the neighbor/boundary term:

\[ \mathbf{g}_{F,\mathrm{ext}}^K = - B_f^T \, W_f \Lambda_f^{-} \, \mathbf{f}_F^{+}. \]

This makes the structure clear: - the interior contribution resembles a face “mass-like” matrix \(B_f^T W_f \Lambda_f^{+} B_f\), - the exterior contribution is a projection of external trace values back into the element DOFs.

Tensor-product factorization for quadrilateral elements

On quadrilateral elements, the volume basis is tensor-product, and face interpolation is also tensor-product.

Let the 1D basis size be \(p = N_{1D}\), and the 1D face quadrature size be \(q_f = Q_{1D}\). For 2D quads:

element DOFs can be indexed as \(f_{a,b}\) with \(a,b \in \{1,\dots,p\}\),
face quadrature points are indexed by a single index \(\ell \in \{1,\dots,q_f\}\).

1D interpolation and restriction

Let \(B^{(1D)} \in \mathbb{R}^{q_f \times p}\) be the 1D interpolation matrix from nodal DOFs to quadrature points.

A face is obtained by fixing one coordinate at an endpoint of the reference element, e.g.:

left face: \(\xi = -1\),
right face: \(\xi = +1\),
bottom face: \(\eta = -1\),
top face: \(\eta = +1\).

This introduces a restriction in the fixed direction, which can be represented by a 1D “endpoint evaluation” row vector \(e_\pm^T \in \mathbb{R}^{1\times p}\) such that

\[ (e_\pm^T)_a = \ell_a(\pm 1), \]

where \(\ell_a\) are the 1D basis functions.

For nodal Gauss–Lobatto bases, this restriction is particularly simple: \(e_-^T\) selects the first node and \(e_+^T\) selects the last node.

Example: left/right faces (normal in \(x\))

Consider the left face (\(\xi=-1\)) and right face (\(\xi=+1\)). The trace on such a face depends only on the \(y\)-index and is obtained by restricting in \(x\) and interpolating in \(y\):

Restrict coefficients to the face line:
\( \tilde{f}_b^{(\pm)} = \sum_{a=1}^{p} (e_\pm^T)_a \, f_{a,b}, \qquad b=1,\dots,p. \)
Interpolate to face quadrature points in \(y\):
\( (f_F^{(\pm)})_\ell = \sum_{b=1}^{p} B_y(\ell,b)\, \tilde{f}_b^{(\pm)}, \qquad \ell=1,\dots,q_f, \) where \(B_y = B^{(1D)}\).

In matrix form, the face interpolation operator for left/right faces can be written as a tensor product

\[ B_{f,\text{LR}} = B_y \otimes e_\pm^T, \]

but in practice it is applied via the two-step procedure above (sum factorization), avoiding explicit Kronecker products.

The transpose action \(B_f^T\) (projection from face quadrature back to element DOFs) similarly factorizes into:

Apply \(B_y^T\) to accumulate a \(p\)-vector along the face line,
Scatter back into the element DOFs via \(e_\pm\) in the fixed direction.

Example: bottom/top faces (normal in \(y\))

For bottom/top faces (\(\eta = \mp 1\)), the roles of \(x\) and \(y\) swap:

Restrict in \(y\):
\( \tilde{f}_a^{(\pm)} = \sum_{b=1}^{p} (e_\pm^T)_b \, f_{a,b}, \qquad a=1,\dots,p. \)
Interpolate in \(x\):
\( (f_F^{(\pm)})_\ell = \sum_{a=1}^{p} B_x(\ell,a)\, \tilde{f}_a^{(\pm)}, \qquad \ell=1,\dots,q_f, \) where \(B_x = B^{(1D)}\).

Upwinding, orientations, and matching quadrature points

The face formula

\[ \mathbf{g}_F^K = - B_f^T \, W_f \left( \Lambda_f^{+}\, \mathbf{f}_F^{-} + \Lambda_f^{-}\, \mathbf{f}_F^{+} \right) \]

assumes that interior and exterior traces are evaluated at the same physical face quadrature points and in a consistent order.

On structured tensor-product faces, the primary subtlety is orientation: the neighbor element may parameterize the shared face in the opposite direction. Numerically, this means:

\(\mathbf{f}_F^{+}\) may need to be permuted (reordered) to align with the interior ordering before applying the upwind selection and weight matrices.

This is an implementation detail but is critical for correctness; it is discussed in the DG_Advection2D operator documentation, with checkers against reference formulations.

Boundary faces

On a boundary face, the “exterior” state \(\mathbf{f}_F^{+}\) is replaced with a boundary state \(\mathbf{f}_{F,\mathrm{bc}}\) determined by the selected boundary condition.

The same algebra applies:

\[ \mathbf{g}_F^K = - B_f^T \, W_f \left( \Lambda_f^{+}\, \mathbf{f}_F^{-} + \Lambda_f^{-}\, \mathbf{f}_{F,\mathrm{bc}} \right). \]

This makes it explicit that:

outflow (\(\lambda \ge 0\)) depends only on the interior state,
inflow (\(\lambda < 0\)) requires a boundary-prescribed value.

Implementation in GUERNICA

2D face operators and sum factorization

DG_Advection2D.hxx / DG_Advection2D.cxx

The operator evaluates face traces using tensor-product face interpolation, applies upwinding through \((\mathbf{v}\cdot \mathbf{n})^\pm\), and projects face contributions back to element DOFs using the transpose sequence. Face orientation handling ensures exterior traces are aligned with interior quadrature ordering.

Mass matrix and inverse application

After DG discretization, the semi-discrete system on each element \(K\) can be written as

\[ \mathbf{M}^K \frac{d \mathbf{f}^K}{dt} = \mathbf{R}^K(\mathbf{f}), \]

where \(\mathbf{M}^K\) is the element-local mass matrix and \(\mathbf{R}^K\) contains volume, face, source, and collision contributions.

Exact integration and mass-matrix structure

In GUERNICA, all volume integrals are evaluated using quadrature rules of sufficient order to integrate the polynomial basis exactly. As a result, the mass matrix

\[ M^K_{ab} = \int_K \phi_a \phi_b \, d\mathbf{x} \]

is neither diagonal nor block diagonal, even when nodal bases are used.

This choice avoids mass lumping and preserves the accuracy and symmetry of the DG formulation, at the cost of a more expensive inverse mass application.

Inverse mass via iterative solve

Rather than forming or inverting \(\mathbf{M}^K\) explicitly, GUERNICA applies the inverse mass operator by solving

\[ \mathbf{M}^K \mathbf{y}^K = \mathbf{R}^K \]

using a conjugate-gradient (CG) method at each time step (or stage).

The mass matrix is symmetric positive definite and element-local, making CG a robust and efficient choice. Because the same mass matrix is reused across time steps, this solve can be performed efficiently without global coupling.

Role in the overall operator application

From the perspective of the DG transport operator:

Volume and face terms are evaluated in a matrix-free, quadrature-based fashion.
The resulting residual is converted into time derivatives by applying \((\mathbf{M}^K)^{-1}\) through an iterative solve.
No mass matrix is assembled globally.

This structure preserves the exact DG weak form while remaining compatible with matrix-free and sum-factorized operator evaluation.

Why this page exists

This page defines the algebraic building blocks used by GUERNICA’s 2D transport operator:

tensor-product interpolation \(B_x, B_y\)
tensor-product derivatives \(D_x, D_y\)
quadrature-point geometric factors \(G\) and weights \(WJ\)
face interpolation/projection for flux terms

The operator implementation (DG_Advection2D) should be read as a concrete realization of these factorizations, with additional details about:

face orientations and boundary handling,
memory layout (element-major vs velocity-major),
device execution (CPU/GPU), and
correctness checkers against assembled or reference formulations.

Implementation in GUERNICA

DG transport operators

DG_Advection.hxx / DG_Advection.cxx — 1D DG advection
DG_Advection2D.hxx / DG_Advection2D.cxx — 2D DG advection (PA + sum factorization)

The 2D operator is a matrix-free application of the weak form described above, using tensor-product \(B\) and \(D\) operators and quadrature-point geometry factors.