SIMD Optimization Analysis for Projection Method Solver

Overview

This document analyzes the SIMD (AVX2) optimization of the projection method solver for incompressible Navier-Stokes equations, explaining which parts can be vectorized and why the Poisson solver remains a bottleneck.

Projection Method Steps

The projection method (Chorin's method) consists of three main steps:

Predictor Step: Compute intermediate velocity u* ignoring pressure
Poisson Solve: Solve pressure equation to enforce incompressibility
Corrector Step: Project velocity to be divergence-free

SIMD Optimization Status

Step	Vectorized	Reason
Predictor	❌ No	Scalar due to complex stencil access pattern
RHS Computation	❌ No	Scalar divergence calculation (simple, not performance-critical)
Poisson Solve	✅ Yes	Uses CG SIMD solver (Conjugate Gradient)
Corrector	✅ Yes	Independent gradient subtraction (AVX2)

Note: As of December 2024, the SIMD projection solver uses the CG SIMD Poisson solver (POISSON_SOLVER_CG_SIMD) which provides robust convergence in ~150 iterations.

Why the Poisson Solver Cannot Be Easily Vectorized

The SOR Algorithm's Data Dependencies

The Successive Over-Relaxation (SOR) method uses a Gauss-Seidel update pattern:

for (size_t i = 1; i < nx - 1; i++) {
    size_t idx = (j * nx) + i;

    // p[idx-1] was JUST updated in the previous iteration
    double p_new = (rhs[idx] - (p[idx + 1] + p[idx - 1]) / dx2
                             - (p[idx + nx] + p[idx - nx]) / dy2) * (-inv_factor);

    // Immediate in-place update
    p[idx] = p[idx] + POISSON_OMEGA * (p_new - p[idx]);
}

The Problem: Read-After-Write Dependency

When computing p[i], the algorithm needs p[i-1] which was updated in the same sweep. This creates a sequential dependency chain:

p[1] → p[2] → p[3] → p[4] → ...

SIMD (AVX2) processes 4 doubles simultaneously, but:

Value at index 2 needs the result from index 1
Value at index 3 needs the result from index 2
etc.

This read-after-write hazard prevents parallel execution.

Why This Matters for Convergence

The SOR method's fast convergence relies on using the most recent values. This is what makes Gauss-Seidel converge ~2x faster than Jacobi iteration. Breaking this dependency would require switching to a different algorithm.

Update (March 2025): A Block SOR SIMD implementation has been added that partially breaks this dependency at SIMD block boundaries. Within each block of 4 (AVX2) or 2 (NEON) consecutive cells, the left-neighbor uses stale values, but between blocks the dependency is satisfied. This is a well-known HPC technique that trades slightly slower convergence per iteration for significantly higher throughput. See Block SOR technical note for details.

Performance Impact

Typical Cost Distribution

For a projection method time step:

Component	% of Total Time	SIMD Speedup
Predictor	~15%	1x (scalar)
RHS Computation	~5%	1x (scalar)
Poisson Solve	~70-80%	2-3x (CG SIMD with AVX2 primitives)
Corrector	~5%	2-4x (AVX2)
Boundary Conditions	~5%	1x

Amdahl's Law Analysis (Updated December 2024)

With SIMD Poisson solver now integrated:

Speedup = 1 / ((1 - P) + P/S)

Where:
- P = parallelizable fraction ≈ 0.80 (Poisson + Corrector)
- S = speedup factor ≈ 2.5x (CG SIMD with optimized BLAS primitives)

Speedup = 1 / (0.20 + 0.80/2.5) ≈ 1 / (0.20 + 0.32) ≈ 1.9x

Expected speedup: ~1.3-1.5x with current SIMD implementation. Further gains possible with:

OpenMP parallelization combined with SIMD
Multigrid preconditioner for faster convergence

Alternative Approaches

1. Jacobi Iteration (GPU Solver Approach)

The GPU solver uses Jacobi iteration which IS parallelizable:

// All reads from OLD array - no dependencies!
p_new[idx] = (rhs[idx] - (p_old[idx+1] + p_old[idx-1]) / dx2
                       - (p_old[idx+nx] + p_old[idx-nx]) / dy2) * (-inv_factor);

Trade-off: Jacobi converges ~2x slower, requiring more iterations.

2. Red-Black SOR (Partial Vectorization)

Within each "color" of a red-black ordering, updates are independent:

Red cells:   (i+j) % 2 == 0  → Can be vectorized together
Black cells: (i+j) % 2 == 1  → Can be vectorized together

This allows SIMD within each color sweep while maintaining SOR convergence.

3. Multigrid Methods

Multigrid solvers offer:

O(N) complexity vs O(N²) for iterative methods
Natural parallelism at each grid level
Can be combined with SIMD at each level

Implementation Details

Vectorized Predictor Step (AVX2)

// Load 4 velocity values at once
__m256d u = _mm256_loadu_pd(&field->u[idx]);
__m256d v = _mm256_loadu_pd(&field->v[idx]);

// Load neighbors
__m256d u_xp = _mm256_loadu_pd(&field->u[idx + 1]);
__m256d u_xm = _mm256_loadu_pd(&field->u[idx - 1]);

// Compute derivatives (4 cells simultaneously)
__m256d du_dx = _mm256_mul_pd(_mm256_sub_pd(u_xp, u_xm), dx_inv);

// Convection, diffusion, time integration...
__m256d u_star = _mm256_add_pd(u,
    _mm256_mul_pd(dt_vec, _mm256_sub_pd(visc_u, conv_u)));

// Store 4 results
_mm256_storeu_pd(&ctx->u_star[idx], u_star);

Scalar Remainder Handling

For grid sizes not divisible by 4:

// SIMD loop processes 4 at a time
for (i = 1; i + 4 <= nx - 1; i += 4) {
    // AVX2 code...
}

// Scalar cleanup for remaining 0-3 elements
for (; i < nx - 1; i++) {
    // Scalar code...
}

MSVC Compiler Note

MSVC does not define __AVX2__ even when using /arch:AVX2. The CMake configuration manually defines it:

if(MSVC)
    if(CMAKE_SYSTEM_PROCESSOR MATCHES "x86_64|AMD64|i686")
        target_compile_options(${target_name} PRIVATE /arch:AVX2)
        target_compile_definitions(${target_name} PRIVATE __AVX2__=1)
    endif()
endif()

Recommendations

For CPU-bound workloads: The current SIMD implementation provides modest speedup (~1.2x) but the Poisson solver remains the bottleneck.
For maximum CPU performance: Consider implementing Red-Black SOR with SIMD for the Poisson solver.
For large grids: Use the GPU solver (SOLVER_TYPE_PROJECTION_JACOBI_GPU) which parallelizes the entire algorithm including the Poisson solve.
For production use: Profile your specific workload to determine if the Poisson solver is indeed the bottleneck, as the ratio varies with grid size and iteration counts.

Files Modified

lib/src/solvers/navier_stokes/avx2/solver_projection_avx2.c - AVX2 projection method implementation
lib/src/solvers/linear/avx2/linear_solver_cg_avx2.c - CG SIMD Poisson solver
lib/CMakeLists.txt - Added __AVX2__ define for MSVC

Test Results

After SIMD optimization (December 2024):

All 6 SIMD projection tests pass
L2 difference vs scalar: ~10⁻¹⁹ (essentially machine precision - results match exactly)
Energy decay ratio: 0.95 (proper physical behavior)
Divergence norm: within tolerance
Non-aligned grid sizes (33x35) handled correctly

Poisson Solver SIMD Implementation Status

Implemented Solvers

Two SIMD Poisson solvers have been implemented in separate files:

Solver	File	Convergence	SIMD Efficiency
Jacobi SIMD	`poisson_jacobi_simd.c`	Slow (~2x more iterations)	High (full vectorization)
Red-Black SIMD	`poisson_redblack_simd.c`	Fast (SOR convergence)	Medium (gather/scatter overhead)

Convergence Analysis

Test Problem: Sinusoidal RHS with p = sin(πx)sin(πy) analytical solution

Solver	Max Iterations	Tolerance	Converges on 16x16?
Jacobi SIMD	2000	1e-6	❌ No (reaches max iter)
Red-Black SIMD	1000	1e-6	❌ No (reaches max iter)
Both with Zero RHS	-	1e-6	✅ Yes (immediate)

Key Finding: The SIMD Poisson solvers do NOT fully converge to 1e-6 tolerance on the sinusoidal test problem within the iteration limits. However:

They produce valid results (no NaN/Inf)
They do converge with simpler problems (zero RHS)
The residual decreases but doesn't reach the strict tolerance
This is expected behavior for iterative Poisson solvers without preconditioning

Test Coverage

Both Jacobi and Red-Black SIMD solvers have dedicated test files:

test_poisson_jacobi_simd.c - 8 tests
test_poisson_redblack_simd.c - 10 tests

Tests verify:

✅ Valid output (no NaN/Inf) on challenging problems
✅ Actual convergence with zero RHS
✅ Deterministic results
✅ Non-aligned grid sizes (SIMD correctness)
✅ Boundary conditions preserved
✅ Uniform RHS handling

TODO: Remaining Work

Note: These TODOs have been migrated to ROADMAP.md Section 4.7 for centralized project tracking. The ROADMAP provides current status, priority justifications, and detailed implementation plans. This section remains for historical context and technical reference.

High Priority

~~Integrate SIMD Poisson into Projection Solver~~ ✅ COMPLETED (December 2024)
- solver_projection_avx2.c now uses CG SIMD Poisson solver (POISSON_SOLVER_CG_SIMD)
- Full projection method implemented with AVX2-optimized corrector step
- All tests pass with results matching scalar implementation
- CG provides reliable convergence in ~150 iterations vs. thousands for Jacobi
Improve Convergence for Non-Trivial Problems
- Increase POISSON_MAX_ITER (currently 1000/2000)
- Consider adaptive tolerance based on problem scale
- Add optional multigrid preconditioner
Performance Benchmarking in Release Mode
- Current tests run in Debug mode where SIMD is slower
- Need Release mode benchmarks to measure actual speedup

Medium Priority

~~Add p_temp Buffer to Projection Context~~ ✅ COMPLETED
- Projection context now uses u_new buffer as temp for Poisson solver
- Red-Black SIMD works in-place; Jacobi would need dedicated buffer if used
OpenMP + SIMD Hybrid
- Combine thread parallelism with SIMD
- Jacobi allows full parallelization across rows

Low Priority

Multigrid SIMD Implementation
- Would provide O(N) complexity
- Each level can use SIMD Jacobi/Red-Black smoothing

Files Added

lib/src/solvers/linear/avx2/linear_solver_cg_avx2.c - AVX2 CG implementation
lib/src/solvers/linear/avx2/linear_solver_jacobi_avx2.c - AVX2 Jacobi implementation
lib/src/solvers/linear/avx2/linear_solver_redblack_avx2.c - AVX2 Red-Black SOR
lib/src/solvers/navier_stokes/avx2/solver_projection_avx2.c - AVX2 projection method
tests/solvers/linear/avx2/test_linear_solver_cg_avx2.c - CG tests
tests/solvers/linear/avx2/test_linear_solver_jacobi_avx2.c - Jacobi tests
tests/solvers/linear/avx2/test_linear_solver_redblack_avx2.c - Red-Black tests

FilesExpand file tree

simd-optimization-analysis.md

Latest commit

History