---
title: Split-form DG spectral elements
summary: Space is discretized by a nodal discontinuous Galerkin spectral element method on Gauss–Lobatto points. Its derivative operator satisfies summation by parts, which lets the volume terms be written as differences of symmetric two-point fluxes and gives discrete conservation and entropy statements on curved meshes.
section: theory
order: 2
level: 2
keywords: [DGSEM, summation by parts, flux differencing, collocation, Gauss-Lobatto]
---
<!-- gyre docs: expanded Markdown of theory/dgsem.md (generated tables, verification numbers and chart data are written out). Base level of this page: L2 Practice; passages at another level are bracketed by 'begin level N' / 'end level N' comments and start with [LN Name]. Built from revision ceb1f4b0be. -->

The domain is partitioned into curved quadrilateral or hexahedral elements. Each element is the
image of the reference cube $[-1,1]^d$ under a polynomial mapping $\mathbf x(\boldsymbol\xi)$, and
the solution is a tensor-product polynomial of degree $p$ in $\boldsymbol\xi$, stored by its
values at the $(p+1)^d$ Legendre–Gauss–Lobatto (LGL) nodes. Neighbouring elements are coupled only
through numerical fluxes on shared faces — the solution may jump across them.

<!-- begin level 2 (Practice) -->
**[L2 Practice] Why Gauss–Lobatto collocation**

Using the same LGL points for interpolation and quadrature makes three things true at once:

* the mass matrix is diagonal, $\mathbf M = \operatorname{diag}(w_0,\dots,w_p)$, so no linear
  system is needed to evaluate a time derivative;
* face values are simply nodal values (the end points are nodes), so no interpolation to faces;
* the derivative matrix $\mathbf D$ with $D_{ij} = \ell_j'(\xi_i)$ satisfies **summation by
  parts** (SBP), the discrete integration-by-parts identity below.

The price is that LGL quadrature with $p+1$ points is exact only to degree $2p-1$, one short of
what the mass matrix of degree-$p$ polynomials needs. The split form turns this weakness into a
controlled one.

<!-- end level 2 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] Summation by parts**

$$
\mathbf Q = \mathbf M \mathbf D, \qquad \mathbf Q + \mathbf Q^{\mathsf T} = \mathbf B = \operatorname{diag}(-1, 0, \dots, 0, +1).
$$

For nodal vectors $\mathbf u, \mathbf v$ this is $\mathbf u^{\mathsf T}\mathbf M\mathbf D\mathbf v +
(\mathbf D\mathbf u)^{\mathsf T}\mathbf M\mathbf v = u_p v_p - u_0 v_0$: integration by parts, exactly,
for the *discrete* inner product. Every stability and conservation property of the scheme is a
consequence of this identity; none relies on exact integration. The verification study
"Summation-by-parts operators" checks it to round-off for all degrees in use.

**Verification status: Summation-by-parts operators** (study `operators`): **PASSED**, 1 of 1 checks passed, 0.29 s.

Gauss–Lobatto collocation operators: the SBP identity M D + (M D)^T = B, quadrature exactness to degree 2p-1, exact differentiation of degree-p polynomials, skew-symmetry of the flux-differencing matrix M S and consistency of the p-multigrid transfers.

| Check | Value | Expected | Result |
|---|---|---|---|
| largest defect over all identities and degrees | 6.2172e-15 | < 1e-11 | pass |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T10:22:15-0500; device cpu; quick suite. **Revision check: Same code as this build.** Produced by 39979fcf902ad6116eb550b881b4278d4599c8e2. Every commit between it and the commit this site was built from (ceb1f4b0bed41de53e586df10858079f5a2bcdd8) changes only documentation text, evidence records or bookkeeping: the code is the same. Full record: [operators.json](../data/verification/operators.json).

<!-- end level 3 -->

## From the strong form to flux differencing

In reference coordinates a conservation law reads $J\,\partial_t \mathbf U + \sum_a \partial_{\xi_a}
\tilde{\mathbf F}_a = 0$ with contravariant fluxes $\tilde{\mathbf F}_a = \sum_b (J a^a_b)\,\mathbf F_b$
(see [curvilinear metrics](curvilinear-metrics.md)). The standard strong-form DGSEM applies
$\mathbf D$ to the nodal fluxes and penalises the difference between the numerical flux
$\mathbf F^*$ and the interior flux at the two end nodes of each line of nodes.

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] The split form**

Gyre replaces the derivative of the nodal flux by a *flux-differencing* operator. Along one line
of nodes in direction $a$, at node $i$:

$$
J_i\,\partial_t \mathbf U_i
 + \sum_{m=0}^{p} S_{im}\; \tilde{\mathbf F}^{\#}_a(\mathbf U_i, \mathbf U_m)
 + \frac{\delta_{i\in\partial}}{w_i}\; \mathbf F^*(\mathbf U^-, \mathbf U^+;\mathbf n)\, s_J = 0,
\qquad \mathbf S = 2\mathbf D - \mathbf M^{-1}\mathbf B ,
$$

where $\mathbf F^{\#}$ is a symmetric, consistent **two-point flux** and the contravariant
two-point flux uses the *average* of the metric terms at the two nodes,
$\tilde{\mathbf F}^{\#}_a(\mathbf U_i,\mathbf U_m) = \sum_b \tfrac12\big[(Ja^a_b)_i + (Ja^a_b)_m\big]\,
\mathbf F^{\#}_b(\mathbf U_i,\mathbf U_m)$. The last term lifts the interface flux onto the two
face nodes with the factor $1/w_i$ ($s_J$ is the surface Jacobian, $\mathbf n$ the outward
normal).

Two remarks make this concrete:

* $\mathbf M\mathbf S = \mathbf Q - \mathbf Q^{\mathsf T}$ is skew-symmetric and $\mathbf S$ has a
  zero diagonal. The interior-flux boundary correction of the strong form has been absorbed
  into $\mathbf S$ using consistency, $\mathbf F^{\#}(\mathbf U,\mathbf U) = \mathbf F(\mathbf U)$.
* With the arithmetic-mean flux $\mathbf F^{\#} = \tfrac12(\mathbf F_i + \mathbf F_m)$ (option
  `central`) the scheme is *algebraically identical* to the standard strong-form DGSEM **on
  elements whose metric terms are constant along each node line** (affine elements). On a
  curved element it is not: because the metric terms are averaged as well, the volume term is
  $\sum_m S_{im}\,\tfrac14\big[(Ja)_i + (Ja)_m\big](\mathbf F_i + \mathbf F_m)$, which differs from
  the standard $\sum_m S_{im}\,\tfrac12\big[(Ja\,\mathbf F)_i + (Ja\,\mathbf F)_m\big]$ by
  $-\tfrac14\sum_m S_{im}\big[(Ja)_i - (Ja)_m\big](\mathbf F_i - \mathbf F_m)$. There `central` is
  the metric-split form of Gassner, Winters and Kopriva (2016), which is still conservative
  and free-stream preserving. Other choices of $\mathbf F^{\#}$ correspond to split
  (skew-symmetric) forms of the nonlinear terms, which is where the name comes from.

<!-- end level 3 -->

<!-- begin level 4 (Theory) -->
**[L4 Theory] Conservation and entropy conservation in three lines**

*Conservation.* Multiply the volume term by $w_i$ and sum over $i$: because $\mathbf M\mathbf S$ is
skew-symmetric and $\mathbf F^{\#}$ is symmetric, $\sum_{i,m} (MS)_{im}\,\mathbf F^{\#}(\mathbf U_i,
\mathbf U_m) = 0$. Only the interface fluxes remain, and they cancel between neighbours for any
conservative $\mathbf F^*$.

*Entropy.* Contract with the entropy variables $\mathbf w_i$ instead. Skew-symmetry gives
$$
\sum_{i,m} (MS)_{im}\, \mathbf w_i^{\mathsf T}\mathbf F^{\#}_{im}
 = \tfrac12 \sum_{i,m} (MS)_{im}\, (\mathbf w_i - \mathbf w_m)^{\mathsf T}\mathbf F^{\#}_{im} .
$$
If the two-point flux satisfies Tadmor's condition $(\mathbf w_i - \mathbf w_m)^{\mathsf T}\mathbf
F^{\#}_{im} = \psi_i - \psi_m$, the right-hand side telescopes, by SBP, to boundary values of the
entropy potential — exactly the discrete counterpart of $\int \mathbf w^{\mathsf T}\partial_\xi
\mathbf F = [\psi]$. No quadrature exactness was used, so the statement survives the aliasing of
LGL collocation and, with metric averaging and discretely divergence-free metrics, holds on
curved elements. With an entropy-conservative interface flux the semi-discrete scheme conserves
$\sum w J\,\eta$ to round-off.

*Entropy decay needs a condition on the interface flux.* What remains after the telescoping is,
per interior face, $(\mathbf w^+ - \mathbf w^-)^{\mathsf T}\mathbf F^* - (\psi^+ - \psi^-)\cdot
\mathbf n$. The total entropy decays if this is non-positive for every pair of states, i.e. if
$\mathbf F^*$ is itself *entropy stable*. The `llf` flux is: it is the entropy-conservative
flux minus $\tfrac12\lambda(\mathbf U^+ - \mathbf U^-)$, and $(\mathbf w^+ - \mathbf
w^-)^{\mathsf T}(\mathbf U^+ - \mathbf U^-) \ge 0$ because the entropy is convex. So is
`es_roe`, which subtracts a symmetric positive semi-definite matrix applied to
$\mathbf w^+ - \mathbf w^-$ and keeps the wave-by-wave dissipation of Roe's solver. The default
`roe` flux and `hllc` are not built on the entropy-conservative flux and no such inequality is
proven for them; that their entropy rate is negative in the chart below is a measurement on
the states of the study. The argument also covers only the inviscid terms on interior or
periodic faces, with positive density and pressure; see
[which combinations are provably entropy stable](two-point-fluxes.md#interface-fluxes).

**Chart: Entropy production rate by flux pairing** (Entropy conservation and stability). Bar chart of |d/dt ∫ η| against volume flux / interface flux, x axis category, y axis log; 2 series. The table holds the plotted data. Measured here: with an entropy-conservative volume flux and the dissipation-free interface flux the entropy rate is at most 8.1e-15 in magnitude; with the upwind interface fluxes it lies between -5.5e-02 and -1.8e-02.

| volume flux / interface flux | 2-D | 3-D |
|---|---|---|
| ranocha / ec | 5.9781e-15 | 3.6880e-15 |
| chandrashekar / ec | 8.0969e-15 | 3.5232e-15 |
| kennedy_gruber / ec | 5.9810e-04 | 0.00119758 |
| central / ec | 4.5667e-04 | 0.00744448 |
| ranocha / llf | 0.0550953 | 0.0543174 |
| ranocha / es_roe | 0.0183493 | 0.0177799 |
| chandrashekar / es_roe | 0.0183493 | 0.0177799 |
| ranocha / roe | 0.0183487 | 0.0178301 |
| ranocha / hllc | 0.0183448 | 0.0178184 |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T10:22:32-0500; device cpu; quick suite. **Revision check: Same code as this build.** Produced by 39979fcf902ad6116eb550b881b4278d4599c8e2. Every commit between it and the commit this site was built from (ceb1f4b0bed41de53e586df10858079f5a2bcdd8) changes only documentation text, evidence records or bookkeeping: the code is the same. Full record: [entropy_conservation.json](../data/verification/entropy_conservation.json).

<!-- end level 4 -->

<!-- begin level 3 (Mechanics) -->
**[L3 Mechanics] How the volume term is evaluated**

A naive implementation evaluates $\mathbf F^{\#}$ for all $(p+1)^2$ ordered node pairs per line.
Since $\mathbf F^{\#}$ is symmetric and $S_{ii} = 0$, only the $p(p+1)/2$ unordered pairs are
needed. `Basis1D.pair_tables()` returns index vectors `iL`, `iR` and a matrix `C` so that the
flux-differencing sum is one gather, one batched flux evaluation and one small matrix product
per direction:

```python
f = phys.flux_volume(jnp.take(U, iL, axis=ax), jnp.take(U, iR, axis=ax),
                     0.5 * (jnp.take(Ja[a], iL, axis=ax) + jnp.take(Ja[a], iR, axis=ax)), th)
JR = JR + apply_1d(f, C, ax)          # C[i, q] = S[i, m] for pair q = {i, m}
```

The same code runs for $d = 2$ and $3$: the array has $d$ nodal axes and the loop runs over them.
For $p = 3$ this more than halves the cost of the dominant kernel (6 pairs instead of 16).

<!-- end level 3 -->

<!-- begin level 2 (Practice) -->
**[L2 Practice] Accuracy you should expect**

For smooth solutions the design order of accuracy is $p + 1$. Gyre's verification studies show:

* $p+1$ for linear-in-effect problems and for steady inviscid manufactured solutions on
  Cartesian meshes;
* between $p$ and $p+1$ for the strongly nonlinear unsteady isentropic vortex, and somewhat
  below $p+1$ on strongly warped coarse meshes — reduced rates of this kind are a known
  characteristic of Gauss–Lobatto collocation;
* in viscous, diffusion-dominated flow, reduced rates between $p$ and $p+1$ with pressure the
  slowest variable (see [order of accuracy](../verification/order-of-accuracy.md)).

**Chart: Density error after convection** (Advected density wave). Convergence chart of L2 error of density against h, x axis log, y axis log; 4 series. The table holds the plotted data. Reference slopes drawn: 2 (h^2), 3 (h^3), 4 (h^4), 5 (h^5).

| h | p = 1 | p = 2 | p = 3 | p = 4 |
|---|---|---|---|---|
| 0.125 | 0.0244511 | 0.00125525 | 5.2276e-05 | 2.0736e-06 |
| 0.0625 | 0.00716847 | 1.6201e-04 | 3.4487e-06 | 6.8806e-08 |
| 0.03125 | 0.00188612 | 2.0333e-05 | 2.1592e-07 | 2.1334e-09 |

Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T10:22:56-0500; device gpu · NVIDIA GeForce RTX 3060 Ti; quick suite. **Revision check: Same code as this build.** Produced by 39979fcf902ad6116eb550b881b4278d4599c8e2. Every commit between it and the commit this site was built from (ceb1f4b0bed41de53e586df10858079f5a2bcdd8) changes only documentation text, evidence records or bookkeeping: the code is the same. Full record: [density_wave.json](../data/verification/density_wave.json).

<!-- end level 2 -->