Discrete Leray Projection

• Discrete Leray projection is a finite-dimensional analogue of the continuous Helmholtz–Hodge decomposition, separating divergence-free velocity fields from pressure gradients.
• It is implemented using finite differences, finite elements, or kernel-based methods to enforce incompressibility and accurately simulate fluid dynamics.
• Its effectiveness is demonstrated through rigorous stability, convergence, and error analyses, making it essential for pressure corrections and robust flow computations.

Discrete Leray projection is the finite-dimensional analogue of the Leray–Helmholtz projector that extracts the solenoidal part of a vector field from a decomposition of the form

$$\mathbf{u}=\mathbf{w}+\nabla\phi,\qquad \mathbb{P}\mathbf{u}=\mathbf{w}.$$

In numerical incompressible flow, it is the operator that enforces the divergence-free constraint, separates velocity from pressure-generated gradient fields, and transfers the continuous Helmholtz–Hodge structure to grids, finite element spaces, or meshfree point clouds. In the literature represented here, discrete Leray projections appear as finite-difference Helmholtz–Hodge decompositions on bounded Lipschitz domains, as kernel-based RBF decompositions into divergence-free and curl-free parts, and as finite element orthogonal projections onto discrete divergence-free subspaces built from discrete gradient, divergence, and Neumann operators (Fuselier et al., 2015, Kuroki et al., 2018, Kaltenbach et al., 20 Jul 2025).

1. Continuous operator and discrete counterpart

The continuous starting point is the Leray–Helmholtz decomposition of a sufficiently smooth vector field,

$\mathbf{u}=\mathbf{w}+\nabla\phi,$

where $\mathbf{w}$ is divergence-free and $\nabla\phi$ is curl-free. In this setting the Leray projector satisfies

$$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$

and may be written formally as

$\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$

with the boundary conditions chosen so that the projection is compatible with the physical domain (Fuselier et al., 2015).

In the finite element setting for the unsteady $p$-Stokes equations, the continuous projector is described as the orthogonal projection onto

$$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$

with complementary projection $\mathcal P^\perp=\mathrm{Id}-\mathcal P$, and the corresponding classical decomposition

$$\mathbf L^2(\Omega)=\mathbf W^2_0(\operatorname{div}^0;\Omega)\oplus \nabla W^{1,2}(\Omega)$$

provides the reference structure that the discrete operator is designed to approximate (Kaltenbach et al., 20 Jul 2025).

A discrete Leray projection reproduces this splitting in finite-dimensional form. In finite differences, the discrete decomposition is

$\mathbf{u}=\mathbf{w}+\nabla\phi,$0

or, equivalently, $\mathbf{u}=\mathbf{w}+\nabla\phi,$1 with $\mathbf{u}=\mathbf{w}+\nabla\phi,$2, so that $\mathbf{u}=\mathbf{w}+\nabla\phi,$3 (Kuroki et al., 2018, Maeda et al., 2020). In finite elements, the backbone is the discrete Helmholtz decomposition

$\mathbf{u}=\mathbf{w}+\nabla\phi,$4

and the discrete Leray projector $\mathbf{u}=\mathbf{w}+\nabla\phi,$5 is defined as the orthogonal projection onto $\mathbf{u}=\mathbf{w}+\nabla\phi,$6 (Kaltenbach et al., 20 Jul 2025). In the RBF setting, the decomposition is embedded in a matrix-valued kernel split into divergence-free and curl-free components, and the projection is obtained by extracting the divergence-free part of the interpolant (Fuselier et al., 2015).

This shared structure is what makes the operator central in incompressible discretizations: it enforces a kinematic constraint through a decomposition rather than through a standalone pressure variable.

2. Principal constructions

A first class of constructions is finite-difference based. On a bounded, connected, Lipschitz domain $\mathbf{u}=\mathbf{w}+\nabla\phi,$7, Kuroki–Soga define the grid

$\mathbf{u}=\mathbf{w}+\nabla\phi,$8

prove a discrete Helmholtz–Hodge decomposition with

$\mathbf{u}=\mathbf{w}+\nabla\phi,$9

and impose

$\mathbf{w}$0

The discrete Leray projection is then $\mathbf{w}$1. Their analysis also yields a projection stability bound and the estimate

$\mathbf{w}$2

when $\mathbf{w}$3 (Kuroki et al., 2018).

A related finite-difference construction based on central differences appears in the continuation paper on Chorin’s method. There the discrete Helmholtz–Hodge theorem uses zero Dirichlet data on the divergence-free part together with a zero mean condition on the scalar potential on each parity class,

$\mathbf{w}$4

and again defines the projection by $\mathbf{w}$5. The parity decomposition is necessitated by the fact that the central difference preserves the parity of grid points (Maeda et al., 2020).

A second class is meshfree and kernel-based. In the RBF method for the incompressible unsteady Stokes equations, a scalar radial basis function $\mathbf{w}$6 generates matrix-valued kernels

$\mathbf{w}$7

$\mathbf{w}$8

with

$\mathbf{w}$9

Generalized interpolation then constructs

$\nabla\phi$0

and the discrete Leray projection is

$\nabla\phi$1

The projection is therefore embedded directly in the kernel decomposition rather than recovered from a separate pressure Poisson problem (Fuselier et al., 2015).

A third class is finite element based. For the fully discrete $\nabla\phi$2-Stokes problem, the discrete divergence-free space is

$\nabla\phi$3

and the projection $\nabla\phi$4 is the orthogonal projection onto this space. On $\nabla\phi$5, the paper obtains the representation formulas

$\nabla\phi$6

These formulas make the discrete projector the exact analogue of the continuous Helmholtz splitting at the level of discrete operators (Kaltenbach et al., 20 Jul 2025).

3. Boundary conditions and trace structure

Boundary treatment is one of the main points of divergence between formulations. The RBF-based method emphasizes that the generalized interpolation framework can match both tangential and normal components of divergence-free vector fields on the boundary. If

$\nabla\phi$7

the projected field is constructed so that the correct boundary behavior is enforced in both components. This is the mechanism by which the method enforces incompressibility without any time-splitting or pressure boundary conditions (Fuselier et al., 2015).

Finite-difference formulations on bounded domains encode boundary information differently. In the forward/backward-difference construction the divergence-free part and the potential satisfy homogeneous boundary conditions,

$\nabla\phi$8

which is the discrete counterpart of no-slip plus a homogeneous potential normalization (Kuroki et al., 2018). In the central-difference construction, the divergence-free field still satisfies zero boundary values, but the potential is constrained by zero mean on each parity class rather than by a boundary condition. This formulation also yields a boundary-sensitive estimate showing that, for smooth divergence-free $\nabla\phi$9 with $$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$0,

$$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$1

a fact used later in the error analysis on bounded domains (Maeda et al., 2020).

For the $$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$2-Stokes equations with slip boundary conditions, the discrete projection must coexist with impermeability and perfect Navier slip. These are incorporated either weakly via Lagrange multipliers or strongly in the discrete velocity space. The paper explicitly ties the pressure analysis to the way the discrete Leray projection approximates its continuous counterpart and states that the interplay of boundary conditions and projection stability governs the accuracy of pressure approximations (Kaltenbach et al., 20 Jul 2025).

In a more general FEEC setting, trace preservation becomes an explicit design criterion. The projections $$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$3 constructed for the three-dimensional de Rham complex preserve discrete trace data and commute with the exterior derivative: $$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$4

$$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$5

The paper states that this is not a pure Helmholtz/Leray projector onto divergence-free fields, but on the $$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$6 level it is Leray-compatible in the sense that

$$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$7

so that divergence-free fields are mapped to discrete divergence-free fields while polynomial normal traces are preserved exactly (Ern et al., 30 Apr 2026).

4. Role in projection methods and pressure elimination

In incompressible Stokes and Navier–Stokes discretizations, the discrete Leray projection is typically the operator that removes the pressure-generated gradient component from the momentum equation. For the unsteady Stokes system

$$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$8

application of the continuous projector gives

$$\mathbb{P}\mathbf{u}=\mathbf{w},\qquad \mathbb{P}(\nabla\phi)=0,$$9

because $\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$0. The RBF method realizes this idea directly and discretely through generalized interpolation, so the projection is not obtained by solving a separate Poisson equation for pressure (Fuselier et al., 2015).

In Chorin-type methods, the projection is the second half of a split update. In the finite-difference analysis on bounded Lipschitz domains, the scheme first computes an intermediate velocity $\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$1 and then sets

$\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$2

The paper identifies this as the exact discrete counterpart of the projected Navier–Stokes system

$\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$3

and the projection is the ingredient that makes $\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$4 exactly divergence-free on the grid (Kuroki et al., 2018). The continuation paper retains the same logic in a central-difference framework and enforces

$\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$5

at every time step (Maeda et al., 2020).

For second-order projection methods, the pressure-correction step can be read directly as a discrete Leray projection. In the BDF2 finite element scheme, the correction step is written as

$\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$6

so the projected velocity lies in the kernel of the discrete divergence operator $\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$7. The paper states that the projection step is literally a discrete Leray projection because it removes the discrete gradient part from the intermediate velocity and leaves a velocity in the discrete divergence-free subspace (Weber, 15 Sep 2025).

In a priori error analysis for the $\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$8-Stokes equations, the role of the projector is more explicitly quantitative. The kinematic pressure estimate uses the decomposition

$\mathbb{P}=I-\nabla\Delta^{-1}\nabla\cdot,$9

with the first term controlled by projection approximation and the second by the discrete error equation. This makes the discrete Leray projection part of the pressure norm itself rather than only a mechanism for incompressibility enforcement (Kaltenbach et al., 20 Jul 2025).

5. Stability, convergence, and accuracy

For RBF discretizations of the incompressible unsteady Stokes equations, the discrete projection is reported to deliver high-order convergence in both space and time. The numerical results show convergence between 5th and 6th order in space and up to 4th order in time for model problems in two dimensional irregular geometries (Fuselier et al., 2015).

The finite-difference analysis of Chorin’s method on arbitrary bounded Lipschitz domains establishes unconditional solvability of the implicit intermediate-velocity step for any $p$0, weak convergence of step-function interpolants, and strong $p$1-convergence under the scaling condition

$p$2

The limit is proved to be a Leray–Hopf weak solution of the incompressible Navier–Stokes equations with no-slip boundary condition and forcing $p$3 (Kuroki et al., 2018). The same paper emphasizes that the projection is central to divergence control, energy estimates, and compactness, even though $p$4 and discrete differentiation do not commute.

The continuation study of the fully discrete Chorin method strengthens the long-time picture. It proves time-global solvability and convergence, derives exponential damping estimates, constructs time-periodic discrete solutions by Brouwer’s fixed point theorem under time-periodic forcing, and proves convergence of such periodic discrete solutions to time-periodic Leray–Hopf weak solutions. For smooth exact solutions on bounded domains with Condition A and scaling

$p$5

it obtains

$p$6

The paper attributes this reduced rate to the mismatch between $p$7 and the grid boundary $p$8, together with the fact that the discrete projection error $p$9 is only $$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$0 in the bounded-domain Dirichlet case (Maeda et al., 2020).

For the fully discrete $$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$1-Stokes problem, the main quantitative statement is that the discrete projector approximates the continuous projector with first-order accuracy: $$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$2 A dual estimate of the same order is also proved. These bounds feed directly into the quasi-best approximation and error-decay analysis for the velocity vector field and kinematic pressure, and the paper states that the projection approximation is first-order optimal under the stated stability and regularity assumptions (Kaltenbach et al., 20 Jul 2025).

Higher-order time discretization does not remove the need for a projection argument. The BDF2 projection method on bounded Lipschitz domains is shown to converge, up to subsequence, to a Leray–Hopf weak solution under the minimal assumptions

$$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$3

The analysis uses a discrete energy inequality and a compactness argument based on Simon’s theorem together with refined time-continuity estimates. The result is described as the first rigorous convergence proof for a higher-order projection method under no additional assumptions on the solution beyond those following from the standard a priori energy estimate (Weber, 15 Sep 2025).

6. Scope, related notions, and conceptual extensions

The term “discrete Leray projection” is used most directly for numerical operators that project onto divergence-free subspaces in incompressible flow. Not every discrete construction carrying the name “Leray” is such an operator in the Helmholtz–Hodge sense.

One distinct usage appears in topological data analysis. Decorated mapper graphs are presented as discrete approximations of the cellular Leray cosheaf over the Reeb graph. The paper explicitly states that the analogy with a projection is conceptual rather than an explicitly named projection operator: the continuous side assigns

$$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$4

whereas the discrete side assigns

$$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$5

The resulting approximation satisfies

$$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$6

so the discrete model converges to the true Leray cosheaf as the cover resolution tends to zero (Curry et al., 2023). This is a Leray-type discretization, but not a solenoidal projection for fluid variables.

A different related usage appears in reduced-order modeling. The Leray reduced order model regularizes the Navier–Stokes convective term by replacing the advecting velocity with a ROM differential filter output $$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$7, defined variationally on the POD space, and the paper also defines the ROM $$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$8-projection

$$\mathbf{H}:=\mathbf{W}^2_0(\operatorname{div}^0;\Omega),$$9

Here the central objects are an explicit ROM spatial filter and a reduced-space projection, not a Helmholtz–Hodge projector onto a divergence-free complement, even though the reduced space is built from solenoidal POD modes (Xie et al., 2017).

The FEEC trace-preserving projections occupy an intermediate position. They are not presented as classical Leray projectors, but they preserve differential constraints and boundary traces simultaneously. This suggests a broader viewpoint in which a discrete Leray projection may be regarded as one instance of a more general family of constraint-preserving projections: in incompressible flow the preserved constraint is $\mathcal P^\perp=\mathrm{Id}-\mathcal P$0, while in FEEC the preserved structure is the commuting de Rham complex together with discrete trace data (Ern et al., 30 Apr 2026).

Within numerical PDEs, however, the classical meaning remains precise. A discrete Leray projection is the operator that realizes a Helmholtz-type splitting at the discrete level, enforces incompressibility exactly or approximately within the chosen discretization, and provides the analytical bridge from pressure-corrected or provisional velocities to divergence-free discrete states and, ultimately, to Leray–Hopf weak solutions (Kuroki et al., 2018, Weber, 15 Sep 2025).