# Cut-Cell Geometry Support for GPU Multigrid

> Architecture update (2026-09-16): Fine metrics are borrowed through the solver geometry contract; coarse geometry owns distinct SDF samples. Pressure halos and immutable boundary descriptors are separated from velocity state/BC uploads. Geometry changes invalidate weights, diagonals and connectivity. See [runtime contracts](/docs/RUNTIME_CONTRACTS) for the supported API; older sketches below retain their historical context.

**Status:** ✅ **Implemented and Tested**  
**Date Created:** November 17, 2025  
**Last Reviewed:** 2026-09-16; original implementation record dated November 17, 2025

**Related:** `GPU_MULTIGRID_DESIGN.md`, `MOVING_GEOMETRY_DESIGN.md`

**Current contract (2026-09-26):** the implemented level plan, coarse geometry rebuilds, transfers (volume-weighted
restriction, injection prolongation; no trilinear interpolation) and tests are in the
[Multigrid README](/src/PolyCfd.Gpu/Multigrid/README). `PcgMgSolverGpu` now lives in `src/PolyCfd.Gpu/Multigrid/`
and the V-cycle in `MultigridCycle`; the file names, the `VCycle()` task and the `--mg` flag below are historical.

## Current implementation notes (2026-09-16)

- Face apertures use corner-based fractions and cell volumes use a fitted plane,
  with a centre-based fallback for sub-cell features; the original centre-sample
  sketches below are historical.
- `SdfGeometryDescriptor.FromStl` supports stationary STL geometry. Rotating STL
  cases use `CachedSdf` / `RotatedCachedSdf` and rebuild coarse metrics from that
  SDF on each update; they do not automatically fall back to Jacobi-PCG.
- `CoarseThickeningMode.ConstantPhysical` keeps the finest physical offset at
  each coarse level and is now the default; `SameCellFraction` remains available.
- `MultigridSymmetryTests` check the weighted V-cycle. `FanHierarchyTests` cover
  six frozen blade angles, including the first two adaptive steps. Its compatible
  test RHS does not establish moving-source compatibility: the evolving fan also
  needs per-region RHS compatibility and the rigid wall-flux source. See
  [MOVING_FAN_VALIDATION.md](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/MOVING_FAN_VALIDATION.md).
- (2026-09-25, SOL-01) The finest-level PCG projects every connected pressure
  region without a fixed-pressure boundary (`PressureComponentsGpu`, owned by the
  solver) out of the recursive residual after each update; coarse levels project
  nothing. Fine cut-cell geometry commonly has more regions than its coarse levels
  (a 12-cell fan pocket at 128³ is merged at 64³); the fine projection handles them.
  See [pressure nullspaces](/docs/PHYSICS_GUIDE#pcg-and-multigrid).
- Timing gates live in [PERFORMANCE_VALIDATION.md](/docs/PERFORMANCE_VALIDATION).
  The test counts, speedup ratios, and task estimates below are dated records,
  not current suite counts or performance guarantees.

## Original implementation record

**Implementation Progress:**
- ✅ Phase 1: Core abstraction (`IGeometryDescriptor`, `SdfGeometryDescriptor`) complete
- ✅ Phase 2: Coarse geometry construction in MG hierarchy complete
- ✅ Phase 3: Geometry-aware transfer operators (restriction/prolongation) complete
- ✅ Phase 3: Explicit solid cell masking in `PressureOperatorGpu` complete
- ✅ Phase 4: Integration testing complete - all tests passing

**Test Results:**
- 6/6 unit tests passing (`GeometryDescriptorTests`)
- 3/3 integration tests passing (`CutCellMultigridTests`)
  - Sphere obstacle: MG converges successfully
  - Box obstacle: MG 1.5-3x faster than Jacobi-PCG  
  - Multi-level geometry: Maintains geometric fidelity across resolutions

---

## 1. Overview

This document tracks the design and implementation of cut-cell geometry support for the GPU geometric multigrid (GMG) solver. The goal is to enable the `PcgMgSolverGpu` to handle interior obstacles (cut-cell geometries) while maintaining geometric fidelity and multigrid efficiency across coarse levels.

### Current State

✅ **What Works:**
- GPU multigrid solver for uniform grids (no cut-cells)
- Cut-cell geometry on CPU via `CutCellBuilder` (SDF-based)
- GPU operators (`PressureOperatorGpu`, `DivGradGpu`) respect geometry metrics at finest level
- V-cycle with restriction, prolongation, Jacobi smoothing

❌ **What's Missing:**
- Coarse-level geometry construction for cut-cells
- Geometry-aware restriction/prolongation (skip solid cells)
- Architecture to rebuild geometry at multiple resolutions

---

## 2. Design Principles

### 2.1 Geometric Fidelity via SDF Re-evaluation

**Decision:** Use the **original signed distance function (SDF)** to rebuild geometry at each multigrid level, rather than averaging fine-level metrics.

**Rationale:**
- **Correctness:** A sphere obstacle remains a sphere at all resolutions; averaging volumes/apertures distorts geometry
- **Accurate apertures:** Face apertures computed from actual SDF values at coarse face centers, not averaged from fine faces
- **Proper normals:** Interface normals computed from SDF gradients at coarse resolution
- **Standard practice:** Rediscretization is the canonical approach for embedded-boundary multigrid (e.g., Chombo, EB-AMR)

**Alternative rejected:**
- Volume/area/aperture averaging: geometrically incorrect, causes convergence issues

### 2.2 Geometry Abstraction Layer

**Decision:** Introduce `IGeometryDescriptor` interface to decouple geometry definition from grid resolution.

**Purpose:**
- Enable geometry rebuilding at arbitrary resolutions (for MG levels)
- Support multiple geometry sources (SDF primitives, composed shapes, future: STL)
- Cleanly separate "geometry definition" from "discretized geometry" (`CellGeometry`)

**Interface:**
```csharp
public interface IGeometryDescriptor
{
    /// <summary>Build discretized geometry for given grid resolution</summary>
    CellGeometry Build(GridSize size, float dx, float dy, float dz);
    
    /// <summary>Optional: get underlying SDF if available</summary>
    Func<float, float, float, float>? GetSdf();
}
```

**Implementations:**
- `SdfGeometryDescriptor`: Wraps a primitive SDF (sphere, box, plane, etc.)
- Planned, not yet implemented: `CompositeGeometryDescriptor` (combine multiple SDFs via union/intersection/difference)
- Future: `StlGeometryDescriptor` (mesh-based, requires SDF reconstruction or alternate coarsening strategy). STL meshes currently enter through `MeshSdfGpu` / `SdfBuilder` rather than a dedicated descriptor.

### 2.3 Multigrid Scope Limitation

**Decision:** Only support MG+cut-cells for geometries backed by `IGeometryDescriptor` (SDF-based).

**For geometries without descriptor:**
- Fall back to finest-level-only (no coarsening), or
- Disable MG entirely, use Jacobi-PCG

**Rationale:**
- Keeps implementation honest and maintainable
- STL/precomputed geometries require different coarsening strategies (deferred to future work)
- SDF-based geometries cover most validation cases (sphere, box, channel obstacles)

### 2.4 Rediscretization Strategy

**Decision:** Use operator rediscretization at each level (not Galerkin coarsening).

**Approach:**
- Each coarse level: call `CutCellBuilder.Build()` with original SDF
- Coarse operator computed from coarse geometry metrics (matrix-free)
- No need to assemble `R·A·P` or store coarse matrices

**Benefits:**
- Consistent with matrix-free design
- Preserves geometric semantics
- Simpler implementation than Galerkin

### 2.5 Solid Cell Handling in Operators

**Critical detail:** The pressure operator and divergence operator must explicitly handle solid cells.

**CPU Implementation:**
- `DivGrad.Divergence()`: Checks `CellClassification[idx] == CellClass.Solid`, sets `div = 0` and skips computation
- `PressureOperator.Apply()`: Checks `CellClassification[idx] == CellClass.Solid`, sets `result = 0` to zero out solid rows
- `PressureOperator.GetDiagonal()`: Returns `0.0f` for solid cells (masked from system)

**GPU Implementation:**
- `DivGradGpu.DivergenceKernel()`: Currently does NOT check for solid cells (relies on zero volume/aperture)
- `PressureOperatorGpu.LaplacianKernel()`: Currently does NOT check for solid cells

**Consequence:** Solid cells naturally get zero contribution via `Vol=0` and `Alpha=0`, but explicit masking is more robust.

**Design decision:** Add explicit solid cell masking to GPU kernels for consistency and numerical safety.

### 2.6 Cut-Face Boundary Conditions

**Key insight:** Cut-face projection is **separate** from the pressure solve and multigrid.

**Where it happens:**
- In `TimeIntegrator.AdvanceOneStep()`:
  1. **After diffusion/sources** (line 127): `CutFaceProjector.ProjectCutFaces()` enforces interior obstacle BCs on velocity
  2. Velocity field enters pressure projection with BCs already satisfied
  3. Pressure solve does NOT need to know about cut-face BCs (already baked into velocity)

**What it does:**
- **NoSlip mode**: Sets velocity to zero at faces with `alpha < 0.999` or adjacent to solid cells
- **FreeSlip mode**: Removes normal component (u·n = 0), preserves tangential

**Multigrid implication:**
- MG operates only on pressure (scalar field)
- No velocity projection needed in MG cycle
- Solid cells handled via zero volume/area/aperture in operator
- Cut-face BCs are "pre-applied" to velocity before pressure solve

---

## 3. Architecture Changes

### 3.1 New Components  

#### `IGeometryDescriptor` (Core)
**Location:** `src/PolyCfd.Core/Geometry/IGeometryDescriptor.cs`

```csharp
public interface IGeometryDescriptor
{
    CellGeometry Build(GridSize size, float dx, float dy, float dz);
    Func<float, float, float, float>? GetSdf();
}
```

#### `SdfGeometryDescriptor` (Core)
**Location:** `src/PolyCfd.Core/Geometry/SdfGeometryDescriptor.cs`

```csharp
public sealed class SdfGeometryDescriptor : IGeometryDescriptor
{
    private readonly Func<float, float, float, float> _sdf;
    private readonly CutCellBuilder.BuildOptions _options;
    
    public SdfGeometryDescriptor(
        Func<float, float, float, float> sdf,
        CutCellBuilder.BuildOptions? options = null)
    {
        _sdf = sdf;
        _options = options ?? new CutCellBuilder.BuildOptions();
    }
    
    public CellGeometry Build(GridSize size, float dx, float dy, float dz)
    {
        return CutCellBuilder.Build(size, dx, dy, dz, _sdf, _options);
    }
    
    public Func<float, float, float, float>? GetSdf() => _sdf;
}
```

#### `CompositeGeometryDescriptor` (Core) — planned, not yet implemented
**Intended location:** `src/PolyCfd.Core/Geometry/CompositeGeometryDescriptor.cs`

Combines multiple SDFs using CSG operations (union, intersection, difference).

### 3.2 Modified Components

#### `GridBlock` (Core)
**Changes:**
- Add optional `IGeometryDescriptor? GeometryDescriptor { get; }` property
- Constructor accepts descriptor instead of pre-built `CellGeometry`
- Build geometry lazily or on-demand

**Signature:**
```csharp
public GridBlock(
    GridSize size, 
    float dx, float dy, float dz, 
    IGeometryDescriptor? geometryDescriptor = null)
{
    // Build geometry if descriptor provided
    Geometry = geometryDescriptor?.Build(size, dx, dy, dz);
    GeometryDescriptor = geometryDescriptor;
    
    // ... rest of initialization
}
```

**Backward compatibility:** Existing code that creates `GridBlock` with no geometry continues to work (trivial uniform geometry).

#### `GpuMgHierarchy` (GPU)
**Changes:**
- Constructor accepts `IGeometryDescriptor?` from finest block
- `BuildCoarseBlock()` uses descriptor to rebuild geometry at coarse resolution

**Modified logic:**
```csharp
private GridBlock BuildCoarseBlock(
    GridBlock fineBlock, 
    IGeometryDescriptor? geometryDescriptor,
    int coarseNx, int coarseNy, int coarseNz)
{
    var coarseSize = new GridSize(coarseNx, coarseNy, coarseNz);
    float coarseDx = fineBlock.Metrics.Dx * fineBlock.Size.Nx / coarseNx;
    float coarseDy = fineBlock.Metrics.Dy * fineBlock.Size.Ny / coarseNy;
    float coarseDz = fineBlock.Metrics.Dz * fineBlock.Size.Nz / coarseNz;
    
    // Rebuild geometry at coarse resolution if descriptor available
    return new GridBlock(coarseSize, coarseDx, coarseDy, coarseDz, geometryDescriptor);
}
```

#### `GpuMgLevel` (GPU)
**Changes:**
- Store `DeviceBuffer<byte> CellClass` for cell classifications
- Used by restriction/prolongation to skip solid cells

#### `MultigridKernels` (GPU)
**Changes:**
- Add `cellClass` parameter to `Restrict()` and `Prolongate()`
- Skip solid cells in restriction averaging
- Skip solid cells in prolongation injection

#### `DivGradGpu` (GPU) - Solid Cell Masking
**Changes:**
- Add explicit solid cell check in `DivergenceKernel()`
- Set `div = 0` for solid cells
- Requires passing `CellClass` buffer to kernel

#### `PressureOperatorGpu` (GPU) - Solid Cell Masking
**Changes:**
- Add explicit solid cell check in `LaplacianKernel()`
- Set `ApInterior = 0` for solid cells
- Requires passing `CellClass` buffer to kernel

---

## 4. Implementation Plan

### Phase 1: Core Abstraction (Foundation)
**Status:** ✅ **Completed**

**Tasks:**
1. ✅ Design doc created
2. ✅ Create `IGeometryDescriptor` interface
3. ✅ Implement `SdfGeometryDescriptor`
4. ✅ Add `GeometryDescriptor` property to `GridBlock`
5. ✅ Update `GridBlock` constructor to accept descriptor

**Estimated effort:** 2-3 hours

**Validation:**
- Unit test: `SdfGeometryDescriptor` builds geometry identical to direct `CutCellBuilder.Build()`
- Ensure `GridBlock` backward compatibility (no descriptor = trivial geometry)

### Phase 2: Multigrid Coarse Geometry (Core Feature)
**Status:** ✅ **Completed**

**Tasks:**
1. ✅ Modify `GpuMgHierarchy.BuildCoarseBlock()` to use `IGeometryDescriptor`
2. ✅ Thread descriptor through hierarchy constructor
3. ✅ Add `DeviceBuffer<byte> CellClass` to `GpuMgLevel`
4. ✅ Upload cell classifications to device for each level

**Estimated effort:** 3-4 hours

**Validation:**
- Unit test: 16³ → 8³ → 4³ hierarchy with sphere obstacle
  - Verify each level has correct geometry (cell classifications, apertures)
  - Check coarse sphere remains centered and circular (not averaged/distorted)

### Phase 3: Geometry-Aware Transfer Ops (Kernel Updates)
**Status:** ✅ **Completed** (Solid cell masking in operators deferred)

**Tasks:**
1. ✅ Add `cellClass` parameter to `RestrictionKernel`
2. ✅ Skip solid cells in restriction averaging
3. ✅ Add `cellClass` parameter to `ProlongationKernel`
4. ✅ Skip solid cells in prolongation injection
5. ✅ Update `MultigridKernels.Restrict()` and `Prolongate()` signatures
6. ⏳ Add solid cell masking to `DivGradGpu.DivergenceKernel()` (deferred - implicit masking sufficient for now)
7. ⏳ Add solid cell masking to `PressureOperatorGpu.LaplacianKernel()` (deferred - implicit masking sufficient for now)

**Estimated effort:** 3-4 hours (updated for operator masking)

**Notes:** 
- Explicit solid cell masking in `DivGradGpu` and `PressureOperatorGpu` deferred to Phase 4 after initial testing
- Current implicit masking via zero volume/aperture should be sufficient for initial validation
- If numerical issues arise during testing, explicit masking can be added

**Validation:**
- Kernel test: Restrict residual from 8³ to 4³ with central solid block
  - Verify solid cells ignored in averaging
  - Check coarse residual only reflects fluid cells
- Kernel test: Prolongate correction from 4³ to 8³
  - Verify solid cells receive no correction

### Phase 4: Integration & Testing (Full MG Cycle)
**Historical planning status:** originally not started; implementation and tests
subsequently completed (see the current notes and dated summary). The checklist
below preserves the original proposed tasks.

**Tasks:**
1. [ ] Update `PcgMgSolverGpu.VCycle()` to pass cell classifications
2. [ ] Create validation case: channel with planar obstacle
3. [ ] Run GPU MG solver with cut-cell geometry
4. [ ] Compare convergence vs Jacobi-PCG

**Estimated effort:** 4-5 hours

**Validation:**
- Integration test: 16³ channel with thin plane at x=0.5
  - MG should converge (no crash, finite iterations)
  - Solution respects geometry (zero pressure gradient across solid)
- Performance test: 32³ sphere in cavity
  - MG converges faster than Jacobi-PCG
  - Final divergence ≈ same as Jacobi-PCG (correctness check)

### Phase 5: Composite Geometries (Enhancement)
**Status:** Not Started (Optional)

**Tasks:**
1. [ ] Implement `CompositeGeometryDescriptor` with CSG operations
2. [ ] Add SDF union/intersection/difference helpers
3. [ ] Test with complex geometries (multiple spheres, box + sphere, etc.)

**Estimated effort:** 3-4 hours

**Validation:**
- Validation case: two spheres at different locations
- Validation case: hollow box (difference of two boxes)

---

## 5. Design Decisions Log

### Decision 1: SDF Re-evaluation vs Averaging (Nov 17, 2025)
**Question:** Should coarse geometry be built by averaging fine-level metrics or re-evaluating the SDF?

**Chosen:** SDF re-evaluation (rediscretization)

**Rationale:**
- Averaging distorts geometry (sphere becomes blocky aggregate)
- Apertures computed from averaged face SDFs ≠ apertures from actual SDF at coarse faces
- Rediscretization is standard for embedded-boundary MG
- Simpler implementation (reuse `CutCellBuilder.Build()`)

**Alternatives considered:**
- Volume-weighted averaging: rejected due to geometric incorrectness
- Galerkin coarsening (`R·A·P`): rejected due to matrix-free design, complexity

### Decision 2: Geometry Abstraction (Nov 17, 2025)
**Question:** How to enable geometry rebuilding at multiple resolutions?

**Chosen:** `IGeometryDescriptor` interface

**Rationale:**
- Decouples geometry definition from discretization
- Enables multi-resolution support for MG
- Extensible to non-SDF geometries (future STL, meshes)
- Clean separation of concerns

**Alternatives considered:**
- Store `Func<float,float,float,float>` on `GridBlock`: too rigid, couples to SDF
- Store pre-built `CellGeometry` only: prevents rebuilding at coarse resolutions
- Pass SDF through MG hierarchy directly: tight coupling, not extensible

### Decision 3: MG Scope for Non-SDF Geometries (Nov 17, 2025)
**Question:** What to do with STL/precomputed geometries that lack an SDF?

**Chosen:** Disable MG coarsening for geometries without `IGeometryDescriptor`

**Options:**
1. Finest-level-only MG (no coarsening, just smoothing)
2. Fall back to Jacobi-PCG
3. Require explicit coarse geometry stack (user-provided)

**Chosen:** Option 2 (fall back to Jacobi-PCG) for simplicity

**Future work:** SDF reconstruction from STL meshes, or mesh coarsening strategies

---

## 6. Testing Strategy

### Unit Tests
**Location:** `test/PolyCfd.Tests/Unit/`

1. **SdfGeometryDescriptor**:
   - Build geometry at 8³, verify vs direct `CutCellBuilder.Build()`
   - Rebuild same descriptor at 16³, verify resolution independence

2. **Coarse Geometry Correctness**:
   - Sphere at center: 16³ → 8³ → 4³
   - Verify cell classifications match actual SDF
   - Verify apertures at coarse faces ≈ analytical sphere apertures

3. **Restriction with Solid Cells**:
   - 8³ grid with 2³ solid block at center
   - Random residual in fluid cells, zero in solid
   - Restrict to 4³, verify only fluid cells contribute

4. **Prolongation with Solid Cells**:
   - 4³ coarse correction (non-zero in fluid)
   - Prolongate to 8³, verify solid cells unchanged

### Integration Tests
**Location:** `test/PolyCfd.Gpu.Tests/Multigrid/`

1. **MG Convergence with Cut-Cells**:
   - 16³ grid with sphere (radius 0.15)
   - Random RHS, zero initial guess
   - PCG-MG should converge to same tolerance as Jacobi-PCG
   - Check iteration count (MG should be fewer)

2. **Channel with Obstacle**:
   - 24×16×16 channel with vertical plane at x=0.5
   - Pressure driven flow setup
   - Verify pressure field is smooth, no penetration

### Validation Tests
**Location:** `validation/PolyCfd.Validation/`

1. **Cavity with Sphere**:
   - 32³ lid-driven cavity with sphere obstacle at center
   - GPU backend with `--mg` flag
   - Compare vs CPU solver, check divergence and velocity profiles

2. **Taylor-Green with Cut-Cells** (future):
   - Add obstacles to TG vortex test
   - Verify vortex evolution with embedded boundaries

---

## 7. Performance Expectations

### Baseline (Uniform Grid, no cut-cells)
- 16³: MG converges in ~7 iterations vs Jacobi ~39 iterations (5.5× speedup)
- 32³: MG converges in ~8 iterations vs Jacobi ~60+ iterations (7.5× speedup)

### Expected with Cut-Cells
- **Iteration count:** Similar to uniform grid (7-10 iterations for MG)
- **Per-iteration cost:** Slightly higher due to:
  - Coarse geometry build (one-time, negligible amortized)
  - Cell classification checks in restriction/prolongation (minimal branching)
- **Overall speedup vs Jacobi-PCG:** 4-6× (slightly lower than uniform due to geometry complexity)

### Degradation Scenarios
- Very fine cut-cells (many slivers): coarsening may create more cut cells at coarse levels
- High solid fraction (>50%): fewer fluid cells at coarse levels, potential for premature coarsening termination

**Mitigation:**
- Sliver agglomeration at finest level reduces pathological coarse cells
- Stop coarsening when fluid cell count drops below threshold (e.g., 64 cells)

---

## 8. Known Limitations

### 8.1 STL/Mesh Geometries
STL meshes are supported through SDF reconstruction: `SdfGeometryDescriptor.FromStl`
for static geometry and GPU cached SDFs for rotation. A dedicated
`StlGeometryDescriptor` class is unnecessary for those paths. This does not provide
arbitrary non-SDF mesh coarsening.

### 8.2 Moving Geometries
Rotating geometry is implemented. The integrator updates projection and multigrid
metrics using the end-of-step rotated SDF; coarse geometry and diagonals are
rebuilt for that shape. Connectivity may change, so pressure-component
compatibility is refreshed as part of moving projection. Rebuild cost is included
in solver-step performance measurements. Moving viscous-wall coefficients remain
separate work.

### 8.3 Anisotropic Coarsening
**Issue:** 2× coarsening in all directions may be suboptimal for thin channels with obstacles.

**Current behavior:** Semi-coarsening not implemented.

**Future solution:** Directional coarsening based on geometry aspect ratio.

---

## 9. Critical Implementation Details

### 9.1 Solid Cell Masking in Operators

**CPU Behavior (reference implementation):**

```csharp
// DivGrad.Divergence()
if (geom != null && geom.CellClassification[cellIdx] == CellClass.Solid)
{
    div[cellIdx] = 0.0f;
    continue;  // Skip computation entirely
}

// PressureOperator.Apply()
if (_block.Geometry != null && _block.Geometry.CellClassification[idx] == Geometry.CellClass.Solid)
{
    v = 0.0f; // solid rows inert
}
yInterior[idx] = v;
```

**GPU Implementation Required:**

```csharp
// In DivGradGpu.DivergenceKernel
if (cellClass != null && cellClass[cellIdx] == 1) // 1 = CellClass.Solid
{
    divInterior[cellIdx] = 0.0f;
    return;
}

// In PressureOperatorGpu.LaplacianKernel
float divergence = /* ... compute ... */;
if (cellClass != null && cellClass[cellIdx] == 1)
{
    ApInterior[cellIdx] = 0.0f;
}
else
{
    ApInterior[cellIdx] = -divergence;
}
```

**Rationale:**
- Explicit masking is more robust than relying on zero volume/aperture
- Prevents numerical issues from near-zero denominators
- Matches CPU reference implementation behavior
- Required for multigrid correctness (solid cells must not participate in residual)

### 9.2 Cut-Face Projection is Not Part of MG

**Common misconception:** Cut-face projection happens during pressure solve.

**Actual flow:**
1. Time integrator advances velocity via advection/diffusion
2. **Cut-face projection applied to velocity** (sets u=0 at interior obstacles)
3. Velocity enters pressure projection step
4. Pressure solve (with MG) operates on scalar pressure field
5. Gradient of pressure corrects velocity
6. Velocity BCs re-applied (domain boundaries only)

**Implication for MG:**
- MG only needs to handle pressure field (no velocity)
- Solid cells are passive (zero contribution via operator masking)
- No special "cut-face BC" in MG cycle itself
- Cut-face effects enter via divergence of velocity (already projected)

### 9.3 Volume/Aperture Zero Handling

**Current approach (implicit):**
- Solid cells have `Vol = 0`, `Alpha = 0`, `Area = 0`
- Operators naturally get zero contribution from solid cells
- Works but can create numerical instabilities (0/0 divisions)

**Enhanced approach (explicit + implicit):**
- Keep zero geometry metrics (implicit masking)
- **Add explicit classification checks** (robust masking)
- Best of both worlds: geometrically correct + numerically safe

### 9.4 Coarse Geometry Consistency

**Critical requirement:** Coarse geometry must come from **same SDF** as fine geometry.

**Example - 16³ → 8³ coarsening with sphere (radius 0.15):**

**Correct:**
```csharp
// Fine level
var fineGeom = CutCellBuilder.Build(fineSize, dx, dy, dz, sphereSDF, options);

// Coarse level
var coarseGeom = CutCellBuilder.Build(coarseSize, 2*dx, 2*dy, 2*dz, sphereSDF, options);
// ^ Same SDF function, different grid spacing
```

**Incorrect (don't do this):**
```csharp
// Averaging fine geometry - WRONG!
var coarseGeom = AverageFineCellVolumes(fineGeom, ...);  // Distorts sphere
```

**Why correct approach works:**
- Sphere equation unchanged: `|r| - 0.15`
- Coarse grid samples sphere at different locations
- Apertures computed from actual SDF at coarse face centers
- Normals computed from actual SDF gradients at coarse cells
- Geometry remains faithful to analytical obstacle

---

## 10. Open Questions

### Q1: Coarsening Strategy for Complex Geometries
**Context:** Multiple disconnected obstacles, thin features.

**Options:**
1. Always coarsen 2× in all directions (current)
2. Adaptive coarsening: stop if geometry becomes degenerate
3. Selective coarsening: coarsen only in directions without thin features

**Decision needed:** After testing with validation cases.

### Q2: Nullspace Handling on Coarse Levels
**Context:** Finest level projects out mean for periodic BCs. What about coarse levels?

**Resolved (2026-09-25):** only the finest-level PCG residual is projected, per connected
region and after every update; coarse levels project nothing. See the current
implementation notes above.

### Q3: Should GPU Operators Always Check CellClass?
**Context:** CPU operators always check solid cells. GPU could rely on zero metrics alone.

**Options:**
1. Always check `CellClass` buffer (safe, explicit)
2. Only check when geometry present (optimization)
3. Never check, rely on zero metrics (current, risky)

**Recommendation:** Option 2 - check when geometry present, skip check for uniform grids.

### Q4: CellClass Buffer Size at Each Level
**Context:** Each MG level needs cell classifications on device.

**Options:**
1. Upload full array (Nx×Ny×Nz bytes per level)
2. Use bitpacked representation (Nx×Ny×Nz bits / 8)
3. Sparse indexing (only store cut/solid cell indices)

**Recommendation:** Option 1 for simplicity (memory cost is minimal: 64³ = 256KB).

---

## 11. Next Steps (Priority Order)

1. ✅ **Complete design doc** (this document)
2. ✅ **Review CPU implementation** (cut-face projection, solid cell masking)
3. ✅ **Implement Phase 1** (Core Abstraction):
   - ✅ `IGeometryDescriptor` interface
   - ✅ `SdfGeometryDescriptor` implementation
   - ✅ Update `GridBlock` constructor
4. ✅ **Implement Phase 2** (Coarse Geometry):
   - ✅ Modify `GpuMgHierarchy.BuildCoarseBlock()`
   - ✅ Add cell classification device buffers
5. ✅ **Implement Phase 3** (Kernel Updates):
   - ✅ Geometry-aware restriction/prolongation
   - ⏳ Operator solid cell masking (deferred to Phase 4)
6. [ ] **Implement Phase 4** (Integration & Testing):
   - [ ] Modify `ChannelWithObstacle.cs` to use `SdfGeometryDescriptor` when creating GridBlock
   - [ ] Run validation case with GPU + MG + cut-cell geometry
   - [ ] Compare convergence vs Jacobi-PCG
   - [ ] If numerical issues occur, add explicit solid cell masking to operators
7. [ ] **Update Documentation**:
   - [ ] Add usage examples to `API.md`
   - [ ] Update `GEOMETRY_PROGRESS.md` with MG status

**Implementation Notes (November 17, 2025):**
- Phases 1-3 completed successfully
- All code compiles without errors
- Solid cell masking in `DivGradGpu` and `PressureOperatorGpu` deferred to Phase 4
  - Current implicit masking via zero volume/aperture should be sufficient
  - If numerical stability issues arise during testing, explicit masking can be added
- Next step: Create test case that uses `SdfGeometryDescriptor` to enable MG with cut-cells

---

## 11. References

- **Embedded Boundary Methods:**
  - Johansen & Colella (1998), "A Cartesian Grid Embedded Boundary Method"
  - Schwartz et al. (2006), "A Cartesian Grid Multigrid Method for the Poisson Equation with Variable Coefficients"
  - Chombo library documentation (LBL)

- **Related PolyCFD Docs:**
  - `GPU_MULTIGRID_DESIGN.md` - GPU MG implementation for uniform grids
  - `GEOMETRY_PROGRESS.md` - Cut-cell geometry implementation status
  - `design.md` - Overall architecture and solver design

---

**Document Version:** 2.0  
**Last Updated:** November 17, 2025  
**Status:** ✅ **Implementation Complete - All Tests Passing**

## Implementation Summary

### Completed Work (November 17, 2025)

**Phase 1-3 Implementation:**
- ✅ Core geometry descriptor abstraction (`IGeometryDescriptor`, `SdfGeometryDescriptor`)
- ✅ Multi-resolution geometry support in `GridBlock` and `GpuMgHierarchy`
- ✅ Geometry-aware multigrid transfer operators (restriction/prolongation)
- ✅ Explicit solid cell masking in `PressureOperatorGpu.LaplacianKernel()`
- ✅ Cell classification upload to GPU for all MG levels

**Phase 4 Testing:**
- ✅ 6 unit tests created and passing (`GeometryDescriptorTests`)
  - Geometry descriptor correctness
  - Multi-resolution rebuilding
  - GridBlock integration
- ✅ 3 integration tests created and passing (`CutCellMultigridTests`)
  - Sphere obstacle convergence
  - Box obstacle: 1.5-3x MG speedup vs Jacobi
  - Multi-level geometry preservation (32³→16³→8³→4³)

**Key Findings:**
1. **Explicit operator masking required**: Initial tests failed with NaN due to implicit masking. Adding explicit `cellClass` check in `PressureOperatorGpu.LaplacianKernel()` resolved the issue.
2. **Transfer operators work correctly**: Restriction and prolongation properly skip solid cells.
3. **MG efficiency maintained**: Even with cut-cells, MG provides 1.5-3x iteration reduction vs Jacobi-PCG.
4. **Geometric fidelity preserved**: SDF re-evaluation at coarse levels maintains obstacle shape correctly.

**Files Created/Modified:**
- `src/PolyCfd.Core/Geometry/IGeometryDescriptor.cs` (new)
- `src/PolyCfd.Core/Geometry/SdfGeometryDescriptor.cs` (new)
- `src/PolyCfd.Core/Core/GridBlock.cs` (modified)
- `src/PolyCfd.Gpu/Multigrid/GpuMgHierarchy.cs` (modified)
- `src/PolyCfd.Gpu/Multigrid/GpuMgLevel.cs` (modified)
- `src/PolyCfd.Gpu/Multigrid/MultigridKernels.cs` (modified)
- `src/PolyCfd.Gpu/PressureOperatorGpu.cs` (modified)
- `src/PolyCfd.Gpu/PcgMgSolverGpu.cs` (modified)
- `test/PolyCfd.Tests/Unit/GeometryDescriptorTests.cs` (new)
- `test/PolyCfd.Gpu.Tests/Multigrid/CutCellMultigridTests.cs` (new)

**Next Steps for Production Use:**
1. Update validation cases (e.g., `ChannelWithObstacle.cs`) to use `SdfGeometryDescriptor` for MG support
2. Add performance benchmarks comparing CPU vs GPU MG with cut-cells
3. Test with more complex geometries (multiple obstacles, thin features)
4. Consider adding explicit masking to `DivGradGpu` if needed for time integration stability