Comparison: This Implementation vs. Steger & Sorenson (1979)

A detailed comparison of the ellipticSS solver in GridGeneration.jl against the original Steger & Sorenson method and its production implementation in the GRAPE code (Sorenson, NASA TM-81198, 1980).

Overview

The ellipticSS implementation is a faithful simplified version of the Steger & Sorenson approach. It captures the core mathematical ideas – elliptic PDE smoothing with automatic wall-forcing for orthogonality – but omits several algorithmic accelerations that GRAPE uses for production efficiency.

1. PDE Formulation: Essentially Identical

Both solve the same transformed Poisson system:

alpha * x_{xi,xi} - 2*beta * x_{xi,eta} + gamma * x_{eta,eta} = RHS_x
alpha * y_{xi,xi} - 2*beta * y_{xi,eta} + gamma * y_{eta,eta} = RHS_y

with identical metric coefficient definitions. The discretization uses the same 9-point stencil with central differences. No difference here.

2. Orthogonality Enforcement: Same Core Idea, Simplified Execution

What Steger & Sorenson Do

The original method computes control functions P(xi,eta) and Q(xi,eta) by evaluating the full Poisson equations at the boundary with desired derivatives imposed:

alpha * x_{xi,xi} - 2*beta * x_{xi,eta} + gamma * x_{eta,eta} + J^2*(P*x_xi + Q*x_eta) = 0

At the wall, with beta_b = 0 (orthogonality enforced) and known desired derivatives, this yields explicit values of P and Q at every boundary point. These P, Q values are then decayed into the interior and appear as source terms throughout the domain, coupled through the J^2 factor.

What This Implementation Does

The code computes the boundary residual directly:

RHS_x = -(alpha_b * x_{xi,xi} + gamma_b * x_{eta,eta})

This is equivalent to computing J^2*(P*x_xi + Q*x_eta) at the boundary, but skips the explicit decomposition into P and Q. The forcing is applied as a single combined RHS term rather than separate P, Q control functions multiplied by first derivatives.

Practical Consequence

The difference is subtle but meaningful:

• Steger & Sorenson: P and Q are scalar functions. When propagated into the interior, they multiply the local first derivatives (xxi, xeta, etc.) at each interior point. This means the forcing adapts to the local grid geometry as it moves away from the wall.

• This implementation: The RHS values computed at the wall are propagated as fixed vectors that don't interact with interior metrics. The forcing at an interior point (i,j) is simply the wall value times exp(-decay * distance), regardless of what the local grid looks like at (i,j).

For grids where the interior geometry doesn't differ dramatically from the boundary geometry, this is a reasonable simplification. For highly curved or stretched grids, the Steger & Sorenson approach provides more physically consistent forcing in the interior.

3. Wall Spacing Control: Simplified

Steger & Sorenson / GRAPE

GRAPE allows the user to prescribe a full spacing distribution s(xi) along each boundary. This can be:

• Constant spacing (uniform boundary layer)
• Geometric stretching (exponential growth away from corners)
• Arbitrary user-specified distribution
• Matched to a flow solution (adaptive refinement)

This Implementation

The spacing is determined by linear interpolation between the actual spacings at the two corner points of each boundary:

s1 = ||point[1, wall] - point[1, wall+1]||
s2 = ||point[Ni, wall] - point[Ni, wall+1]||
s_vec = LinRange(s1, s2, Ni)

This means:

• The user doesn't directly control wall spacing – it's inherited from the initial grid
• The spacing varies linearly along the wall, which may not match the desired distribution
• There's no option for geometric stretching or custom distributions

Impact: For grids where the initial TFI/edge-solver already provides good boundary spacing, this works fine. For boundary layer grids requiring precise wall-normal spacing control, the GRAPE approach is more capable.

4. Linear Solver: The Biggest Algorithmic Difference

Steger & Sorenson / GRAPE: Line Gauss-Seidel

GRAPE reorganizes the discretized equations at each j-level into a tridiagonal system in the i-direction:

a_i * x_{i-1,j}^{n+1} + b_i * x_{i,j}^{n+1} + c_i * x_{i+1,j}^{n+1} = d_i

where:

• a_i = alpha[i,j]
• b_i = -2*(alpha[i,j] + gamma[i,j])
• c_i = alpha[i,j]
• d_i = -gamma[i,j](x[i,j+1] + x[i,j-1]) + beta[i,j]/2(cross terms using old values) - RHS[i,j]

This tridiagonal system is solved exactly in O(Ni) operations using the Thomas algorithm. The key advantage: all points along a j-line are updated simultaneously and implicitly coupled, leading to much faster information propagation across the grid.

After sweeping all j-levels, the same process is applied along i-levels for the eta-direction (alternating direction).

This Implementation: Point SOR

Each interior point is updated individually using only its current neighbors:

x_new = (alpha*(x[i+1,j]+x[i-1,j]) + gamma*(x[i,j+1]+x[i,j-1])
         - beta/2*cross_terms + RHS) / (2*(alpha + gamma))
x[i,j] = (1 - omega)*x[i,j] + omega*x_new

Performance Comparison

The line solver is dramatically faster because information propagates across an entire grid line in a single step, whereas point SOR can only move information one grid point per iteration.

Example: For a 100 x 100 grid:

• Line Gauss-Seidel: ~100-500 iterations
• Point SOR (omega=0.2): ~5000-50000 iterations

This is the single biggest performance gap between the implementation and GRAPE.

5. Acceleration Techniques: Not Implemented

GRAPE's Coarse-Fine Sequencing

GRAPE uses a two-level coarse-fine strategy (a simplified multigrid):

1. Extract every 2nd (or 4th) point to create a coarse grid
2. Solve the elliptic system on the coarse grid to convergence
3. Interpolate the coarse solution onto the fine grid
4. Use this as the initial guess for the fine-grid solve
5. The fine-grid solve converges much faster with a good initial guess

This can reduce total iteration count by a factor of 4-10x.

This Implementation

No multigrid or coarse-fine acceleration. The solver starts from the TFI grid directly and iterates at full resolution. This is simpler but slower for large grids.

6. Relaxation Strategy: Opposite Approaches

GRAPE: Over-Relaxation (omega > 1)

Because GRAPE uses line Gauss-Seidel (which is inherently more stable due to implicit coupling), it can afford over-relaxation with omega in the range 1.0-1.8. The optimal omega for the line method is typically around 1.4-1.6, providing significant acceleration.

GRAPE also uses variable relaxation: different omega values for different iteration phases, and sometimes different omega for x vs. y updates.

This Implementation: Under-Relaxation (omega = 0.2)

The point SOR method with recomputed nonlinear metrics requires under-relaxation for stability. The default omega = 0.2 means each update step is only 20% of the full Gauss-Seidel correction.

Why the difference?

• Line methods implicitly solve for all points on a line simultaneously, averaging out local errors. This makes them tolerant of aggressive relaxation.
• Point methods update each point in isolation using neighbors that may be a mix of old and new values. The nonlinear metric feedback can amplify errors if steps are too large.

7. Convergence Monitoring

GRAPE

GRAPE monitors multiple convergence metrics:

• Maximum absolute change in x and y coordinates (L-infinity norm)
• RMS change across the grid
• Maximum Jacobian ratio (measure of cell quality)
• Orthogonality angle at boundaries

It also uses adaptive iteration control: if convergence stalls, it can adjust omega or switch between solver phases.

This Implementation

Monitors the larger of the L2 norms of the x- and y-displacements, warns when max_iter is reached, and errors on a non-finite update:

err = max(norm(x - x_old), norm(y - y_old))

This is adequate for detecting convergence but doesn't distinguish between:

• Grid still moving significantly (not converged)
• Grid quality improving but positions are stable
• Grid quality degrading (rare, but possible with bad parameters)

y-displacement is not independently monitored.

8. Boundary Derivative Computation

Steger & Sorenson

The original paper uses second-order one-sided differences for normal derivatives at the wall, with careful treatment of the ghost-point approach to impose desired first derivatives.

This Implementation

Uses a specific modified stencil (line 86 of EllipticSolver.jl):

x_etaeta = 0.5*(-7*x[i,wall] + 8*x[i,wall+1] - x[i,wall+2]) - dir*3*x_eta_desired

This formula is consistent with the Steger & Sorenson approach. It can be derived by:

1. Writing the Taylor expansion at the wall with the desired first derivative:

   x[wall+1] = x[wall] + x_eta_desired + (1/2)*x_etaeta + ...
   x[wall+2] = x[wall] + 2*x_eta_desired + 2*x_etaeta + ...

2. Solving these two equations for x_etaeta in terms of known values.

The coefficients (-7, 8, -1, -3) arise from this derivation. This matches the original method's approach of using the desired derivative to compute a consistent second derivative.

However, let us verify the stencil. From a second-order one-sided formula for f''(0) using points at 0, h, 2h with a known f'(0):

f(h)  = f(0) + h*f'(0) + (h^2/2)*f''(0) + (h^3/6)*f'''(0) + ...
f(2h) = f(0) + 2h*f'(0) + 2h^2*f''(0) + (4h^3/3)*f'''(0) + ...

With h = 1 (unit computational spacing):

f(1) = f(0) + f'_desired + (1/2)*f''
f(2) = f(0) + 2*f'_desired + 2*f''

From the second equation: f'' = (f(2) - f(0) - 2*f'_desired) / 2

But the code uses a different linear combination, suggesting a higher-order or modified stencil. The exact derivation likely involves ensuring consistency with the specific discretization used in the elliptic equation at the boundary.

9. Multi-Block Support

GRAPE

The original GRAPE is a single-block solver. Multi-block extensions were developed later (e.g., GRAPE3D by Sorenson & Steger, 1984), where interface conditions between blocks must be carefully handled to ensure continuity of the grid and its derivatives.

This Implementation

Has native multi-block support through SmoothBlocks(), which iterates over each block independently:

for i in eachindex(blocks)
    xr, yr, finalError, finalIter = EllipticSolver(blocks[i][1,:,:], blocks[i][2,:,:], ...)
end

However, the blocks are smoothed independently – there is no inter-block coupling during the elliptic solve. This means:

• Block interfaces maintain their positions (treated as fixed boundaries)
• No guarantee of derivative continuity across block interfaces
• Each block converges independently

A fully coupled multi-block elliptic solver would iterate across all blocks simultaneously, updating interface positions to ensure C1 continuity (matching first derivatives) across block boundaries.

10. Summary: What's Preserved and What's Simplified

Faithfully Preserved from Steger & Sorenson

1. The elliptic PDE formulation in computational coordinates
2. The 9-point stencil discretization with alpha, beta, gamma metrics
3. Automatic orthogonality enforcement via boundary residual computation
4. Exponential decay of forcing into the domain interior
5. The overall iterate-until-convergence framework
6. Independent wall control (enable/disable per boundary)
7. The modified one-sided derivative stencil at boundaries

Simplified or Omitted

1. Point SOR instead of line Gauss-Seidel – the largest performance difference
2. No coarse-fine sequencing – no multigrid acceleration
3. Under-relaxation instead of over-relaxation – necessary due to point SOR
4. Direct RHS forcing instead of separate P,Q control functions – minor accuracy impact
5. Linear spacing interpolation instead of prescribed s(xi) – less control over wall spacing
6. Single displacement-based convergence metric (max of x and y L2 norms) – no grid-quality diagnostics
7. Independent block smoothing – no inter-block coupling for derivative continuity

Suggested Improvements (in priority order)

1. Line Gauss-Seidel: Replace point SOR with a Thomas-algorithm-based line solver. This would allow over-relaxation and reduce iteration counts by 10-100x. Partly done: EllipticParams(sweep = :line) adds alternating line Gauss-Seidel for the current equations (100x fewer iterations without forcing, ~6x with wall forcing, where omega must stay <= 0.3). Proper P,Q control functions (item 5) are likely needed before over-relaxation works with forcing.

2. Prescribed wall spacing: Allow user to specify s(xi) along each boundary rather than interpolating from corner spacings.

3. ~~Monitor both x and y convergence: Track max(norm(x-xold), norm(y-yold)).~~ Done.

4. Coarse-fine sequencing: Add a simple two-level multigrid for large grids.

5. P,Q control functions: Decompose the boundary forcing into proper P,Q functions that multiply local first derivatives when propagated into the interior.

References

1. Steger & Sorenson (1979). "Automatic mesh-point clustering near a boundary in grid generation with elliptic partial differential equations." J. Comp. Phys. 33, 405-410.

2. Sorenson (1980). "A computer program to generate two-dimensional grids about airfoils and other shapes by the use of Poisson's equation." NASA TM-81198 (the GRAPE code).

3. Sorenson & Steger (1984). "Grid generation in three dimensions by Poisson equations with control of cell size and skewness at boundary surfaces." NASA TM-86301.

4. White (1990). "Two-dimensional grid generation with derivatives of the physical coordinates on boundary orthogonal grids." NASA CR-4348. (Improved forcing propagation using transfinite interpolation.)

5. Hilgenstock (1988). "A fast method for the elliptic generation of three- dimensional grids with full boundary control." In Numerical Grid Generation in Computational Fluid Mechanics. (Efficient computation of control functions.)

6. Thompson, Warsi & Mastin (1985). Numerical Grid Generation: Foundations and Applications. (Comprehensive textbook covering the full theory.)