Convex Limiting in High-Order Discretizations

• Convex limiting is an algebraic stabilization technique that constrains high-order corrections within convex admissible sets to preserve conservation laws and positivity.
• It decomposes the numerical update into a robust low-order scheme and an antidiffusive high-order correction, ensuring invariant-domain preservation through symmetric limiting.
• This approach is applied in CFD, radiative transfer, and shallow-water models, delivering entropy stability, well-balanced treatments, and robust performance across discretization frameworks.

Convex limiting is an algebraic stabilization strategy for high-order discretizations of conservation laws and related PDEs in which the numerical update is constrained so that discrete states remain in a convex admissible set or invariant domain. In the formulations developed by Guermond, Popov, Tomas, Kuzmin, and collaborators, the method starts from a low-order invariant-domain-preserving approximation, represents the difference to a high-order target as conservative antidiffusive corrections, and then limits those corrections by symmetric coefficients or inequality-constrained flux modifications so that positivity, local maximum principles, entropy bounds, or realizability conditions are retained while recovering as much high-order content as possible (Guermond et al., 2018, Kuzmin et al., 2020, Kuzmin, 2 Feb 2026)

1. Invariant sets, admissible states, and bar states

The fundamental object in convex limiting is a convex set of admissible states. For a generic hyperbolic system with source,

$\partial_t u + \nabla\!\cdot f(u)=S(u),$

Guermond–Popov–Tomas define an invariant set $\mathcal B\subset\mathcal E$ so that all relevant Riemann solutions remain in $\mathcal B$, and a numerical mapping is invariant-domain preserving if nodal states in $\mathcal B$ are mapped back into $\mathcal B$ at each step. In scalar conservation laws this set is typically an interval $[u_{\min},u_{\max}]$; for the compressible Euler equations it is built from positivity of density, positivity of internal energy, and lower bounds on specific entropy; for the $M_1$ model of radiative transfer it is the realizable cone

$$R_1=\{u=[\psi^{(0)},\psi^{(1)}]^T:\psi^{(0)}>0,\ |\psi^{(1)}|<\psi^{(0)}\}.$$

These admissible sets are convex, which is the structural reason convex-combination arguments can be used to prove preservation properties (Guermond et al., 2018, Guermond et al., 2017, Moujaes et al., 9 Sep 2025)

A recurring algebraic device is the bar state. In low-order graph-viscosity or LLF/Rusanov-type schemes, one rewrites the update in terms of auxiliary states such as

$$\overline U_{ij} =\tfrac12(U_i+U_j) -\frac{(f(U_j)-f(U_i))\cdot c_{ij}}{2\,d_{ij}},$$

or, in the scalar CG formulation of Kuzmin–Quezada de Luna,

$$\bar u_{ij} =\tfrac12(u_i+u_j)-\frac{c_{ij}\cdot(f_j-f_i)}{2\,d_{ij}^{\max}}.$$

Under the corresponding CFL restriction, the low-order update is a convex combination of nodal states and bar states. The theoretical burden is therefore transferred to proving that bar states themselves belong to the desired invariant set (Guermond et al., 2017, Kuzmin et al., 2020)

2. Low-order/high-order decomposition

Convex limiting is built on a low-to-high decomposition. The low-order scheme is chosen for robustness: Guermond–Popov–Tomas use graph viscosity with coefficients

$\mathcal B\subset\mathcal E$0

while Kuzmin–Quezada de Luna use a low-order LLF flux

$\mathcal B\subset\mathcal E$1

Both constructions yield conservative, first-order, invariant-domain-preserving updates when the time step satisfies the corresponding CFL condition (Guermond et al., 2018, Kuzmin et al., 2020)

The high-order method is then expressed as a correction to the low-order backbone. In Guermond–Popov–Tomas one writes

$\mathcal B\subset\mathcal E$2

where $\mathcal B\subset\mathcal E$3 is the antidiffusive correction. In Kuzmin–Quezada de Luna one decomposes any target flux as

$\mathcal B\subset\mathcal E$4

with $\mathcal B\subset\mathcal E$5 the antidiffusive remainder. In the element-based residual-distribution interpretation, the same role is played by element contributions $\mathcal B\subset\mathcal E$6 satisfying $\mathcal B\subset\mathcal E$7. Antisymmetry or zero-sum structure is essential: it is the discrete mechanism that preserves global conservation after limiting (Kuzmin et al., 2020, Guermond et al., 2018, Kuzmin, 2 Feb 2026)

This decomposition separates accuracy from admissibility. The high-order target may be entropy-consistent, residual-viscosity stabilized, WENO-sensor stabilized, or based on the consistent Galerkin residual, but it is not assumed to preserve invariant domains. Convex limiting is the algebraic stage that restores those properties without discarding the entire high-order update (Kuzmin, 2 Feb 2026)

3. Limiting strategies

Two main algorithmic families appear in the literature. The first is flux-corrected transport in which one computes a low-order predictor and then adds limited antidiffusion. The second is monolithic convex limiting, in which limiting is built directly into the semi-discrete or fully discrete algebraic form. Kuzmin’s 2026 review makes this distinction explicit and notes that both formulations enforce inequality constraints on scalar functions of intermediate states that are required to stay in convex invariant sets (Kuzmin, 2 Feb 2026)

In the pairwise flux formulation of Guermond–Popov–Tomas, one introduces symmetric coefficients $\mathcal B\subset\mathcal E$8 and defines

$\mathcal B\subset\mathcal E$9

The computation of $\mathcal B$0 is reduced to one-dimensional line searches along admissible segments. The framework is particularly general because the constraints are encoded by quasiconcave functionals $\mathcal B$1; examples listed in the paper include density, internal energy, specific entropy, and negative kinetic energy. A key novelty stated in the paper is that the bounds enforced at each time step are necessarily satisfied by the low-order approximation (Guermond et al., 2018)

In edge-based monolithic convex limiting for scalar conservation laws, the corrected bar state

$\mathcal B$2

is required to remain inside neighbor-based bounds. For $\mathcal B$3, Kuzmin–Quezada de Luna impose

$\mathcal B$4

with the corresponding max-formula for $\mathcal B$5. They then add an optional entropy-fix limiter $\mathcal B$6 so that the final flux

$\mathcal B$7

is simultaneously conservative, bound-preserving, and entropy stable (Kuzmin et al., 2020)

More recent high-order continuous finite element variants reconstruct and limit antidiffusion at the element level. In the graph-Poisson formulation, the difference between high-order and low-order updates is represented as minimum-energy antisymmetric graph fluxes $\mathcal B$8 obtained from a graph-Laplacian solve $\mathcal B$9, after which a first limiting stage produces a bounded candidate state $\mathcal B$0, and a second capacity-based redistribution stage restores mass while preserving local bounds whenever sufficient admissible capacity is available (Nazarov, 22 Jun 2026)

4. Preserved properties and analytical guarantees

The primary analytical guarantee is invariant-domain preservation. In the scalar setting, each stage of an SSP Runge–Kutta method applied to the limited system can be shown to be a convex combination of states lying in the local interval $\mathcal B$1, which yields a local maximum principle and prevents new overshoots or undershoots. In the abstract Guermond–Popov–Tomas framework, Theorem 7.24 states that if the low-order update lies in the intersection of admissible sets defined by finitely many quasiconcave functionals, then convex limiting preserves all user-imposed constraints while maintaining conservation (Kuzmin et al., 2020, Guermond et al., 2018)

Entropy stability is a second major theme. Kuzmin–Quezada de Luna derive a sufficient semi-discrete entropy condition,

$\mathcal B$2

and combine minimal entropy diffusion $\mathcal B$3, entropy viscosity $\mathcal B$4, and an optional entropy-fix limiter to obtain a scheme that is simultaneously invariant-domain preserving and entropy stable. For the Euler equations, Guermond et al. use convex limiting to preserve positivity of density, positivity of internal energy, and the local minimum principle on the specific entropy in a second-order method (Kuzmin et al., 2020, Guermond et al., 2017)

Several extensions add problem-specific structure. In steady convection–diffusion–reaction equations, Knobloch, Kuzmin, and Jha modify bar states and raw fluxes by incorporating source-term components so that certain linear equilibria are preserved exactly; for the linear equilibrium $\mathcal B$5, one obtains $\mathcal B$6, and the discrete system is satisfied exactly by the nodal values of $\mathcal B$7. This is the well-balanced property in the sense used in the paper (Knobloch et al., 2024)

For the $\mathcal B$8 model of radiative transfer, the relevant admissible set is realizability. The tailor-made two-step MCL algorithm first imposes local bounds on each conserved component, ensuring positivity of the zeroth moment, and then applies a scalar scaling $\mathcal B$9 so that the corrected bar state satisfies $\mathcal B$0. The resulting fully discrete scheme is proved invariant-domain preserving, including the diagonally implicit treatment of reactive terms (Moujaes et al., 9 Sep 2025)

Implicit formulations require additional fixed-point analysis. For backward Euler discretizations of the compressible Euler equations, monolithic convex limiting is combined with a fixed-point framework satisfying a Krasnoselskii-type theorem. The existence of an invariant-domain-preserving limit is part of the analysis, and positivity preservation of intermediate states is used as a stopping criterion for nonlinear iterations (Moujaes et al., 2024)

5. Discretization frameworks and computational realizations

Convex limiting is not tied to a single discretization family. The published formulations span continuous Galerkin, high-order Bernstein-basis finite elements, residual-distribution methods, nodal DG with subcell corrections, and DGSEM on adaptive meshes.

Framework	Low-order backbone	Limiting object
Continuous $\mathcal B$1 CG for scalar conservation laws	LLF / graph-Laplacian dissipation	Edge fluxes $\mathcal B$2
Steady CDR finite elements	M-matrix low-order operator with source balancing	Edge fluxes and source-modified bar states
Bernstein-basis CG with nonlinear stabilization	Element-average LLF states on macrocells	Facet antidiffusive fluxes and element slopes
High-order CFEM on fine $\mathcal B$3 submesh	Fine-submesh graph-viscosity scheme	Graph-Poisson edge fluxes plus mass redistribution
Sparse DG with subcell FCT	Sparse-stencil IDP DG low-order method	Subcell fluxes
LGL-DGSEM with AMR mortars	Graph-viscosity volume and mortar fluxes	Antidiffusive corrections blended with low-order states

The continuous finite element line is particularly rich. Kuzmin–Quezada de Luna work with standard continuous $\mathcal B$4 or $\mathcal B$5 bases on unstructured meshes and show $\mathcal B$6 complexity per time step. Kuzmin–Hajduk–Vedral use Bernstein polynomials as local basis functions because they are nonnegative and form a partition of unity; this makes element averages and slope-limited nodal values genuine convex combinations of admissible states, and it supports a matrix-free, hardware-aware implementation with explicit SSP-RK updates (Kuzmin et al., 2020, Kuzmin et al., 4 Sep 2025)

The high-order CFEM formulation of 2026 separates low-order admissibility from high-order accuracy by constructing the low-order scheme on a fine $\mathcal B$7 submesh induced by the high-order nodes. The graph-Poisson reconstruction then connects the high-order and low-order updates without requiring direct limitation on wide high-order stencils. Numerical tests show that residual-viscosity+CIP with graph-Poisson convex limiting converges to the correct entropy solution in Burgers’ and KPP benchmarks, whereas a CIP-only high-order variant does not (Nazarov, 22 Jun 2026)

On the DG side, Pazner’s subcell convex limiting uses sparse low-order stencils whose size does not grow with polynomial degree, then limits subcell fluxes in a way that preserves conservation and invariant-domain properties. In nonconforming LGL-DGSEM with AMR, convex limiting requires an invariant-domain-preserving mortar construction; the 2026 AMR work introduces sparsified mortar couplings based on LGL subcell characteristic functions so that nonconforming interfaces fit into graph-viscosity-based low-order schemes used for convex limiting (Pazner, 2020, Bolm et al., 7 Jul 2026)

6. Applications, extensions, and adjacent directions

Convex limiting has been applied to both scalar and system problems. For the compressible Euler equations, Guermond et al. build a second-order method by combining a first-order invariant-domain-preserving GMS-GV1 scheme with a higher-order entropy-consistent method and then convex-limiting the corrections. For steady-state Euler calculations, implicit pseudo-time stepping with monolithic convex limiting and adaptive explicit underrelaxation has been used to obtain high-resolution, non-oscillatory solutions while maintaining positivity of density and internal energy (Guermond et al., 2017, Moujaes et al., 2024)

In machine-learning-assisted shallow-water models, convex limiting functions as a property-preserving wrapper around learned subgrid fluxes. In the 2024 finite-volume study, neural-network flux corrections are limited in fluctuation form so that positivity of water depth and local maximum principles for height and velocity are preserved; when the training set is small, the unconstrained NN-reduced model develops spurious oscillations and slight negative depths, whereas the MCL-NN-reduced model recovers stability, positivity, and a much closer match to DNS. In the 2026 long-term shallow-water parametrization study, flux limiting is combined with a four-point-stencil neural network and is reported to reduce oscillations near shocks while leaving energy spectra intact in the numerical experiments summarized in the paper (Timofeyev et al., 2024, Mojamder et al., 30 Jan 2026)

Problem-specific adaptations continue to widen the scope of the methodology. The realizability-preserving MCL scheme for the $\mathcal B$8 radiative transfer model is motivated by hyperbolicity of the entropy-based closure and by radiotherapy dose calculation. The well-balanced steady CDR scheme shows how source terms can be built directly into intermediate states and fluxes. Kuzmin’s 2026 review also emphasizes that multidimensional convex limiting can be interpreted in residual-distribution terms when antidiffusive element contributions, rather than pairwise fluxes, are constrained (Moujaes et al., 9 Sep 2025, Knobloch et al., 2024, Kuzmin, 2 Feb 2026)

Taken together, these formulations present convex limiting not as a single limiter formula but as a general algebraic principle: construct a robust low-order update whose intermediate states are known to lie in a convex admissible set; represent the high-order correction conservatively; and admit only that portion of the correction that preserves the desired convex constraints. The invariant-domain argument, the use of bar states, and the preservation of conservation through antisymmetric or zero-sum corrections are the common structural features across the continuous Galerkin, DG, AMR, residual-distribution, and machine-learning-enhanced variants surveyed here (Guermond et al., 2018, Nazarov, 22 Jun 2026)