Monolithic Convex Limiting (MCL)

• Monolithic convex limiting (MCL) is a family of methods that decomposes high-order updates into a robust low-order invariant-domain-preserving scheme plus antidiffusive corrections.
• It enforces prescribed admissibility constraints through convex blending coefficients and localized optimization, ensuring conservation and bound preservation.
• MCL is applied in various discretization frameworks—including finite elements, DG, and DGSEM—to achieve high-order accuracy while controlling entropy, positivity, and physical bounds.

Monolithic convex limiting (MCL) is a family of property-preserving limiting procedures for high-order discretizations of conservation laws and related PDE systems. Its defining pattern is the construction of a robust low-order invariant-domain-preserving or bound-preserving scheme, the representation of the high-order target scheme as that low-order core plus antidiffusive corrections or subcell flux differences, and the selection of convex blending coefficients or limited corrections so that prescribed admissibility constraints remain satisfied while conservation is retained. In the finite-element and discontinuous-Galerkin literature, the term monolithic distinguishes a single nonlinear residual or semi-discrete ODE from predictor-corrector flux-corrected transport workflows; in later variants, it also denotes the solution of a single convex optimization problem enforcing bounds and conservation (Guermond et al., 2018, Hajduk, 2020, Rueda-Ramírez et al., 2023, Lin et al., 2023).

1. Foundational framework

The modern MCL literature is anchored in invariant-domain preservation. In the discretization-independent formulation of Guermond, Nazarov, Popov, and Tomas, a low-order approximation is designed so that discrete updates are convex combinations of admissible states produced by Riemann averages and source-term contributions; higher-order methods are then corrected by convex limiting when they violate the prescribed invariant-domain properties. A central point of that framework is that the bounds enforced at each time step are necessarily satisfied by the low-order approximation, and after limiting the resulting methods satisfy all invariant-domain properties imposed by the user (Guermond et al., 2018).

This framework is especially explicit for compressible Euler flows. The admissible set used in the second-order convex-limiting method of 2017 is

$$\mathcal B := \{U : \rho>0, e>0, s(\rho,e)\ge s_{\min}\},$$

so the enforced constraints are positivity of density, positivity of internal energy, and a local minimum principle on specific entropy. The limiting step acts on skew-symmetric flux differences between low-order and high-order updates, and conservation follows from that skew symmetry together with symmetric limiter coefficients (Guermond et al., 2017).

A recurrent structural ingredient is the use of bar states or intermediate states. In low-order graph-viscosity or Rusanov-type schemes, these states are constructed so that they lie in a convex invariant set whenever the neighboring nodal or cell states do. MCL then restricts the antidiffusive part of the high-order method so that the corrected intermediate states remain in the same admissible set. This bar-state interpretation appears in continuous finite elements, discontinuous Galerkin methods, DGSEM subcell formulations, shallow-water closures, and realizability-preserving radiative-transfer discretizations (Hajduk, 2020, Rueda-Ramírez et al., 2023, Moujaes et al., 9 Sep 2025).

2. Algebraic structure and the meaning of “monolithic”

In its standard finite-element or DG form, MCL starts from a decomposition of the high-order residual into a low-order invariant-domain-preserving residual plus antidiffusive fluxes. For the discretization-independent and finite-element formulations, this leads to updates of the type

$$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$

where $A_{ij}^n=-A_{ji}^n$ are flux differences between low-order and high-order schemes. The limiter coefficients are chosen by checking quasiconcave constraints or convex invariant-set conditions along one-dimensional segments from the low-order state toward the unlimited high-order candidate (Guermond et al., 2018, Guermond et al., 2017).

For DG discretizations, Hajduk’s formulation makes the same idea explicit at the level of volume and face terms. Raw antidiffusive corrections $f_{ij}^e$ and $f_{i,k}^e$ are defined so that they recover the target DG scheme from the low-order mass-lumped, graph-viscosity method. Limiting is then applied directly to these semi-discrete corrections, yielding a flux-corrected residual

$$R_i^{\mathrm{MCL}} = R_i^{\mathrm{LO}} + \sum_{j\in\mathcal N_i} f_{ij}^{e,*} + \sum_{k\in\mathcal F_i} f_{i,k}^{e,*}.$$

Because the limiting is imposed on the semi-discrete ODE rather than on a separate predictor state, steady states of the DG scheme remain fixed points of the MCL-modified system, and convergence to steady-state solutions is not inhibited by the algorithm (Hajduk, 2020).

The same distinction is stressed in the subcell Bernstein finite-element formulation. There, the final scheme is written as a single residual

$$m_i\dot u_i = L_i^{\rm LO}(u_h) + \sum_{K\ni i}\sum_{j\neq i}f_{ij}^{K,*} + \tilde b_i(u_h,u_{\rm in}),$$

and the paper emphasizes that, unlike predictor-corrector AFC, there is only one nonlinear residual $R_i(u)$; no separate predictor step is needed (Kuzmin et al., 2019).

This monolithic character has two practical consequences. First, conservation is preserved by construction because pairwise or elementwise antidiffusive terms are skew-symmetric or zero-sum. Second, the limiter can be analyzed at the level of the semi-discrete dynamics, which is why many MCL schemes combine naturally with SSP Runge-Kutta methods and preserve invariant-domain properties stage by stage (Hajduk, 2020, Rueda-Ramírez et al., 2023).

3. DGSEM and the entropy-stable subcell formulation

The extension of MCL to Legendre-Gauss-Lobatto DG spectral element methods exploits the natural subcell-flux decomposition induced by collocation and the SBP structure. In the DGSEM formulation of 2023, subcell fluxes extracted from the flux-differencing representation are blended with compatible low-order Rusanov fluxes. Interface-wise limiting of the antidiffusive component yields blended subcell fluxes

$$f^{\mathrm{MCL}}_{(i,j)} = \alpha_{(i,j)} f^{hi}_{(i,j)} + (1-\alpha_{(i,j)}) f^{lo}_{(i,j)}, \qquad 0\le \alpha_{(i,j)}\le 1,$$

and the corresponding limited bar states are constrained to satisfy positivity, local bounds on scalar quantities, and optionally Tadmor-type entropy conditions (Rueda-Ramírez et al., 2023).

The 2023 entropy-stable subcell paper strengthens this construction by formulating the subcell limiting factors themselves as the solution of an optimization problem. On each element,

$$M\,\dot u^H = \Delta\,\bar f^H, \qquad M\,\dot u^L = \Delta\,\bar f^L,$$

and the monolithic subcell-limited residual is

$$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$0

where $$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$1 collects subcell limiting factors $$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$2. Because the high-order and low-order methods share the same interface numerical fluxes, this construction is locally conservative and reduces to the high-order method when all $$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$3 (Lin et al., 2023).

For a convex entropy $$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$4 with entropy variables $$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$5 and entropy potential $$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$6, the cell entropy balance reduces to a single linear inequality on the volume contributions,

$$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$7

with

$$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$8

The limiter then chooses the least restrictive admissible subcell blending by solving

$$U_i^{n+1} = U_i^{L,n+1} + \frac{1}{m_i}\sum_{j}\ell_{ij}A_{ij}^n, \qquad \ell_{ij}=\ell_{ji}\in[0,1],$$9

subject to the entropy inequality and the bounds $A_{ij}^n=-A_{ji}^n$0, where $A_{ij}^n=-A_{ji}^n$1 are upper bounds enforcing other convex constraints such as positivity or maximum principles (Lin et al., 2023).

In multiple dimensions, the same construction is applied dimension by dimension, solving one linear program per coordinate direction. The paper identifies each LP as a continuous knapsack problem and gives a deterministic greedy algorithm that is provably optimal and runs in $A_{ij}^n=-A_{ji}^n$2 time, dominated by sorting. Feasibility is immediate because choosing $A_{ij}^n=-A_{ji}^n$3 recovers the low-order scheme, which already satisfies the semi-discrete cell-entropy inequality. The resulting framework enforces a discrete cell-entropy inequality, preserves arbitrary convex constraints through the bounds $A_{ij}^n=-A_{ji}^n$4, remains locally conservative, and maximizes proximity to the unlimited high-order discretization (Lin et al., 2023).

4. Constraint classes and admissibility notions

The earliest MCL formulations focused on discrete maximum principles for scalar transport and invariant domains for systems, but the constraint set has broadened substantially. In the DG and continuous finite-element literature, admissibility conditions include positivity of density and energy, local minimum principles on entropy, local maximum principles on scalar fields, and generalized convex invariant sets described by quasiconcave functionals (Guermond et al., 2017, Hajduk, 2020, Guermond et al., 2018).

In the entropy-stable DGSEM subcell framework, the same linear program can incorporate general convex constraints before entropy is enforced. The paper lists positivity of density and pressure, maximum principles on modified entropy quantities such as $A_{ij}^n=-A_{ji}^n$5, and TVD-like bounds on primary fields. These are encoded through interface-wise upper bounds $A_{ij}^n=-A_{ji}^n$6, after which the entropy LP simply replaces the unconstrained upper bounds by $A_{ij}^n=-A_{ji}^n$7. The monolithic nature of the limiter means that entropy, positivity, and maximum-principle constraints are satisfied simultaneously (Lin et al., 2023).

A distinct but related admissibility notion appears in continuous Galerkin discretizations of the $A_{ij}^n=-A_{ji}^n$8 model of radiative transfer. There, the realizability set is the convex cone

$A_{ij}^n=-A_{ji}^n$9

and the MCL procedure is tailored to keep the moments inside that set. The algorithm is explicitly two-step: first, local bounds are imposed on each conserved variable to avoid spurious oscillations and maintain positivity of the scalar-valued zeroth moment; second, a scalar factor $f_{ij}^e$0 is chosen so that the vector-valued first moment remains realizable. The resulting flux-corrected finite-element scheme is provably invariant-domain preserving and hyperbolic (Moujaes et al., 9 Sep 2025).

For scalar continuous Galerkin discretizations with entropy constraints, MCL can also be combined with algebraic entropy fixes. In that setting, the limited fluxes are required to satisfy both local bounds and a Tadmor-type discrete entropy inequality, and the procedure is described as involving no free parameters. This line of work shows that entropy consistency alone is generally insufficient to prevent violations of local bounds in shock regions, motivating the coupled use of entropy constraints and convex limiting (Kuzmin et al., 2020).

5. Algorithmic realizations and computational complexity

MCL is not a single algorithmic template but a class of closely related constructions. The dominant implementations are local edge- or element-based flux limiting, subcell interface blending, and global convex post-processing. The following variants are representative.

Formulation	Core unknowns	Representative papers
Edge/element flux limiting in FE or DG	$f_{ij}^e$1, $f_{ij}^e$2, or element antidiffusion	(Guermond et al., 2018, Hajduk, 2020, Kuzmin et al., 2019)
DGSEM subcell entropy limiting	subcell factors $f_{ij}^e$3 or $f_{ij}^e$4 from a small LP	(Rueda-Ramírez et al., 2023, Lin et al., 2023)
Global convex post-processing	corrected averages $f_{ij}^e$5 from a constrained minimization	(Liu et al., 2023)

In local FE and DG schemes, limiter coefficients are usually obtained in closed form or by small one-dimensional searches. Examples include clamp-type formulas for scalar bounds, quadratic constraints for pressure or internal energy, and pairwise symmetrization procedures that retain conservation. These schemes are naturally coupled with forward Euler or SSP Runge-Kutta stages, for which invariant-domain preservation follows from convex-combination arguments under CFL-like conditions (Hajduk, 2020, Rueda-Ramírez et al., 2023).

The DGSEM entropy-stable variant replaces ad hoc local formulas by a small linear program per element and coordinate direction. Its deterministic greedy solver is optimal for the continuous-knapsack structure and runs in $f_{ij}^e$6 time. This is a local optimization problem embedded in each element, not a global nonlinear program (Lin et al., 2023).

A different interpretation of MCL appears in the Cahn-Hilliard-Navier-Stokes work, where MCL is a post-processing strategy for raw cell averages. The corrected state is defined as the minimizer of

$f_{ij}^e$7

subject to box constraints $f_{ij}^e$8 and the conservation constraint $f_{ij}^e$9. The nonsmooth convex problem is solved by generalized Douglas-Rachford splitting. Because the proximal maps are available in closed form, each iteration requires one clipping operation, one inner product, and a few saxpy operations, so the cost is $f_{i,k}^e$0 per iteration. The reported performance is at most $f_{i,k}^e$1 iterations per time step and cost at most $f_{i,k}^e$2, with bounds and conservation enforced up to round-off error (Liu et al., 2023).

Implicit realizations also exist. For steady Euler calculations with backward Euler pseudo-time stepping, MCL is combined with a fixed-point iteration satisfying a Krasnoselskii-type theorem. The first iteration gives a linearized semi-implicit solution that is conservative but generally not invariant-domain preserving; further iterations are performed only if non-IDP states are detected, and positivity preservation is used as a stopping criterion (Moujaes et al., 2024).

6. Numerical behavior, applications, and recurring tensions

Across the literature, MCL is used where high-order resolution must coexist with hard admissibility constraints. In DG and DGSEM simulations of the compressible Euler equations, reported applications include shock tubes, double-Mach reflection, shallow-water dam break, Sedov blast waves, supersonic bow shocks, Kelvin-Helmholtz instability, and hypersonic jets (Hajduk, 2020, Rueda-Ramírez et al., 2023, Lin et al., 2023).

The entropy-stable DGSEM subcell study gives especially sharp evidence for the method’s role in entropy control. In the 1D modified Sod problem and the 2D KPP problem, the pure DGSEM produces non-entropic artifacts such as discontinuous rarefactions and wrong jets, whereas the MCL limiter restores the correct smooth rarefaction and the physically admissible 2D KPP solution. For the 2D isentropic vortex, the MCL-limited DGSEM converges at rates between $f_{i,k}^e$3 and $f_{i,k}^e$4, while simpler limiters based on a minimum-entropy principle stall at $f_{i,k}^e$5–$f_{i,k}^e$6. In Kelvin-Helmholtz and astrophysical-jet simulations, the limiter imposes both a cell-entropy bound and relaxed positivity/TVD bounds, yielding robust non-oscillatory results even at Mach $f_{i,k}^e$7 (Lin et al., 2023).

The DGSEM literature also emphasizes a subtle point about dissipation. Because MCL is imposed at the semi-discrete level, the nonlinear flux formula is independent of the time-step size; numerically, this leads to dissipation that is nearly insensitive to the CFL number, in contrast to FCT/IDP procedures whose dissipation can be highly CFL-dependent (Rueda-Ramírez et al., 2023).

Beyond gas dynamics, continuous-Galerkin MCL has been used for the $f_{i,k}^e$8 model of radiative transfer, where strict enforcement of $f_{i,k}^e$9 is necessary for physical admissibility and hyperbolicity. Numerical tests reported rotational symmetry preservation, behavior extremely close to the realizability boundary without breakdown, and robust transient and steady solutions for lattice-type problems (Moujaes et al., 9 Sep 2025). In the energy-dependent $$R_i^{\mathrm{MCL}} = R_i^{\mathrm{LO}} + \sum_{j\in\mathcal N_i} f_{ij}^{e,*} + \sum_{k\in\mathcal F_i} f_{i,k}^{e,*}.$$0 model for proton therapy dose calculation, MCL is combined with Strang-type operator splitting and SSP-RK transport updates so that all nodal states remain in the interior of the realizable set during backward-in-energy evolution (Moujaes et al., 11 Mar 2026).

The literature also records an important limitation of entropy constraints taken in isolation. In the local subcell DG/FV framework on unstructured grids, fully discrete subcell-wise entropy dissipation for all entropies forces $$R_i^{\mathrm{MCL}} = R_i^{\mathrm{LO}} + \sum_{j\in\mathcal N_i} f_{ij}^{e,*} + \sum_{k\in\mathcal F_i} f_{i,k}^{e,*}.$$1 and reverts to first-order accuracy, while two-point Tadmor entropy stability for a single entropy pair generally reduces the scheme to second order. The same study argues that a semi-discrete, cell-wise entropy inequality enforced through a small knapsack-type constraint is the only entropy constraint in that framework that preserves arbitrary high-order accuracy. It also reports that, for nonconvex problems, entropy-only constraints for one chosen entropy may still pick non-physical solutions, whereas positivity and local maximum-principle conditions are needed to capture the correct entropy solution on coarse meshes (Vilar, 2024).

Taken together, these results characterize MCL less as a single limiter than as a general methodology for embedding convex admissibility into high-order discretizations. The common invariant is the same: a conservative high-order method is decomposed relative to a guaranteed-safe low-order scheme, and the admissible portion of the high-order correction is selected by convexity arguments, local algebra, or small optimization problems so that accuracy is recovered where possible and failure modes are confined by construction (Guermond et al., 2018, Lin et al., 2023).