Schur-complement domain decomposition solver for large-scale FD (memory-bounded, parallel substructuring)

[milanofthe]

Summary

Add a Schur-complement (non-overlapping substructuring) domain-decomposition solver to the frequency-domain backend so we can run larger FD problems without the direct-solver memory wall, with a path to parallel scaling. The goal is purely solver scalability for a single monolithic FD solve, not coupling separate models (that is the substructuring-for-co-simulation track, see "Related" below).

Today the FD backend is purely direct: assemble_and_solve builds one global complex-symmetric system and factorizes it with PARDISO (mtype=6, METIS-ordered, upper-CSR) or faer LU via the SparseSolver trait. The fill-in of that single global factorization is the hard memory/time wall for large RFIC-scale meshes. Schur DD never factorizes the full matrix: it factorizes each subdomain's interior block independently and reduces the problem to a much smaller interface ("skeleton") system. Peak memory is then bounded by the largest subdomain factor plus the interface system, not by the whole.

Why Schur DD (and how it relates to #2)

• Iterative solver: AMS-preconditioned Krylov for H(curl) systems #2 (AMS-preconditioned Krylov) is the iterative route: best asymptotic O(N) memory, but FD curl-curl is indefinite with a large gradient null-space, and high-frequency/indefinite Helmholtz convergence is genuinely hard. High ceiling, high risk.
• Schur DD is the direct/hybrid route: each subdomain is solved by a direct factorization, so there is no global convergence fragility; only the (much smaller) interface problem is solved directly or iteratively. This is the approach production FD-FEM tools (e.g. HFSS DDM) use to fit large problems into bounded memory, and it is robust at high frequency because the local solves are direct. The two are complementary: AMS could later precondition the interface or the interior solves.

Math sketch (complex-symmetric, mtype=6)

Partition the free DOFs into per-subdomain interior sets I_s and the shared interface set Γ. With interior DOFs ordered first, the reduced system is

[ K_II   K_IΓ ] [ x_I ]   [ b_I ]
[ K_ΓI   K_ΓΓ ] [ x_Γ ] = [ b_Γ ]

K_II is block-diagonal across subdomains (no interior-interior coupling between subdomains). Per subdomain s:

S_s = K_ΓΓ^(s) − K_ΓI^(s) (K_II^(s))^{-1} K_IΓ^(s)        (local Schur complement)
g_s = b_Γ^(s)   − K_ΓI^(s) (K_II^(s))^{-1} b_I^(s)         (condensed RHS)

Assemble S = Σ_s S_s, g = Σ_s g_s, solve the interface system S x_Γ = g, then back-substitute per subdomain: x_I^(s) = (K_II^(s))^{-1} (b_I^(s) − K_IΓ^(s) x_Γ).

• S is complex-symmetric (inherits mtype=6), so the same solver family applies.
• Multi-RHS maps cleanly onto the existing "factor once, solve many RHS" port loop: b_I, b_Γ become matrices with one column per driven port; interior back-solves and the interface solve are multi-column.
• Memory: each (K_II^(s))^{-1} is a small factor; the global factor is never formed. The remaining cost is S (dense N_Γ × N_Γ if formed explicitly) — see phasing.

Where it plugs into the code (grounded)

• Solver seam: crates/rapidfem-fd/src/solver/mod.rs already defines SparseSolver (factorize/refactorize/solve) with PARDISO/faer/Accelerate backends and SolverChoice::from_env(). Schur DD is a layer that holds one SparseSolver per subdomain interior block plus a solver for the interface — it reuses, not replaces, the existing backends.
• Assembly: crates/rapidfem-fd/src/assembly.rs assemble_and_solve_with_pml builds COO (rows, cols, data_e, data_b) via assemble_global_matrices, forms K = E − k0²·B + Robin, eliminates PEC DOFs into free_dofs / dof_to_free, then factorizes + solves per-port RHS (port_vectors). The DD path forks here: element-filtered assembly per subdomain into the I/Γ blocks instead of one global COO.
• DOF model: crates/rapidfem-fd/src/basis.rs Nedelec2Basis — n_field = 2·n_edges + 2·n_tris, edge_to_field (2 DOFs/edge), tri_to_field (8 DOFs/tri). DOFs live on edges and faces, which determines interface-DOF identification.
• Mesh adjacency for partitioning + interface detection: crates/rapidfem-core/src/mesh.rs — tri_to_tet (face→ up to 2 tets), tet_to_tri, tet_to_edge, vtag_to_tet. A tet dual graph (tets joined across shared faces via tri_to_tet) feeds the partitioner; a DOF is an interface DOF iff its supporting edge/face is touched by tets in ≥2 subdomains.

Proposed design and phasing

Phase 0 — correctness spike (single process, explicit S).

• Partition tets into k subdomains; classify DOFs into I_s / Γ.
• Element-filtered assembly of the four blocks per subdomain.
• Factorize each K_II^(s) (reuse SparseSolver); form S and g explicitly; solve S x_Γ = g (dense or sparse direct); back-substitute.
• Gate: bit-for-bit agreement with the monolithic solve on WR-90 + coax (same S-parameters to solver precision). This is the non-negotiable correctness check.

Phase 1 — integration + memory relief.

• Expose as SolverChoice::Schur (RAPIDFEM_SOLVER=schur) and/or ProblemFD(..., n_subdomains=k).
• Instrument peak RAM + factorization time; produce the N vs peak-RAM curve vs monolithic on a mesh-refinement series (report-harness plot).
• Target: clearly lower peak memory than monolithic at equal accuracy for k ≥ 2.

Phase 2 — scalable interface solve.

• Replace explicit S with a matrix-free interface GMRES (matvec S·v = Σ_s [K_ΓΓ^(s) v − K_ΓI^(s)(K_II^(s))^{-1}(K_IΓ^(s) v)] via cached local factorizations) so the dense N_Γ × N_Γ is never stored — this is the real scaling lever for many subdomains/large interfaces.
• Add an interface preconditioner (lumped / Neumann-Neumann, later a coarse space for subdomain-count-independent convergence).

Phase 3 — parallelism + robustness.

• Parallelize subdomain factorizations and back-solves with rayon (already a dep behind the parallel feature).
• Optional 2-level coarse space.

Partitioner

• MVP: recursive coordinate bisection (RCB) on tet centroids — no new dependency, deterministic, good enough to validate and to relieve memory. Aligns with the pure-Rust direction in Pure-Rust geometry kernel + tet mesher (replace gmsh/OCC dep) #3.
• Later: graph partitioning (METIS-style multilevel) on the tet dual graph for better edge-cut / interface-size minimization. Evaluate a Rust METIS binding vs a native multilevel partitioner.

Key risk: interior resonance

A non-overlapping interior block K_II^(s) can be singular or ill-conditioned at the solve frequency if the subdomain has a spurious interior resonance (a known DD pitfall for indefinite Maxwell). Mitigations to evaluate: a small complex frequency shift on interior blocks, or Robin/impedance transmission conditions on interior interfaces (which is exactly the optimized-Schwarz variant and reuses the existing port-Robin assembly). Phase 0/1 should detect near-singular K_II^(s) (PARDISO already uses pivot perturbation, iparm[9]=13) and flag it.

Validation plan

• Correctness: S-parameters identical to monolithic (to solver precision) on WR-90, coax, and at least one multi-mode / reflective case; reciprocity and passivity gates from the existing FD suite must still pass.
• Memory/scaling: peak-RAM and wall-time vs N and vs k (subdomain count), monolithic vs Schur, plotted in the report harness. Local testing only (never CI — FEM solves do not belong in runners).

Scope / non-goals

• In scope: a memory-bounded, eventually-parallel solver for one FD problem; conforming interfaces (already produced by g.fragment() within one gmsh model).
• Out of scope here: non-conforming/mortar interfaces; coupling separately authored meshes; ROM-level substructuring for transient co-simulation (shares the subdomain machinery but is tracked separately).

Related

• Iterative solver: AMS-preconditioned Krylov for H(curl) systems #2 Iterative AMS-preconditioned Krylov for H(curl) — complementary; AMS could precondition the Schur interior/interface solves.
• Pure-Rust geometry kernel + tet mesher (replace gmsh/OCC dep) #3 Pure-Rust geometry kernel + tet mesher — an RCB/native partitioner aligns with dropping external deps.
• The same subdomain-factorization + interface machinery is the foundation for substructuring multiport ROMs (passive blocks reduced and connected at ports for transient co-simulation with embedded transistors).

References

• A Parallel Direct Domain Decomposition Solver Based on Schur Complement for Electromagnetic Finite Element Analysis (IEEE TAP, 2021)
• Direct Domain Decomposition Method (D3M) for Finite Element Electromagnetic Computations (arXiv:2002.05018)
• Direct Schur complement method by domain decomposition based on H-matrix approximation (Comput. Vis. Sci.)
• A preconditioner for the Schur complement domain decomposition method (Cros, DD14)
• Non-Overlapping Domain Decomposition via BURA Preconditioning of the Schur Complement (Mathematics, 2022)

[milanofthe]

Implementation-fit review (grounded in the current FD code)

Read the real assembly/solve path (assembly.rs::assemble_and_solve_with_pml, tet_assembly.rs::assemble_global_matrices, basis.rs::Nedelec2Basis, mesh.rs). Verdict: primal Schur fits cleanly, and is in fact the cleanest DD variant for this codebase — cleaner than an iterative Robin/FETI-2LM approach — for three concrete reasons found in the code.

Why primal Schur is the clean fit here

1. Global assembly is already per-tet additive with implicit tet ownership. assemble_global_matrices builds nnz = n_tets*400 triplets, one 20×20 element block per tet, written in 400-entry chunks (tet_assembly.rs:399-449). The owning tet of triplet i is exactly i/400. So we can route every triplet into a subdomain block with zero new element-assembly code — just tet_subdomain[i/400].
2. It stays in the existing global DOF / free_dofs space — no submesh extraction. Primal Schur partitions DOFs and triplets in place. The iterative Robin/optimized-Schwarz alternative would instead need real submesh objects (subset of tets, renumbered nodes/edges/tris, interface faces re-tagged as Robin) feeding assemble_and_solve per subdomain — that submesh-slicing + re-tag layer is the new infra we avoid by going primal. (We keep Robin transmission in our pocket as the phase-2 robustness upgrade; see interior resonance.)
3. Interface DOFs are directly classifiable from existing maps. Face DOFs: mesh.tris already contains all interior faces, and tri_to_tet gives the (≤2) adjacent tets — a face DOF is interface iff its two tets are in different subdomains. Edge DOFs: invert tet_to_edge once to get edge→tets, then an edge DOF is interface iff its incident tets span ≥2 subdomains. Map edge/tri → field DOFs via basis.edge_to_field / basis.tri_to_field. Interior DOFs belong to exactly one subdomain, so K_II is genuinely block-diagonal and K_II^(s) extraction has no double-counting; K_ΓΓ is just the global interface block (no per-subdomain split needed because the correction term Σ_s K_ΓI^(s)(K_II^(s))^{-1}K_IΓ^(s) is naturally additive).

Net: we assemble the global COO once (existing code, O(nnz)≈O(N), cheap), then extract blocks by DOF set. We never factorize the global matrix — only the per-subdomain K_II^(s). The global COO in RAM is not the memory wall; the global factor fill-in is, and that is exactly what we skip.

Precise build list (what we actually need to add)

1. Tet partitioner → tet_subdomain: Vec<usize>. MVP: recursive coordinate bisection on tet centroids (no dep, deterministic; aligns with Pure-Rust geometry kernel + tet mesher (replace gmsh/OCC dep) #3). Later: multilevel graph partition on the tet dual graph (tri_to_tet) to minimise interface size.
2. Free-DOF classification interior_s / Γ: build edge_to_tet (invert tet_to_edge), then classify each free DOF via the face/edge straddle test above.
3. COO block router: from the global (rows, cols, data_e, data_b) with tet = i/400 and K = E − k0²B, emit per-subdomain K_II^(s), K_IΓ^(s) and the global K_ΓΓ in local interior/interface numbering. Robin entries (bempty) only sit on exterior port faces, which are never interface faces, so they never touch Γ — they land in their owning subdomain's interior block.
4. Schur driver over SparseSolver: factor each K_II^(s) (one trait instance per subdomain); form S = K_ΓΓ − Σ_s K_ΓI^(s)(K_II^(s))^{-1}K_IΓ^(s) and g = b_Γ − Σ_s K_ΓI^(s)(K_II^(s))^{-1}b_I^(s); solve S xΓ = g; back-substitute x_I^(s) = (K_II^(s))^{-1}(b_I^(s) − K_IΓ^(s) xΓ). The per-port RHS loop becomes the multi-column RHS.
5. Interior-resonance guard (see risk): detect near-singular K_II^(s) (factor failure / tiny pivots — PARDISO already perturbs, iparm[9]=13) and fall back to a small complex frequency shift on the interior block, or to Robin transmission on Γ (phase 2).
6. Interface solve: phase 0 forms S explicitly (dense/sparse direct) — fine for moderate subdomain counts and the correctness gate; phase 2 makes it matrix-free GMRES (S·v via cached K_II^(s) back-solves) so the dense N_Γ × N_Γ is never stored.

Sharpened caveats (honest)

• PARDISO/faer already do nested-dissection ordering = implicit substructuring. So the memory win over monolithic-PARDISO is not automatic. It comes specifically from (a) bounding peak RAM by the largest K_II^(s) using separate solver instances (free each after its Schur contribution when S is matrix-free), (b) parallel subdomain factorizations, (c) reusing the subdomain factorizations across the many port RHS and the whole frequency sweep. Phase 1 must measure this against monolithic, not assume it.
• Explicit S is itself O(N_Γ²) memory — for many subdomains / large interfaces it just relocates the wall. The matrix-free interface solve (item 6) is the real scaling lever, not an optional extra.
• Interior resonance is inherent to primal Schur (the interior blocks have a homogeneous-Dirichlet artificial boundary on Γ, so a subdomain can be singular at a solve frequency). The Robin/FETI-2LM variant removes it but reintroduces the submesh/transmission machinery we deliberately skipped — hence: ship primal first (cleanest integration), keep Robin as the phase-2 robustness path.

Revised phase 0 (unchanged in spirit, sharper in mechanism)

Partition (RCB) → classify free DOFs → route the existing global COO into blocks → factor K_II^(s) → form S explicitly → solve → back-substitute, gated by bit-for-bit agreement with the monolithic solve on WR-90 + coax. No new assembler, no submesh objects, pure DOF-space surgery on top of the current pipeline.

[milanofthe]

Phase-0 implemented and validated (branch feature/schur-dd)

Primal Schur-complement DD is built, wired end-to-end, and validated. Commits: 91834b6 (core solver + assembly fork into ReducedSystem), a3caef0 (sweep wiring), 858284b (memory proxy + validation script).

What landed

• crates/rapidfem-fd/src/solver/schur.rs: RCB tet partitioner, free-DOF interior/interface classification, COO block router, and the Schur driver (factor each K_II^(s) via the existing SparseSolver, form/solve the interface system, back-substitute). Reuses the global assembly — no new element assembler, no submesh extraction (works in global DOF space, exactly as the fit-review predicted).
• n_subdomains threads SolverConfig → frequency_sweep_with_pml (Schur branch) → Python ProblemFD(geometry, n_subdomains=k).

Key finding (changed the validation gate)

Bit-equality with the monolithic solve is the wrong criterion: this FD curl-curl system is ill-conditioned (cf. #2), so the monolithic direct factorization itself leaves an O(1e-3) residual. The honest gate is the residual ‖Kx−b‖/‖b‖. With faer, Schur DD reaches 1e-14 while the monolithic solve sits at 2e-3 — DD is not just equivalent, it is more accurate here (it factorizes smaller, better-conditioned blocks). S-parameters (the physical observable) are robust and agree to ~1e-4.

Validation (Schur k=4 vs monolithic, max|ΔS| over band+ports, gate < 5e-3)

Geometry	class	max|ΔS|
WR-90	hollow, modal port	1.7e-4
coax	TEM port	8.4e-5
iris filter	internal PEC, reflective	2.0e-4
microstrip (60k DOFs)	inhomogeneous wave port + open ABC	3.7e-4

(Tier-1 + Tier-2 + Tier-4 essence. RFIC multilayer cases deferred per scope.)

Memory proxy (the actual goal)

Direct-LU fill-in scales super-linearly, so the bound is the largest interior block we factorize, never the global matrix:

k	largest interior block	interface
1	100%	0%
2	48%	4%
4	23%	9%
8	12%	16%

On the big microstrip case (353k free DOFs), k=4 leaves only 1.3% interface — the win grows with problem size, as intended.

Remaining (phase-1/2)

• Matrix-free interface GMRES (avoid the explicit dense S) — the real scaling lever for many subdomains.
• Interior-resonance robustness: hard-singular detection is in place (factor failure → clear error); the complex-shift / Robin-transmission fallback is phase-2.
• True peak-RAM/wall-time instrumentation + parallel subdomain factorization (rayon).
• Graph partitioner (METIS-style) to shrink the interface vs RCB.