Solving Wall Distance in O(N log N): Eikonal Equation and Fast Marching

In 1996, James Sethian of UC Berkeley was refining algorithms to track vessel boundaries in medical imaging. The resulting Fast Marching Method is now used in unexpected places. Every time a RANS solver runs around an aircraft wing, this very algorithm generates the wall distance for each grid point. The Spalart-Allmaras and k-ω SST models require the value of "where is the nearest wall" as an input for every cell.

Start on the single-cylinder pattern. Raising the speed $F$ sends the yellow band farther in the same instant. Switching the obstacle to a slit changes which wall the band leaves from.

The yellow band is the Considered set being updated, which is effectively the "wavefront" itself. Fronts start simultaneously from the wall and cylinder boundaries, and distances are finalized where they meet. In the two-cylinder pattern, the collision line of the two fronts forms near the medial axis of the cylinders. The slit wall pattern shows the wavefront rounding the gap to reach the back side. This "rounding" is where it diverges from naive Euclidean distance. In real aircraft or turbine geometries, turbulence model inputs for regions behind walls require this rounded distance.

Why Turbulence Models Need Wall Distance#

Statistics of boundary layer turbulence depend strongly on the distance from the wall, $y$. RANS (Reynolds-Averaged Navier-Stokes) models incorporate this dependency via algebraic expressions. Spalart-Allmaras uses $d_w$ (wall distance) directly in its destruction term. In k-ω SST, $d_w^2$ appears within the blending functions $F_1$ and $F_2$. Wall distance also acts as a switch in hybrid LES models (DES, DDES, IDDES).

The problem is that $d_w$ is a geometric property of the mesh. If the mesh moves or deforms, it must be recalculated at every time step. This is the case for rotating machinery (compressors, pumps) where the rotor-stator boundary involves relative motion.

The Pitfalls of the Naive Approach#

The most immediate method is to "calculate the distance between every wall face and every cell and pick the minimum." This is actually used in textbook examples.

The cost is $O(N_{\text{cell}} \cdot N_{\text{wall}})$. Imagine solving for an entire aircraft. With 10M cells and 0.5M wall faces, that would be 5 trillion distance calculations. If this runs every time step, the simulation will never finish. While spatial data structures like KD-trees can lower this to $O(N_{\text{cell}} \log N_{\text{wall}})$, near-wall queries for complex geometries remain expensive. The same all-pairs habit shows up in particle contact. A DEM sector sized from a single particle cuts that list. The objects differ. The reason to clip candidates by geometry does not.

Changing Perspective: The Eikonal Equation#

Let's view the wall distance as the solution to a PDE. A function that is "zero at the wall and has a gradient magnitude of 1 in all directions" is precisely the distance field.

Here, $d$ is the "arrival time" to the wall, $F$ is the propagation speed (defaulting to 1), and $\Gamma_w$ is the set of wall boundaries. If $F=1$, $d$ is simply the Euclidean distance.

This equation is hyperbolic. Information flows in only one direction along characteristic lines (straight lines extending outward from the wall). This unidirectionality is the key to an efficient algorithm. By exploiting the property that "values only enter from points where the value is already determined," visiting each point once is sufficient.

Fast Marching: The Continuous Dijkstra#

Sethian's idea can be summarized in two steps:

1. Divide all grid points into three sets: Accepted, Considered, and Far.
2. From the Considered set, identify the point with the smallest arrival time, move it to Accepted, and update its neighbors.

To quickly extract the smallest value, a min-heap is natural. The structure is identical to Dijkstra's shortest path algorithm on graphs. The difference lies in how neighbors are updated. In a graph, it is simple addition; in Eikonal, we solve an upwind quadratic equation.

The update formula on a 2D structured grid is as follows:

Where $T_x$ and $T_y$ are the upwind neighbor values in the horizontal and vertical directions (the smaller of the two neighbors), and $h$ is the grid spacing. Among the two roots, only the one satisfying $T \ge \max(T_x, T_y)$ is accepted. This condition is called causality. It enforces the Eikonal unidirectionality, ensuring that information only comes from "already determined, smaller" neighbors.

If information comes from only one direction (e.g., if other neighbors are Far), the equation simplifies:

Since min-heap insertion and deletion are $O(\log N)$ and each point is processed a constant number of times on average, the total cost is $O(N \log N)$.

Python Implementation#

import heapq
import numpy as np
 
INF = np.inf
FAR, CONSIDERED, ACCEPTED = 0, 1, 2
 
 
def solve_eikonal_neighbor(T, j, i, h, F):
    """Obtains the quadratic upwind solution at a single point."""
    ny, nx = T.shape
    Tx = min(
        T[j, i - 1] if i > 0 else INF,
        T[j, i + 1] if i < nx - 1 else INF,
    )
    Ty = min(
        T[j - 1, i] if j > 0 else INF,
        T[j + 1, i] if j < ny - 1 else INF,
    )
    a, b = sorted([Tx, Ty])
    inv_f = h / F
    if not np.isfinite(b) or b >= a + inv_f:
        return a + inv_f
    diff = b - a
    disc = 2 * inv_f * inv_f - diff * diff
    if disc < 0:
        return a + inv_f
    return 0.5 * (a + b + np.sqrt(disc))
 
 
def march_wall_distance(walls, h=1.0, F=1.0):
    """Calculates wall distance using Fast Marching on a 2D structured grid."""
    ny, nx = walls.shape
    T = np.full_like(walls, INF, dtype=np.float64)
    status = np.zeros_like(walls, dtype=np.uint8)
    heap = []
 
    T[walls] = 0.0
    status[walls] = ACCEPTED
 
    # Initialize neighbors of walls as Considered
    for j in range(ny):
        for i in range(nx):
            if not walls[j, i]:
                continue
            for dj, di in [(-1, 0), (1, 0), (0, -1), (0, 1)]:
                nj, ni = j + dj, i + di
                if 0 <= nj < ny and 0 <= ni < nx and status[nj, ni] != ACCEPTED:
                    new_t = solve_eikonal_neighbor(T, nj, ni, h, F)
                    if new_t < T[nj, ni]:
                        T[nj, ni] = new_t
                        status[nj, ni] = CONSIDERED
                        heapq.heappush(heap, (new_t, nj, ni))
 
    while heap:
        t, j, i = heapq.heappop(heap)
        if status[j, i] == ACCEPTED:
            continue  # Old entry
        status[j, i] = ACCEPTED
        for dj, di in [(-1, 0), (1, 0), (0, -1), (0, 1)]:
            nj, ni = j + dj, i + di
            if not (0 <= nj < ny and 0 <= ni < nx):
                continue
            if status[nj, ni] == ACCEPTED:
                continue
            new_t = solve_eikonal_neighbor(T, nj, ni, h, F)
            if new_t < T[nj, ni]:
                T[nj, ni] = new_t
                status[nj, ni] = CONSIDERED
                heapq.heappush(heap, (new_t, nj, ni))
    return T
 
 
# Verification: Flat channel
walls = np.zeros((61, 121), dtype=bool)
walls[0, :] = True
walls[-1, :] = True
 
d = march_wall_distance(walls, h=1.0)
 
# The wall distance at the centerline (y=30) should be 30
print(f"Centerline wall distance = {d[30, 60]:.4f} (Exact: 30.0)")
print(f"Max relative error = {np.max(np.abs(d[30, :] - 30)):.4f}")

In a flat channel, the analytical solution $d=\min(y, H-y)$ emerges cleanly. The centerline output of the code above matches 30.0 almost exactly. While the 2D first-order upwind approximation has slight errors near corners, it is sufficient for practical applications.

One accepted neighbor adds 1 at every cell#

The channel check uses a 61 by 121 grid, spacing $h=1$, speed $F=1$. The bottom wall is index 0, the top wall is 60, the centerline is 30. If only one neighbor is Accepted, the update collapses to $T = T_x + h/F$. With $h/F=1$, each cell away from the wall adds 1 to the arrival time. Thirty cells give 30. Coming down from the top wall also gives 30. The minimum of the two fronts is the centerline value 30.0. That is the exact number written next to the print. The corner residual is where the discriminant goes negative and the update falls back to the one-sided branch.

Ten million cells times half a million walls is five trillion. Fast Marching touches the heap at $O(\log N)$ per point and promotes each point to Accepted once. At $N=10^{7}$, $\log N$ is 7 in base 10 and about 24 in base 2. The operation count has a hard time clearing a few hundred million. That is four orders under five trillion. A KD-tree still multiplies the cell count by a log of $N_{\text{wall}}$. The wall query remains heavier than the march.

If the neighbor gap $b-a$ is large enough that $2(h/F)^2-(b-a)^2$ is negative, the root returns to the smaller neighbor plus $h/F$. In the channel one side is still Far, so that branch is chosen from the start. On a flat wall the distance field is the cell count itself.

Extension to Unstructured Meshes#

While formulas for structured grids are elegant, half of CFD happens on unstructured meshes. Sethian & Vladimirsky (2000) solved this. The value at a point $\mathbf{x}_0$ is obtained from the values at other vertices of the simplex (triangle, tetrahedron) containing that point.

By solving the normal equations from a matrix $P$ of direction vectors $\mathbf{p}_r = \mathbf{x}_r - \mathbf{x}_0$ and directional derivatives $v_r = (T_0 - T_r)/|\mathbf{p}_r|$, the Eikonal equation $|\nabla T|^2 = 1/F^2$ becomes a quadratic equation for $T_0$. Only the root where the approximate gradient points inside the simplex is accepted. This simplex causality generalizes the simple upwind condition of structured grids.

Practical Checklist#

• Moving Mesh: Instead of full recalculation, lazy updates for regions with significant deformation are useful. Fast Marching is sequential and thus difficult to parallelize, but narrow-band variants are easier to implement.
• Periodic BC: The graph must be connected such that values propagate from the periodic side back toward the wall. Otherwise, wall distances near periodic boundaries will be overestimated.
• Multiple Materials: In rotor-stator analysis, include both sets of moving walls in $\Gamma_w$. Regions with sharp distance gradients, like turbine tip gaps, require sufficiently fine meshes.
• P-Poisson Alternative: Solving $|\nabla d|^p = 1$ ($p=4 \sim 10$) implicitly is slower than Fast Marching but easier to parallelize. This is often preferred in large-scale parallel solvers.
• Local Negative Errors: Depending on the numerical discretization, $d_w$ may result in negative values in extremely thin boundary layers. Always clamp to $\max(d_w, 0)$.

March the whole mesh every step, or hand it to P-Poisson#

On a moving rotor mesh, running Fast Marching over every cell at every step repeats the heap work for 10M cells on every time advance. Skip the narrow band of cells that actually moved and the $O(N \log N)$ saving gets multiplied by the step count again, the wrong way.

If parallelism comes first, switch to $|\nabla d|^p = 1$, with $p$ between 4 and 10. How fast the iteration moves is set by the linear solver, from Jacobi to BiCGSTAB. Make that swap only if the march's causality can be given up.

References#

• J. A. Sethian, A. Vladimirsky, "Fast methods for the Eikonal and related Hamilton–Jacobi equations on unstructured meshes", PNAS 97(11), 5699–5703, 2000.
• P. Tucker, "Assessment of geometric multilevel convergence robustness and a wall distance method for flows with multiple internal boundaries", Applied Mathematical Modelling 22, 1998.