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. This post summarizes the cost of calculating this value naively, the idea of transforming the problem into the Eikonal equation (a hyperbolic PDE of the form |∇d|=1/F), and how to obtain the solution in $O(N \log N)$ using Fast Marching, the continuous version of Dijkstra's algorithm.

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.

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.

Visualizing Wavefront Propagation#

Try interacting with the simulation below. Start by clicking "Single Cylinder" and try increasing the speed $F$ or changing the obstacles.

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.

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)$.

Key Takeaways#

• Wall distance can be viewed as the solution to the Eikonal equation $|\nabla d|=1/F$, opening the door to $O(N \log N)$ algorithms.
• Fast Marching manages the Considered set using a min-heap and finalizes the wavefront by solving upwind quadratic equations point by point.
• Choose between Fast Marching or P-Poisson relaxation based on moving mesh, periodic boundary, and parallelization requirements.

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.