Sharp Collocated Projection Method

• The method introduces a nodal projection operator that enforces divergence-free velocity fields in incompressible flows using adaptive octree grids.
• It leverages supra-convergent finite differences and a hybrid finite difference–finite volume scheme for sharp boundary and interfacial jump treatment.
• Numerical tests confirm stable, second-order accurate results without extra stabilization, even under complex geometries and adaptive refinement.

Searching arXiv for papers on the Sharp Collocated Projection Method and closely related work. The Sharp Collocated Projection Method is a family of collocated, or “nodal,” projection methods for incompressible flow in which all scalar and vector unknowns are stored at the same set of grid nodes, typically on non-graded adaptive quadtree or octree meshes. In the formulation introduced in “Stable nodal projection method on octree grids,” the method targets the incompressible Navier–Stokes equations with arbitrary boundaries, combines supra-convergent finite differences with sharp boundary treatments, and defines a nodal projection operator $P_N = I - G_N (L_N)^{-1} D_N$ for the pressure-correction step (Blomquist et al., 2023). A later extension applies the same collocated quadtree/octree framework to immiscible two-phase Navier–Stokes flow, where interfacial jump conditions are treated sharply by a hybrid finite difference-finite volume methodology (Binswanger et al., 14 Aug 2025). A related antecedent is the variational collocated projection method of Gagniere et al., which enforces incompressibility on collocated velocity grids over regular meshes with cut-cell geometry (Gagniere et al., 2020).

1. Conceptual definition and lineage

In the octree-based formulation, “nodal” and “collocated” are synonymous: all variables live at the same set of nodes, and on an adaptive non-graded octree grid all scalar and vector unknowns $(u,v,(w),p)$ are stored at the mesh nodes (Blomquist et al., 2023). The principal motivation stated for this layout is that it “reduces the overhead in code development through data collocation” while retaining second-order accuracy and supporting dynamic grid adaptivity with arbitrary geometries (Blomquist et al., 2023).

The method is “projection” based in the sense of Chorin splitting: an intermediate velocity is computed first, and then a pressure- or Hodge-variable correction is applied so that the updated velocity is approximately divergence free. In the nodal octree formulation this is expressed by the operator

$P_N = I - G_N (L_N)^{-1} D_N,$

together with the discrete Poisson problem

$D_N[G_N \phi] = D_N u^*.$

Applying $P_N$ to $u^*$ enforces approximately $\nabla \cdot u^{n+1}=0$ (Blomquist et al., 2023).

A related collocated projection construction appeared earlier in a hybrid particle/grid setting on regular grids. There, the projection is formulated variationally over a collocated multiquadratic-B-spline velocity grid and a multilinear-pressure grid, with cut-cell geometry for irregular flow domains (Gagniere et al., 2020). This suggests that “sharp collocated projection” is best understood not as a single stencil, but as a broader design pattern: collocated storage, projection-based incompressibility enforcement, and sharp treatment of embedded boundaries or interfaces.

2. Governing equations and collocated discretization framework

For the single-phase incompressible case, the governing equations in the fluid region $\Omega^- \subset \mathbb{R}^{2,3}$ are

$$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$

solved on adaptive non-graded quadtree or octree meshes (Blomquist et al., 2023).

For the immiscible two-phase extension, two incompressible Newtonian fluids occupy $\Omega = \Omega^+ \cup \Omega^-$, separated by a sharp interface $(u,v,(w),p)$0. In each phase the density and viscosity are constant, with $(u,v,(w),p)$1 in $(u,v,(w),p)$2 and $(u,v,(w),p)$3 in $(u,v,(w),p)$4. The strong form is

$(u,v,(w),p)$5

$(u,v,(w),p)$6

with surface tension represented as

$(u,v,(w),p)$7

where $(u,v,(w),p)$8 is the surface tension coefficient, $(u,v,(w),p)$9, $P_N = I - G_N (L_N)^{-1} D_N,$0, and $P_N = I - G_N (L_N)^{-1} D_N,$1 is the Dirac delta supported on $P_N = I - G_N (L_N)^{-1} D_N,$2 (Binswanger et al., 14 Aug 2025).

The canonical jump conditions across $P_N = I - G_N (L_N)^{-1} D_N,$3 are

$P_N = I - G_N (L_N)^{-1} D_N,$4

$P_N = I - G_N (L_N)^{-1} D_N,$5

equivalently

$P_N = I - G_N (L_N)^{-1} D_N,$6

for the tangential shear component (Binswanger et al., 14 Aug 2025). The interface is represented implicitly by a level-set function $P_N = I - G_N (L_N)^{-1} D_N,$7, reinitialized as a signed distance so that $P_N = I - G_N (L_N)^{-1} D_N,$8 (Binswanger et al., 14 Aug 2025).

In both the single-phase and two-phase formulations, the computational framework is a non-graded adaptive quadtree in two dimensions or octree in three dimensions. Each cell has four or eight children when refined, and refinement depth is bounded by user-specified minimum and maximum levels (Binswanger et al., 14 Aug 2025). At T-junctions, where a node lacks a direct neighbor, ghost-node values are created by a third-order interpolation formula of the form

$P_N = I - G_N (L_N)^{-1} D_N,$9

which is used in both the nodal octree method and the two-phase extension (Blomquist et al., 2023, Binswanger et al., 14 Aug 2025).

3. Discrete differential operators and the projection mechanism

The collocated octree method defines nodal analogues of the Laplacian, divergence, and gradient. At node $D_N[G_N \phi] = D_N u^*.$0, with neighbors left $D_N[G_N \phi] = D_N u^*.$1, right $D_N[G_N \phi] = D_N u^*.$2, bottom $D_N[G_N \phi] = D_N u^*.$3, and top $D_N[G_N \phi] = D_N u^*.$4, the nodal Laplacian is

$D_N[G_N \phi] = D_N u^*.$5

On uniform grids this reduces to the 5-point second-order stencil, and on adaptive grids with ghosts it retains second-order accuracy of both $D_N[G_N \phi] = D_N u^*.$6 and $D_N[G_N \phi] = D_N u^*.$7, described as supra-convergence (Blomquist et al., 2023).

The divergence operator is

$D_N[G_N \phi] = D_N u^*.$8

and is stated to be formally only first-order accurate on adaptive meshes, chosen for stability (Blomquist et al., 2023). The gradient operator is

$D_N[G_N \phi] = D_N u^*.$9

The same nodal stencil family is reused in the two-phase method as part of its sharp hybrid FD–FV discretization (Binswanger et al., 14 Aug 2025).

The central projection step is the discrete Hodge projection. In the single-phase method, the scalar Hodge variable $P_N$0 is obtained from

$P_N$1

after which the updated field is obtained by repeated application of $P_N$2 (Blomquist et al., 2023). In the two-phase method, the projection step becomes a variable-density Poisson problem: $P_N$3 with

$P_N$4

and Dirichlet or Neumann conditions on $P_N$5 (Binswanger et al., 14 Aug 2025).

No extra pressure-velocity correction terms are needed beyond the iterative projection in the nodal octree formulation (Blomquist et al., 2023). This is directly relevant to the standard collocated-grid concern about pressure–velocity decoupling.

4. Sharp treatment of boundaries and interfaces

A defining feature of the method is its insistence on sharp, rather than smeared, treatment of geometric constraints. In the single-phase octree method, a level-set $P_N$6 defines $P_N$7 versus solids $P_N$8. On cut cells the hybrid finite-volume/finite-difference discretization “imposes no-slip $P_N$9 exactly on $u^*$0” and homogeneous Neumann $u^*$1 for the Hodge variable (Blomquist et al., 2023). Curved and moving interfaces are represented by a level set on the octree, and cut-cells near $u^*$2 are handled sharply via hybrid FV/FD stencils that respect the exact geometry and impose $u^*$3 on $u^*$4 without smearing (Blomquist et al., 2023).

In the two-phase extension, the implicit viscosity step requires solving a generalized Poisson-type jump problem

$u^*$5

with prescribed interfacial Dirichlet and flux jumps (Binswanger et al., 14 Aug 2025). The inhomogeneous Dirichlet jump is removed by constructing a smooth extension $u^*$6, after which the modified system is discretized by finite volumes on nodal-centered control volumes. Volumetric integrals use the midpoint rule, uncut faces use the same finite-difference stencils as the nodal operators, cut faces are split into $u^*$7 and $u^*$8 portions, and mixed-derivative terms are computed by averaging adjacent finite-difference quotients on the face, weighted by neighboring cell areas (Binswanger et al., 14 Aug 2025). The resulting globally coupled linear system is solved by preconditioned BiCGSTAB (Binswanger et al., 14 Aug 2025).

The 2020 variational method uses a different discrete apparatus but the same sharp-geometry principle. It stores a level set of the solid boundary, polygonizes it by marching cubes, intersects each computational cell with the fluid region, and converts all volume and boundary integrals into integrals over clipped polyhedra and cut faces (Gagniere et al., 2020). Because the true geometry is carried into the variational integrals, Dirichlet and free-surface conditions are enforced “to sub-cell resolution—no ‘voxelized’ smearing” (Gagniere et al., 2020).

A plausible implication is that “sharp collocated projection” denotes not only nodal storage, but a broader commitment to embedding geometry and jump conditions directly into the discrete operators rather than regularizing them over several cells.

5. Temporal integration, iterative correction, and adaptivity

In the single-phase octree method, time advancement uses “a two-step fractional-step SLBDF2 scheme with implicit viscosity, explicit advection on characteristics” (Blomquist et al., 2023). The predictor computes $u^*$9 via semi-Lagrangian BDF2,

$\nabla \cdot u^{n+1}=0$0

with departure points found by RK2; an improved version extrapolates $\nabla \cdot u^{n+1}=0$1 from $\nabla \cdot u^{n+1}=0$2 to keep second-order at large CFL (Blomquist et al., 2023). The corrector iterates

$\nabla \cdot u^{n+1}=0$3

until

$\nabla \cdot u^{n+1}=0$4

up to $\nabla \cdot u^{n+1}=0$5, and a boundary correction then adjusts a correction $\nabla \cdot u^{n+1}=0$6 on $\nabla \cdot u^{n+1}=0$7 so that $\nabla \cdot u^{n+1}=0$8 approaches exact no-slip within $\nabla \cdot u^{n+1}=0$9 (Blomquist et al., 2023).

The two-phase method organizes the advance from $\Omega^- \subset \mathbb{R}^{2,3}$0 into four steps. It begins with a pressure-guess solve,

$\Omega^- \subset \mathbb{R}^{2,3}$1

followed by an iteration over viscosity and projection substeps until jump residuals are small (Binswanger et al., 14 Aug 2025). In the semi-implicit viscosity step, the departure values are found by RK2 backward-tracking of characteristics; if a characteristic leaves its original phase, the original-phase velocity field is used for interpolation; and if the SLBDF2 weights produce overshoot, the method falls back locally to first-order SL backward Euler at that node (Binswanger et al., 14 Aug 2025). The projection solve is then followed by updates to interfacial correction terms $\Omega^- \subset \mathbb{R}^{2,3}$2 and related quantities (Binswanger et al., 14 Aug 2025).

Interface advection in the two-phase case uses a Volume-Preserving Reference Map (VPRM),

$\Omega^- \subset \mathbb{R}^{2,3}$3

with a volume-preserving Poisson correction for $\Omega^- \subset \mathbb{R}^{2,3}$4 in a shell around $\Omega^- \subset \mathbb{R}^{2,3}$5 (Binswanger et al., 14 Aug 2025). The reported extension is that this produces $\Omega^- \subset \mathbb{R}^{2,3}$6 relative mass loss (Binswanger et al., 14 Aug 2025).

Dynamic adaptivity is integral to both octree formulations. In the single-phase case, a cell $\Omega^- \subset \mathbb{R}^{2,3}$7 is refined if either $\Omega^- \subset \mathbb{R}^{2,3}$8, to capture the interface band, or $\Omega^- \subset \mathbb{R}^{2,3}$9, to resolve vortical regions (Blomquist et al., 2023). In the two-phase case, refinement uses an interface-band criterion $$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$0, with $$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$1, together with a vorticity or velocity-gradient criterion and bounds $$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$2 (Binswanger et al., 14 Aug 2025).

6. Stability, checkerboard suppression, and numerical evidence

The single-phase nodal octree paper establishes a stability result on uniform periodic grids by relating nodal operators to their staggered counterparts through simple linear interpolations $$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$3 and $$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$4. Under this construction, $$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$5 is contracting: $$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$6 and its eigenvalues lie in $$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$7 (Blomquist et al., 2023). On adaptive octrees, a sufficient condition is given for convergence of the iterated projection $$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$8: if $$\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \mu \Delta u + f, \qquad \nabla \cdot u = 0,$$9 measures interpolation error from nodes to ghost to edges to nodes, and $\Omega = \Omega^+ \cup \Omega^-$0 is the exact staggered projection, then convergence follows when $\Omega = \Omega^+ \cup \Omega^-$1 (Blomquist et al., 2023). The paper further states that in practice $\Omega = \Omega^+ \cup \Omega^-$2 and rapid convergence is observed (Blomquist et al., 2023).

Because all variables are collocated, pressure–velocity decoupling or “checkerboarding” is an obvious concern. The formulation addresses this directly: “one might fear pressure–velocity decoupling (‘checkerboarding’). In our method the iterative projection restores coupling: by repeatedly applying $\Omega = \Omega^+ \cup \Omega^-$3 until $\Omega = \Omega^+ \cup \Omega^-$4 is tiny, one systematically suppresses spurious divergence modes. No ad-hoc stabilization (Rhie–Chow, etc.) is needed” (Blomquist et al., 2023). This is a central methodological claim rather than an incidental implementation detail.

The reported numerical validation for the single-phase method includes a projection-only supra-convergence test on a random non-graded quadtree, where both $\Omega = \Omega^+ \cup \Omega^-$5 and $\Omega = \Omega^+ \cup \Omega^-$6 errors in the divergence and in $\Omega = \Omega^+ \cup \Omega^-$7 converge at second order as the grid refines (Blomquist et al., 2023). Fully coupled Navier–Stokes tests include an analytic unsteady vortex, a lid-driven cavity at $\Omega = \Omega^+ \cup \Omega^-$8, an oscillating cylinder at $\Omega = \Omega^+ \cup \Omega^-$9, Kármán vortex street simulations around a cylinder and a sphere, and flow past a $(u,v,(w),p)$00 sculpture at $(u,v,(w),p)$01 on an 11-level octree with $(u,v,(w),p)$02 nodes in $(u,v,(w),p)$03 on 40 cores (Blomquist et al., 2023). In all cases the solver is reported stable at CFL up to $(u,v,(w),p)$04 when using the improved departure-point reconstruction, and no spurious checkerboard modes are observed (Blomquist et al., 2023).

For the two-phase extension, the reported benchmarks include projection stability under Dirichlet, Neumann, and mixed boundary conditions; an analytic vortex with observed $(u,v,(w),p)$05 convergence in $(u,v,(w),p)$06 and first-order in $(u,v,(w),p)$07; a static circular droplet with parasitic currents $(u,v,(w),p)$08; an oscillating bubble whose measured period converges to Lamb’s prediction in 3D $(u,v,(w),p)$09; rising bubbles over the Bhaga–Weber cases $(u,v,(w),p)$10; and multi-bubble interactions including bubbles rising through a solid funnel (Binswanger et al., 14 Aug 2025). The paper attributes less than $(u,v,(w),p)$11 volume loss in the rising-bubble tests to VPRM, while also noting that case $(u,v,(w),p)$12 $(u,v,(w),p)$13 exhibits tip-rupture under under-resolution but remains stable (Binswanger et al., 14 Aug 2025).

The earlier variational collocated projection method reports exact preservation of a hydrostatic pool with arbitrary cut-cell geometry, second-order convergence for its BSLQB advection in a 2D Burgers’ study, and demonstrations involving vortex shedding, dam break, smoke in a bunny, and narrow-band water, emphasizing sharp enforcement of free-surface and wall-boundary conditions without ghost currents (Gagniere et al., 2020). This suggests a continuity of concerns across the collocated-projection literature: incompressibility enforcement, geometric sharpness, and avoidance of spurious pressure–velocity artifacts.

7. Scope, advantages, and extensions

Across the octree-based papers, the advantages stated for the method are consistent. The single-phase formulation is described as second-order accurate overall, capable of dynamic grid adaptivity with arbitrary geometries, and simpler in data layout because all variables are collocated (Blomquist et al., 2023). The two-phase formulation states analogous benefits: “All variables collocated $(u,v,(w),p)$14 simpler data structures on AMR trees and unified stencils,” “Sharp treatment of interfacial jump conditions without smearing or continuum surface force oscillations,” and “Monolithic discretization of full stress tensor $(u,v,(w),p)$15 no iterative velocity-component splitting, reduced spurious currents” (Binswanger et al., 14 Aug 2025).

The two-phase paper also states that the framework is “easily extended to variable viscosity/density, Marangoni forces $(u,v,(w),p)$16, non-Newtonian rheologies, and wall contact-angle dynamics by incorporating additional jump terms” (Binswanger et al., 14 Aug 2025). Since these are presented as extensions rather than implemented cases, the appropriate interpretation is prospective rather than demonstrative.

A common misconception is that collocated layouts are intrinsically unsuitable for sharp-interface incompressible flow because they inevitably require Rhie–Chow-type stabilization or suffer from checkerboarding. The cited work does not support that blanket view. Instead, it supports a narrower statement: with the particular iterative projection, sharp FD–FV boundary treatment, and octree ghost-node construction described above, stable collocated projection can be achieved on non-graded adaptive trees without ad hoc stabilization (Blomquist et al., 2023). Another plausible implication is that the decisive issue is not collocation alone, but the compatibility of the discrete divergence, gradient, Laplacian, and boundary/interface closures.

Within that scope, the Sharp Collocated Projection Method denotes a technically specific synthesis: collocated unknown placement, adaptive quadtree or octree discretization, projection-based incompressibility enforcement, and sharp treatment of embedded boundaries or phase interfaces. In the cited literature, this synthesis is used for arbitrary-boundary incompressible Navier–Stokes flow (Blomquist et al., 2023), extended to immiscible two-phase flow with surface tension and interfacial jumps (Binswanger et al., 14 Aug 2025), and related to variational collocated projection on cut-cell regular grids (Gagniere et al., 2020).