Motivation
The direct sparse-LU solver (spsolve / UMFpack) scales roughly as O(N^1.5) for square N×N grids and somewhat better for elongated domains (see benchmarks in benchmarks/bench_solvers.py). For very large grids this becomes the bottleneck, and an iterative solver has not proven practical: the biharmonic FD matrix is not an M-matrix, so standard smoothers and AMG fail to converge for geologically realistic (heterogeneous) Te fields. See Issue #57 for the full investigation.
The most viable path to large-grid scaling is domain decomposition with direct subdomain solvers.
Approach
Decompose the domain into overlapping or non-overlapping strips along the longer axis. Each strip is solved directly (UMFpack). The strips are coupled iteratively (e.g., Schwarz alternating method) until the interface residual falls below a tolerance.
Why strips?
For a domain with ny rows and nx columns stored in row-major order:
- The sparse-LU fill-in scales as O(min(nx, ny)² × max(nx, ny)) — the short dimension sets the separator size under nested dissection.
- A strip of height
ny_s ≪ ny (or width nx_s ≪ nx) replaces the global min with ny_s, dramatically reducing fill.
- This is consistent with the empirical benchmark result that elongated (e.g. 4:1) domains solve faster than square domains of the same total cell count.
Coupling options
| Scheme |
Interface condition |
Convergence |
| Non-overlapping (Schur complement) |
Continuity of w and ∂w/∂n |
Direct (one solve per strip) — no iteration needed if the Schur complement is formed explicitly |
| Overlapping Schwarz |
Exchange of w in overlap zone |
Iterative — simple to implement, converges geometrically |
| FETI / dual primal |
Lagrange multipliers at interfaces |
Parallel-friendly, robust |
The overlapping Schwarz approach is probably the right starting point: it is straightforward, leverages the existing direct solver on each strip, and the biharmonic equation is elliptic so convergence is guaranteed.
Expected benefit
For an N×N domain decomposed into K strips of height N/K:
- Per-strip fill: O((N/K)² × N)
- Total: O(N³ / K²) vs O(N^1.5) for the global solve (which itself beats O(N^1.5) only for square domains)
- Break-even over global direct solve: around K ≈ √(N / min_dim) — meaningful for grids larger than ~500×500
Relationship to other issues
Scope
Not assigned to any release milestone. This is a significant algorithmic addition that warrants its own branch and careful benchmarking before it touches master.
Motivation
The direct sparse-LU solver (
spsolve/ UMFpack) scales roughly as O(N^1.5) for square N×N grids and somewhat better for elongated domains (see benchmarks inbenchmarks/bench_solvers.py). For very large grids this becomes the bottleneck, and an iterative solver has not proven practical: the biharmonic FD matrix is not an M-matrix, so standard smoothers and AMG fail to converge for geologically realistic (heterogeneous) Te fields. See Issue #57 for the full investigation.The most viable path to large-grid scaling is domain decomposition with direct subdomain solvers.
Approach
Decompose the domain into overlapping or non-overlapping strips along the longer axis. Each strip is solved directly (UMFpack). The strips are coupled iteratively (e.g., Schwarz alternating method) until the interface residual falls below a tolerance.
Why strips?
For a domain with
nyrows andnxcolumns stored in row-major order:ny_s ≪ ny(or widthnx_s ≪ nx) replaces the globalminwithny_s, dramatically reducing fill.Coupling options
wand∂w/∂nwin overlap zoneThe overlapping Schwarz approach is probably the right starting point: it is straightforward, leverages the existing direct solver on each strip, and the biharmonic equation is elliptic so convergence is guaranteed.
Expected benefit
For an N×N domain decomposed into K strips of height N/K:
Relationship to other issues
Scope
Not assigned to any release milestone. This is a significant algorithmic addition that warrants its own branch and careful benchmarking before it touches
master.