Add ColumnwiseTridiagonalPreconditioner for the CG Poisson solver

[xkykai]

Summary

Adds a block-diagonal preconditioner for ConjugateGradientPoissonSolver based on Marshall et al. (1997, §4.2.2). For each horizontal column (i, j) it solves the vertical (k-direction) sub-system of the symmetric volume-weighted Laplacian V∇² exactly while discarding horizontal couplings — i.e. M = Lz⁻¹, applied as one batched Thomas sweep per column (reusing the existing BatchedTridiagonalSolver over ZDirection).

For ocean-like problems (large Nz, thin layers, stretched vertical grids) this is a much stronger preconditioner than DiagonallyDominantPreconditioner. On the staircase_convection configuration (16×16×128, Bounded³ + slope, Δz/Δx = 0.125) it gives a 4–5× reduction in CG iterations (478 → 85 in Float64, 184 → 28 in Float32).

What's in the kernel

A naïve transcription of K₃D fails even in its target regime; the implementation includes three corrections:

1. BC / immersed-face masking. The preconditioner builds Lz directly where Bounded-domain and immersed face couplings are masked with inactive_cell(i, j, k±1, grid).
2. Neumann–Neumann regularization. Each column's Lz is singular (its rows sum to zero), and its null space is Nx·Ny-dimensional vs. the operator's single global constant — so the global gauge fix cannot cover it. A multiplicative diagonal shift (1 + ε) with ε = 1/100 lifts the column-constant modes off the kernel (a Helmholtz-style screening) without disturbing the rough vertical modes K₃D needs. ε = 1/100 is empirically found to be useful for both Float32 and Float64.
3. Vertically-isolated active cells. A PartialCellBottom thin surface cell with both vertical neighbors inactive has no vertical sub-system; the diagonal would collapse to exactly zero and stall CG. These are treated as identity (b = 1), like inactive cells.

Regime of effectiveness

κ(MA) ≈ 1 + (Δz/Δx)²·Nz²/π²:

• Oceanic / hydrostatic ((Δz/Δx)²Nz² ≪ 1): near-exact, CG terminates in O(1) iterations.
• Isotropic (Δz/Δx ≈ 1): conditioning is worse than no preconditioner — use DiagonallyDominantPreconditioner or the FFT solver instead. This is expected, not a bug.

Also in this PR

• Fix DiagonallyDominantPreconditioner: remove the spurious V⁻¹ᶜᶜᶜ factor from the Ax±/Ay±/Az± coefficients so the preconditioner matches the V∇² operator on non-uniform grids. heuristic_residual is unchanged.

[Copilot]

Pull request overview

This PR adds a new block-diagonal preconditioner for ConjugateGradientPoissonSolver that inverts the vertical (columnwise) tridiagonal subsystem of the symmetric volume-weighted Laplacian V∇², intended to significantly reduce CG iterations in strongly anisotropic (ocean-like) grids. It also fixes the existing DiagonallyDominantPreconditioner so its coefficients match the V∇² operator on non-uniform grids.

Changes:

• Implement ColumnwiseTridiagonalPreconditioner, backed by a BatchedTridiagonalSolver in ZDirection, with masking and regularization to handle immersed/BC effects and per-column nullspaces.
• Fix DiagonallyDominantPreconditioner coefficient definitions by removing an incorrect V⁻¹ᶜᶜᶜ factor.
• Extend CG Poisson solver tests to cover anisotropic grids and PartialCellBottom cases with multiple preconditioners.

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated no comments.

File	Description
src/Solvers/conjugate_gradient_poisson_solver.jl	Adds the new ColumnwiseTridiagonalPreconditioner and corrects coefficient definitions for DiagonallyDominantPreconditioner to align with V∇².
test/test_conjugate_gradient_poisson_solver.jl	Adds new testsets for anisotropic grids and immersed PartialCellBottom to validate CG solver + preconditioner behavior.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[tomchor]

How hard would it be to create a test that actually tests the result?