Split-form DG spectral elements

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

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.

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.

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:

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).

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).