# GPU Geometric Multigrid for Pressure PCG

> Architecture update (2026-09-16): Pressure operators use halo-only `PressureWorkspace` objects borrowing pressure interiors. Immutable pressure descriptors are shared across levels. `PrepareGeometry` binds borrowed fine metrics while preserving separate coarse samples and retained scratch. See [runtime contracts](/docs/RUNTIME_CONTRACTS) for the supported API; older sketches below retain their historical context.

**Status:** Implemented in `src/PolyCfd.Gpu/Multigrid/` (`PcgMgSolverGpu`, `MultigridCycle`, `MultigridSmoother`,
`GpuMgHierarchy`, `GpuMgLevel`, `MultigridKernels`); scenarios select it with `numerics.pressure.method: "pcgMultigrid"`
(the `--mg` validation flag in the sections below no longer exists). The implemented contract (level plan, buffer roles,
V-cycle, smoother seam, transfers, footprints, tests) is the [Multigrid README](/src/PolyCfd.Gpu/Multigrid/README).
Cut-cell support is described in `CUT_CELL_MULTIGRID_DESIGN.md`. The rest of this document is the original design and is
kept for the reasoning behind it.

**Loop update (2026-09-26, SOL-04e):** both pressure solvers run one PCG loop, `PcgIteration`, with results
unchanged; `PcgMgSolverGpu` is a shell that lends its hierarchy and `MultigridCycle` (the multigrid
`IPressurePreconditioner`), and `PcgSolverGpu` lends its operator and `JacobiPreconditioner`. The descriptions of
per-solver loops, buffers and metrics below are historical; the loop contract is the
[Projection README](/src/PolyCfd.Gpu/Projection/README#one-pcg-loop).

**Structure update (2026-09-26, SOL-04d):** the V-cycle moved into `MultigridCycle` with unchanged arithmetic and the damped Jacobi
smoother (ω 0.7, 3/3/10 sweeps) into `MultigridSmoother`, the only owner of the smoothing method; the solver constructor
no longer takes smoother parameters, and its hierarchy diagnostics are internal. `GpuMgHierarchy.PlanLevels` is the one
coarsening rule, shared with the memory estimator.

**Interface update (2026-09-26, SOL-04b):** every `IPcgSolverGpu` member is required; `Solve`,
`PreparePressureRegions` and `GaugePressure` take `PressureBoundaryConditions` instead of a whole
`DevicePatchBcSet`, as does the multigrid `ApplyPreconditioner` diagnostic; `PrepareGeometry` is the only solver
geometry entry. `PcgSolverMetrics`, the legacy and concrete `UpdateGeometry` overloads and `IGridBlasGpu`
were removed, so the sketches below that show them are historical. The current contract is in the
[Projection README](/src/PolyCfd.Gpu/Projection/README#solver-interface).

**Current update (2026-09-16):** Cut-cell and rotating STL geometry are supported.
Coarse geometry is rebuilt from the SDF; `CoarseThickeningMode.ConstantPhysical`
is the default. Transfers use volume-weighted restriction and injection
prolongation (trilinear interpolation was dropped for cut cells); weighted-symmetry and frozen fan hierarchy tests cover this path.
The evolving fan additionally requires per-component RHS compatibility in
`ProjectionGpu`. See [MOVING_FAN_VALIDATION.md](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/MOVING_FAN_VALIDATION.md). Current
timing regression checks are described in [PERFORMANCE_VALIDATION.md](/docs/PERFORMANCE_VALIDATION).

**Nullspace update (2026-09-25, SOL-01):** the global-only, `HasPressureNullspace`-gated
mean projection described in sections 1.2, 1.4, 2.2, 5.4 and 7 is superseded. Both PCG
solvers own a `PressureComponentsGpu` and project every connected region without a
fixed-pressure boundary out of the recursive residual after each update (the multigrid
loop previously projected only the initial residual, which stalled tight periodic solves);
see [pressure nullspaces](/docs/PHYSICS_GUIDE#pcg-and-multigrid). The V-cycle itself projects
nothing; its coarse right-hand-side aliasing is documented in `MultigridCycle` and the
[Multigrid README](/src/PolyCfd.Gpu/Multigrid/README#v-cycle).

**Scope:** Add a *new* GPU pressure solver that uses a geometric multigrid (GMG) preconditioner underneath PCG, without modifying the existing `PcgSolverGpu` behavior. The goal is to enable apples-to-apples comparisons between Jacobi-PCG and MG-preconditioned PCG on the GPU, using the same matrix-free pressure operator and dual-buffer MAC layout.

---

## 1. Current GPU Pressure Solver Overview

This section summarizes the current design (as of 2025-11-16) to clarify how multigrid will plug in.

### 1.1 Core types and layout

Relevant types:

- `GridBlock` (CPU, `PolyCfd.Core.Core`)
  - Encapsulates `GridSize` and `GridMetrics`.
  - `GridMetrics` holds interior-only metrics and geometry: `Vol`, `Area{X,Y,Z}`, `Alpha{X,Y,Z}`, `Dx/Dy/Dz`.
- `GridDescriptor` (GPU, `PolyCfd.Gpu.Core`)
  - Derived from `GridBlock`; describes `(Nx, Ny, Nz)` and halo geometry for device kernels.
- `DeviceMacState` (GPU, `PolyCfd.Gpu.Core`)
  - Device-resident MAC state using **dual buffers**: interior arrays + packed face halos (no halo-inclusive unified arrays).
  - For pressure, uses interior `P` buffer plus packed `PHalo` with 6-face layout and in-plane indexing consistent with `HaloOffsets`/`HaloIndexers`.
- `VolumeWeights` (GPU, `PolyCfd.Gpu.Numerics`)
  - Owns an interior-only volume buffer on the device and the global volume sum.

### 1.2 `PcgSolverGpu`

`PcgSolverGpu` implements `IPcgSolverGpu` and is the current GPU pressure solver:

- Owns device buffers for interior-only vectors: `_d_x`, `_d_r`, `_d_z`, `_d_p`, `_d_Ap`, `_d_Minv`.
- Builds a Jacobi preconditioner on the CPU via `BuildJacobiPreconditioner()`, using interior metrics only.
- Uses `PressureOperatorGpu` for `A·v` and `GridBlasGpu` for vector ops and weighted reductions.
- Uses `BoundaryKernels` + `DeviceMacState` for GPU-side ghost refreshes before each operator application.
- Supports `Solve(DeviceBuffer<float> d_rhs, DeviceBuffer<float> d_pressure, SolverOptions options, DevicePatchBcSet patchBcSet)` for GPU-resident data.
- Handles pressure null space by volume-weighted zero-mean projection of the initial guess and residual **only when** `bcSet.HasPressureNullspace()` is true.
- Tracks optional performance metrics via `PcgSolverMetrics` and `GpuTimingHelper`.

### 1.3 Operator and BC integration

- Operator: `PressureOperatorGpu.Apply(DeviceBuffer<float> pInterior, DeviceBuffer<float> pHalo, DeviceBuffer<float> result)`
  - Matrix-free Laplacian using interior metrics and halo pressure values.
  - Assumes halos are pre-refreshed and uses dual-buffer MAC layout.
- BCs: `BoundaryKernels.ApplyPressureBcs(DeviceMacState state, BcSet bcSet, DeviceBuffer<float> pInterior)`
  - Populates `state.PHalo` from interior `p` given `bcSet`.
  - Updates `DeviceMacState.HaloVersion` to detect stale halos.

### 1.4 Nullspace & stopping criteria

- Nullspace:
  - For all-Neumann/periodic pressure BCs, the system is singular up to a constant.
  - `PcgSolverGpu` projects the mean out of `x` and `r` using `GridBlasGpu.ProjectOutMean` with volume weights.
- Stopping criteria:
  - Weighted residual norm `||r||_w` via `Nrm2Weighted`.
  - Dual criteria: `||r||_w < AbsTolerance` **OR** `||r||_w / ||r0||_w < RelTolerance`.
  - Optional breakdown tracking (`PcgBreakdownReason`) for `pᵀAp ≈ 0`, `rᵀz ≈ 0`, inconsistent RHS mean, or max iterations.

**Key point for multigrid:** Multigrid must respect the same weighted inner product (volume weights) and BC treatment so that the overall PCG+MG scheme remains SPD/symmetric.

---

## 2. High-Level Multigrid Design

### 2.1 Goals

- Provide a **geometric multigrid preconditioner** for the GPU pressure operator.
- Implemented as a *separate* solver class (e.g., `PcgMgSolverGpu`) to avoid changing `PcgSolverGpu` behavior.
- Support uniform Cartesian blocks initially (no cut cells / `CellGeometry` variations beyond trivial uniform geometry).
- Use the same matrix-free `PressureOperatorGpu` at each level.
- Respect dual-buffer layout and GPU boundary kernels at all levels.
- Reuse as much existing infrastructure as possible (`GridBlasGpu`, `BoundaryKernels`, `VolumeWeights`).

### 2.2 Multigrid strategy

- Type: **V-cycle** multigrid used as a **right preconditioner** inside PCG.
- Coarsening: standard 2× coarsening in each dimension:
  - Level `ℓ` has size `(Nxℓ, Nyℓ, Nzℓ)` with `Nxℓ+1 = ceil((Nxℓ)/2)`, etc., while respecting halo face counts.
  - Stop when any dimension drops below a minimum threshold (e.g., 4 or 2 cells).
- Discretization:
  - At each level, reuse same form of pressure operator:
    - Use level-specific `GridBlock` and `GridMetrics` with recomputed `Vol` and `Area` (for uniform grids these are analytic).
  - Coarse operator is **not** assembled; multigrid remains matrix-free.
- Smoother:
  - Damped Jacobi (few iterations per pre- and post-smoothing, e.g., 3–5) on the GPU.
  - Uses diagonal approximations similar to `PcgSolverGpu.BuildJacobiPreconditioner()` but per-level.
- Restriction & prolongation:
  - **Restriction:** volume-weighted for cell-centered fields (residual).
  - **Prolongation:** trilinear interpolation for cell-centered corrections.
- Nullspace handling:
  - At the finest level, PCG applies the same zero-mean projection as `PcgSolverGpu` for nullspace configurations.
  - Coarse levels should maintain zero-mean via restriction/prolongation that preserve sums under volume weights.

### 2.3 Interaction with PCG

We introduce a new solver class, e.g., `PcgMgSolverGpu`, with the following characteristics:

- Implements `IPcgSolverGpu` but internally:
  - Uses PCG on the finest level.
  - Uses a multigrid V-cycle as the preconditioning operation `M^{-1}`.
- From the PCG perspective:
  - Replace Jacobi preconditioner `z = M^{-1} r` with `z = VCycle(r)`.
  - V-cycle must be **symmetric** (or at least sufficiently close in practice) to maintain PCG convergence; the classic symmetric V-cycle (same smoother forward/backward) is appropriate.

---

## 3. Level Hierarchy and Data Structures

### 3.1 Multigrid level descriptor

Introduce a new GPU-side data structure representing a single multigrid level:

```csharp
internal sealed class GpuMgLevel : IDisposable
{
    public GridSize Size { get; }
    public GridBlock Block { get; }
    public GridDescriptor Descriptor { get; }
    public DeviceMacState PressureState { get; }
    public PressureOperatorGpu Operator { get; }
    public GridBlasGpu Blas { get; }
    public VolumeWeights VolumeWeights { get; }

    // Smoother: Jacobi diagonal on this level
    public DeviceBuffer<float> DInv { get; }

    // Scratch vectors: r, e, Ae
    public DeviceBuffer<float> Residual { get; }
    public DeviceBuffer<float> Correction { get; }
    public DeviceBuffer<float> Temp { get; }
}
```

Notes:

- `Block` holds the per-level `GridMetrics` (dx, dy, dz, Vol, Area, Alpha), recomputed for coarse grids.
- `PressureState` is a minimal `DeviceMacState` variant focused on pressure and its halo; velocity components are not needed for MG.
- `VolumeWeights` is per-level (volume sum differs across levels).
- `DInv` holds per-level Jacobi diagonals for the smoother.

### 3.2 Multigrid hierarchy container

Add a hierarchy class managing all levels:

```csharp
internal sealed class GpuMgHierarchy : IDisposable
{
    public IReadOnlyList<GpuMgLevel> Levels { get; }

    // Construction from finest GridBlock
    public GpuMgHierarchy(IlgpuBackend backend, GridBlock finestBlock, BcSet finestBc, int maxLevels); 

    public int LevelCount => Levels.Count;
}
```

Responsibility:

- Build coarsened `GridSize` / `GridBlock` sequences from finest to coarsest.
- Construct `GpuMgLevel` for each level, including `GridDescriptor`, `PressureOperatorGpu`, `GridBlasGpu`, `VolumeWeights` and `DInv`.
- Handle BC propagation to coarse levels (see §4.3).

---

## 4. Operators, Smoothers, and Transfers

### 4.1 Smoother (Jacobi)

For each level `ℓ`:

- Precompute diagonal `Dℓ` on the CPU using level metrics, mirroring `BuildJacobiPreconditioner()` in `PcgSolverGpu`.
- Transfer `DInvℓ = 1 / Dℓ` to the GPU as `DeviceBuffer<float>`.
- Jacobi smoother kernel (device):
  - `x ← x + ω * DInvℓ ∘ (r - A x)` (with `ω` damping, e.g., 0.7).
  - Use `BoundaryKernels` to refresh halos for `x` before applying `A` at each smoothing sweep.
- Wrap smoothing steps in a helper on `GpuMgLevel`, e.g., `Smooth(x, b, int iterations, BcSet bc)`.

### 4.2 Restriction and prolongation

**Restriction (fine → coarse):**

- For uniform grids, use volume-weighted averaging:

  - Conceptually: for each coarse cell `(I,J,K)`, average the residuals of the `2×2×2` fine cells it covers, weighted by cell volumes:

    $$ r_c(I,J,K) = \frac{\sum_{i,j,k \in \text{children}} V_f(i,j,k) r_f(i,j,k)}{\sum_{i,j,k \in \text{children}} V_f(i,j,k)} $$

- Implemented as a GPU kernel operating on interior-only arrays; no halo data required for residuals.

**Prolongation (coarse → fine):**

- Use trilinear interpolation of coarse corrections to fine grid points:

  - For each fine cell center, identify the owning coarse cell and its neighbors and compute trilinear weights.
  - For even/odd indices, we can use a standard stencil (identity on coarse nodes, averages on midpoints, etc.).

- Prolongation is applied to **corrections** `e_c` to obtain `e_f`; then `x_f ← x_f + e_f`.

### 4.3 Boundary conditions on coarse levels

BC strategy per level:

- The **type** of BC (Dirichlet, Neumann, Periodic) remains the same across levels.
- `BcSet` is reused at all levels without modification to type.
- Geometry and metrics at coarse levels are recomputed assuming uniform spacing; ghost treatment uses the same `BoundaryKernels` logic.
- This is sufficient for the current M1 uniform-domain focus. For cut-cell geometries (future M3), the coarse-level α/area/Vol would need careful restriction rules.

---

## 5. New Solver: `PcgMgSolverGpu`

### 5.1 API

Introduce a new GPU solver class alongside `PcgSolverGpu`:

```csharp
public sealed class PcgMgSolverGpu : IPcgSolverGpu, IDisposable
{
    public PcgSolverMetrics? Metrics { get; set; }

    public PcgMgSolverGpu(
        IlgpuBackend backend,
        GridBlock block,
        int maxLevels,
        bool projectOutMean = true,
        ILogger? logger = null);

    public SolveResult Solve(DeviceBuffer<float> d_rhs, DeviceBuffer<float> d_pressure, SolverOptions options, DevicePatchBcSet patchBcSet);
}
```

Notes:

- `SolveDual` can mirror the `PcgSolverGpu` behavior by using periodic defaults or, in future, hooking into a GPU-aware `BcRegistry`.
- The existing `PcgSolverGpu` remains unchanged; callers can select either solver.

### 5.2 Internal structure

`PcgMgSolverGpu` will:

- Own:
  - A `GpuMgHierarchy` built from the finest `GridBlock`.
  - Finest-level device vectors for PCG: `_d_x`, `_d_r`, `_d_z`, `_d_p`, `_d_Ap`.
  - Finest-level `VolumeWeights` shared with the finest `GpuMgLevel`.
- Use the same nullspace and convergence logic as `PcgSolverGpu` at the finest level.
- Replace Jacobi preconditioning step with a V-cycle call, e.g.:

```csharp
// Instead of z = Minv * r
MgVCycle(level = 0, z = _d_z, r = _d_r, bcSet);
```

### 5.3 V-cycle pseudocode

For level `ℓ` with solution `xℓ`, RHS `bℓ`, and residual `rℓ`:

```text
VCycle(ℓ, xℓ, bℓ):
  if ℓ == Lmax (coarsest level):
      // Coarse solve: small fixed iteration count of Jacobi or mini-PCG
      Smooth(ℓ, xℓ, bℓ, iterations = N_coarse)
      return

  // Pre-smoothing
  Smooth(ℓ, xℓ, bℓ, iterations = N_pre)

  // Compute residual rℓ = bℓ - Aℓ xℓ
  ApplyPressureBcs(ℓ, xℓ)
  rℓ = bℓ - Aℓ xℓ

  // Restrict residual to coarse level: bℓ+1 = Rℓ→ℓ+1 rℓ
  bℓ+1 = Restrict(rℓ)

  // Initialize coarse correction eℓ+1 = 0
  eℓ+1 = 0

  // Recursive V-cycle on coarse level
  VCycle(ℓ+1, eℓ+1, bℓ+1)

  // Prolongate correction and update xℓ
  xℓ += Pℓ+1→ℓ eℓ+1

  // Post-smoothing
  Smooth(ℓ, xℓ, bℓ, iterations = N_post)
```

For PCG preconditioning we will:

- Allocate temporary vectors for `xℓ` and `bℓ` at each level.
- On finest level, for each preconditioning call `z = M^{-1} r`:
  - Initialize solution on finest: `x0 = 0`, `b0 = r`.
  - Invoke `VCycle(0, x0, b0)`.
  - Return `z = x0`.

### 5.4 Nullspace consistency

To preserve compatibility with the finest-level nullspace handling:

- Ensure restriction and prolongation operators preserve zero-mean subspace under volume weights.
- Optionally apply a mean-projection at the end of each V-cycle on the finest level when `bcSet.HasPressureNullspace()` and `_projectOutMean` is true.
- This keeps the MG preconditioner aligned with the SPD structure expected by PCG.

---

## 6. Integration and Testing Strategy

### 6.1 Validation targets

- **Unit tests (GPU):**
  - New tests in `PolyCfd.Gpu.Tests` for:
    - Restriction/prolongation correctness on simple analytic fields (constants, linear, quadratic).
    - Smoother behavior (residual reduction for a single application on synthetic problems).
    - Multigrid V-cycle reducing residual for a random RHS on small grids (e.g., 8³, 16³).
  - Compare `PcgMgSolverGpu` vs `PcgSolverGpu` on uniform grids:
    - Check they converge to similar residual norms / solutions for the same RHS and BCs.
    - Check iteration counts and wall-clock timings.

- **Integration tests:**
  - Hook `PcgMgSolverGpu` into `PolyCfd.Validation` optionally (e.g., via command-line flag) to run lid-driven cavity and Taylor–Green on GPU with MG preconditioner.
  - Verify final max|div| and global metrics against the existing Jacobi-PCG GPU solver.

### 6.2 Performance metrics

- Extend or reuse `PcgSolverMetrics` to:
  - Track number of V-cycles per solve and per iteration.
  - Track operator applications per level.
  - Track relative time spent in smoothing vs restriction/prolongation vs top-level PCG operations.

### 6.3 Configuration and feature gating

- The validation CLI selects the solver with `--mg` (alias `--use-multigrid`); without it the Jacobi-preconditioned `PcgSolverGpu` is used. (The originally proposed `--pressure-solver` option was not implemented.)

---

## 7. Task List

This section captures the concrete engineering tasks needed to implement the above design.

### 7.1 Infrastructure and data structures

- [x] **T1:** Add new design doc (`GPU_MULTIGRID_DESIGN.md`) – *this file*.
- [x] **T2:** Define `GpuMgLevel` class (internal) in `PolyCfd.Gpu`:
  - [x] Holds per-level `GridBlock`, `GridDescriptor`, `DeviceMacState`, `PressureOperatorGpu`, `GridBlasGpu`, `VolumeWeights`.
  - [x] Allocates per-level Jacobi diagonal buffer (`DInv`) and scratch vectors (`Residual`, `Correction`, `Temp`).
  - [x] Implements `Dispose()` to free device resources.
- [x] **T3:** Define `GpuMgHierarchy` class (internal) responsible for:
  - [x] Building level sizes from finest grid using 2× coarsening.
  - [x] Creating per-level `GridBlock` with recomputed metrics for uniform grids.
  - [x] Constructing `GpuMgLevel` instances for each level.

### 7.2 Level operators and smoother

- [x] **T4:** Implement per-level Jacobi diagonal computation (CPU-side):
  - [x] Mirror `PcgSolverGpu.BuildJacobiPreconditioner()` but operate on a general `GridBlock`.
  - [x] Transfer `DInv` to device for each level.
- [x] **T5:** Implement Jacobi smoother kernel and wrapper on `GpuMgLevel`:
  - [x] GPU kernel performing `x ← x + ω * DInv ∘ (b - A x)`.
  - [x] Use `BoundaryKernels.ApplyPressureBcs` on each smoothing sweep.
  - [x] Respect level-specific `BcSet` (mirroring finest level BCs).

### 7.3 Transfer operators (restriction & prolongation)

- [x] **T6:** Implement restriction kernel for cell-centered residuals:
  - [x] Volume-weighted `2×2×2` aggregation from fine to coarse grid.
  - [x] Handle non-divisible grid sizes via appropriate mapping.
- [x] **T7:** Implement prolongation kernel for cell-centered corrections:
  - [x] Trilinear interpolation from coarse to fine grid.
  - [x] Validate against known analytic fields.

### 7.4 V-cycle implementation

- [x] **T8:** Implement recursive or iterative V-cycle routine using `GpuMgHierarchy`:
  - [x] Pre-smoothing and post-smoothing counts configurable (e.g., via `SolverOptions` or MG-specific options).
  - [x] Coarse solve implemented as a small fixed-count Jacobi or mini-PCG on coarsest level.
  - [x] Optionally enforce zero-mean at finest level for nullspace BCs.

### 7.5 New solver `PcgMgSolverGpu`

- [x] **T9:** Add `PcgMgSolverGpu` class in `PolyCfd.Gpu`:
  - [x] Construct `GpuMgHierarchy` in constructor.
  - [x] Implement PCG logic, mirroring `PcgSolverGpu` but replacing Jacobi preconditioner with V-cycle.
  - [x] Integrate nullspace handling: use volume-weighted mean projection of `x` and `r` when `bcSet.HasPressureNullspace()`.
  - [x] Hook in `PcgSolverMetrics` / `GpuTimingHelper` for instrumentation.
- [x] **T10:** Implement `IPcgSolverGpu` interface and legacy overloads for backwards compatibility.

### 7.6 Tests and validation

- [x] **T11:** Add GPU unit tests (in `PolyCfd.Gpu.Tests`) for:
  - [x] Restriction/prolongation consistency on constants and linear fields.
  - [x] Smoother reducing residual on synthetic Laplacian problems.
  - [x] Single V-cycle residual reduction on 8³ and 16³ grids.
- [x] **T12:** Add comparison tests between `PcgSolverGpu` and `PcgMgSolverGpu`:
  - [x] Same RHS and BCs, assert similar final residual and solution on simple problems.
  - [x] Assert MG-preconditioned PCG converges in fewer iterations than pure Jacobi-PCG.

### 7.7 Integration and benchmarking

- [x] **T13:** Expose solver choice in GPU backend and/or validation runner:
  - [x] Add flag or configuration to select `PcgMgSolverGpu`.
  - [x] Ensure default still uses `PcgSolverGpu` to keep current behavior.
  - Implementation notes:
    - Added `UseMultigrid` flag to `ValidationParametersBase`
    - Modified `TimeIntegratorGpu` to accept `useMultigrid` parameter
    - Updated `ProjectionGpu` to dynamically select solver based on flag
    - Added `DeviceBuffer` overload to `PcgMgSolverGpu.Solve()` for consistency
    - Modified `LidDrivenCavity` and `TaylorGreenVortex` validation cases to support multigrid
    - Backend tags: "gpu" (Jacobi), "gpu-mg" (Multigrid), "gpu-cublas" (Jacobi+cuBLAS)
- [ ] **T14:** Add simple benchmark scenarios:
  - [ ] 64³, 128³ lid-driven cavity GPU runs comparing iteration counts and wall time.
  - [ ] Document results in `PERF_NOTES.md` once stable.

---

## 9. Historical Implementation Status

**Date:** January 2025

**Summary:** Core multigrid implementation complete and integrated into validation runner.

### Completed

1. **Infrastructure (T1-T3):** ✅
   - Created `GpuMgLevel` class with per-level operators, BLAS, and scratch vectors
   - Created `GpuMgHierarchy` with automatic 2× coarsening
   - Jacobi diagonal computation implemented for all levels

2. **Smoother and Transfer Operators (T4-T7):** ✅
   - Damped Jacobi smoother kernel implemented
   - Volume-weighted restriction (2×2×2 averaging)
   - Injection-based prolongation (nearest neighbor)
   - All kernels working on GPU

3. **V-Cycle and Solver (T8-T10):** ✅
   - Recursive V-cycle implementation with configurable pre/post smoothing
   - `PcgMgSolverGpu` class implementing `IPcgSolverGpu`
   - Nullspace handling integrated (zero-mean projection for periodic BCs)

4. **Tests (T11-T12):** ✅
   - Unit tests in `PolyCfd.Gpu.Tests/Multigrid/`
   - Zero RHS test (sanity check)
   - Convergence speed comparison: MG converges 5-6× faster than Jacobi-PCG
   - Example: 16³ grid with random RHS: MG = 7 iterations, Jacobi = 39 iterations

5. **Integration (T13):** ✅
   - Added `UseMultigrid` flag to validation runner parameters
   - Modified `TimeIntegratorGpu` and `ProjectionGpu` to support solver selection
   - Both `LidDrivenCavity` and `TaylorGreenVortex` cases support multigrid
   - Backend tagging: "gpu-mg" for multigrid runs, "gpu" for Jacobi

### Remaining

- **Performance benchmarking (T14):** Large-scale validation runs needed
  - Planned: 64³ and 128³ lid-driven cavity comparisons
  - Document iteration counts, wall-clock times, and speedup factors
  - Record results in `PERF_NOTES.md` or baseline documentation

### Limitations at that snapshot (superseded where noted above)

- Prolongation uses injection (nearest neighbor) rather than full trilinear interpolation
  - Sufficient for uniform grids and provides good convergence
  - Could be enhanced for better high-frequency error reduction
- No cut-cell/non-uniform geometry support yet (M3 milestone)
- MG hierarchy always built with 2× coarsening; no adaptive coarsening

---

## 8. Future Extensions

- **Cut-cell support:** extend restriction/prolongation to handle `CellGeometry` with non-trivial `Alpha`/`Vol`/`Area` fields, respecting `MinVolumeFraction` and agglomeration.
- **Alternative smoothers:** Chebyshev(2) smoother, red-black Gauss–Seidel on GPU.
- **Adaptive coarsening:** stop coarsening earlier in highly anisotropic grids, or perform semi-coarsening when one dimension is much smaller.
- **Backend abstraction:** mirror MG design on CPU side and unify via a backend-agnostic `IMgHierarchy` interface when multiblock/GPU-portable backends (M6) are implemented.