---
title: Global Flux Quadrature for Shallow Water DG-SEM
url: https://www.emergentmind.com/topics/global-flux-quadrature
type: topic
---

# Global Flux Quadrature for Shallow Water DG-SEM

Searching arXiv for the primary paper and closely related works on Global Flux Quadrature.
First, retrieving the main paper by arXiv ID 2212.11931.
Global Flux Quadrature is a discretization strategy for hyperbolic balance laws in which the source term is absorbed into an additional flux primitive, so that the governing equation is rewritten in a quasi-conservative form and the discrete steady states are characterized by constant global flux. In the discontinuous Galerkin spectral element method (DG-SEM) for the shallow water equations, this construction replaces the source integral by a local quadrature based on Gauss–Lobatto nodes, yields a discrete well-balanced property without any a-priori knowledge of the steady equilibrium, does not involve the explicit solution of any local auxiliary problem, and is equivalent at steady state to LobattoIIIA collocation of the steady flux ODE [2212.11931].

## 1. Governing equations and discrete setting

The primary setting is the one-dimensional shallow water system with bathymetry \(b(x)\),
\[
\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),
\]
with
\[
U = \begin{bmatrix} h \\ hu \end{bmatrix},\qquad
F(U) = \begin{bmatrix} hu \\ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad
S(U,b,\partial_x b) = \begin{bmatrix} 0 \\ - g h \,\partial_x b \end{bmatrix}.
\]
For the pseudo-1D rotating shallow water model with friction, the source becomes
\[
S = -h \begin{bmatrix} 0 \\ \partial_x \varphi + c_f u + \omega v \\ -\omega u \end{bmatrix},\qquad \varphi = g b.
\]
The classical lake-at-rest equilibrium is
\[
u=0,\qquad h+b=H_0,
\]
so that the free surface \(\zeta=h+b\) is constant [2212.11931].

In the DG-SEM formulation, each element is mapped from a reference interval and discretized with nodal Gauss–Lobatto points \(\{\xi_i\}_{i=0}^p\) and Lagrange basis functions \(\{\phi_i\}\) of degree \(p\). The nodal interpolation is
\[
u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.
\]
The mass matrix is diagonal,
\[
M=\operatorname{diag}(w_i),
\]
with Gauss–Lobatto weights \(w_i\), and the derivative operator is represented by \(D_x\), together with the SBP-scaled form
\[
\widetilde D_x = M D_x M^{-1},
\]
and the boundary sampling matrix
\[
B=\operatorname{diag}(-1,0,\dots,0,+1).
\]

The standard semi-discrete DG-SEM strong residual form is
\[
\dfrac{d \mathbf{U}}{dt} + D_x \mathbf{F} + M^{-1} B ( \widehat{\mathbf{F}} - \mathbf{F}) = 0,
\]
or, with explicit numerical flux
\[
\widehat{F} = \alpha F^+ + (1-\alpha) F^- + \mathcal{D}(U^+ - U^-),
\]
equivalently
\[
\dfrac{d \mathbf{U}}{dt} + \widetilde{D}_x \mathbf{F} + M^{-1} B(\alpha [\![\mathbf{F}]\!]) + M^{-1} B(\mathcal{D} [\![\mathbf{U}]\!]) = 0.
\]

## 2. Global flux construction and quadrature formula

The central idea is the global flux approach: construct an additional flux \(R(U;x)\) as a primitive of the source so that steady states satisfy
\[
F(U)=R(U;x)+\text{const}.
\]
With
\[
R(U;x)=r_0+\int_{x_0}^{x} S(U(s),\varphi(s))\,ds,
\]
the equation is rewritten as
\[
\partial_t U + \partial_x G(U;x)=0,\qquad G=F-R.
\]
At steady state, \(\partial_x F(U)=S(U;x)\) is equivalent to
\[
F(U(x))=F(x_0)+\int_{x_0}^{x} S(U(s);\varphi(s))\,ds,
\]
so constant global flux is the discrete analogue of equilibrium [2212.11931].

In nodal spectral elements, the source primitive is represented through a local cumulative quadrature. If
\[
S_h \approx \sum_{k=0}^{p}\phi_k(x)S_k,
\]
the integration matrix is
\[
\mathcal{I}_{jk}:=\int_{0}^{\xi_j}\phi_k(\xi)\,d\xi.
\]
On an element of size \(h\), the nodal values of the primitive are
\[
\mathbf{R}=\mathbf{R}^- - h\,\mathcal{I}\,\mathbf{S},
\]
where \(\mathbf{R}^-\) is a constant vector, taken as the last value of the left neighbor to avoid jumps for continuous data.

The DG-SEM written in terms of the global flux becomes
\[
\dfrac{d \mathbf{U}}{dt} + \widetilde{D}_x \mathbf{G} + M^{-1} B( \widehat{\mathbf{G}} - \mathbf{G} ) = 0,\qquad \mathbf{G}=\mathbf{F}-\mathbf{R}.
\]
Using the explicit numerical flux and the definition of \(\mathbf{R}\), the semi-discrete form is
\[
\dfrac{d \mathbf{U}}{dt} + \widetilde{D}_x \mathbf{F} + M^{-1} B( \alpha [\![\mathbf{F}]\!] ) + M^{-1} B( \mathcal{D} [\![\mathbf{U}]\!] ) = h \,\widetilde{D}_x \mathcal{I}\,\mathbf{S}.
\]

This yields the global flux quadrature formula:
\[
\int_K \phi_i f_h \, dx \approx h\,(\widetilde{D}_x \mathcal{I}\,\mathbf{f})_i.
\]
Hence the source contribution reduces to the local matrix–vector product
\[
h\,\widetilde{D}_x \mathcal{I}\,\mathbf{S},
\]
which requires only the nodal source values and the precomputed integration tableau \(\mathcal{I}\) defined by the Gauss–Lobatto nodes [2212.11931].

## 3. Collocation equivalence, superconvergence, and well-balancedness

At steady state, the method is governed by the flux ODE
\[
F'(U^*)(x)=S(U^*(x);\varphi(x)),
\]
with initial condition \(F(U(x_0))=F_0\). Provided \(F(U)\) is invertible and the inverse \(U(F)\) is uniquely defined and bounded away from singular or critical points, the global flux Gauss–Lobatto DG-SEM admits a nodally continuous discrete steady state obtained by integrating this ODE with the fully implicit LobattoIIIA collocation method having \(p+1\) stages at the Gauss–Lobatto nodes. On each element,
\[
\mathbf{F}=h\,\mathcal{I}\,\mathbf{S},
\]
hence
\[
F_i-R_i=F_0,\qquad \forall i,
\]
and continuity of \([\![\mathbf{F}]\!]\) implies continuity of \(U^*\) [2212.11931].

The collocation system is
\[
F(\xi_i)=F(\xi_0)+h\sum_{k=0}^{p}\mathcal{I}_{ik}S(U(F(\xi_k)),\varphi(\xi_k)),\qquad i=0,\dots,p.
\]
This is exactly the discrete equation enforced by global flux quadrature at steady state. Under the additional assumption that there exists a flux linearization
\[
F(U)-F(V)=A(U,V)(U-V),
\]
with \(A\) diagonalizable and eigenvalues bounded away from zero,
\[
\min_j |\lambda_j(A)| \ge C_A>0,
\]
the discrete steady state inherits the collocation superconvergence:
- at element endpoints, order \(2p\) in \(h\);
- at internal Gauss–Lobatto nodes, order \(\min(p+2,2p)\).

The well-balanced property is formulated directly at the semi-discrete level. The residual is
\[
\dfrac{d\mathbf{U}}{dt} + \widetilde{D}_x \mathbf{F} + M^{-1} B(\alpha [\![\mathbf{F}]\!]) + M^{-1} B(\mathcal{D} [\![\mathbf{U}]\!]) - h \widetilde{D}_x \mathcal{I} \mathbf{S} = 0.
\]
If
\[
F = h \mathcal{I} S + \text{const}
\]
at the Gauss–Lobatto nodes and
\[
[\![\mathbf{F}]\!] = [\![\mathbf{U}]\!] = 0,
\]
then
\[
\widetilde{D}_x \mathbf{F} - h \widetilde{D}_x \mathcal{I} \mathbf{S}=0,
\]
the interface terms vanish, and the discrete residual is zero. In this sense, any smooth steady state driven by \(F'=S\) is preserved as a collocation steady state [2212.11931].

Exact preservation of lake-at-rest requires a modified nodal source,
\[
S_l = -\begin{bmatrix}
0 \\
g \zeta_l \,\partial_x b_h(x_l) - \partial_x p_h(b)(x_l) + c_f h_l u_l + \omega h_l v_l \\
-\omega h_l u_l
\end{bmatrix},
\]
with \(\zeta=h+b\) and \(p(h)=gh^2/2\). For \(\zeta\) constant and \(u=v=0\),
\[
(h \mathcal{I}\mathbf{S})_i
= g\zeta^*(b_i-b_0) - g\left(\frac{b_i^2}{2}-\frac{b_0^2}{2}\right)
= p(h_i)-p(h_0),
\]
so
\[
\widetilde D_x(F-F_0-h\mathcal I S)=0,
\]
and the lake-at-rest state is preserved exactly.

This construction differs from hydrostatic reconstruction and path-conservative formulations. Hydrostatic reconstruction modifies local reconstructions to enforce exact balance at prescribed equilibria and often requires a priori knowledge of the equilibrium form. Path-conservative formulations redesign flux/source discretizations for nonconservative systems relying on prescribed paths or invariants. Global flux quadrature instead embeds the source integral through a collocation-based primitive \(R\), requires no auxiliary local steady solves, and no explicit reconstruction of the equilibrium [2212.11931].

## 4. Entropy control and cell entropy correction

The method is paired with an entropy framework. For shallow water without bathymetry in the flux, the entropy pair is
\[
\eta(U)=p(h)+hk,\qquad k=\frac{u^2+v^2}{2},
\]
with entropy flux
\[
q(U)=F_\eta(U)=hu(gh+k).
\]
With spatially constant potential, the total entropy and flux are
\[
\eta_{\varphi}(U)=p(h)+hk+h\varphi,\qquad
q_{\varphi}(U)=hu(gh+k+\varphi)=hu(g\zeta+k).
\]
In frictionless cases,
\[
\partial_t \eta_\varphi + \partial_x q_\varphi = 0,
\]
whereas with friction,
\[
\partial_t \eta_\varphi + \partial_x q_\varphi = -2hc_fk \le 0
\]
[2212.11931].

To control entropy production, the scheme introduces a symmetric positive definite artificial viscosity at the cell level:
\[
w_i \dfrac{d U_i}{dt} + \Phi_i + \Psi_i^{L} + \Psi_i^{R} + \alpha_K \, \mathcal{D}_i^K = 0,
\]
with
\[
\mathcal{D}_i^K := \int_K \partial_x \phi_i \, A_0 \, \partial_x W_h \, dx,\qquad
W = \partial_U \eta,\qquad A_0^{-1}=\partial_{UU}\eta.
\]
Dotting by \(W_i^T\) and summing yields the cell entropy equation
\[
|K| \dfrac{d \bar{\eta}_K}{dt} + \Phi_\eta^K + \alpha_K \| \partial_x W_h \|^2_{L^2_{A_0}(K)} = 0.
\]
The correction parameter is chosen as
\[
\alpha_K=\dfrac{\Psi_\eta^K-\Phi_\eta^K}{\| \partial_x W_h \|^2_{L^2_{A_0}(K)}}.
\]

Two numerical entropy fluxes are used. The analytical flux-based option is
\[
\hat q(U^+,U^-)=\lambda q(U^+) + (1-\lambda)q(U^-),
\]
for example \(\lambda=1/2\), with a local entropy balance consistent with analytically balanced steady states. The global-flux-consistent option is
\[
\Psi_\eta^K = \lambda [\![ q ]\!]^R + \int_K W_h^T \partial_x G_h + (1-\lambda)[\![ q ]\!]^L,
\]
which aligns the entropy control with the global flux quadrature and preserves discrete collocation steady states exactly under the correction.

The resulting entropy balance is element-wise:
\[
|K| \dfrac{d \bar{\eta}_K}{dt} + \Psi_\eta^K = 0 \qquad \text{(frictionless)},
\]
and globally dissipative in the frictional case,
\[
\sum_K |K| \dfrac{d \bar{\eta}_K}{dt} = - \sum_K \mathcal{D}_f \le 0.
\]

## 5. DG-SEM implementation and numerical behavior

Implementation requires only a minimal modification of standard DG-SEM. One chooses the polynomial degree \(p\), precomputes Gauss–Lobatto nodes, weights, basis functions, the mass and differentiation matrices,
\[
M=\operatorname{diag}(w_i),\qquad D_x,\qquad \widetilde D_x = M D_x M^{-1},\qquad B=\operatorname{diag}(-1,\dots,0,\dots,+1),
\]
and a surface numerical flux
\[
\hat F = \alpha F^+ + (1-\alpha)F^- + \mathcal D(U^+-U^-).
\]
The volume source contribution is then assembled by precomputing
\[
\mathcal I_{jk}=\int_0^{\xi_j}\phi_k(\xi)\,d\xi
\]
and evaluating the nodal source values \(S_i\), followed by
\[
\text{RHS}_{\text{source}} = h\,\widetilde D_x \mathcal I\,\mathbf S.
\]
The semi-discrete update is
\[
\dfrac{d\mathbf U}{dt}
= -\widetilde D_x \mathbf F
- M^{-1}B(\alpha[\![\mathbf F]\!])
- M^{-1}B(\mathcal D[\![\mathbf U]\!])
+ h\,\widetilde D_x \mathcal I\,\mathbf S,
\]
optionally augmented by the local entropy correction. Explicit Runge–Kutta of order \(p+1\) may be used, with CFL controlled by the eigenvalues \(|u|\pm\sqrt{gh}\). Relative to standard DG-SEM, the additional operations are multiplication by the dense \((p+1)\times(p+1)\) matrix \(\mathcal I\), one extra derivative matvec, and optional local cell correction integrals, so overall complexity and memory footprint remain comparable to standard DG-SEM [2212.11931].

The numerical tests establish several characteristic properties. For steady moving equilibria, the global flux solution is superconvergent, with order \(2p\) at endpoints and order \(\min(p+2,2p)\) at internal nodes. In transient tests, global flux quadrature reduces steady-state error by orders of magnitude versus classical DG-SEM, thereby enabling accurate resolution of small perturbations. For lake-at-rest with nontrivial bathymetry, exact preservation is obtained with the modified source, and perturbations of amplitude down to \(10^{-5}\,\mathrm m\) are resolved without spurious oscillations, whereas non-well-balanced DG-SEM produces errors larger than the perturbation amplitude on the same mesh. Trans-critical flows also perform well numerically despite the fact that the formal superconvergence proof does not cover critical points [2212.11931].

The paper further reports two-dimensional tests, including perturbations of 1D equilibria on 2D domains, stationary vortices with bathymetry, anticyclonic vortex propagation on a \(\beta\)-plane, geostrophic adjustment, and an equatorial Kelvin front. The well-balanced formulation suppresses spurious oscillations present in non-well-balanced schemes, is less diffusive in several of these tests, and exhibits robust entropy behavior.

## 6. Related formulations, extensions, and limitations

Global Flux Quadrature is part of a broader flux globalization program. In finite volume WENO, the same principle is used by reconstructing a global flux \(\mathcal G\) rather than the conservative variables, with a tailored quadrature for the source primitive \(\mathcal R\); exact discrete steady states are then characterized by constant global fluxes [2205.13315]. A finite-difference variant computes the source primitive from multi-step ODE weights, so that discrete steady states align with the underlying Adams–Bashforth or Adams–Moulton integrator and the steady-state accuracy is determined solely by the ODE method order [2501.06155]. For shallow water moment equations with non-conservative products, flux globalization integrates both source terms and non-conservative products into the divergence term, again reconstructing the global flux to preserve steady states without prior analytical knowledge of them [2507.00573].

The same viewpoint has been extended beyond DG and finite volumes. For nodal continuous finite elements on Cartesian grids, Global Flux quadrature replaces element-wise volume operators by line- or surface-integrated operators, yielding constraint-compatible stabilization and vorticity-preserving schemes for linear acoustics [2407.10579]. A closely related multidimensional formulation rewrites the entire spatial operator as a mixed derivative of a single global flux and proves stationarity preservation and steady-state superconvergence for stabilized nodal finite element methods [2510.02928]. A positivity-preserving PAMPA formulation uses Gauss–Lobatto global flux quadrature together with a first-order local Lax–Friedrichs blending to preserve still-water equilibria, positivity of water height, and wet–dry fronts in one-dimensional shallow water models [2510.26862].

Within the DG-SEM shallow water setting, the main assumptions and limitations are explicit. The formal superconvergence proof assumes continuous bathymetry \(b(x)\) and potential \(\varphi(x)\), smooth steady states, and invertibility of \(U(F)\) with bounded inverse away from critical points. Exact lake-at-rest preservation requires the modified nodal source; the standard nodal source yields only a high-order approximation. The two-dimensional extension treated in the paper is dimensionally split, so genuinely multidimensional well-balancing remains future work. Dry/wet fronts are not treated. The paper also identifies several extensions: multidimensional tensor DG-SEM and truly multidimensional global flux tensors, application to other balance laws such as Euler with gravity or MHD, use in finite volume and residual distribution methods, and the incorporation of interface jump terms in \([\![\mathbf R]\!]\) for discontinuous bathymetry [2212.11931].

A plausible implication is that Global Flux Quadrature is best understood not as a single algorithm, but as a structural principle: the source term is embedded into a globally reconstructed primitive so that the discrete operator acts on a quantity that is constant at equilibrium. In the DG-SEM shallow water formulation, this principle is realized by the local Gauss–Lobatto quadrature
\[
h\,\widetilde D_x \mathcal I\,\mathbf S,
\]
the steady-state LobattoIIIA equivalence, and the accompanying cell entropy correction, which together provide high-order accuracy, discrete well-balancing, and controlled entropy production [2212.11931].

Source: https://www.emergentmind.com/topics/global-flux-quadrature