Skip to content

Kronecker preconditioners: StiffnessPreconditioner, refactor of the mass preconditioners; bump feectools - #731

Merged
spossann merged 20 commits into
develfrom
kron-stencil-matmul
Oct 9, 2026
Merged

spossann merged 20 commits into
develfrom
kron-stencil-matmul

Conversation

@spossann

@spossann spossann commented Oct 8, 2026 •

Copy link
Copy Markdown
Member

Summary

Points the feectools submodule to the branch of struphy-hub/feectools#100, so that the struphy test suite runs against it before that PR is merged.

struphy-hub/feectools#100 changes KroneckerStencilMatrix and KroneckerLinearSolver:

  • factors on groups of axes (e.g. 2d x 1d);
  • A @ B returns a new ComposedKroneckerStencilMatrix (operands in multiplicands, exact Kronecker factors of the product in mats);
  • ComposedLinearOperator.multiplicants → multiplicands, with a deprecated alias;
  • KroneckerSumSolver, the exact inverse of sums of Kronecker products by fast diagonalization;
  • factor_ndims in KroneckerLinearSolver;
  • docstrings and type annotations; feectools 0.7.0.

Existing calls are unchanged, so struphy needs no code changes. struphy builds these objects in feec/preconditioner.py (MassMatrixPreconditioner) and feec/projectors.py.

multiplicants → multiplicands

  • All struphy uses (preconditioner, projectors, multigrid coarsening) now use multiplicands.
  • linop_nbytes in feec/memory.py adds up all child attributes, so a ComposedKroneckerStencilMatrix counts both its operands and the factors of the product.

Refactor of MassMatrixPreconditioner

  • The construction of the Kronecker approximation is split into module-level helpers in feec/preconditioner.py. They cover the boundary conditions, gathering array weights over MPI, the 1d weight, the 1d mass matrix and its BCs, the 1d solver, the process-local 1d factor, the block-diagonal assembly and the operator to invert.
  • MassMatrixDiagonalPreconditioner now calls these helpers instead of repeating the same ~200 lines.
  • Values and the order of MPI calls are unchanged. The copy into the process-local 1d stencil matrix is now vectorized, and the block-diagonal operators are no longer hard-coded to 3 components.
  • Docstrings, comments and type annotations.

MassMatrixPreconditioner: new option and fixes

  • New option weight_reduction="midpoint" | "average" sets how the weight is reduced to 1d in the directions other than dim_reduce. The default "midpoint" gives the same results as before, verified for M0, M1, M2, Mv, M1n and array weights, serially and on 2 ranks. "average" uses the mean over those directions: Gauss quadrature of the Derham for array weights, 16-point Gauss–Legendre for callables.
  • Array weights are reduced with a single Allreduce, instead of a subcommunicator plus Bcast.
  • The mass matrix is located in the composed operator by identity, and the approximate inverse is applied exactly there. Both preconditioners share this code (_apply_composed). Before, a mass matrix in the left-most position was applied instead of inverted.
  • FFTSolver copies its column, which was a view of the dense matrix also used for the stencil factor. It now stabilises a singular matrix once, at construction. is_circulant is vectorized.

PCG iterations (tolerance 1e-8) on HollowTorus, grid 12×16×6, degree (2,3,2), as a reference for choosing weight_reduction:

M0 M1 M2 M3 Mv M1n M2n Mvn
midpoint 13 22 22 13 31 673 600 1000 (not converged)
average 13 15 15 13 30 798 727 1000 (not converged)

On HollowCylinder both need 2 iterations for every matrix.

Kronecker preconditioners: base class, stiffness preconditioner

KroneckerPreconditioner is a new base class for preconditioners of the form P = B E · S Ã⁻¹ S · Eᵀ Bᵀ:

  • Ã (matrix) approximates the core operator, and its exact inverse is solver. Subclasses build the approximation.
  • S is an optional diagonal scaling D̂^{1/2} D^{-1/2} (Loli–Sangalli–Tani), with D the diagonal of the core operator and D̂ that of the approximation.
  • The base class holds the common interface (solve, dot, the composition with B E, the properties).

Mass preconditioners

  • MassMatrixPreconditioner is a subclass. It gets the options diagonal_scaling (default False) and dim_reduce=None (unit weights in all directions), and the method update_mass_operator. With diagonal_scaling=True, the PCG iterations for M1 on a HollowTorus drop from 20 to 8.
  • MassMatrixDiagonalPreconditioner is now MassMatrixPreconditioner(dim_reduce=None, diagonal_scaling=True), kept for its name in the solver options. It no longer assembles a separate 3d logical mass matrix; D̂ comes from the Kronecker approximation. Results are unchanged up to rounding (1e-13), checked on HollowTorus and IGAPolarCylinder, including update_mass_operator, serially and on 2 ranks.

StiffnessPreconditioner (new) preconditions the stabilized stiffness operators Gᵀ M1 G + σ M0, Cᵀ M2 C + σ M1 and Dᵀ M3 D + σ M2:

  • With Kronecker mass matrices, each diagonal block is a sum of Kronecker products of 1d stiffness and mass matrices. It is inverted exactly with KroneckerSumSolver, one solver per component (the components have different spline types, e.g. DNN/NDN/NND for M1). For grad, this is exact on the logical cube.
  • curl and div neglect the off-diagonal blocks (block Jacobi). Block Jacobi alone overestimates the operator on the kernel of the derivative, so a kernel correction σ⁻¹ d₋ P₋ d₋ᵀ is added (default, needs σ > 0): with d₋ = G and the grad preconditioner for curl, and d₋ = C and the curl block Jacobi for div.
  • Weights: weights="average" (default) approximates the mass-matrix weights by one common separable shape per direction plus a mean per term, which keeps the fast diagonalization exact. weights="unit" uses the logical cube. Diagonal scaling (default on) adds the geometry pointwise; for it, the diagonal of the stiffness operator is computed by probing.
  • Polar splines are not supported yet.

PCG iterations (tolerance 1e-8), HollowTorus, grid 12×16×6, degree (2,3,2):

no preconditioner mass preconditioner StiffnessPreconditioner (defaults)
grad, σ = 0 187 – 41
curl, σ = 1 > 3000 716 46
div, σ = 1 > 3000 636 48

On a stretched cuboid (constant, anisotropic weights), grad needs 2 iterations. For large σ (1e4), the mass preconditioner is as good or slightly better (about 13 vs 25 iterations).

Tests

Run locally with the submodule branch: all preconditioner tests (test_preconditioner_transpose.py, the new test_stiffness_preconditioner.py, test_mass_matrices.py::test_mass_preconditioner*, test_reduced_weight_1d) pass serially and on 2 ranks (52 passed each). The new stiffness tests cover the approximation against the operator on the unit cube, PCG iteration counts for grad/curl/div, the kernel correction, weights="average", and independence of the MPI decomposition. The rest is left to CI.

Don't raise the feectools pin in pyproject.toml until 0.7.0 is on PyPI. Before merging, the submodule pointer should move to the merged devel-tiny commit.

🤖 Generated with Claude Code

Points the feectools submodule to the kron-stencil-matmul branch
(KroneckerStencilMatrix @, factors on groups of axes, factor_ndims in
KroneckerLinearSolver; feectools 0.6.0).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
spossann and others added 12 commits October 8, 2026 10:16
- Split the construction of the Kronecker approximation into module-level
  helpers (boundary conditions, MPI weight gathering, 1d weight, 1d mass
  matrix + bcs, 1d solver, process-local 1d factor, block-diagonal
  assembly, operator to invert); MassMatrixDiagonalPreconditioner reuses
  them instead of a copy of the same code.
- Vectorized copy into the process-local 1d stencil matrix; block-diagonal
  operators no longer hard-coded to 3 components.
- Docstrings, comments and type annotations.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
- New option weight_reduction="midpoint" (default, unchanged results) or
  "average" (mean of the weight over the other two directions; Gauss
  quadrature of the derham for array weights, 16-point Gauss-Legendre for
  callables). Array weights are reduced with a single Allreduce instead of
  a subcommunicator + Bcast.
- Locate the mass matrix in the composed operator by identity and apply
  the approximate inverse exactly there (shared _apply_composed for both
  preconditioners); previously a left-most mass matrix was applied instead
  of inverted.
- FFTSolver: copy the column (it was a view of the matrix used for the
  stencil factor) and stabilize a singular matrix once at construction.
- Vectorized is_circulant.
- Tests: weight_reduction in the transpose and MPI array-weight tests, new
  test of the reduced 1d weights against analytic values.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
- Rename uses of ComposedLinearOperator.multiplicants to multiplicands
  (preconditioner, projectors, multigrid coarsening, comments).
- linop_nbytes: sum over all child attributes, so that the factors of a
  ComposedKroneckerStencilMatrix (multiplicands and mats) are counted.
- Bump feectools to the commit with ComposedKroneckerStencilMatrix.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
- New base class KroneckerPreconditioner: P = B E S A_c~^{-1} S E^T B^T,
  where A_c~ approximates the core operator of the composition of a mass
  operator and S is an optional diagonal scaling (D_hat^{1/2} D^{-1/2},
  Loli-Sangalli-Tani). It holds the common interface (matrix, solver, solve,
  dot, properties); subclasses build the approximation.
- MassMatrixPreconditioner and MassMatrixDiagonalPreconditioner are
  subclasses; their results are unchanged (checked serial and on 2 ranks).
- MassMatrixPreconditioner gets the option diagonal_scaling (default False),
  e.g. on HollowTorus the PCG iterations of M1 drop from 20 to 8.
- Helpers for local diagonals of Kronecker/sum/block operators.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
New subclass of KroneckerPreconditioner for the stabilized stiffness
operators G^T M1 G + sigma M0 (and C^T M2 C + sigma M1, D^T M3 D + sigma M2).
With unit weights each diagonal block is a sum of Kronecker products of 1d
stiffness and mass matrices, inverted exactly by fast diagonalization
(KroneckerSumSolver). For grad this is exact on the logical cube; curl and
div are block Jacobi for now (kernel correction follows). The geometry is
included by diagonal scaling (default), with the diagonal of the stiffness
operator computed by probing. Polar splines are not supported yet.

Tests: approximation against the operator on the unit cube, PCG for grad
(2 iterations on the unit cube, HollowCylinder 89 -> 48), independence of
the MPI decomposition, options.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Block Jacobi alone overestimates C^T M2 C + sigma M1 on gradients (where it
equals sigma M1) and was never better than MassMatrixPreconditioner. With
kernel_correction (default, needs sigma > 0) the preconditioner is
P_BJ + sigma^{-1} G P_G G^T, with P_G the preconditioner of G^T M1 G, which
is exact for gradients on the logical cube.

PCG iterations (12x16x6, degree (2,3,2), sigma=1): Cuboid 432 -> 14
(mass preconditioner 144), HollowTorus 1947 -> 150 (mass preconditioner 716).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Same idea as for curl: P_BJ + sigma^{-1} C P_C C^T, where P_C is the curl
block Jacobi preconditioner (sigma = 0 is fine there: on the range of C^T
the gradients do not matter, G^T C^T = 0).

PCG iterations (12x16x6, degree (2,3,2), sigma=1): Cuboid 411 -> 18
(mass preconditioner 142), HollowTorus 2734 -> 175 (mass preconditioner 636).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Fast diagonalization needs one 1d mass matrix per direction for all terms,
so the weights are approximated by a common separable shape prod_d phi_d
(mean over the other directions, averaged over the diagonal blocks of
M_{k+1}) times the mean weight of each term. Exact for constant
anisotropic weights. Default stays weights="unit" with diagonal scaling.

PCG iterations (12x16x6, degree (2,3,2); unit+scaling -> average+scaling):
stretched Cuboid grad 56 -> 2, HollowTorus grad 58 -> 41,
curl (sigma=1) 150 -> 46, div (sigma=1) 175 -> 48.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Average + diagonal scaling gave the fewest PCG iterations in all cases
tested (e.g. HollowTorus curl 150 -> 46, div 175 -> 48).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
- MassMatrixPreconditioner: dim_reduce=None for unit weights in all
  directions; update_mass_operator (rebuilds the approximation only if it
  depends on the weights, updates the diagonal scaling).
- MassMatrixDiagonalPreconditioner is now
  MassMatrixPreconditioner(dim_reduce=None, diagonal_scaling=True), kept for
  its name (solver options). The diagonal of the logical mass matrix is
  taken from the Kronecker approximation instead of assembling the 3d
  logical mass matrix; results are unchanged up to rounding (checked on
  HollowTorus and IGAPolarCylinder, incl. update_mass_operator, serial and
  on 2 ranks).
- Test: only 1d mass matrices are assembled.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@spossann spossann changed the title Bump feectools: KroneckerStencilMatrix matmul and grouped factors Kronecker preconditioners: StiffnessPreconditioner, refactor of the mass preconditioners; bump feectools Oct 8, 2026
WeightedMassOperator sets blocks with zero weight to None, also after
reassembly (e.g. Mrho in VariationalDensityEvolve). _local_diagonal raised
for such blocks, so update_mass_operator of MassMatrixDiagonalPreconditioner
failed in the VariationalBarotropicFluid and VariationalPressurelessFluid
model tests. A missing diagonal block now gives a zero diagonal (scaling 1
for that block); before the refactor, the stale diagonal of the previous
assembly was kept.

Checked: all models with VariationalDensityEvolve/Viscosity/Resistivity
(11 models) and the variational propagator tests pass.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🟡 Changes recommended

The code requires feectools 0.6 APIs while the package metadata still restricts installations to feectools 0.5.0.

1 open finding
What changed in this PR

Adds reusable Kronecker-based mass and stiffness preconditioners and adopts feectools 0.6 APIs.

Changes:

  • Refactors mass preconditioners and adds weight reduction and diagonal scaling.
  • Adds StiffnessPreconditioner with MPI-aware tests.
  • Migrates composed operators to multiplicands.
File Description
feectools Points to the paired feectools development commit.
src/​struphy/​feec/​preconditioner.py Implements Kronecker preconditioners and helpers.
src/​struphy/​feec/​memory.py Counts all composite operator children.
src/​struphy/​feec/​projectors.py Uses the renamed composition API.
src/​struphy/​linear_algebra/​multigrid/​coarsen.py Updates composition coarsening.
src/​struphy/​propagators/​variational_viscosity.py Updates a commented attribute reference.
src/​struphy/​propagators/​variational_resistivity.py Updates a commented attribute reference.
src/​struphy/​feec/​tests/​test_stiffness_preconditioner.py Tests stiffness preconditioning behavior.
src/​struphy/​feec/​tests/​test_preconditioner_transpose.py Expands mass-preconditioner tests.
src/​struphy/​feec/​tests/​test_mass_matrices.py Tests weight reduction and MPI consistency.

🧠 Review effort: Balanced


Give feedback about Copilot approvals in this survey to enter a drawing for a $150 gift card.

Comment thread src/struphy/feec/preconditioner.py
spossann and others added 4 commits October 8, 2026 13:56
Local diagonal of a KroneckerStencilMatrix: use the global row offset, so
that factors may own more rows than the codomain on this process.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
_process_local_matrix_1d stored entry (i, j) at diagonal (j - i + p) mod n.
For periodic directions with 2p + 1 > n (e.g. 1 element, degree 1) this put
the main diagonal at another equivalent offset; dot was unaffected, but the
diagonal of the Kronecker approximation (used for diagonal scaling since
MassMatrixDiagonalPreconditioner no longer assembles the 3d logical mass
matrix) was zero, so PCG returned NaN in
test_variational_mag_field_evolve. Periodic offsets are now taken nearest
to 0; non-periodic ones are not wrapped (the modulo was also wrong there
for 2p >= n).

Regression test with a (4, 4, 1) grid. CI shard 2 passes locally.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@spossann
spossann requested a review from max-models October 8, 2026 16:31
@spossann

spossann commented Oct 9, 2026

Copy link
Copy Markdown
Member Author

@max-models this is ready now!

@spossann
spossann merged commit 0432020 into devel Oct 9, 2026
28 checks passed
@spossann
spossann deleted the kron-stencil-matmul branch October 9, 2026 08:00
max-models added a commit that referenced this pull request Oct 9, 2026
…ss-preconditioner-cupy

Conflict in src/struphy/feec/preconditioner.py: #731's refactor (KroneckerPreconditioner base class,
module-level helpers, StiffnessPreconditioner) is kept; the CuPy changes of this branch are
re-applied on top of it:

- _kronecker_approximation, _solver_1d, is_circulant, FFTSolver: the dense 1d matrices and the 1d
  solvers are host (NumPy) setup data on every backend. KroneckerLinearSolver (feectools#96, in
  0.7.0) solves device data with dense inverses built from these host solvers, so FFTSolver no
  longer needs the host round trip of the previous version of this branch.
- _process_local_matrix_1d: band filled on the host, copied to the factor's data in one step;
  StiffnessPreconditioner passes its host 1d matrices directly.
- _local_diagonal: outer product by broadcasting instead of ufunc.outer (not available on every
  backend), so diagonal scaling (MassMatrixDiagonalPreconditioner) works on CuPy.

test_preconditioner_cupy.py also covers weight_reduction="average" with diagonal_scaling and
checks that applying the preconditioners makes no host/device transfers.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants