DomainDecomposition.coarsen, refine fix, ghost-sync fixes in axpy and essential BCs - #91
Merged
Merged
Conversation
coarsen(factors) builds an aligned coarse decomposition sharing the process topology of self (needed for geometric multigrid). refine() now shares the topology too instead of re-running compute_dims/Create_cart on the old ncells without mpi_dims_mask, and computes local_ncells from the new starts/ends. Bump version to 0.3.0. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Zeroing boundary coefficients left ghost copies on neighbouring processes stale while ghost_regions_in_sync stayed True, so a following StencilMatrix or derivative dot used outdated values. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
y += a*x left y marked in sync when x was not (stale ghost data in y) and wrongly changed the flag of the input x. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Member
Author
|
Paired struphy PR: struphy-hub/struphy#91 (the struphy CI runs against this branch via the submodule pointer). |
max-models
pushed a commit
to struphy-hub/struphy
that referenced
this pull request
Oct 5, 2026
…iscretization) (#91) Replaces the previous multigrid attempt (#60, and the earlier state of this PR) with a new geometric multigrid that runs under MPI. Its V-cycle is used as a preconditioner for CG. **Paired feectools PR:** struphy-hub/feectools#91 (into `devel-tiny`). The submodule points to its branch `ddm-coarsen`; the `feectools` pin in `pyproject.toml` can be raised once 0.3.0 is on PyPI. ## What it does `MultiGridSolver(A, derham, domain, MultiGridOptions(...))` solves `A x = b` for a symmetric positive (semi-)definite `A` on one Derham space. Typical example: `sigma * M0 + grad.T @ M1 @ grad`. - **Grid hierarchy** (`multigrid/hierarchy.py`): halves the number of elements in every direction where possible, and stops in a direction that gets too small (semi-coarsening). The MPI decomposition of each coarse level is aligned with the fine one, via `DomainDecomposition.coarsen` and the new `Derham(..., domain_decomposition=...)` argument. - **Grid transfer** (`multigrid/transfer.py`): - `SplineProlongation` is the exact embedding of the nested spline spaces, for any degree, periodic or clamped, B- and D-splines, with homogeneous Dirichlet BCs. - It is applied one direction at a time on the local ghosted arrays. - The restriction is its transpose. - **Coarse operators** (`multigrid/coarsen.py`): the fine operator is walked as an expression tree and rebuilt on each coarse grid. No string recipes are needed. - Sums, compositions, scalings and block operators are rebuilt from their children; derivative, boundary and identity operators are recreated on the coarse spaces. - Mass and basis projection operators are recreated from their new `to_dict()`/`from_dict()`. - Coarse pieces are cached, so changing a scalar re-assembles nothing. - **Smoothers** (`multigrid/smoothers.py`), selected via `MultiGridOptions`: - Chebyshev (default) with an approximate mass inverse, the exact Jacobi diagonal, or the identity inside; - damped Jacobi; - a fixed number of PCG steps. The exact diagonals of composite operators are computed by colored probing. - **V-cycle** (`multigrid/preconditioner.py`): symmetric pre- and post-smoothing. The coarsest level is solved directly on every rank (or by CG). There is an optional `nullspace="constants"` for singular problems. - **Propagators:** `ImplicitDiffusion` and `PoissonSolve` get `precond="MultiGrid"` and a `multigrid: MultiGridOptions` field. The preconditioner is only updated when `sigma_1` (e.g. `sigma_1/dt`) changes. The other `precond` values keep their previous behavior (`pc=None`). - **Removed:** the old `multigrid_solver.py` (eval-string operators, serial-only transfer, 2D coarsening) and the 3,246-line research file `feec/tests/test_multigrid.py`. ## Results Poisson, CG iterations to a relative tolerance of 1e-8, default options: | Case | Multigrid + CG | Plain CG | |---|---|---| | 2D, 16² → 128², p = 2, 3, Dirichlet or periodic, 1 and 4 ranks | 7 at every size | 21 → 149 | | 3D, 8³ → 32³, p = 2, 3, 4 ranks | 6 – 9 | 28 – 135 | ## Notes and known limitations - **Which smoother.** The default mass-inverse smoother is robust in the spline degree but degrades on strongly curved mappings (Colella α = 0.1: 35 iterations). The Jacobi-based smoother is robust to the mapping (10 iterations there) but slower at degree 3 in 3D (24–31 iterations). It also costs a one-time setup of up to (2p+1)³ operator applications per level. - **1-forms and 2-forms.** Problems such as curl-curl will need dedicated smoothers (Hiptmair / Arnold–Falk–Winther type). The coarsening and grid transfer already support them. - Polar splines are not supported (clear error). The transfer operators are NumPy-only. - **Coarsest level.** The coarse matrix is assembled by applying the operator to every unit vector, which is fine for small coarse grids. - **Performance with Dirichlet BCs.** With several ranks, the existing `MassMatrixPreconditioner` is the main cost: it calls SuperLU on many right-hand sides. A banded solver in feectools would help. ## Tests - New in `linear_algebra/tests/`, run serially and on 4 ranks: - `test_multigrid_transfer.py`: exactness of the transfer, R = Pᵀ, R M_h P = M_H. - `test_multigrid_coarsen.py`: `to_dict` round trips, R A P = A_H for `GᵀM₁G` and `CᵀM₂C`. - `test_multigrid_solver.py`: smoother symmetry, V-cycle symmetric positive definite and contracting, iteration counts that don't grow with the grid, `update`. - New in `propagators/tests/test_poisson.py`: `PoissonSolve` with multigrid on the Colella mapping (periodic, Dirichlet, Neumann), and `ImplicitDiffusion` with a changing `dt`. - No new test folders, so no CI shard changes. - Existing tests re-run: mass / basis-op tests (89 passed), Poisson / gyrokinetic Poisson (92 passed). 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Needed for geometric multigrid in struphy (paired struphy PR: struphy-hub/struphy#91), plus two ghost-region bug fixes found while testing it under MPI.
Changes
DomainDecomposition.coarsen(factors)(new). Builds the decomposition of a coarsened grid that is aligned withself: every process owns exactly the coarse cells covering its fine cells. It shares the process topology and communicators, so it is not collective. RaisesValueErrorif starts/ends are not divisible by the factor.DomainDecomposition.refinefixed.compute_dims/Create_carton the oldncellswithoutmpi_dims_mask, so it could pick a different process grid. It also calledCreate_cartagain, which every rank has to join.local_ncellsfrom the old starts/ends.Both methods now use one helper,
_with_element_partition, which validates the partition.apply_essential_bc_stencilnow setsghost_regions_in_sync = Falseon vectors. Boundary coefficients may be ghost entries of neighbouring processes. Those kept stale values while still being marked in sync, so the next matrix-vector product used outdated data.StencilVectorSpace.axpy(behindmul_iadd) updated the sync flag of the inputxinstead of the outputy.y += a*xtherefore stayed marked in sync whenxwas not.Version bumped to 0.3.0.
Tests
ddm/tests/test_coarsen.py(serial +mpi),api/tests/test_essential_bc_ghosts.py(serial +mpi),linalg/tests/test_axpy_ghost_sync.py.linalgserial tests: 7544 passed.femMPI tests on 4 ranks: 144 passed.test_cart_3d.py,test_block.pyandtest_toarray.pyfail to collect locally because theparallelmarker is not registered. That was already the case ondevel-tiny.🤖 Generated with Claude Code