← Explainer · Markdown · Source on GitHub
# Boundary Conditions and Halo/Ghost Cells
This document describes how boundary conditions interact with halo/ghost cells in the PolyCfd staggered-grid MAC discretization, and when BCs are applied in the overall time-stepping cycle.
**Note:** As of December 2025, the solver uses ILGPU (CUDA or the supported CPU accelerator). Velocity halos have been removed from `DeviceMacState`; pressure halos remain for pressure operators, gradients, and multigrid. A halo-only `PressureWorkspace` borrows the pressure interior, and immutable pressure BC descriptors are shared across levels. The integrator owns stage BC application. References to `BcRegistry`, `BcType`, and `BcDescriptor` below describe the removed CPU/legacy path and are kept for context; the current types are `PatchBcSet` (host) and `DevicePatchBcSet` (device).
## 1. Grid Layout and Halos
### 1.1 Staggered MAC grid
We use a standard MAC (Marker-And-Cell) staggered layout on a logical grid of size $(N_x, N_y, N_z)$ cells:
- **Pressure / scalar cells `p`**: cell-centered, indices $(i, j, k)$ with
- $i \in [0, N_x-1]$, $j \in [0, N_y-1]$, $k \in [0, N_z-1]$.
- **Velocity components `(u, v, w)`**:
- `u` at **x-faces** (normal to x): $(i+\tfrac12, j, k)$
- stored as `UInterior`: size $(N_x+1) \times N_y \times N_z$, indexed by `iHalf` in $[0, N_x]$.
- `v` at **y-faces** (normal to y): $(i, j+\tfrac12, k)$
- stored as `VInterior`: size $N_x \times (N_y+1) \times N_z$, indexed by `jHalf` in $[0, N_y]$.
- `w` at **z-faces** (normal to z): $(i, j, k+\tfrac12)$
- stored as `WInterior`: size $N_x \times N_y \times (N_z+1)$, indexed by `kHalf` in $[0, N_z]$.
### 1.2 Halo/ghost cells (Historical - Velocity Halos Removed)
> **IMPORTANT:** As of December 2025, velocity halos (`UHalo`, `VHalo`, `WHalo`) have been **removed** from `DeviceMacState`. Only pressure halos (`PHalo`) remain on the GPU, for pressure operators, gradients, and multigrid. Host `MacState` still retains its velocity halo arrays. The section below describes the historical design.
For pressure fields we maintain a **packed halo buffer** that stores ghost values on the six boundary faces:
- Faces: `XMin`, `XMax`, `YMin`, `YMax`, `ZMin`, `ZMax`.
- For each field, the halo faces have natural 2D extents matching that field’s **face plane** at the domain boundary, e.g.
- `U` halos on `XMin`/`XMax`: size `Ny × Nz` (one value per `(j, k)` at a fixed `iHalf`).
- `U` halos on `YMin`/`YMax`: size `(Nx+1) × Nz` (one value per `(iHalf, k)` at a fixed `j`).
- `U` halos on `ZMin`/`ZMax`: size `(Nx+1) × Ny` (one value per `(iHalf, j)` at a fixed `k`).
- analogous layouts for `V`, `W`, and `P`.
Halo buffers are packed linearly using the `HaloOffsets` helper. Specific face accessors like `GetXMinHaloIndex`, `GetYMinHaloIndex` etc. map logical halo coordinates to packed indices.
**Key principle**: interior arrays represent physical DOFs; halo arrays are purely auxiliary storage for stencil operations and BC enforcement.
## 2. Boundary Condition Types
We support several BC types via the `BcType` / `BcDescriptor` abstraction. Conceptually:
- **Periodic** (fully periodic domain)
- **NoSlipVelocity** (no-slip wall)
- **MovingWall** (prescribed wall velocity, generalizing inflow/outflow walls)
- **FreeSlipVelocity** (slip / symmetry wall)
- **NeumannZero** (zero normal gradient for scalar/velocity component)
- **Dirichlet** (fixed scalar value, e.g. pressure reference)
The CPU BC handlers (`NoSlipWallBc`, `FreeSlipBc`, `InflowBc`, `OutflowBc`, `PeriodicBc`) are the reference model; GPU kernels in `BoundaryKernels` must match their semantics.
## 3. What Each BC Does (Conceptually)
### 3.1 Periodic
- **Pressure** `p` (cell-centered):
- values on one side of the domain wrap to the opposite side.
- Halo on `XMin` is filled from interior `XMax-1` cell, `XMax` halo from `0`, similarly in `Y` and `Z`.
- **Velocities**:
- For each component, the ghost face at one side wraps to the interior face at the opposite side in the normal direction.
- Interior boundary values are usually also adjusted to enforce exact periodicity (e.g. copy `u` at `iHalf=0` from `iHalf=Nx`, etc.), so that the velocity field is globally periodic.
### 3.2 No-slip wall
- Represents a solid wall with **zero velocity at the wall**.
- **Normal velocity** at the wall is zero.
- **Tangential velocity** at the wall is also zero; ghosts mirror the interior with opposite sign to support centered stencils.
Consequences:
- Interior velocity on the wall faces is set to zero for the component normal to the wall.
- Ghosts for tangential components reflect the interior with a sign flip.
- Pressure typically uses a **Neumann** condition (zero normal gradient) on solid walls.
### 3.3 Moving wall / inflow
- Prescribes a specific velocity $(u_w, v_w, w_w)$ at a boundary (e.g., lid-driven cavity or inflow boundary).
- **Normal component** equals the prescribed wall/inflow normal velocity.
- **Tangential components** equal the prescribed tangential components.
- Ghosts are set so that the wall value is enforced in the centered stencil (usually by reflecting around the wall value).
Pressure typically uses **Neumann** (zero normal gradient) for inflow or solid moving walls; **Dirichlet** may be used on outflow faces.
### 3.4 Free-slip wall / symmetry
- Models a symmetry plane or inviscid wall.
- **Normal velocity** at the wall is zero (no penetration).
- **Tangential velocity** has **zero normal derivative** at the wall (no shear stress in inviscid case, or symmetry plane).
Implementation practice:
- Ghost normal velocity is the negative of the interior value (antisymmetric), which, together with the pressure projection, drives normal velocity to zero at the wall.
- Ghost tangential velocity is copied from the interior (symmetric), giving zero normal derivative.
- Pressure uses a **Neumann** BC (zero normal gradient).
### 3.5 Outflow
- Represents a convective / open boundary where fluid exits the domain.
- **Velocity**: frequently modeled as zero normal gradient (Neumann) for all components, or as a weak condition combined with pressure.
- **Pressure**: often uses a **Dirichlet** reference value (e.g., $p = 0$) at the outflow, which closes the Poisson problem.
In this codebase:
- `OutflowBc` uses **Neumann** for velocity (`NeumannZero`) and **Dirichlet(pRef)`** for pressure.
## 4. Per-face Behavior on the Staggered Grid
The following table summarizes the idealized behavior of BCs for each field and face. Signs and indices are expressed conceptually; CPU implementations in `FreeSlipBc`, `NoSlipWallBc`, `InflowBc`, `OutflowBc`, and `PeriodicBc` are the source of truth.
### 4.1 Legend
- `u_n` / `u_t`: normal / tangential velocity components.
- `ghost`: halo cell on a face, e.g. `u(i=-1, j, k)` at `XMin`.
- `interior`: nearest interior face or cell.
- `\partial_n`: derivative in normal direction.
### 4.2 Pressure BCs
| BC type | Face | Condition on `p` | Halo fill (conceptual) |
|--------------|--------------|------------------------------------|----------------------------------------------------------------|
| Periodic | all | $p$ wraps across domain | Ghost from opposite-side interior cell |
| NeumannZero | wall, inflow | $\partial_n p = 0$ | Ghost equals nearest interior cell |
| Dirichlet | outflow | $p = p_\text{ref}$ at boundary | Ghost set so that interpolation yields $p_\text{ref}$ (often just constant in halo) |
In practice in this project:
- `PeriodicP`: halos are filled from opposite-side interior cells.
- `NeumannP`: halos copy adjacent interior cell.
- `DirichletP`: halos are set to the specified value.
For this project, **pressure BC handlers only ever write to halos**. Interior pressure cells, including those adjacent to boundaries, are treated as unknowns of the Poisson solve and are not explicitly clamped by BC handlers.
### 4.3 Velocity BCs (per face)
Below we describe the intended behavior for the component normal to the face and the two tangential components.
#### 4.3.1 X faces (`XMin`, `XMax`)
Normal component: `u` (x-face velocities)
| BC type | Face | Normal `u` at wall | `u` ghost | Tangential `v, w` ghosts |
|------------------|---------|----------------------------------------|-----------------------------------|--------------------------------------------------------|
| Periodic | XMin/XMax | Periodic | Wrap from opposite side | Wrap from opposite side |
| NoSlipVelocity | XMin | `u(0, j, k) = 0` | `u(-1, j, k) = -u(0, j, k)` | `v(-1, j+1/2, k) = -v(0, j+1/2, k)` and similarly for `w` |
| NoSlipVelocity | XMax | `u(Nx, j, k) = 0` | `u(Nx+1, j, k) = -u(Nx, j, k)` | tangential mirrored with sign flip |
| MovingWall | XMin | `u(0, j, k) = u_w` | `u(-1, j, k) = 2u_w - u(0, j, k)` | tangential ghosts reflect around wall velocity |
| FreeSlipVelocity | XMin | no-penetration (enforced by projection) | `u(-1, j, k) = -u(0, j, k)` | `v, w` ghosts copy interior (zero normal gradient) |
| NeumannZero | XMin/XMax | $\partial_n u = 0$ (e.g. outflow) | `u_\text{ghost} = u_\text{interior}` | tangential usually also zero-gradient or left unchanged |
#### 4.3.2 Y faces (`YMin`, `YMax`)
Normal component: `v` (y-face velocities)
| BC type | Face | Normal `v` at wall | `v` ghost | Tangential `u, w` ghosts |
|------------------|---------|----------------------------------------|-----------------------------------|--------------------------------------------------------|
| Periodic | YMin/YMax | Periodic | Wrap from opposite side | Wrap from opposite side |
| NoSlipVelocity | YMin | `v(i, 0, k) = 0` | `v(i, -1, k) = -v(i, 0, k)` | `u(i+1/2, -1, k)` and `w(i, -1, k+1/2)` mirrored with sign flip |
| NoSlipVelocity | YMax | `v(i, Ny, k) = 0` | `v(i, Ny+1, k) = -v(i, Ny, k)` | tangential mirrored with sign flip |
| MovingWall | YMin/YMax | `v` set to wall normal velocity | `v_\text{ghost}` reflected around wall velocity | tangential extrapolated or mirrored |
| FreeSlipVelocity | YMin | no-penetration | `v(i, -1, k) = -v(i, 0, k)` | `u` and `w` halos copy interior (zero normal gradient) |
| NeumannZero | YMin/YMax | $\partial_n v = 0$ | `v_\text{ghost} = v_\text{interior}` | tangential components typically also Neumann |
#### 4.3.3 Z faces (`ZMin`, `ZMax`)
Normal component: `w` (z-face velocities)
| BC type | Face | Normal `w` at wall | `w` ghost | Tangential `u, v` ghosts |
|------------------|---------|----------------------------------------|-----------------------------------|--------------------------------------------------------|
| Periodic | ZMin/ZMax | Periodic | Wrap from opposite side | Wrap from opposite side |
| NoSlipVelocity | ZMin | `w(i, j, 0) = 0` | `w(i, j, -1) = -w(i, j, 0)` | `u(i+1/2, j, -1)` and `v(i, j+1/2, -1)` mirrored |
| NoSlipVelocity | ZMax | `w(i, j, Nz) = 0` | `w(i, j, Nz+1) = -w(i, j, Nz)` | tangential mirrored with sign flip |
| MovingWall | ZMin/ZMax | `w` set to wall normal velocity | `w_\text{ghost}` reflected around wall velocity | tangential extrapolated or mirrored |
| FreeSlipVelocity | ZMin/ZMax | no-penetration | `w_\text{ghost} = -w_\text{interior}` | tangential components extrapolated |
| NeumannZero | ZMin/ZMax | $\partial_n w = 0$ | `w_\text{ghost} = w_\text{interior}` | tangential components typically also Neumann |
> Note: The exact loops and index ranges are implemented in `FreeSlipBc`, `NoSlipWallBc`, `PeriodicBc`, etc. The GPU `BoundaryKernels` must match those semantics when filling packed halos.
## 5. When BCs Are Applied in the Time Stepping
At a high level, a typical time step for incompressible flow looks like:
1. **Compute advection/diffusion** of velocities and scalars in the interior.
2. **Apply velocity BCs**:
- Use `BcRegistry.ApplyVelocityGhosts` (CPU) or `BoundaryKernels.ApplyVelocityBcs` (GPU) to fill velocity halos based on the current interior fields.
- This ensures that all stencils used by the pressure solve / projection have consistent boundary data.
3. **Assemble and solve pressure Poisson equation**:
- Use divergence of provisional velocity and density, with pressure BCs (Dirichlet on outflow, Neumann on walls/inflow, periodic as appropriate).
- Apply pressure BCs: `BcRegistry.ApplyPressureGhosts` (CPU) or `BoundaryKernels.ApplyPressureBcs` (GPU) to fill pressure halos.
4. **Projection / velocity correction**:
- Update velocities using the pressure gradient to enforce incompressibility.
- After projection, velocities satisfy the normal BCs (no-penetration, outflow, etc.) in a discrete sense.
Practical notes for this codebase:
- **Velocity BCs before Poisson**: halos are filled so that the discretized divergence operator sees the correct boundary fluxes.
- **Pressure BCs during Poisson**: **only pressure halos are written by BC handlers**; interior pressure cells (including boundary-adjacent cells) are left to the Poisson solve.
- **No-slip interior zeroing** (for normal components) is currently done in the BC handlers themselves on the CPU, and mirrored in the GPU kernels. Free-slip BCs, by contrast, only modify halos and rely on the projection step to enforce zero normal velocity.
## 6. Implementation Guidelines
To keep CPU and GPU behavior consistent and easier to reason about:
1. **Halos first-class, interiors stable**
- Prefer BC handlers that treat interior arrays as read-only and write only to halo buffers, except where interior modification is strictly necessary (e.g. no-slip walls for normal components).
- For GPU kernels in `BoundaryKernels`, aim for an API where interior views are logically read-only.
2. **Mirror CPU behavior exactly on GPU**
- For each BC type, follow the CPU handler’s loops and signs closely.
- Pay special attention to:
- Which indices are included (e.g. `iHalf` in `[0..Nx]` vs `[0..Nx-1]`).
- How corners and edges (intersection of two faces) are handled.
- Whether some halo entries are intentionally left at the default value.
3. **Ordering matters when interiors are modified**
- If a BC modifies interior boundary faces (e.g. NoSlip sets `u=0` on `XMin`), order of application relative to other BCs that read those values (e.g. FreeSlip on `YMin`) must be consistent between CPU and GPU.
- In tests that compare CPU vs GPU, register handlers in **face index order** (`XMin..ZMax`) so that `BcRegistry` iteration order matches `BoundaryKernels`.
4. **Tests as specification**
- GPU tests in `PolyCfd.Gpu.Tests` (`test/PolyCfd.Gpu.Tests/BC/PatchBoundaryKernelsTests.cs`, `test/PolyCfd.Gpu.Tests/BC/Kernels/BoundaryKernelComparisonTests.cs`) exercise the patch kernels for a variety of BCs and grid sizes.
- When changing BC logic, update both CPU and GPU consistently, then extend tests to cover the new behavior.
5. **Division of responsibilities (project policy)**
- **Pressure BCs**: handlers populate `P` halos only; they must **not** write interior `P` cells. The discrete operator and Poisson solver enforce pressure conditions using these halo values.
- **Velocity BCs**: handlers may modify interior velocity faces where required by physics (e.g. zeroing normal components for no-slip, prescribing inflow/outflow/moving wall velocities), and must also fill corresponding halos consistently with the CPU reference.
- **Free-slip BCs**: velocity handlers for free-slip/symmetry walls **do not** modify interior velocities; they only write halos for normal and tangential components as described above.
- Any new BC type should explicitly document whether it is allowed to modify interior DOFs or is halo-only.
## 7. Summary Table
The table below summarizes intended behavior at a high level. `u_n` is the component normal to the face, `u_t` are tangential components.
| BC type | Quantity | Wall type / face | Condition at wall | Ghost value (typical) |
|------------------|----------|------------------------|-------------------------------------------|---------------------------------------------------|
| Periodic | p | any | $p$ periodic | ghost = opposite-side interior |
| Periodic | u_n, u_t | any | velocities periodic | ghost = opposite-side interior |
| NoSlipVelocity | u_n | solid wall | $u_n = 0$ | ghost = $-u_n$ (interior) |
| NoSlipVelocity | u_t | solid wall | $u_t = 0$ | ghost = $-u_t$ (interior) |
| MovingWall | u_n,u_t | moving wall/inflow | $u = u_w$ (prescribed) | ghost = $2u_w - u_\text{interior}$ |
| FreeSlipVelocity | u_n | slip/symmetry wall | no penetration (via projection) | ghost = $-u_n$ |
| FreeSlipVelocity | u_t | slip/symmetry wall | $\partial_n u_t = 0$ | ghost = $u_t$ |
| NeumannZero | p | solid, inflow, slip | $\partial_n p = 0$ | ghost = interior |
| Dirichlet | p | outflow / reference | $p = p_\text{ref}$ | ghost = constant value (or mirrored accordingly) |
This document should be kept in sync with:
- `src/PolyCfd.Core/BC/Patches/*.cs` for the patch definitions (`PatchBcSet`, `PatchBuilder`, `PatchBcSetFactory`).
- `src/PolyCfd.Gpu/BC/BoundaryKernels.cs` and `src/PolyCfd.Gpu/BC/Patches/*.cs` for GPU kernels.
- `test/PolyCfd.Gpu.Tests/BC/PatchBoundaryKernelsTests.cs` and `test/PolyCfd.Gpu.Tests/BC/Kernels/BoundaryKernelComparisonTests.cs` for tests.
## 8. BC-Aware Diffusion (Halo-Free Stencil Computation)
The GPU diffusion solver (`DiffusionGpu`) computes boundary neighbor contributions inline in its
one row evaluator, so no velocity halo arrays are needed during diffusion iterations. This section
derives the ghost relations; the implemented operator, method and tests are in the
[Diffusion README](/src/PolyCfd.Gpu/Diffusion/README).
### 8.1 Mathematical Background
For the implicit diffusion equation (backward Euler):
$$
(I - \alpha \nabla^2) u^{n+1} = u^n
$$
where $\alpha = \nu \cdot dt$, the row at a grid point $(i,j,k)$ is:
$$
D\,u^{n+1}_{i,j,k} - \sum_{\text{neighbors}} \frac{\alpha}{h_e^2}\,u^{n+1}_{\text{neighbor}} = u^{n}_{i,j,k}
$$
where $D = 1 + 2\alpha(dx^{-2} + dy^{-2} + dz^{-2})$ is the diagonal coefficient.
At boundaries, the "neighbor" value is a ghost cell that must be computed from the BC.
### 8.2 BC Type Handling
The `VelocityDiffusionBcs` struct captures BC types and values for all 6 domain faces. The row
evaluator (`EvaluateRow`, shared by every diffusion kernel and all three components) branches on
the face type inline. The base diagonal $D$ already counts every neighbor, including the ghost, so a
boundary row only *adjusts* $D$ and the neighbor sum.
The rules below are for **tangential** components, whose DOF sits half a cell ($h/2$) from the
wall. The **normal** component of a holding face (Dirichlet types, slip, symmetry) sits on the wall
itself: it is a prescribed row (residual p − u) whose audited result holds its patch value; a partial
velocity mask holds only its selected faces. How each phase of the method updates it is in the Diffusion README.
#### 8.2.1 Dirichlet BCs (NoSlip, FixedValue, MovingWall)
The wall value is enforced by reflecting around it (linear extrapolation to the wall):
$$
u_{\text{ghost}} = 2\,u_{\text{wall}} - u_{\text{interior}}
$$
NoSlip is the case $u_{\text{wall}} = 0$. Substituting into the stencil moves the
$-u_{\text{interior}}$ term onto the diagonal:
- the neighbor sum gains $2\alpha\,u_{\text{wall}}/h^2$
- the diagonal gains $\alpha/h^2$ (one extra copy of the center; the base $D$ already holds one)
The diagonal adjustment is an $O(1)$ matter, not a discretization detail: with $2\alpha/h^2$
instead of $\alpha/h^2$ (as an earlier version had) the wall-adjacent row converges to
$\tfrac{2}{3}u_{\text{wall}}$ as $h \to 0$, which made the lid-driven cavity about 50% too weak.
`DiffusionWallBcTests` pins the correct coefficient with a Couette profile.
#### 8.2.2 Neumann BCs (ZeroGradient, Slip, SymmetryPlane)
$u_{\text{ghost}} = u_{\text{interior}}$, so the neighbor sum gains nothing and the diagonal
loses $\alpha/h^2$. An earlier version left the diagonal untouched, which made the missing
neighbor behave like a zero-valued ghost: a plug flow next to a slip wall or an outflow face
decayed towards the wall by about $\tfrac{\alpha/h^2}{1 + 6\alpha/h^2}$ per step.
`DiffusionNeumannBcTests` pins the correct behaviour (a uniform field next to a zero-gradient
face is a steady solution).
#### 8.2.3 Periodic BCs
The kernel reads the neighbor from the opposite side of the domain directly; no halo is needed.
#### 8.2.4 Immersed walls (cut-cell geometry)
With `ImmersedWallCoefficients` attached (`--immersed-shear` in the validation suite), the
kernel also puts the no-slip condition at the *true* immersed surface. For a fluid unknown
$u_f$ whose stencil neighbour at distance $h$ lies in or behind the solid, the wall intercept
$d$ along that connection is found by bisection on the SDF and the neighbour is replaced by
the linearly reconstructed ghost
$$u_g = u_w\,\frac{h}{d} - u_f\,\frac{h-d}{d},$$
which is the domain-wall reflection above when $d = h/2$. The row's diagonal gains
$\alpha (h-d)/(h^2 d)$ and its right-hand side $\alpha\,u_w\,h/(h^2 d)$; both are geometry-only
and are built once per geometry. Closed faces (zero aperture) are held at their wall velocity
for the whole solve; a face with a small positive aperture stays an ordinary unknown, because
the projection corrects it every step. Intercepts are searched up to $2h$ so a wall lying
just beyond a blocked neighbour is placed where it is; without a crossing the link falls back
to $d = h$ (the neighbour itself is the wall), and an unknown whose own centre is inside the
solid uses $d = 0.1h$. Tests: `ImmersedWallDiffusionTests` (translated wall exact, tilted wall
exact on interior rows, circular Couette to 0.2 %).
### 8.3 Implementation Details
Files:
- `src/PolyCfd.Gpu/Diffusion/DiffusionBcDescriptor.cs`: `FaceBc` and `VelocityDiffusionBcs` structs
- `src/PolyCfd.Gpu/Diffusion/DiffusionGpu.cs`: the row evaluator and the kernels over it
- `test/PolyCfd.Gpu.Tests/Diffusion/DiffusionWallBcTests.cs`: wall-coefficient regression tests
- `test/PolyCfd.Gpu.Tests/Diffusion/DiffusionOperatorTests.cs`: every fold against an independent host row
Key methods:
- `VelocityDiffusionBcs.ForU/V/W(DevicePatchBcSet)`: Extract BC descriptors from patch BC set
- `EvaluateRow`: the one row evaluator for all three components (`MacComponent` carries the
extents and layout); `Neighbour` folds the domain BC of each stencil link, `ImmersedWallViews`
the immersed-wall coefficients
### 8.4 Advantages
1. **Eliminates halo updates during diffusion iterations**: No need to call `ApplyVelocityBcsHaloOnly`
2. **Reduced memory traffic**: Ghost values computed on-the-fly instead of read from halo arrays
3. **Simpler data flow**: Only interior arrays needed during diffusion solve
### 8.5 Limitations
1. **Moving immersed walls**: the immersed-wall coefficients are built once for static
geometry; the moving-geometry path does not rebuild them per step yet.
2. **Non-uniform BCs per face**: Each face has a single BC type. Complex geometries may need halo approach.