Boundary conditions as data, not methods

Proposal generated as part of full AI review of SolidStateDetectors codebase:

Written up from the 2026-07 review
on the `improvements` branch; line numbers refer to that branch.

## Summary

The SOR boundary-condition (BC) code is ~25 nearly identical Julia methods spread
over three files, one method per combination of axis role and boundary types. Each
method is a one- or two-line ghost-row copy, and the only information that varies is
*which source row* is read and *which factor* is applied. During the review, this
duplication produced one latent bug, one coverage hole, and one unexplained
inconsistency (details below). The GPU port added on the branch already encodes each
face as a small data tuple `(active, source_row, factor)`; the proposal is to make
those tables the single source of truth and have the CPU path consume them too. This
deletes most of the three files, makes the BC rules unit-testable in one place, and
forces the open questions to be answered exactly once.

## Background: how BCs are applied

The SOR solver stores each red-black color as a 4-array with one ghost row on every
face; dimension 1 is the parity-compressed dimension (cylindrical `z`, Cartesian
`x`). After each half-iteration, `apply_boundary_conditions!` rewrites the ghost
rows from interior rows:

* fixed: ghost stays untouched,
* reflecting: ghost mirrors the interior neighbor,
* periodic: ghost wraps to the opposite interior row,
* infinite: ghost = decay factor × outermost interior row, with the factor
precomputed in `grid_boundary_factors`.

Two layout facts drive every index in these methods:

* Non-compressed dimensions have real cells in rows `2 … end−1`; the mirror
neighbor for a reflecting face is row `3` / `end−2`.
* The compressed dimension holds only one parity, so the "one row inside" of the
boundary already is the ±1-neighbor of the other color: reflecting/infinite reads
row `2` / `end−1` there.

## Current state

| file | methods | role |
|---|---|---|
| `BoundaryConditionsCartesian.jl` | 8 z-axis + 9 x-axis + 9 y-axis | broadcast ghost-row copies per (left,right) boundary-type pair |
| `BoundaryConditionsCylindrical.jl` | 2 φ + 3 r + r0 special case | ditto, plus the φ-averaged `r = 0` axis update |
| `BoundaryConditionsGPU.jl` (new, on the branch) | 2 op tables + 1 kernel | data encoding of the same rules |

Notes on the structure:

* The cylindrical `z` axis is handled by calling
`apply_boundary_conditions_on_x_axis!` (`BoundaryConditionsCylindrical.jl:129`) —
the name says "x axis", the semantics are "compressed dimension".
* The compressed-dimension methods carry comments like *"hmm, this is probably not
fully correct since this is the red-black dimension"* and *"anyhow its just an
approximation"* (`BoundaryConditionsCartesian.jl:51–52`) — the duplication makes
it genuinely hard to reason about, and the uncertainty was committed instead of
resolved.
* The `r = 0` axis update (`apply_boundary_conditions_at_r0!`) is not a ghost-row
copy but a φ-weighted average over the first real ring; it is genuinely different
math and should stay a separate function.

## Evidence that the duplication bites

**1. The compressed-dimension periodic wrap was wrong** (fixed in `14ff8e52`).
`(:periodic, :periodic)` on the compressed dimension copied ghost rows into ghost
rows: row 1 read row `end` (a ghost) and row `end` read the just-overwritten row 1.
The correct wrap — matching every other method of the same file — is
`ghost[1] ← interior[end−1]`, `ghost[end] ← interior[2]`. The method sat next to
eight siblings that all did it right for their cases; a shared implementation would
have made the error impossible to write. (Latent in practice because periodic
boundaries on the compressed dimension are rare; now pinned by a red-black layout
unit test in `test/test_electric_field.jl`.)

**2. The Cartesian z-axis method set is incomplete.** `…_on_x_axis!` and
`…_on_y_axis!` each cover all nine (left,right) combinations plus periodic;
`…_on_z_axis!` is missing `(:infinite, :reflecting)` and `(:reflecting, :infinite)`
(`BoundaryConditionsCartesian.jl:1–42`). A Cartesian config using such a z-axis
combination fails with a `MethodError` deep in the solver. Combinatorial method
explosion invites exactly this kind of hole; a table covers the product space by
construction.

**3. The r-axis outer face reads a different row than every other axis.**
`(:r0, :infinite)` and `(:r0, :reflecting)` read `size − 2`
(`BoundaryConditionsCylindrical.jl:20,25`), i.e. the *second-to-last* real row.
Every other non-compressed axis reads `size − 1` for infinite faces and `size − 2`
only for reflecting ones. The decay factor itself is defined as a tick ratio of the
last real tick vs. the ghost tick
(`PotentialCalculationSetupCylindrical.jl:144`: `r_ext[end−1] / r_ext[end]`), which
is only consistent with reading row `end−1`. So the infinite-r ghost currently
extrapolates from one row further in than its factor assumes. The effect is small
(infinite BCs are approximations, and the outer rows differ little on padded
grids), but nobody can tell from the code whether it is intended. The GPU tables
transcribe it 1:1 on purpose (bug-compatible, `BoundaryConditionsGPU.jl:58–63`)
so that CPU and GPU results stay bit-identical until this is decided.

## What the GPU path already does

`BoundaryConditionsGPU.jl` reduces each face to
`_GhostFaceOp{T} = Tuple{Bool, Int, T}` — `(active, source_row, factor)` with
`ghost = factor * pot[source_row]`. Two ten-entry tables
(`_ghost_ops_compressed`, `_ghost_ops`) map axis boundary types to op pairs, e.g.

```julia
_ghost_ops(::DiscreteAxis{T, :infinite, :reflecting}, n::Int, gbf) where {T} =
((true, 2, gbf[1]), (true, n - 2, one(T)))
```

and one kernel applies all six faces per color. The tables are pure functions of
the axis type, the array size and the boundary factors — trivially unit-testable
without building a solver setup.

## Proposal

1. Move the op tables out of the GPU file into a backend-neutral
`BoundaryConditions.jl`.
2. Replace all ~25 broadcast methods with one CPU applier (~15 lines): for each
dimension, fetch the op pair and execute the two `ghost .= factor .* view(...)`
broadcasts when active. The GPU applier keeps consuming the same tables via the
fused kernel.
3. Keep `apply_boundary_conditions_at_r0!` as the one genuine special case (it
now also has a KernelAbstractions variant on the branch).
4. While unifying, decide the two open questions explicitly, with tests:
* infinite-r source row: `end−1` (consistent with the factor definition) or
`end−2` (status quo);
* complete the boundary-type product space (the table form gives the two
missing z-axis combinations for free).
5. Port the bit-identity harness used during the review: run the example configs
(BEGe, InvertedCoax, CGD, 2D IVC; with and without depletion handling) before
and after and require `data_before == data_after` exactly — except where 4.
changes behavior deliberately.

Net effect: ~250 lines deleted, one place to read the BC rules, both backends
provably identical, and the remaining physics questions surfaced instead of
duplicated.

Start by reading BoundaryConditionsGPU.jl, then compare its operation tables with the methods in BoundaryConditionsCartesian.jl and BoundaryConditionsCylindrical.jl. Run the red-black layout coverage in test/test_electric_field.jl and review the cited factor definition in PotentialCalculationSetupCylindrical.jl. Done means the CPU and GPU paths share tested boundary rules, the missing combinations are covered, and the example configurations remain bit-identical except for deliberate decisions.

Written by the indexing model from the issue text.