Characteristic Bending (CB) Method

• Characteristic Bending (CB) method is an advanced framework for advecting scalars and interfaces in incompressible flows by constructing a reference map with backward characteristic tracing.
• It employs a second-order RK2 integration and quadratic cell-wise interpolation combined with a minimal L2 volume-preserving projection to systematically correct compressibility errors.
• Numerical benchmarks show that CB maintains second-order mass conservation and interface accuracy even when standard methods like SL and VPRM struggle with divergence challenges.

The Characteristic Bending (CB) method is a framework for advecting scalar and interface quantities under incompressible velocity fields, designed to systematically enforce incompressibility in numerical advection schemes. Building on semi-Lagrangian (SL) methods, CB interprets backward-in-time characteristic tracing as the construction of a reference map—an explicit diffeomorphism between the current and initial geometries of the advected domain. The core of the method is a volume-preserving projection on this map, which corrects the spurious compressibility introduced by time integration, interpolation, and discretization errors, ensuring mass and geometric feature preservation even in highly advective incompressible flows (Blomquist et al., 6 Dec 2025).

1. Mathematical Foundations

The CB method considers scalar transport under an incompressible velocity field $u(x,t)$ with $\nabla\cdot u=0$, governed by

$$\frac{\partial\phi}{\partial t} + u\cdot\nabla\phi = 0,\quad \phi(x,0) = \phi_0(x).$$

Advection is formulated through Lagrangian trajectories $x(t)$, for which $\phi(x(t),t)=\phi_0(x_0)$. The backward-in-time construction yields a reference map $\xi(x^{n+1},t^{n+1}) = x_d^n$, where

$$x_d^n = x^{n+1} - \Delta t\,u(\hat{x}, t^n+\tfrac{1}{2}\Delta t) + \mathcal{O}(\Delta t^3),$$

with $\hat{x}$ as the intermediate point for RK2 integration. The single-step reference map $$\hat{\phi}(x,t+\Delta t) = x - \Delta t\,v(\hat{\phi}(x,t), t)$$ is computed for each grid node.

Despite the divergence-free field, $\hat{\phi}$ typically violates the volume-preservation constraint $\det\nabla\phi=1$, accumulating compressibility errors $E = \Delta t^p + \Delta x^q + \epsilon$ from integration, interpolation (orders $p,q$), and inexact divergence ($\epsilon$). The remedy is to apply a minimal $L^2$ volume-preserving projection. Introducing an adjoint variable $\lambda(x)$ via the Poisson equation,

$$-\Delta\lambda = 1 - \det\nabla\hat{\phi}, \quad \lambda|_{\partial\Omega}=0,$$

the corrected map is generated by the near-identity correction $\gamma^{-1}(x)=x-\nabla\lambda(x)$ such that the projected map $\phi(x) = \hat{\phi}(\gamma^{-1}(x))$ strictly satisfies $\det\nabla\phi=1$ to leading order.

2. Algorithmic Structure and Implementation

The single-step CB workflow, and its integration into long-term reference map formulations (RMCB), consist of:

• Advection (unprojected map $\hat{\phi}$): At each Eulerian grid node $x^{n+1}$, compute the RK2 backward departure point $x_d^n$ using velocity snapshots $u^n$, $u^{n-1}$. Store $\hat{\phi}(x^{n+1}) = x_d^n$.
• Volume-preserving projection: Assemble the right-hand side $r(x) = 1 - \det\nabla\hat{\phi}(x)$, solve $-\Delta\lambda = r$ with Dirichlet boundary conditions, compute $\gamma^{-1}(x) = x - \nabla\lambda(x)$, then evaluate the projected map $\phi(x) = \hat{\phi}(\gamma^{-1}(x))$.
• Scalar reconstruction: Update the advected scalar by sampling $\phi^{n+1}(x) = \phi^n(\phi(x))$.

The following pseudocode captures the main steps:

for n=0…N−1: for each grid-node xi: # Advection step (φ̂) xi_d = RK2_departure(xi, u^n, u^{n−1}) φ̂[i] = xi_d # Projection step r[i] = 1 − det(grad φ̂)[i] λ = solve_Poisson(−Δλ = r, λ|_{∂Ω}=0) γinv[i] = xi − grad λ[i] for each node i: φ_proj[i] = φ̂(γinv[i]) # Reconstruction for each node i: φ^{n+1}[i] = interp(φ^n, φ_proj[i]) endfor

Data structures are based on non-graded quad/octrees with ghost-node support via third-order interpolation. Characteristic tracing employs two-level RK2 backward integration, and quadratic cell-wise interpolation (quadENO) with WENO-limited second derivatives is used for spatial interpolation.

3. Error Sources, Projection Properties, and Convergence

Spurious compressibility in $\hat{\phi}$ arises from:

• Time integration error: $\mathcal{O}(\Delta t^p)$ with RK2 ($p=2$).
• Interpolation error: $\mathcal{O}(\Delta x^q)$, typically $q=2$ for quadratic interpolation.
• Inexact divergence: $\|\nabla\cdot u\| = \mathcal{O}(\epsilon)$ from approximate projection.

Volume error scales as $\mathcal{O}(E)$ before projection. The linearized projection mapping reduces volume error to $\mathcal{O}(E^2)$ in one step, with

$$\begin{align*} \|\hat{\phi} - \mathrm{Id}\|_{L^\infty} &= \mathcal{O}(\Delta t^{p} + \Delta x^{q})\ \|1-\det\nabla\hat{\phi}\|_{L^\infty} &= \mathcal{O}(\Delta t^{p} + \Delta x^{q} + \epsilon)\ \|\phi - \hat{\phi}\|_{L^\infty} &= \mathcal{O}(\Delta t^{p} + \Delta x^{q} + \epsilon)\ \|1-\det\nabla\phi\|_{L^\infty} &= \mathcal{O}\big((\Delta t^{p} + \Delta x^{q} + \epsilon)^2\big) \end{align*}$$

provided smoothness of $u$ and $\phi$, and $\epsilon = \mathcal{O}(\Delta t^p + \Delta x^q)$.

4. Numerical Benchmarks and Performance Evaluation

Numerical benchmarks evaluate CB against standard SL, CLSRM, and VPRM methods:

Benchmark	Conventional SL/CLSRM/VPRM	CB/RMCB Outcome
2D Gaussian rotation	SL, CB: both 2nd-order $L^\infty$ for $\epsilon=0$. SL drops to 1st-order mass-conservation under artificial divergence.	CB retains 2nd-order mass-conservation even for $\mathcal{O}(\Delta x)$ divergence.
Slotted cylinder	SL: 1st-order interface, 1–2nd-order volume loss. CLSRM accurate for $\epsilon=0$; degrades otherwise. VPRM reduces mass loss but perturbs interface.	CB: Matches SL interface accuracy for $\epsilon=0$; maintains volume and interface accuracy for $\epsilon>0$.
Single vortex deformation	SL/CLSRM/VPRM degrade for $\epsilon>0$;	CB/RMCB invariant across divergence errors (2D/3D).
Incompressible Euler	SL: velocity error 1–2nd-order; divergence stalls.	CB: 2nd-order velocity, nearly vanishing divergence.
Analytic vortex, 2-phase	SL+VPRM vs CB+RMCB: CB+RMCB lower velocity/interface error, slight volume loss improvement.	As described.
Rising bubbles (Navier–Stokes)	VPRM distorts shape; SL/CB similar otherwise.	CB delays pinch-off; RMCB preserves mass/shape longer.

In all cases, the volume-preserving correction of CB leads to improved or preserved mass conservation and interface geometry, especially in the presence of non-negligible velocity divergence.

5. Practical Implementation Guidelines

Key recommendations for robust CB deployment:

• Time integration: RK2 backward characteristics for second-order accuracy; SETTLS optional for strong advection.
• Spatial interpolation: Quadratic cell-wise with WENO-limited derivatives (quadENO).
• Stability: CFL numbers up to 5 observed stable.
• Grid: Non-graded quad/octree, ghost nodes via third-order interpolation.
• Projection solver: Supra-convergent Poisson solver (Min et al.), tolerance $\sim 10^{-10}$ for second-order accuracy in $\lambda$.
• Interface reinitialization: PDE-based method (Sussman '94); for RM-based only at restarts, shell-width $=8\Delta x_{\min}$.
• Map restart: Trigger when columns of $\det(\nabla\xi)$ are nearly collinear (threshold $\sim 10^{-3}$).

6. Limitations and Known Failure Modes

Identified constraints and potential failure scenarios include:

• Projection cost: Global elliptic (Poisson) solve each time-step forms the principal computational overhead.
• Strong compressibility: If compressible errors exceed step/grid scales ($\epsilon\gg\Delta x,\Delta t$), the near-identity projection assumption fails, reducing correction effectiveness.
• Interface perturbation: For exactly divergence-free $u$, projection can slightly perturb interfaces (as in VPRM), but this effect is mitigated in CB by the stepwise, near-identity assumption.
• Highly nonlinear/discontinuous flows: Interpolation limiter adjustments may be necessary.

A plausible implication is that CB is most robust and cost-effective when compressibility errors remain modest and when grid and time step sizes are chosen commensurate with the method's accuracy assumptions (Blomquist et al., 6 Dec 2025).