TheP8I Interface: Patch-Based Tracking

• TheP8I Interface is a patch-based formulation that decomposes a 2D moving interface into single-valued graph segments for direct geometric representation.
• It employs CFL-stable explicit and implicit second-order time-stepping schemes with finite differences, ensuring robust and accurate interface evolution.
• The method integrates with immersed boundary/interface flow solvers and boundary-integral Stokes formulations, enabling precise curvature evaluation and reliable mass conservation.

TheP8I Interface denotes the patch-based interface tracking formulation introduced in “Interface tracking using patches” (Lindbo et al., 2013). It represents a moving interface $\Gamma(t)$ by decomposing it into segments that are single-valued graphs in local coordinate systems, while evolving each segment on a fixed one-dimensional uniform grid. The resulting construction is explicit in the sense that the interface location is directly available as coordinates, yet Eulerian in the sense that the numerical update is performed through PDE discretization on stationary grids. In the same framework, the paper derives CFL-stable and second order accurate explicit and implicit time-stepping schemes, develops a second-order treatment of boundary integral Stokes flow, and embeds the interface representation in immersed boundary and immersed interface formulations for two-phase Navier–Stokes flow (Lindbo et al., 2013).

1. Geometric representation and domain decomposition

The starting point is a moving interface in two dimensions written parametrically as

$\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$

with evolution law

$\partial_t \Gamma = u(\Gamma,t).$

The key representational step is to decompose the interface into a collection of segments,

$\Gamma = \bigcup_i \gamma_i,$

and to require that each segment $\gamma_i$ be representable as a single-valued graph in its own local coordinate system (Lindbo et al., 2013).

For each segment, the paper introduces a rigid mapping between global and local coordinates,

$$M x = T(\theta)x + b, \qquad Lx = M^{-1}x = T(-\theta)(x-b),$$

where $T(\theta)$ is a rotation matrix and $b$ is a translation vector. In local coordinates, the segment is represented as

$M\gamma_i = \{(\xi,f(\xi))\}.$

This domain decomposition is not merely geometric bookkeeping. It is the device that allows the interface to be stored explicitly as piecewise graph patches while still permitting standard finite difference operations on uniform grids local to each patch (Lindbo et al., 2013).

The method’s “explicit and yet Eulerian” characterization follows directly from this construction. It is explicit because the interface location is stored directly through the coordinates $(\xi,f(\xi))$, rather than implicitly through a level set or phase field. It is Eulerian because each patch is discretized on a fixed grid

$\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$0

so the interface evolution is computed by solving a PDE on stationary mesh points. A common misconception is that explicit interface representation must imply a purely Lagrangian marker method; the patch framework shows that direct geometry and Eulerian differencing can coexist within the same formulation (Lindbo et al., 2013).

2. Local evolution equation and moving frames

Differentiating the local graph representation in time yields the patch evolution equation. For a fixed patch mapping $\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$1, the global velocity transformed into local coordinates is written as

$\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$2

and the local interface dynamics reduce to

$\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$3

The paper identifies this as a nonlinear advection equation with source term. The nonlinearity arises because the coefficients $\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$4 and $\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$5 depend on the current interface position through the mapping back to global coordinates (Lindbo et al., 2013).

This reduction is central to the method. Instead of advecting an entire interface in global coordinates, one evolves a set of one-dimensional graph functions, one per patch. That shift makes it possible to reuse standard hyperbolic-PDE machinery on each patch while retaining a direct description of interface geometry. Tangent slope, normals, curvature, and quadrature weights can then be computed within the same one-dimensional framework.

The paper also extends the construction to time-dependent patch mappings by allowing

$\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$6

The velocity field is decomposed as

$\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$7

with $\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$8, and the frame motion is chosen through

$\Gamma(t)=\{(x(s,t),y(s,t)) : s\in[a,b)\},$9

This produces an updated patch PDE of the form

$\partial_t \Gamma = u(\Gamma,t).$0

The practical motivation is that large rigid-body translation or rotation can be absorbed into the patch frame, leaving only the deformational component to be handled by the PDE solver. The paper states that this is combined naturally with the explicit method, together with Strang splitting to preserve second-order accuracy (Lindbo et al., 2013).

3. Patchwise finite differences, explicit stepping, and communication between patches

Because each patch is defined on a fixed one-dimensional grid, the method uses standard finite differences,

$\partial_t \Gamma = u(\Gamma,t).$1

These operators are used both for the evolution equation and for geometric quantities derived from the local graph (Lindbo et al., 2013).

For the explicit update, the paper derives a second-order Lax–Wendroff-type scheme for

$\partial_t \Gamma = u(\Gamma,t).$2

The predictor step is

$\partial_t \Gamma = u(\Gamma,t).$3

followed by a second-order correction involving reevaluation of $\partial_t \Gamma = u(\Gamma,t).$4 and $\partial_t \Gamma = u(\Gamma,t).$5 at the predicted state. The method is described as CFL-stable in the standard Lax–Wendroff sense, with the time step restricted by the maximum characteristic speed on the patch. The explicit scheme is formally second-order accurate in time and space, provided that coefficient evaluation and ghost-point interpolation are performed at compatible order (Lindbo et al., 2013).

Patch coupling is handled through ghost points. When a centered stencil near a patch boundary requires values outside the local patch, those values are taken from a neighboring patch after transforming the neighboring graph into the local frame of the current patch:

$\partial_t \Gamma = u(\Gamma,t).$6

This requires neighboring patches to overlap by a finite amount in global space. If overlap is lost, or if a segment ceases to be single-valued in its local frame, the decomposition must be rebuilt. The practical explicit algorithm therefore proceeds by exchanging ghost points, transforming to global coordinates, evaluating the local coefficients, taking the predictor and corrector steps, checking patch validity, and repartitioning the interface if necessary (Lindbo et al., 2013).

4. Implicit Crank–Nicholson formulation, error structure, and observed convergence

To address stiff coupled interface-flow problems, the paper also derives a patchwise Crank–Nicholson method. For a semi-discrete equation

$\partial_t \Gamma = u(\Gamma,t).$7

the update is

$\partial_t \Gamma = u(\Gamma,t).$8

and for the patch advection equation the operators are

$\partial_t \Gamma = u(\Gamma,t).$9

Because $\Gamma = \bigcup_i \gamma_i,$0 depends on the unknown interface, the resulting system is nonlinear and must be solved iteratively (Lindbo et al., 2013).

The iterative procedure described in the paper initializes from an explicit guess, computes velocity coefficients, exchanges ghost points, evaluates the nonlinear residual

$\Gamma = \bigcup_i \gamma_i,$1

and iterates until convergence. Quasi-Newton updates such as BFGS or SR1 are mentioned, with practical convergence reported in about 2–5 iterations. The implicit method relaxes the CFL restriction, although the paper notes that it cannot take arbitrarily large time steps because patch overlap must still remain adequate for ghost-point exchange (Lindbo et al., 2013).

The truncation error is given as

$\Gamma = \bigcup_i \gamma_i,$2

where $\Gamma = \bigcup_i \gamma_i,$3 is the interpolation order used in ghost-point transfer. Using a partition-of-unity argument, the authors state that the convergence order of the reconstructed interface is the same as that of the individual segments: second order. Numerical results report second-order convergence of interface location in $\Gamma = \bigcup_i \gamma_i,$4, curvature convergence that is first order in $\Gamma = \bigcup_i \gamma_i,$5 and close to second order in $\Gamma = \bigcup_i \gamma_i,$6, and second-order mass conservation behavior for both the explicit and implicit schemes. The method has no built-in exact mass conservation, but the observed relative mass loss,

$\Gamma = \bigcup_i \gamma_i,$7

shows second-order convergence with a very small error constant. The paper emphasizes that many interface methods, especially level set methods, have much poorer mass conservation (Lindbo et al., 2013).

5. Curvature evaluation and boundary-integral Stokes formulation

A notable strength of the patch representation is that curvature can be computed directly through the graph formula

$\Gamma = \bigcup_i \gamma_i,$8

The paper states that this is discretized with centered differences and that curvature is captured without filtering or smoothing, a point presented as important for interface methods (Lindbo et al., 2013).

The same local graph structure is then used to formulate a boundary-integral method for two-dimensional Stokes flow. The velocity at a point $\Gamma = \bigcup_i \gamma_i,$9 is written as

$\gamma_i$0

with the 2D Stokeslet

$\gamma_i$1

For equal viscosities inside and outside, the density is

$\gamma_i$2

and for surface tension force

$\gamma_i$3

the nondimensional form becomes

$\gamma_i$4

The interface integral is decomposed patchwise as

$\gamma_i$5

Because the geometric quantities are already available in local coordinates, the boundary integral can be expressed entirely patch by patch, and the paper reports second-order accuracy from trapezoidal quadrature with one unequal endpoint interval where necessary (Lindbo et al., 2013).

Singular integrals are treated explicitly when the evaluation point lies on the interface. The logarithmic term is split according to

$\gamma_i$6

with

$\gamma_i$7

and the last term is integrated analytically. The paper further gives the limits required to regularize the Stokeslet rational terms, such as

$\gamma_i$8

which together provide a second-order-consistent singular quadrature treatment (Lindbo et al., 2013).

6. Coupling to immersed-boundary and immersed-interface flow solvers, and scope of the method

The patch framework is embedded in an immersed boundary (IB) discretization for two-phase Navier–Stokes flow governed by

$\gamma_i$9

The flow solver is specified as a standard staggered-grid Chorin/pressure-correction scheme with explicit convection and forcing, implicit diffusion, and pressure projection. The interface contributes the singular force

$$M x = T(\theta)x + b, \qquad Lx = M^{-1}x = T(-\theta)(x-b),$$0

and the interface velocity is interpolated from the flow field through

$$M x = T(\theta)x + b, \qquad Lx = M^{-1}x = T(-\theta)(x-b),$$1

For both spreading and interpolation, the paper uses a cubic regularized delta function $$M x = T(\theta)x + b, \qquad Lx = M^{-1}x = T(-\theta)(x-b),$$2 for $$M x = T(\theta)x + b, \qquad Lx = M^{-1}x = T(-\theta)(x-b),$$3, with the piecewise definition of $$M x = T(\theta)x + b, \qquad Lx = M^{-1}x = T(-\theta)(x-b),$$4 given explicitly in the paper (Lindbo et al., 2013).

The immersed interface method (IIM) is also considered. Instead of regularizing the singular force, it imposes jump corrections such as the pressure jump

$$M x = T(\theta)x + b, \qquad Lx = M^{-1}x = T(-\theta)(x-b),$$5

and for surface tension

$$M x = T(\theta)x + b, \qquad Lx = M^{-1}x = T(-\theta)(x-b),$$6

The flow solve is modified by correction terms,

$$M x = T(\theta)x + b, \qquad Lx = M^{-1}x = T(-\theta)(x-b),$$7

On the MAC grid, the correction terms are placed on grid lines cut by the interface. The paper notes that the IIM is generally second-order accurate in space, whereas the IBM with regularized delta functions is typically only first order for nonsmooth interface problems (Lindbo et al., 2013).

The reported Navier–Stokes experiments indicate that both IBM and IIM, when coupled to the patch interface method, give overall first-order convergence because the flow-solver time discretization is first order. The IIM can be more unstable and may require smaller time steps than the IBM in the tests presented. In the Stokes limit $$M x = T(\theta)x + b, \qquad Lx = M^{-1}x = T(-\theta)(x-b),$$8, the Navier–Stokes/IBM results approach those from the boundary-integral Stokes formulation, which the paper interprets as a consistency check (Lindbo et al., 2013).

In summary, the method’s stated advantages are explicit interface geometry, Eulerian finite-difference numerics on uniform grids, second-order accurate advection schemes, direct curvature computation, natural patchwise boundary integrals, and straightforward integration into IB, IIM, and Stokes-flow frameworks. The stated limitations are that the formulation is currently developed primarily in 2D; that a 3D extension is conceptually natural but technically difficult because closed surfaces are harder to patch robustly; that patch overlap must be maintained; that large deformations can invalidate local single-valuedness and force repartitioning; and that implicit coupling can remain expensive and is not necessarily more efficient in practice. A plausible implication is that the method is best understood not as a replacement for all interface-capturing strategies, but as a specialized framework for problems where direct access to geometry, curvature, and patchwise quadrature is central (Lindbo et al., 2013).