---
title: Sharp Collocated Projection Method
url: https://www.emergentmind.com/topics/sharp-collocated-projection-method
type: topic
---

# Sharp Collocated Projection Method

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 [2306.09957]. 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 [2508.11107]. 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 [2003.12227].

## 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 [2306.09957]. 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 [2306.09957].

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\) [2306.09957].

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 [2003.12227]. 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 [2306.09957].

For the immiscible two-phase extension, two incompressible Newtonian fluids occupy \(\Omega = \Omega^+ \cup \Omega^-\), separated by a sharp interface \(\Gamma\). In each phase the density and viscosity are constant, with \((\rho^+,\mu^+)\) in \(\Omega^+\) and \((\rho^-,\mu^-)\) in \(\Omega^-\). The strong form is
\[
\nabla \cdot u = 0 \quad \text{in } \Omega^+ \cup \Omega^-,
\]
\[
\rho \left(\frac{\partial u}{\partial t} + u \cdot \nabla u\right) = -\nabla p + \nabla \cdot \tau + f_{\text{body}},
\quad
\tau = \mu (\nabla u + \nabla u^T),
\]
with surface tension represented as
\[
f_{\text{st}}(x) = \sigma \kappa n \delta_\Gamma(x),
\]
where \(\sigma\) is the surface tension coefficient, \(\kappa = \nabla \cdot n\), \(n = \nabla \phi / |\nabla \phi|\), and \(\delta_\Gamma\) is the Dirac delta supported on \(\Gamma\) [2508.11107].

The canonical jump conditions across \(\Gamma\) are
\[
[u] \equiv u^+ - u^- = 0 \quad \text{on } \Gamma,
\]
\[
[(-pI+\tau)\cdot n] = \sigma \kappa n,
\]
equivalently
\[
[p] = \sigma \kappa,
\qquad
[\mu(\nabla u + \nabla u^T)\cdot n]\cdot t = 0
\]
for the tangential shear component [2508.11107]. The interface is represented implicitly by a level-set function \(\phi(x,t)\), reinitialized as a signed distance so that \(|\nabla \phi| = 1\) [2508.11107].

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 [2508.11107]. At T-junctions, where a node lacks a direct neighbor, ghost-node values are created by a third-order interpolation formula of the form
\[
\phi_r =
\frac{r_b \phi_{r_t} + r_t \phi_{r_b}}{r_t+r_b}
-
\frac{r_t r_b}{t+b}
\left[
\frac{\phi_t-\phi_0}{t} - \frac{\phi_0-\phi_b}{b}
\right],
\]
which is used in both the nodal octree method and the two-phase extension [2306.09957, 2508.11107].

## 3. Discrete differential operators and the projection mechanism

The collocated octree method defines nodal analogues of the Laplacian, divergence, and gradient. At node \(0\), with neighbors left \(\ell\), right \(r\), bottom \(b\), and top \(t\), the nodal Laplacian is
\[
L_N \phi|_0
=
\frac{2}{r+\ell}
\left[
\frac{\phi_r-\phi_0}{r} - \frac{\phi_0-\phi_\ell}{\ell}
\right]
+
\frac{2}{t+b}
\left[
\frac{\phi_t-\phi_0}{t} - \frac{\phi_0-\phi_b}{b}
\right].
\]
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 \(\phi\) and \(\nabla \phi\), described as supra-convergence [2306.09957].

The divergence operator is
\[
D_N(u,v)|_0 = \frac{u_r-u_\ell}{r+\ell} + \frac{v_t-v_b}{t+b},
\]
and is stated to be formally only first-order accurate on adaptive meshes, chosen for stability [2306.09957]. The gradient operator is
\[
G_N \phi|_0 =
\left(
\frac{\ell}{r+\ell}\cdot\frac{\phi_r-\phi_0}{r}
+
\frac{r}{r+\ell}\cdot\frac{\phi_0-\phi_\ell}{\ell},
\;
\frac{b}{t+b}\cdot\frac{\phi_t-\phi_0}{t}
+
\frac{t}{t+b}\cdot\frac{\phi_0-\phi_b}{b}
\right).
\]
The same nodal stencil family is reused in the two-phase method as part of its sharp hybrid FD–FV discretization [2508.11107].

The central projection step is the discrete Hodge projection. In the single-phase method, the scalar Hodge variable \(\phi\) is obtained from
\[
D_N[G_N \phi] = D_N u^*,
\]
after which the updated field is obtained by repeated application of \(P_N\) [2306.09957]. In the two-phase method, the projection step becomes a variable-density Poisson problem:
\[
\nabla \cdot \left(\frac{1}{\rho}\nabla \phi\right) = \nabla \cdot u^*
\quad \text{in } \Omega,
\]
with
\[
[\phi]=0,
\qquad
\left[\left(\frac{1}{\rho}\nabla \phi\right)\cdot n\right]=0
\quad \text{on } \Gamma,
\]
and Dirichlet or Neumann conditions on \(\partial \Omega\) [2508.11107].

No extra pressure-velocity correction terms are needed beyond the iterative projection in the nodal octree formulation [2306.09957]. 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 \(\phi(x)\) defines \(\Omega^-\) versus solids \(\Omega^+\). On cut cells the hybrid finite-volume/finite-difference discretization “imposes no-slip \(u=U_{\text{solid}}\) exactly on \(\Gamma\)” and homogeneous Neumann \(\partial \phi/\partial n = 0\) for the Hodge variable [2306.09957]. Curved and moving interfaces are represented by a level set on the octree, and cut-cells near \(\Gamma\) are handled sharply via hybrid FV/FD stencils that respect the exact geometry and impose \(u=U_{\text{solid}}\) on \(\Gamma\) without smearing [2306.09957].

In the two-phase extension, the implicit viscosity step requires solving a generalized Poisson-type jump problem
\[
\eta v - \nabla \cdot [\mu(\nabla v + \nabla v^T)] = r
\quad \text{in } \Omega^+ \cup \Omega^-,
\]
with prescribed interfacial Dirichlet and flux jumps [2508.11107]. The inhomogeneous Dirichlet jump is removed by constructing a smooth extension \(\psi_{\text{ext}}\), 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 \(\Omega^+\) and \(\Omega^-\) portions, and mixed-derivative terms are computed by averaging adjacent finite-difference quotients on the face, weighted by neighboring cell areas [2508.11107]. The resulting globally coupled linear system is solved by preconditioned BiCGSTAB [2508.11107].

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 [2003.12227]. Because the true geometry is carried into the variational integrals, Dirichlet and free-surface conditions are enforced “to sub-cell resolution—no ‘voxelized’ smearing” [2003.12227].

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” [2306.09957]. The predictor computes \(u^*\) via semi-Lagrangian BDF2,
\[
\rho \left[\alpha\frac{u^*-u_d^n}{\Delta t} + \beta\frac{u_d^n-u_d^{n-1}}{\Delta t_{n-1}}\right]
=
\mu \Delta u^* + f,
\]
with departure points found by RK2; an improved version extrapolates \(u^{n+1}\) from \(u^n, u^{n-1}\) to keep second-order at large CFL [2306.09957]. The corrector iterates
\[
u^{(m+1)} = P_N u^{(m)}
\]
until
\[
\|u^{(m+1)}-u^{(m)}\| \le \epsilon_i \|u^{(m+1)}\|,
\]
up to \(K_{\max} \approx 5\), and a boundary correction then adjusts a correction \(c\) on \(\Gamma\) so that \(u^{n+1}|_\Gamma\) approaches exact no-slip within \(O(\Delta t)\) [2306.09957].

The two-phase method organizes the advance from \(t^n \to t^{n+1}\) into four steps. It begins with a pressure-guess solve,
\[
\nabla^2 p_g = 0 \quad \text{in } \Omega,
\qquad
[p_g] = -\sigma \kappa,
\qquad
\left[\left(\frac{1}{\rho}\nabla p_g\right)\cdot n\right]=0 \text{ on } \Gamma,
\]
followed by an iteration over viscosity and projection substeps until jump residuals are small [2508.11107]. 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 [2508.11107]. The projection solve is then followed by updates to interfacial correction terms \(X_k\) and related quantities [2508.11107].

Interface advection in the two-phase case uses a Volume-Preserving Reference Map (VPRM),
\[
\frac{\partial \xi}{\partial t} + u \cdot \nabla \xi = 0,
\qquad
\phi(x,t)=\phi_0(\xi(x,t)),
\]
with a volume-preserving Poisson correction for \(\xi\) in a shell around \(\Gamma\) [2508.11107]. The reported extension is that this produces \(O(10^{-3}\!-\!10^{-4})\) relative mass loss [2508.11107].

Dynamic adaptivity is integral to both octree formulations. In the single-phase case, a cell \(C\) is refined if either \(|\phi(\text{node})| \le B \cdot \operatorname{Lip}(\phi)\cdot \operatorname{diag}(C)\), to capture the interface band, or \(\operatorname{diag}(C)\cdot |\nabla \times u|/|u_\infty| \ge T_v\), to resolve vortical regions [2306.09957]. In the two-phase case, refinement uses an interface-band criterion \(\min_{n \in C} |\phi(n)| \le B \cdot \Delta x \cdot \operatorname{Lip}(\phi)\), with \(\operatorname{Lip}(\phi)\approx 1\), together with a vorticity or velocity-gradient criterion and bounds \(\ell_{\min} \le \operatorname{level}(C) \le \ell_{\max}\) [2508.11107].

## 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 \(I_E^N\) and \(I_N^E\). Under this construction, \(P_N\) is contracting:
\[
\|P_N X\|_2 \le \|X\|_2,
\]
and its eigenvalues lie in \([0,1]\) [2306.09957]. On adaptive octrees, a sufficient condition is given for convergence of the iterated projection \(P_N^k \to P_\infty\): if \(\epsilon = I-I_p\) measures interpolation error from nodes to ghost to edges to nodes, and \(Q\) is the exact staggered projection, then convergence follows when \(\rho(\epsilon Q) < 1\) [2306.09957]. The paper further states that in practice \(\epsilon \approx O(1)\) and rapid convergence is observed [2306.09957].

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 \(u \leftarrow P_N u\) until \(\|u-P_Nu\|\) is tiny, one systematically suppresses spurious divergence modes. No ad-hoc stabilization (Rhie–Chow, etc.) is needed” [2306.09957]. 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 \(L^1\) and \(L^\infty\) errors in the divergence and in \(u^{n+1}\) converge at second order as the grid refines [2306.09957]. Fully coupled Navier–Stokes tests include an analytic unsteady vortex, a lid-driven cavity at \(Re=1000\), an oscillating cylinder at \(Re=100, KC=5\), Kármán vortex street simulations around a cylinder and a sphere, and flow past a \(12\,\mathrm{m}\) sculpture at \(Re \approx 1.9 \times 10^7\) on an 11-level octree with \(O(10^6)\) nodes in \(O(10^2)\,\mathrm{sec/step}\) on 40 cores [2306.09957]. In all cases the solver is reported stable at CFL up to \(O(10)\) when using the improved departure-point reconstruction, and no spurious checkerboard modes are observed [2306.09957].

For the two-phase extension, the reported benchmarks include projection stability under Dirichlet, Neumann, and mixed boundary conditions; an analytic vortex with observed \(1\!-\!2^\circ\) convergence in \(u\) and first-order in \(\phi\); a static circular droplet with parasitic currents \(\|u\|_\infty \sim O(\Delta x)\); an oscillating bubble whose measured period converges to Lamb’s prediction in 3D \((T \approx 3.63)\); rising bubbles over the Bhaga–Weber cases \(a\!-\!h\); and multi-bubble interactions including bubbles rising through a solid funnel [2508.11107]. The paper attributes less than \(0.2\%\) volume loss in the rising-bubble tests to VPRM, while also noting that case \(f\) \((Re \sim 259)\) exhibits tip-rupture under under-resolution but remains stable [2508.11107].

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 [2003.12227]. 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 [2306.09957]. The two-phase formulation states analogous benefits: “All variables collocated \(\to\) 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 \(\to\) no iterative velocity-component splitting, reduced spurious currents” [2508.11107].

The two-phase paper also states that the framework is “easily extended to variable viscosity/density, Marangoni forces \((\nabla \sigma)\), non-Newtonian rheologies, and wall contact-angle dynamics by incorporating additional jump terms” [2508.11107]. 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 [2306.09957]. 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 [2306.09957], extended to immiscible two-phase flow with surface tension and interfacial jumps [2508.11107], and related to variational collocated projection on cut-cell regular grids [2003.12227].

Source: https://www.emergentmind.com/topics/sharp-collocated-projection-method