PAMPA: Point–Average–Moment Framework

• PAMPA is a high-order framework for hyperbolic conservation laws that integrates point values, cell averages, and additional moments into a locally continuous polynomial approximation.
• It combines conservative cell-average updates with non-conservative point-value evolution using residual distribution techniques to enhance stability and accuracy.
• The method extends to various geometries and systems—including Euler, shallow water, and ideal MHD—by incorporating invariant-domain preservation and well-balanced strategies.

Searching arXiv for recent PAMPA papers and related formulations. Point–Average–Moment PolynomiAl-interpreted, usually abbreviated PAMPA, is a high-order framework for hyperbolic conservation laws in which the numerical solution is represented by point values on element boundaries, cell averages, and, in the general formulation, additional moments, all interpreted as defining a local polynomial or a virtual-element approximation. Across its one-dimensional, triangular, polygonal, well-balanced, and ideal-MHD variants, PAMPA combines a conservative evolution for cell averages with a non-conservative or residual-based evolution for point values, while maintaining a globally continuous representation at the boundary degrees of freedom (Abgrall et al., 2024, Abgrall et al., 2024).

1. Concept and nomenclature

The acronym PAMPA is expanded in the literature as Point–Average–Moment PolynomiAl-interpreted. “Point” refers to pointwise degrees of freedom, typically located at vertices, edge midpoints, or Gauss–Lobatto nodes on edges. “Average” refers to the cell average. “Moment” refers to integral moments of the solution against monomials or related polynomial bases, although in several third-order constructions the only internal moment used is the cell average itself. “PolynomiAl-interpreted” means that these degrees of freedom are read as defining a polynomial approximation, either explicitly through basis functions or implicitly through projection operators (Abgrall et al., 2024, Abgrall et al., 20 Nov 2025, Abgrall et al., 2024).

A recurring structural property is that the discrete solution is sought as globally continuous while remaining locally conservative at the level of cell averages. In the discontinuous-Galerkin reinterpretation on triangles, the scheme is also formulated so that there is no mass matrix to invert in the final algorithm; the mass matrix appears only conceptually and is absorbed into a projection or residual-distribution step (Abgrall et al., 20 Nov 2025). In the virtual-element formulation on polygonal grids, the reconstruction is described as polynomial-free in the sense that gradients and projected polynomials are computed from degrees of freedom without explicitly constructing basis polynomials on each arbitrary polygon (Abgrall et al., 2024).

The name therefore denotes both a data layout and a discretization philosophy. In one dimension it appears as a compact third-order scheme based on interface point values and cell averages; on triangles it is realized by seven degrees of freedom per element; on general polygons it is coupled to Virtual Element Method technology; and in later work it is specialized to bound-preserving, well-balanced, oscillation-controlled, and positivity-preserving settings for systems such as Euler, shallow water, and ideal MHD (Abgrall et al., 2024, Abgrall et al., 14 Feb 2025, Liu et al., 22 Oct 2025).

2. Degrees of freedom and local representation

In the one-dimensional third-order formulation, the mesh consists of intervals $I_{j+\frac12}=[x_j,x_{j+1}]$. The stored degrees of freedom are a point value $u_j$ at each grid point and a cell average $\bar u_{j+\frac12}$ in each interval. These are interpreted as defining a quadratic polynomial on each cell, with basis functions $\ell_0$, $\ell_{1/2}$, and $\ell_1$ chosen so that the polynomial interpolates the two endpoint values and has the prescribed average (Abgrall et al., 2024).

On unstructured triangular meshes, the third-order PAMPA scheme uses seven degrees of freedom per triangle: point values at the three vertices $\sigma_1,\sigma_2,\sigma_3$, point values at the three edge midpoints $\sigma_4,\sigma_5,\sigma_6$, and one cell average $\bar u_K$. The local approximation is written as

$$u_h|_K=\sum_{i=1}^{6}u_{\sigma_i}\,\varphi_{\sigma_i}+\bar u_K\,\bar\varphi,$$

with $u_j$0, vertex basis functions $u_j$1, and midpoint basis functions $u_j$2 in barycentric coordinates (Abgrall et al., 20 Nov 2025). In the shallow-water triangular formulation, an equivalent seven-function basis $u_j$3 is used so that the first six functions interpolate the boundary degrees of freedom and the seventh has unit cell average and vanishes at the six boundary points (Liu, 2024).

For polygonal meshes, the virtual-element formulation defines a local space $u_j$4 with $u_j$5 such that the restriction to each edge is a polynomial of degree $u_j$6, the trace is continuous on $u_j$7, and $u_j$8. The degrees of freedom are point values at vertices, point values at the $u_j$9 internal Gauss–Lobatto nodes on each edge, and internal moments up to degree $\bar u_{j+\frac12}$0,

$\bar u_{j+\frac12}$1

The dimension is

$\bar u_{j+\frac12}$2

and the polynomial interpretation is realized by projectors such as $\bar u_{j+\frac12}$3 and $\bar u_{j+\frac12}$4 (Abgrall et al., 2024).

A useful unifying observation is that the “moment” component is formulation-dependent. In the 1D third-order and triangular third-order schemes, the internal datum is the cell average, i.e. the zeroth moment. In the higher-order and virtual-element formulations, additional cell moments are admitted explicitly, and the approximation space is enlarged accordingly (Abgrall et al., 2024, Abgrall et al., 23 Aug 2025).

3. Semi-discrete evolution of averages and point values

The defining algorithmic split in PAMPA is between the update of cell averages and the update of point values. Cell averages are advanced through a conservative balance. For a polygon $\bar u_{j+\frac12}$5, the semi-discrete average evolution is

$\bar u_{j+\frac12}$6

and on each edge the boundary integral is approximated by Gauss–Lobatto quadrature using the point values on that edge. Because the reconstruction is globally continuous, the physical flux $\bar u_{j+\frac12}$7 can be used directly on interior edges, and no numerical Riemann solver is required there (Abgrall et al., 2024).

Point values are updated by a non-conservative or residual-distribution relation. In the polygonal VEM-active-flux formulation,

$\bar u_{j+\frac12}$8

with high-order residuals built from the Jacobian $\bar u_{j+\frac12}$9, a reconstructed gradient $\ell_0$0, and an upwind matrix splitting $\ell_0$1 (Abgrall et al., 2024). In the one-dimensional third-order construction, the point update is expressed through high-order residuals derived from the derivative of the local quadratic polynomial and a Jacobian splitting $\ell_0$2 or $\ell_0$3 for systems (Abgrall et al., 2024). On triangles in the DG reinterpretation, the elementwise DG load is first computed and then projected back onto shared globally continuous point values with weights $\ell_0$4, yielding a high-order residual-distribution form (Abgrall et al., 20 Nov 2025).

Several variants add a low-order companion scheme. For averages, this is typically a monotone finite-volume update based on a local Lax–Friedrichs or Rusanov flux. For points, it is a low-order residual-distribution scheme constructed on subcells or dual control volumes around the point degrees of freedom (Abgrall et al., 2024, Abgrall et al., 2024, Abgrall et al., 20 Nov 2025). These low-order operators serve as fallback schemes in MOOD, convex blending, or bound-preserving constructions.

Time integration is generally performed by strong-stability-preserving Runge–Kutta methods. The triangle DG paper states that the time integration is done via SSP Runge–Kutta and residual distribution rather than characteristic backtracking, and the 1D and MHD positivity-preserving formulations are explicitly designed to remain compatible with arbitrarily high-order SSP time integration (Abgrall et al., 2024, Liu et al., 22 Oct 2025).

4. Relationship to active flux, VEM, DG, and SBP

PAMPA belongs to the active-flux family in the precise sense that it evolves cell averages together with point values on cell boundaries, and computes fluxes from those boundary point values. The polygonal VEM paper describes it as a VEM- and RD-based generalization of active flux: it retains the active-flux structure of boundary point values, cell averages, and global continuity, but replaces characteristic backtracking by residual distribution and uses VEM projectors to obtain polynomial information on arbitrary polygons (Abgrall et al., 2024).

The Virtual Element Method enters through the projector-based reconstruction. In the polygonal setting, the scheme does not require a reference-element polynomial basis for each cell; instead, the coefficients of the projected polynomial are computed from degrees of freedom by solving local systems involving $\ell_0$5- or $\ell_0$6-projection operators (Abgrall et al., 2024). This is why several PAMPA papers describe the method as especially suited to arbitrary polygonal grids, including dual grids, Voronoi-type grids, and regular hexagon-like polygons (Abgrall et al., 2024, Abgrall et al., 14 Feb 2025).

A second major reinterpretation identifies PAMPA with a projected discontinuous Galerkin method. On unstructured triangles, the DG formulation is written in the local PAMPA basis, the DG residuals are computed elementwise, and then the discontinuous boundary updates are projected to a single globally continuous value at each shared point through upwind-biased weights (Abgrall et al., 20 Nov 2025). Later analysis generalizes this perspective and states that AF/PAMPA schemes can be interpreted so that the discontinuous Galerkin scheme is one of their building blocks, not only in 1D but also for triangles and polygonal/VEM settings (Abgrall et al., 23 Aug 2025).

In one dimension, a further structural result is that PAMPA has the Summation by part property. The 2025 analysis derives an SBP operator for the third-order PAMPA scheme and shows that, after projection to the globally continuous AF/PAMPA space, the resulting global derivative operator satisfies the periodic SBP identity (Abgrall et al., 23 Aug 2025). This links PAMPA to the finite-difference and DG SBP tradition and provides a discrete integration-by-parts structure for the one-dimensional formulation.

One persistent misconception is that PAMPA always requires an explicit high-order reconstruction stage in implementation. The one-dimensional bound-preserving paper states that its construction does not need any explicit reconstruction within each element, and the ideal-MHD paper emphasizes a polynomial-interpreted but implementation-level reconstruction-free design (Abgrall et al., 2024, Liu et al., 22 Oct 2025). The polynomial interpretation is therefore conceptual and algebraic rather than necessarily explicit.

5. Invariant-domain preservation, positivity, oscillation control, and equilibrium preservation

A central development in PAMPA has been the incorporation of invariant-domain and bound-preserving mechanisms. In one dimension, a bound-preserving PAMPA scheme is obtained by blending a first-order BP scheme with the third-order PAMPA update. The resulting coefficients $\ell_0$7 are derived explicitly, and for the Euler equations they are built using the geometric quasi-linearization framework so that both density and internal energy remain positive for cell averages and point values (Abgrall et al., 2024).

A complementary one-dimensional line of work develops a simpler invariant-domain-preserving framework. It first analyzes the role of the midpoint value in the original PAMPA average update, introduces a local scaling limiter for cell midpoint values, replaces the continuous-flux average update by an IDP numerical flux, and then uses an automatic variable reformulation inspired by Softplus and Clipped ReLU so that evolving point values is limiter-free and invariant-domain-preserving by construction (Abgrall et al., 2024). This yields two distinct stabilization philosophies within PAMPA: explicit convex blending between high- and low-order operators, and structural IDP through midpoint control and variable mapping.

On polygonal meshes, the bound-preserving formulation extends the blending strategy to both cell averages and point values. High-order PAMPA residuals and fluxes are convexly combined with local Lax–Friedrichs operators, and for the Euler equations the admissibility constraints are again expressed through GQL-based linear forms so that density and internal energy stay positive on general polygonal meshes (Abgrall et al., 14 Feb 2025). On unstructured triangular meshes in DG form, the robust PAMPA scheme supplements bound preservation with an oscillation-eliminating parameter $\ell_0$8 built from jumps of higher-order derivatives across element interfaces; the effective blending parameters are the minima of the BP and OE parameters, and boundary conditions are handled within the DG weak formulation (Abgrall et al., 20 Nov 2025).

Later analysis isolates intrinsic bound preserving properties of AF/PAMPA averages and gives numerical illustrations in which limiting is applied only to point values (Abgrall et al., 23 Aug 2025). This suggests that the necessity of average limiting is formulation-dependent rather than universal, although the strongest nonlinear guarantees in the current literature still rely on explicit convex limiting, midpoint limiting, or both.

The well-balanced extensions are equally prominent. For the two-dimensional shallow water equations with horizontal temperature gradients on unstructured triangular meshes, PAMPA uses pressure-momentum-temperature variables in the non-conservative update and suitable Gaussian quadrature in the conservative update to prove well-balancedness for both lake at rest and isobaric steady states (Liu, 2024). For one-dimensional shallow water balance laws, a global flux quadrature reformulation yields a well-balanced PAMPA method that preserves a large family of smooth moving equilibria, including supercritical and subcritical flows, preserves still water exactly by an appropriate quadrature strategy for the source, and combines this with positivity-preserving blending for water depth and wet-dry fronts (Abgrall et al., 30 Oct 2025).

The ideal-MHD extension introduces a positivity-preserving PAMPA scheme with a discrete divergence-free magnetic field on Cartesian grids. Its key ingredients are a limiter-free positivity-preserving update of interface point values through a new nonconservative reformulation, a local DDF projection, a positivity-preserving limiter only at the cell center, a PP flux with appropriate wave-speed bounds, and a suitable discretization of the Godunov–Powell source term. The PP proof again uses GQL, and the method is described as the first active-flux-type ideal-MHD method that is rigorously positivity-preserving for both cell averages and interface point values while maintaining DDF throughout (Liu et al., 22 Oct 2025).

6. Applications, accuracy, and scope

The application range of PAMPA is broad. The polygonal VEM-active-flux paper treats scalar conservation laws, the acoustics system, and the Euler equations in two dimensions on general polygonal meshes, with benchmarks involving smooth solutions, shock waves, and other discontinuities (Abgrall et al., 2024). For the scalar convection test on triangular and polygonal meshes, the reported slopes are around $\ell_0$9–$\ell_{1/2}$0 in $\ell_{1/2}$1 and $\ell_{1/2}$2 norms for both averages and point values, confirming third-order behavior in smooth problems (Abgrall et al., 2024). The same work reports long-time preservation of a steady vorticity field for acoustics up to $\ell_{1/2}$3, and Euler results for a moving vortex, Liu–Lax and Kurganov–Tadmor two-dimensional Riemann problems, and Double Mach Reflection (Abgrall et al., 2024).

The one-dimensional BP literature covers linear advection, Burgers, Buckley–Leverett, Jiang–Shu’s non-oscillatory test, Sod’s shock tube, Shu–Osher, LeBlanc, double rarefaction, Sedov blast, and near-vacuum isentropic problems. In these tests the BP PAMPA schemes remain stable in regimes where the unlimited high-order PAMPA scheme can violate density or pressure positivity, and the reported convergence orders for smooth problems are typically around $\ell_{1/2}$4–$\ell_{1/2}$5, with some degradation under very strict global bounds (Abgrall et al., 2024, Abgrall et al., 2024). The robust triangular DG variant reports convergence order approximately $\ell_{1/2}$6–$\ell_{1/2}$7 on linear advection, and adds Zalesak’s rotating disk, a KPP rotational non-convex flux problem, Kurganov–Tadmor, Double Mach Reflection, and shock diffraction at a convex corner, with bound preservation and substantially reduced oscillations near shocks (Abgrall et al., 20 Nov 2025).

For shallow water, PAMPA has developed into a well-balanced balance-law methodology rather than only a conservation-law method. On unstructured triangular meshes with temperature gradients, the method is validated for high-order accuracy, well-balancedness, and robustness (Liu, 2024). In one dimension, the global-flux-quadrature variant is designed for the Saint–Venant system with Manning friction term and rotating shallow water equations; it is stated to preserve hydrostatic and non-hydrostatic equilibria, preserve a large family of smooth moving equilibria in a super-convergent manner, and exactly preserve still water states by an appropriate source quadrature (Abgrall et al., 30 Oct 2025). The detailed numerical section further reports fourth-order behavior for discrete moving equilibria derived from LobattoIIIA or sub-cell LobattoIII collocation (Abgrall et al., 30 Oct 2025).

The MHD extension markedly enlarges the regime of applicability. The positivity-preserving DDF PAMPA scheme is tested on Alfvén waves, smooth vortices, Orszag–Tang, rotor, blast waves, shock–cloud interaction, and astrophysical jets, including a blast wave with plasma $\ell_{1/2}$8 and jets up to Mach $\ell_{1/2}$9 (Liu et al., 22 Oct 2025). The reported outcome is high-order accuracy in smooth cases together with sharp MHD-structure resolution and strong-shock robustness (Liu et al., 22 Oct 2025).

Taken together, these results place PAMPA at the intersection of active flux, residual distribution, VEM, DG, and modern invariant-domain analysis. Its distinctive feature is not a single spatial discretization formula, but a reusable pattern: point values and averages are advanced together, interpreted through a polynomial or projected polynomial structure, and then supplemented as needed by convex limiting, IDP mappings, equilibrium-preserving quadrature, or DDF projection according to the target PDE class (Abgrall et al., 2024, Abgrall et al., 23 Aug 2025).