Function distributed_gs_solve 

pub fn distributed_gs_solve(
    psi: &Array2<f64>,
    source: &Array2<f64>,
    r_axis: &[f64],
    z_axis: &[f64],
    dr: f64,
    dz: f64,
    cfg: &DistributedSolverConfig,
) -> FusionResult<DistributedSolveResult>

Distributed Grad-Shafranov solver using additive Schwarz domain decomposition with Rayon thread-parallelism.

Each Schwarz iteration:

1. Splits the global Ψ into 2D tiles (with halo overlap).
2. Runs inner_sweeps local Red-Black SOR sweeps on each tile in parallel via Rayon.
3. Injects tile cores back into the global array.
4. Exchanges halos (serial reference — MPI-ready interface).
5. Computes global residual; checks convergence.

The 5-point GS stencil: R d/dR(1/R dΨ/dR) + d²Ψ/dZ² = -μ₀ R J_φ

discretises to the same operator as fusion_math::sor::sor_step, but applied independently on each tile with local boundary from halo data.