Stokes_Constrained is slow in parallel despite low KSP iterations

[gthyagi]

Summary

Stokes_Constrained works in parallel, but the 8-rank runtime is much slower than the equivalent Nitsche free-slip solve even though the constrained solve only takes two KSP iterations. This suggests the bottleneck is setup/assembly/field handling for the block-constrained system, not Krylov convergence.

This came up while validating the Zhong et al. spherical-shell internal-boundary benchmark path using SphericalShellInternalBoundary() and exterior free-slip boundaries.

Reproducer Context

External benchmark script:

cd /Users/tgol0006/uw_folder/uw3_git_gthyagi_latest/underworld3/.claude/worktrees/mantle-convection-benchmarks
pixi run -e amr-dev mpirun -np 8 python /Users/tgol0006/uw_folder/uw3-mantle-convection-benchmarks/benchmarks/005_internal_boundary_delta_probe.py

The script uses:

• uw.meshing.SphericalShellInternalBoundary()
• ri = 0.55, rint = 0.775, ro = 1.0
• cellSize = 0.125 under 8 MPI ranks
• velocity P2, pressure P1
• internal radial natural load on the Internal boundary
• exterior free slip on Upper and Lower

The Nitsche comparison is:

pixi run -e amr-dev mpirun -np 8 python /Users/tgol0006/uw_folder/uw3-mantle-convection-benchmarks/benchmarks/005_internal_boundary_delta_probe.py -uw_freeslip_type nitsche

Timing Evidence

Measured with /usr/bin/time -p on the same machine and mesh resolution:

Method	Wall time	KSP iterations	Result
Nitsche free slip	27.93s	1	pass
Stokes_Constrained	434.92s	2	pass
constrained with temporary degree-1 multiplier test	409.62s	2	pass

The temporary degree-1 multiplier test did not materially improve runtime, so multiplier polynomial degree is unlikely to be the main bottleneck.

Observed Metrics

8-rank constrained run:

{
  "cellsize": 0.125,
  "ksp_iterations": 2,
  "ksp_reason": 2,
  "l": 2,
  "max_boundary_area_relative_error": 0.009919085347124204,
  "max_y_l0_norm_error": 0.009849659575170255,
  "mpi_size": 8,
  "passed": true,
  "snes_reason": 5,
  "stokes_tolerance": 1e-05,
  "upper_characteristic_velocity": 0.012565603878091148,
  "upper_normal_velocity_rms": 1.631821482362003e-05
}

8-rank Nitsche comparison:

{
  "cellsize": 0.125,
  "ksp_iterations": 1,
  "ksp_reason": 2,
  "l": 2,
  "max_boundary_area_relative_error": 0.009919085347124204,
  "max_y_l0_norm_error": 0.009849659575170255,
  "mpi_size": 8,
  "passed": true,
  "snes_reason": 5,
  "stokes_tolerance": 1e-05,
  "upper_characteristic_velocity": 0.010110860387364697,
  "upper_normal_velocity_rms": 2.2416675430165106e-05
}

Current Diagnosis

The linear solve is not the issue: constrained free slip reports only 2 KSP iterations. The expensive part is likely one or more of:

• setup/assembly of the extra Lagrange-multiplier fields for Upper and Lower;
• boundary residual/Jacobian registration for the multiplier coupling;
• grouping pressure plus multipliers into the Schur block;
• _constrain_interior_multipliers_in_section() section work in parallel;
• creation or handling of full-domain multiplier fields when only the boundary trace is physical;
• repeated DM/section/fieldsplit setup that could be cached when mesh and constraints are unchanged.

Candidate Fix Directions

• Profile Stokes_Constrained setup and assembly with PETSc/UW timing to locate the exact hotspot.
• Optimize _constrain_interior_multipliers_in_section() if Python-side set/section operations dominate.
• Consider boundary-only or submesh multiplier fields instead of full-domain multiplier fields with interior DOFs constrained out.
• Cache constrained section/fieldsplit setup when the mesh, fields, and constraint boundaries are unchanged.
• Check whether grouped [p, lambda] Schur setup rebuilds too much state each solve.

Related Work

PR #242 fixes SphericalShellInternalBoundary() boundary labels so the internal-boundary benchmark can use the built-in mesh path directly. This issue is separate: after that fix, the constrained solve is correct but slow in parallel.

[gthyagi]

Follow-up from the Zhong 2008 benchmark profiling.

I re-ran the smallest strong reproducer, 005_internal_boundary_delta_probe.py, with PETSc -log_view -log_view_memory on the UW3 mantle-convection benchmark worktree.

Benchmark script:
https://github.com/gthyagi/uw3-mantle-convection-benchmarks/blob/main/benchmarks/005_internal_boundary_delta_probe.py

Docs/report:

• https://github.com/gthyagi/uw3-mantle-convection-benchmarks/blob/main/docs/constrained_freeslip_experiments.md
• https://github.com/gthyagi/uw3-mantle-convection-benchmarks/blob/main/docs/constrained_freeslip_performance.md

Key result for 005_internal_boundary_delta_probe.py, 8 ranks, -uw_cellsize 1/8:

case	wall time	PETSc total	SNESSolve	KSPSolve_FS_Schu	KSPGMRESOrthog	_setup_pointwise_functions	DMPlexJacobianFE	key metric
Nitsche	29.69 s	26.07 s	8.72 s	5.89 s	0.52 s	11.05 s	1.72 s	upper_characteristic_velocity = 1.011086e-02
constrained default	466.74 s	464.2 s	454.55 s	447.78 s	37.06 s	5.10 s	2.12 s	upper_characteristic_velocity = 1.256560e-02

Conclusion: the earlier observation that the outer constrained solve reports only 2 KSP iterations was incomplete. The runtime is dominated by the nested grouped [p,h] Schur fieldsplit solve, specifically KSPSolve_FS_Schu, not by setup/JIT or FE assembly.

I also tested several candidate solver options:

candidate	result	timing	metric impact	decision
fieldsplit_1_ksp_type preonly, fieldsplit_1_pc_type gasm	converged	18.60 s for 005 np8	upper_characteristic_velocity = 1.236730e-02, about 1.6% below default constrained	useful diagnostic, but not accepted as default
same preonly/GASM in 001_shell_stokes_smoke.py np8	converged	16.32 s	smoke metrics moved materially	rejected as correctness-preserving default
fieldsplit_1_ksp_max_it 20	converged	45.41 s for 005 np8	response changed by about 16%	rejected
fieldsplit_1_ksp_max_it 50	converged	88.33 s for 005 np8	response changed by about 19.5%	rejected
fieldsplit_1_ksp_type preonly, fieldsplit_1_pc_type lu	failed	abort after 17.76 s	PETSc bus error on rank 4	rejected
preonly/GASM plus outer ksp_norm_type unpreconditioned	converged	20.31 s	wrong velocity scale	rejected
preonly/GASM plus outer ksp_rtol 1e-10	converged	18.45 s	unchanged from preonly/GASM	rejected

So the immediate fast alternatives prove the bottleneck is the grouped [p,h] Schur sub-KSP, but they change benchmark observables. I reverted the attempted UW3 default change and rebuilt the worktree. Validation after the investigation: pixi run -e amr-dev pytest -q tests/test_1061_constrained_freeslip.py passed with 15 passed in 44.04s.

Next safe fix should target a correctness-preserving pressure-plus-multiplier Schur preconditioner or convergence strategy, rather than a simple cap or one-shot Schur PC application.

[lmoresi]

#244 — Free-slip parallel profile at np=8 (post-#306 three-way)

Setup

Identical mesh, variables, load, tolerance across the three free-slip
formulations:

• SphericalShellInternalBoundary(radiusOuter=1.0, radiusInternal=0.775, radiusInner=0.55, cellSize=1/8, qdegree=2, degree=1) — the test_1064 mesh
• velocity P2, pressure P1 (continuous)
• viscosity = 1, bodyforce = 0
• internal-boundary natural BC Y_l=2 · ê_r (Zhong l=2 load)
• free-slip on Upper and Lower
• tolerance 1e-7, snes_type = ksponly
• development @ 6451896a (post-fix(rotated-bc): Schur preconditioning parity + 3D rotation nullspace — Zhong #248 blow-out (44 its -> 1) #306), PETSc 3.25.2 (conda-forge), MPICH 5.0.1, np=8

Three "free-slip on Upper+Lower" variants:

variant	API	wall-normal imposition
nitsche	add_nitsche_bc(bdy, normal, gamma)	weak (penalty-consistent)
rotated	add_rotated_freeslip_bc(bdy, normal) (#298 + #306)	strong, per-node DOF rotation + Dirichlet on v_n — machine-zero
constrained	add_constraint_bc(bdy, g=0, normal, augmentation_base=1e4, degree=2)	in-saddle Lagrange multiplier λ on [p, λ_upper, λ_lower]

Headline numbers (np=8, test_1064 shell)

method	wall-clock .solve()	ratio vs Rotated	outer Krylov iters	s / outer iter	SNES reason
rotated	21.1 s	1.0×	1	—	0 (SNES bypassed — linear)
nitsche	30.7 s	1.5×	2	7.5 s	5
constrained	729.4 s	34.6×	5	142 s	5

Post-#306, rotated is now the fastest of the three — beating Nitsche
by 32 % on the same problem and delivering machine-zero wall-normal flow
(v·n = 0 to round-off, not the ~1e-3 leakage weak-Nitsche has).

Per-stage breakdown from -log_view

Wall-clock seconds, and share of the total run:

stage	Rotated (post-#306)	Nitsche	Constrained
SolverBaseClass._build	3.46 s (14 %)	10.37 s (30 %)	5.10 s (0.7 %)
↳ _setup_pointwise_functions (JIT weak forms)	2.87 s	9.73 s	4.59 s
↳ _setup_discretisation	0.08 s	0.10 s	0.42 s
↳ _setup_solver	0.55 s	0.58 s	0.46 s
SNESSolve	21.12 s (87 %)	20.31 s (60 %)	724.13 s (99 %)
↳ SNESJacobianEval / DMPlexJacobianFE	8.25 s	4.87 s	5.14 s
↳ KSPSolve (total Krylov)	8.74 s	14.98 s	718.27 s
↳ KSPSolve_FS_0 (velocity block)	0.32 s (1 iter)	0.56 s (2)	3.01 s (5)
↳ KSPSolve_FS_Schu (Schur block)	7.95 s (1 iter)	13.77 s (2)	710.71 s (5)
↳ KSPSolve_FS_Low (lower solve)	0.46 s (1 iter)	0.64 s (2)	4.51 s (5)

Reading the numbers

Setup / assembly is not the bottleneck for Constrained. Constrained's
_build is faster than Nitsche's (5.1 s vs 10.4 s), its JIT is faster (4.6
s vs 9.7 s), and its Jacobian assembly is comparable (5.1 s vs 4.9 s). The
earlier #244 hypotheses — "setup / assembly / section surgery /
_constrain_interior_multipliers_in_section() / repeated DM+section setup"
— are not where the ~700 s goes.

The bottleneck is unambiguously KSPSolve_FS_Schu on Constrained.
Attribution per Schur solve:

• Rotated: 1 Schur solve × ~7.95 s = 7.95 s
• Nitsche: 2 Schur solves × ~6.9 s = 13.8 s
• Constrained: 5 Schur solves × ~142 s = 710.7 s (97 % of the run)

Per Schur solve, Constrained's is 18× more expensive than Nitsche's and
Rotated's. That is where the 34× wall-clock gap lives.

Why Constrained's Schur block is expensive

The Stokes_Constrained Schur block is on the grouped field
[p, λ_upper, λ_lower], not just [p]. Each Schur application touches:

• the pressure block (same size as Nitsche's → cheap);
• the two multiplier blocks λ_upper, λ_lower — full-domain fields whose
  interior DOFs are zeroed by _constrain_interior_multipliers_in_section,
  but the block dimension stays full-domain-sized, not boundary-sized;
• the augmentation coupling augmentation_base=1e4 between p and each
  λ.

At 5 outer iters the inner KSPSolve for that combined block spends most of
its time in MatMult on the augmented [p, λ_u, λ_l] operator whose row
count is the full domain. Whether the multiplier physically lives only on the
boundary is not what the solver sees.

Two recommendations for #244

1. Recommend add_rotated_freeslip_bc as the free-slip default

Post-#306, add_rotated_freeslip_bc is the fastest free-slip option
here, and it delivers exact machine-zero v·n = 0 (Constrained's headline
guarantee, without the multiplier/topography machinery). At np=8 it is 32 %
faster than Nitsche and 34× faster than Constrained on this problem.

Users should reach for Stokes_Constrained only when they need
Stokes_Constrained's recoverable topography / Lagrange multiplier field
— not merely for "strong-imposition free-slip", which the rotated path now
covers cleanly. Consider updating the constrained-vs-nitsche how-to and the
Stokes_Constrained docstring accordingly. (CLAUDE.md already records this
preference — #300.)

2. Shrink the Stokes_Constrained Schur block to boundary-sized

The 97 % of run time inside KSPSolve_FS_Schu is dominated by MatMult on
the full-domain-sized augmented [p, λ_upper, λ_lower] block. Cleanest fix:
give the multiplier variable a codim-1 sub-DM so the Schur block dimension
tracks the boundary node count instead of the volume node count — a factor
~10–20 shrink on this shell mesh, which would move the KSPSolve_FS_Schu
column down from 710 s toward the tens of seconds range and bring
Constrained within 3–5× of Rotated/Nitsche.

The interior-multiplier zeroing that
_constrain_interior_multipliers_in_section already does could then be
retired as unnecessary — the sub-DM would guarantee the multiplier only
lives where it is physically defined.

Reference log_view files

• /tmp/uw3-verify/rotated_post306_np8_logview.txt
• /tmp/uw3-verify/nitsche_np8_logview.txt
• /tmp/uw3-verify/constrained_np8_logview.txt

All three built against development @ 6451896a, PETSc 3.25.2, MPICH
5.0.1, run on the same 48-core Linux host with only one mpirun in flight at
a time. Rotated wall-clock was re-measured after PR #306; Nitsche and
Constrained were not re-measured (neither is touched by #306) — earlier
measurements against a7f0b11 still apply.

[lmoresi]

Correction / full context on "recommendation 2" (codim-1 sub-DM) — check the premise first

The codim-1 sub-DM recommendation in my profiling comment above rests on a premise that
contradicts work that has already landed, and needs a direct measurement before anyone
invests in it.

What is already on development

The boundary-only multiplier reduction landed in PR #224 (2026-06-08) and was
parallel-audited in PR #265 (Item 3):

• _constrain_interior_multipliers_in_section()
  (src/underworld3/cython/petsc_generic_snes_solvers.pyx, ~L6925), default ON
  (_reduce_interior_multiplier = True).
• It constrains every interior (off-constraint-boundary) multiplier DOF directly in the
  fine DM's PetscSection, so interior λ DOFs drop out of the GLOBAL system — the
  assembled [p, λ_upper, λ_lower] block should already be pressure + boundary-trace λ
  only, not full-domain λ.
• PR Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265's audit confirmed the knockout is lossless and active in the 3-D np8 shell runs
  (turning it OFF re-admits the ill-conditioned ε·M interior block and degrades
  convergence — that was the "answer moves a few % when off" behaviour).

So my earlier statement that "the block dimension stays full-domain-sized, not
boundary-sized" is unverified and very likely wrong. One of the two is true:

• (a) The reduction is active on the shell at np8 (what the Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265 audit implies). Then a
  codim-1 sub-DM buys assembly/memory tidiness only — it cannot explain the 18×-per-Schur-solve
  gap, because the grouped Schur block contains the full-domain pressure either way, and
  that is the same size as Nitsche's Schur block, which is 18× cheaper per application. The
  cost must then live in inner-Schur iteration count / preconditioner quality: the grouped
  block's true operator is S_h = C·A⁻¹·Cᵀ + εM (the dynamic-topography boundary operator,
  augmentation-coupled at 1e4), which has no spectrally-equivalent PC. Prior measurements
  (PR Stokes_Constrained: in-saddle multiplier free-slip + topography (supersedes #219) #224 development notes): cheap [p,h] PCs are catastrophic — preonly+jacobi 112 s,
  bjacobi 134 s, lu(mumps) 203 s vs gasm-default 8.2 s on a 42k-DOF case, because a one-shot
  PC destroys the Schur approximation and the outer Krylov explodes. @gthyagi's preonly/GASM
  row in the table above is the same phenomenon from the other side.
• (b) The reduction silently no-ops in this configuration (3-D shell / np8 / this mesh
  path). Then the fix is a section/ordering bug in the existing reduction — much cheaper than
  building a sub-DM.

The decisive probe (cheap — please run this before anything else)

No 1/8-np8 run needed. On the shell at cellSize 1/4, serial and np≥2, with
ksp_max_it=1 (assembly + one iteration only):

# after the first solve (or after solver._setup_solver()):
fields, ises = zip(*[(n, i) for n, i in zip(*solver.dm.createFieldDecomposition()[:2])])
for name, iset in zip(fields, ises):
    print(name, iset.getSizes())   # (local, global) per field

Compare the multiplier fields' global sizes against (i) the boundary point count
(dm.getLabel("UW_Boundaries").getStratumIS(value).getSize() closure) and (ii) the same run
with solver._reduce_interior_multiplier = False. Expected if (a): λ global ≈ boundary trace,
and OFF−ON delta ≈ the ~26–45% DOF reductions recorded in the #224 validation. Expected if
(b): λ global ≈ volume node count, unchanged by the flag.

Also worth capturing while there: the inner Schur KSP iteration count
(-fieldsplit_1_ksp_monitor or -fieldsplit_pressure_ksp_monitor) for constrained vs
Nitsche — if (a) holds, that count is the whole story of the 34×.

Re-framing the remaining #244 work

• If (a): retire the codim-1 sub-DM idea (premise false). The open item becomes a real
  preconditioner for the augmented boundary Schur operator S_h — with the caveat that it
  must be an accurate iterative solve (the cheap-PC routes above are measured dead ends).
  A codim-1 sub-DM is anyway blocked today by two documented limitations: the UW3 JIT
  boundary residual/Jacobian is same-DM only (no MATNEST/DMComposite C/Cᵀ assembly), and the
  manifold vector-FE gap (F1.reshape(dim) vs cdim on embedded surfaces, see the Zhong 2008 isoviscous response: constrained vs Nitsche free-slip behavior #248
  investigation).
• If (b): the open item becomes "fix the section reduction on the 3-D shell / parallel
  path" and the sub-DM idea is superseded.

Independent of (a)/(b): recommendation 1 stands — add_rotated_freeslip_bc is the free-slip
default (#300), and σ_nn / dynamic topography is recoverable from the rotated constraint
reaction (solver.boundary_normal_traction()), so Stokes_Constrained performance matters
only for genuinely multiplier-needing use cases.

Underworld development team with AI support from Claude Code

[lmoresi]

Probe ran — premise settled: it's (a). The λ block is ALREADY boundary-sized; retire the codim-1 sub-DM recommendation

I ran the decisive probe from the previous comment. Result: the interior-multiplier section
reduction is active — serial and parallel, 2-D annulus and the actual 3-D
SphericalShellInternalBoundary mesh. The [p, λ] Schur block is already
pressure + boundary-trace multipliers. My earlier "the block dimension stays
full-domain-sized" claim was wrong, and recommendation 2 (codim-1 sub-DM) should be
retired — it would buy essentially nothing the section reduction hasn't already bought.

Measured (probe script below, development-equivalent build, PETSc 3.25)

2-D annulus, cellSize 0.1, free-slip Upper:

	velocity	pressure	multiplier_0	multiplier full-domain nodes
reduce ON (default), np1	2464	348	128	1296
reduce OFF, np1	2464	348	1296	1296
reduce ON, np2	2464	348	128	—

3-D shell (the #244 mesh family), cellSize 0.3, free-slip Upper+Lower, np2:

velocity	pressure	multiplier_0 (Upper)	multiplier_1 (Lower)	multiplier full-domain nodes
7650	385	790	222	~2760

The multipliers are the P2 boundary-surface traces, a small multiple of the pressure count
— and since λ scales as N² per surface while p scales as N³, the λ share of the grouped
block shrinks further at the benchmark's 1/8 resolution. Block size cannot explain an
18×-per-Schur-solve gap vs Nitsche (whose Schur block is the same full-domain pressure).

Where #244 actually stands

1. The remaining cost is the inner grouped-Schur solve quality, not DOF count: the
   grouped block's true operator is the augmentation-coupled boundary Schur complement
   S_h = C·A⁻¹·Cᵀ + εM, which currently has no spectrally-equivalent preconditioner, and
   every inner iteration is expensive. The measured record so far: cheap one-shot PCs are
   catastrophic (preonly+jacobi/bjacobi/lu all 15–25× worse — PR Stokes_Constrained: in-saddle multiplier free-slip + topography (supersedes #219) #224 dev notes, and
   @gthyagi's preonly/GASM row above is the same trade from the other side).
2. Next measurement (for whoever picks this up on the parallel benchmark machine):
   -fieldsplit_1_ksp_monitor inner iteration counts, constrained vs Nitsche, on the 1/8
   np8 benchmark. Expectation: that count IS the 34×. A correctness-preserving fix then
   means a real PC for S_h (e.g. a boundary mass matrix scaled by μ/h per surface — the
   natural spectral equivalent for a single-layer-type boundary operator — kept alongside
   the pressure-mass PC in the grouped split), not block shrinkage.
3. Recommendation 1 stands unchanged: rotated free-slip is the default; this only matters
   for genuinely multiplier-needing workflows.

Probe script (self-contained, validated on the runs above)

Run from any UW3 checkout; no external files needed.

pixi run -e amr-dev python issue244_schur_size_probe.py                      # 2D annulus, ON
pixi run -e amr-dev python issue244_schur_size_probe.py -uw_reduce 0        # OFF
pixi run -e amr-dev mpirun -np 8 python -u issue244_schur_size_probe.py \
    -uw_case shell -uw_cellsize 0.25                                        # the #244 mesh

Note: pypi gmsh 4.15.1 fails to mesh SphericalShellInternalBoundary at cellSize ≤ 0.25
("Identical points in triangulation") — use cellSize 0.3 or pin gmsh==4.13.1.

"""#244 decisive probe: is the [p, lambda] Schur block boundary-sized or full-domain-sized?

Prints the per-field GLOBAL DOF counts of the assembled Stokes_Constrained
field decomposition, with the interior-multiplier section reduction ON or OFF.
"""

import sympy
import underworld3 as uw

params = uw.Params(
    uw_case="annulus",     # "annulus" (2D, seconds) | "shell" (3D SphericalShellInternalBoundary)
    uw_cellsize=0.1,       # 0.1 annulus; 0.3 shell (gmsh 4.15.1 trap below 0.3; 4.13.1 fine)
    uw_reduce=True,        # interior-multiplier section reduction (default ON in the solver)
    uw_tolerance=1e-2,     # loose: we only need setup + a converged-enough solve, not physics
)

MU = 1.0

if params.uw_case == "shell":
    mesh = uw.meshing.SphericalShellInternalBoundary(
        radiusOuter=1.0, radiusInternal=0.775, radiusInner=0.55,
        cellSize=params.uw_cellsize, qdegree=2, degree=1,
    )
    coords = list(mesh.X)
    r = sympy.sqrt(sum(c**2 for c in coords))
    unit_r = sympy.Matrix([[c / r for c in coords]])
    buoy = sympy.cos(2 * sympy.atan2(coords[1], coords[0]))
    constraint_boundaries = ["Upper", "Lower"]
    dirichlet = None  # fully free-slip, like the #244 benchmark
else:
    mesh = uw.meshing.Annulus(radiusInner=0.5, radiusOuter=1.0,
                              cellSize=params.uw_cellsize, qdegree=3)
    x, y = mesh.X
    r = sympy.sqrt(x**2 + y**2)
    unit_r = sympy.Matrix([[x / r, y / r]])
    buoy = 1.0e2 * sympy.cos(3 * sympy.atan2(y, x)) * (r - 0.5) / 0.5
    constraint_boundaries = ["Upper"]
    dirichlet = "Lower"

blk = uw.systems.Stokes_Constrained(mesh)
blk.constitutive_model = uw.constitutive_models.ViscousFlowModel
blk.constitutive_model.Parameters.shear_viscosity_0 = MU
blk.bodyforce = buoy * unit_r
if dirichlet is not None:
    blk.add_dirichlet_bc((0.0,) * mesh.dim, dirichlet)

multipliers = {b: blk.add_constraint_bc(b, g=0.0, normal=unit_r)
               for b in constraint_boundaries}

blk.petsc_use_pressure_nullspace = True
blk.tolerance = params.uw_tolerance
blk._reduce_interior_multiplier = bool(params.uw_reduce)

blk.solve()

# ---- The measurement: per-field GLOBAL sizes of the assembled system ------ #
names, ises, _dms = blk.dm.createFieldDecomposition()

uw.pprint(f"\ncase={params.uw_case} cellsize={params.uw_cellsize} "
          f"np={uw.mpi.size} reduce={bool(params.uw_reduce)}")
sizes = {}
for name, iset in zip(names, ises):
    _local, global_sz = iset.getSizes()
    sizes[name] = global_sz
    uw.pprint(f"  field {name:<14s} global DOFs = {global_sz}")

p_sz = sizes.get("pressure")
for name, gsz in sizes.items():
    if name.startswith("multiplier") and p_sz:
        uw.pprint(f"  {name}: {gsz / p_sz:.3f} x pressure "
                  f"(boundary trace => small multiple; full-domain P2 => ~volume-sized)")

# Full-domain reference: total nodes of one multiplier variable (local sums).
h0 = next(iter(multipliers.values()))
n_local = h0.coords.shape[0]
n_global = uw.mpi.comm.allreduce(n_local)  # over-counts shared nodes slightly; scale reference only
uw.pprint(f"  multiplier variable full-domain nodes (sum of local, incl. shared) ~ {n_global}")

Underworld development team with AI support from Claude Code

[lmoresi]

#244 — Definitive: the grouped [p, λ] Schur PC is not the right PC (inner Schur hits max_it every outer iter)

Following @lmoresi's premise-check probe: the multiplier section reduction is
active on this benchmark, [p, λ] is boundary-sized, and the sub-DM path is
retired. The remaining question was whether the 34× wall-clock gap is
explained by inner Schur iteration count, not DOF count. It is, and more
sharply than expected.

Setup

Same as REPORT_244 (post-#306): SphericalShellInternalBoundary(radiusOuter=1.0, radiusInternal=0.775, radiusInner=0.55, cellSize=1/8, qdegree=2, degree=1),
internal Y_l=2 · ê_r load, free-slip on Upper+Lower, tol 1e-7, np=8,
snes_type = ksponly, PETSc 3.25.2, MPICH 5.0.1.

Both variants use fgmres on the outer Krylov and on every inner Schur block,
with identical GASM preconditioner defaults on the pressure-side block
(fieldsplit_*_ksp_max_it = 200, fieldsplit_*_ksp_rtol = 1e-4,
fieldsplit_*_pc_type = gasm). The only structural difference is the
grouped block: Constrained's Schur is on [p, λ_upper, λ_lower]; Nitsche's is
on [p] alone. The multipliers are already boundary-sized here (@lmoresi's
probe: annulus 128 vs 348 pressure, shell 790+222 vs 385 pressure), so the
block-size ratio is small.

Inner Schur iteration counts

Ran probe_244_inner_iters.py with fieldsplit_*_ksp_converged_reason
on. Distinct outcomes:

variant	inner Schur fgmres/gasm behaviour	outer FGMRES	wall
Nitsche	CONVERGED_RTOL iterations 24	1 iter, CONVERGED_RTOL	32.88 s
Constrained	DIVERGED_ITS iterations 200 (every outer application)	5 iters, CONVERGED_RTOL	742.26 s

Direct from probe_244_inner_iters.py — the final sub[i] summary from
snes.getKSP().getPC().getFieldSplitSubKSP() after solve:

# Nitsche
sub[0]_ksp_type=fgmres pc_type=gamg its_LAST=65   reason_LAST=2    (CONVERGED_RTOL)
sub[1]_ksp_type=fgmres pc_type=gasm its_LAST=24   reason_LAST=2    (CONVERGED_RTOL)   ← [p]

# Constrained
sub[0]_ksp_type=fgmres pc_type=gamg its_LAST=109  reason_LAST=2    (CONVERGED_RTOL)
sub[1]_ksp_type=fgmres pc_type=gasm its_LAST=200  reason_LAST=-3   (DIVERGED_ITS)     ← [p, λ_u, λ_l]

Nitsche's inner Schur converges tightly. Constrained's inner Schur does not
converge at all — it hits the fieldsplit_1_ksp_max_it = 200 ceiling on
every outer iteration. The outer FGMRES tolerates that (it's flexible; a
non-converged inner counts as a "variable preconditioner"), which is why the
solver still lands on a correct answer — but each outer step burns a full 200
inner FGMRES iterations, and the outer needs 5 of them.

That matches the KSPGMRESOrthog / KSPSolve_FS_Schu cross-check from the
existing log_view files:

	Nitsche	Constrained	ratio
KSPGMRESOrthog count	3,395	168,138	49.5×
KSPSolve_FS_Schu wall (s)	13.77	710.71	51.6×

The two ratios agreeing to 4% confirms the 34× overall wall-clock gap is
entirely explained by iteration count on the grouped Schur block, at
essentially identical per-iteration cost. No DOF-size, no assembly, no JIT.

Diagnosis

The GASM preconditioner is spectrally OK for the pure pressure Schur
complement in Nitsche but is a very poor PC for the augmented boundary
Schur complement S_h = C · A⁻¹ · Cᵀ + εM that appears when the multipliers
are grouped with pressure in Constrained. The measured record on cheaper
one-shot PCs is uniformly worse (from #224's dev notes and @gthyagi's
preonly/GASM row earlier in this thread — either much slower or wrong
physics). So the outer solver keeps a pretty good variable-PC-tolerant
FGMRES and pays the ~30-40s per outer iter cost of grinding the inner GMRES
to max_it.

Where the fix actually lives (per @lmoresi's earlier comment)

A correctness-preserving PC for S_h — a boundary mass matrix scaled by μ/h
per constrained surface is the natural spectral equivalent for a
single-layer-type operator, kept alongside the pressure-mass PC in the
grouped split (not replacing it). That means composing per-field PCs inside
the grouped block, not swapping the whole fieldsplit_1_pc_type.

Everything else is a red herring:

• raising fieldsplit_1_ksp_max_it past 200 — does not help (the residual
  stagnates before 200 on this problem);
• dropping fieldsplit_1_ksp_rtol — same;
• one-shot PCs (preonly + jacobi / bjacobi / lu / gasm) — either much slower
  or wrong physics (all measured in Stokes_Constrained: in-saddle multiplier free-slip + topography (supersedes #219) #224 dev notes and above);
• codim-1 λ sub-DM — retired (block is already boundary-sized).

Recommendation 1 from REPORT_244 still stands: add_rotated_freeslip_bc is
the fast default (#300). Stokes_Constrained performance only matters when
you genuinely need the multiplier field / dynamic-topography reaction that
Rotated already exposes.

Reference

• /tmp/uw3-verify/probe_244_inner_iters.py — the monitor script
• /tmp/uw3-verify/probe244_nitsche.log — nitsche summary
• /tmp/uw3-verify/probe244_constrained.log — constrained summary
• Built against development @ 6451896a, PETSc 3.25.2, MPICH 5.0.1, np=8.

Underworld development team with AI support from Claude Code

[lmoresi]

Cross-cutting: 3D null-space attach + preconditioner-parity recurs in three fixes — #244's inner Schur has the same signature

Two recent fixes and the still-open #244 all share a diagnostic
signature: a Krylov solve that "runs to max_it without converging" on
a 3D problem, whose actual root cause is a genuine near-null-space mode
of the operator being iterated on that PETSc has not been told about.
Preconditioner-quality work only helps once the null space is attached.
Recording the pattern so we can spot it earlier next time.

The three cases

fix	operator	null modes	pre-fix behaviour	fix
#265 (Stokes_Constrained parallel correctness)	outer coupled Jacobian on [u, p, λ]	rigid-rotation of u (SO(3) in 3D) + combined-constant (p+1, λ-1) gauge	outer gmres false-converged on preconditioned norm — variable-PC-tolerant fgmres + unpreconditioned norm required; rigid rotations left partition-dependent tangential u	attach null spaces + project rigid rotations out of the converged solution (_remove_velocity_rotation_gauge), use fgmres + ksp_norm_type unpreconditioned
#306 (rotated free-slip)	rotated saddle  = Q A Qᵀ	constant p + rigid rotations Q·(SO(3)) in 3D	outer FGMRES hit ksp_max_it = 300 on a 3D shell, residual stagnated at 2.7e-7; 2D annulus converged in 19 iters (because 2D nullspace IS attached)	attach all rigid-rotation generators (1 in 2D, up to 3 in 3D) + move to native 1/μ mass USER Schur PC
#244 (Stokes_Constrained slow)	inner grouped-Schur block S_h = C·A⁻¹·Cᵀ + εM (fieldsplit_1) on [p, λ_upper, λ_lower]	combined-constant (p+1, λ_u-1, λ_l-1) restricted to the Schur block	inner FGMRES hits fieldsplit_1_ksp_max_it = 200 on every outer iteration (DIVERGED_ITS) at np=8, cellSize=1/8	not fixed yet

The Constrained solver's _build_stokes_nullspace (petsc_generic_snes_solvers.pyx:6022)
correctly constructs the combined mode via
_build_block_gauge_nullspace_vector (line 5989 — p = +1, λ_upper = λ_lower = -1, u = 0 — orthogonal to the u-row functional because
(+1 - 1)·∫n·δu = 0). And _attach_stokes_nullspace (line 6100)
attaches it. But — critically — only to the outer coupled Jacobian
at lines 6132–6133:

operator_matrix.setNullSpace(nullspace)
operator_matrix.setTransposeNullSpace(nullspace)

Nothing then reaches inside the fieldsplit to attach the restriction
of that mode to the fieldsplit_1 sub-block operator. The inner Krylov
running on that sub-block therefore has no way to reduce residual in
its combined-constant direction — which is exactly the signature the
probe shows: DIVERGED_ITS iterations 200 at every outer application.

That is the same "outer Krylov cannot touch un-attached null-space
components" pathology that PR #306 diagnosed and fixed for the rotated
path in 3D. There it manifested at the outer solve; here it
manifests one level down at the inner Schur solve.

Why Nitsche's inner Schur does not show this

Nitsche's inner Schur block is on [p] alone. Its natural null mode is
p = constant — a literal constant vector, which PETSc handles
implicitly when the outer pressure null space is set (or the Schur PC
selects a preconditioner that already respects it). No explicit
setNullSpace on the inner Schur operator is needed for the
plain-pressure case.

Constrained's inner Schur block on [p, λ_upper, λ_lower] has a
null mode whose per-field entries are not all +1 (p = +1,
λ = -1). This is not a "constant null space" in PETSc's implicit
sense; it must be attached explicitly as a vector. Currently it is
not — and GASM, like any other PC, cannot reduce residual in the
un-attached direction, so the inner FGMRES coasts to max_it.

Preconditioner-quality is downstream

The follow-up preconditioner-quality work (e.g. Louis's suggested
per-boundary μ/h mass composed alongside the pressure-mass PC, as
noted in an earlier #244 comment)
is real and independently valuable. But — as the pre-#306 rotated case
showed — no amount of PC-quality work makes a Krylov solver reduce
residual in a null-space direction the operator has not been told
about. The order is: attach the null space first, then measure
preconditioner quality against the residual outside it.

Testable predictions

If the "unattached inner-Schur null space" hypothesis is right:

1. Attaching the combined-constant vector as an explicit
   MatNullSpace on the fieldsplit_1 sub-operator (after pc.setUp)
   should drop the inner FGMRES from DIVERGED_ITS 200 to a few tens
   of iterations under the same GASM PC.
2. The outer FGMRES iteration count should not increase (the outer
   already knows the coupled null space; the fix is purely at the
   inner level).
3. In 2D the effect should be milder — the annulus / box constrained
   tests show 2D convergence patterns closer to Nitsche's already.

I can run experiment 1 on this machine if it would be useful.

Where the fix lives

_attach_stokes_nullspace is the natural place: after attaching the
coupled null space to the outer Jacobian, also — when
self._block_constraint_bcs is non-empty AND the PC is fieldsplit —
attach the block-restriction of the combined-constant vector to the
fieldsplit_1 sub-operator. Timing constraint: this has to happen
after pc.setUp() builds the sub-KSPs but before the next outer
ksp.solve().

If the experiment confirms the prediction, the fix mirrors the
"attach the null spaces properly" half of #306, one level deeper in
the solver stack.

Underworld development team with AI support from Claude Code

[lmoresi]

#244 update — null-space attach on the inner Schur does NOT resolve the

stagnation. Pattern analysis was too quick — the diagnosis is preconditioner
quality, not null-space attach.

What I tested

Predicted in the previous comment that
attaching the combined-constant vector (p=+1, λ_upper=-1, λ_lower=-1) as
an explicit MatNullSpace on the fieldsplit_1 sub-operator would let the
inner FGMRES converge — the same pattern PR #306 used for the pre-fix
Rotated case.

Ran two consecutive solves of the same test_1064 setup (np=8, PETSc
3.25.2, MPICH 5.0.1) in one process:

• SOLVE 1 (baseline) — reproduces the reported behaviour exactly.
• SOLVE 2 (after attach) — monkey-patched to call
  A_S.setNullSpace(nsp) / A_S.setTransposeNullSpace(nsp) on the
  fieldsplit_1 sub-operator between the two solves, where nsp is
  built by scattering _build_block_gauge_nullspace_vector() into the
  fieldsplit_1 sub-space via the composite pressure+multipliers IS.
  Vector norm before normalisation |v| = 89.58, consistent with
  √(N_pressure + N_multiplier_boundary_trace) on the 1/8 mesh — so
  the vector is being constructed correctly.

Result

	SOLVE 1 (BASELINE)	SOLVE 2 (WITH-INNER-NSP-ATTACHED)
wall (s)	788.33	1217.59 (+54 %)
outer_ksp_its	5	8 (+60 %)
outer_ksp_reason	2 (CONVERGED_RTOL)	2
sub[0] velocity (fgmres/gamg) its_LAST	109, reason=2	88, reason=2
sub[1] Schur (fgmres/gasm) its_LAST	200, reason=-3 (DIVERGED_ITS)	200, reason=-3 (DIVERGED_ITS)

The attach ran and MatSetNullSpace accepted the vector, but the inner
FGMRES on fieldsplit_1 behaves identically — same DIVERGED_ITS 200
per outer application. And the outer FGMRES ends up needing more
iterations (5 → 8) to reach the same rtol, making SOLVE 2 measurably
slower than the baseline. No trace of the "inner iter count collapses"
prediction.

Why the pattern didn't transfer

I over-fit the pattern from PR #306. In that case the outer Krylov was
solving the operator we handed it directly, so MatNullSpaceRemove
(which PETSc invokes on user RHS) actually projected the null-space
component out at solve entry.

On the fieldsplit inner solve, the picture differs on two counts:

1. The inner RHS is not a user vector. The fieldsplit PC computes
   the Schur-block residual at each outer FGMRES iteration and hands
   it to the inner KSP. PETSc's MatNullSpaceRemove on the operator
   is invoked at RHS projection time, but the fieldsplit-generated
   RHS is a residual (already in the operator's range in principle)
   whose null-space component is either genuinely zero or already
   handled by the outer coupled nullspace (which we DID attach). So
   the inner setNullSpace has nothing left to do.

2. GASM is not near-null-space aware. Additive Schwarz doesn't
   consume the attached null space the way GAMG (near-null-space
   basis for coarsening) or LU (kernel-of-the-factorisation) would.
   The attach just sits on the operator without changing what GASM
   does per iteration.

3. The attached vector isn't exactly a null mode. The
   _build_block_gauge_nullspace_vector mode (p=+1, λ=-1) is only
   an approximate null mode of the augmented boundary Schur
   S_h = C·A⁻¹·Cᵀ + εM — the augmentation term εM acts on the
   p-block-only part with ε=1e-4·μ and on the λ-block-only part
   with ε=1/augmentation_base = 1e-4. So S_h · v = O(ε)·v ≠ 0
   pointwise; PETSc treats the vector as a genuine null-mode
   direction and its projection instead frustrates a legitimate
   residual component. That is likely the "SOLVE 2 needs more outer
   iterations" pathology.

What this leaves standing

The definitive iteration-count measurement from
the earlier comment
is unchanged: the grouped Schur block fgmres/gasm hits
DIVERGED_ITS 200 at every outer iteration; the plain-pressure Schur
under the identical GASM PC converges in 24 iters (CONVERGED_RTOL).
Same PC, structurally similar block sizes, very different behaviour.

The remaining explanation is preconditioner quality on the augmented
boundary Schur operator S_h = C · A⁻¹ · Cᵀ + εM. GASM does not
capture its spectral structure and the pattern PR #306 used for the
rotated case does not port down into the fieldsplit interior. So:

Recommendation 2 becomes — design an actual PC for S_h. The most
natural candidate is a boundary mass matrix scaled by μ/h per
constrained surface, kept alongside the pressure-mass PC in the
grouped-Schur block via pc_composite — Louis's earlier
comment
sketched this and it is the spectrally-equivalent choice for a
single-layer boundary operator. That is a real piece of solver-design
work, not a session slice.

Recommendation 1 (rotated free-slip is the default; #300) is
unaffected — add_rotated_freeslip_bc remains the fastest free-slip
option and covers the "strong-imposition, no need for the multiplier
field" use cases without touching this bottleneck.

What NOT to spend time on (this session's negative results)

• Attaching the combined-constant nullspace on fieldsplit_1 via
  A_S.setNullSpace — measured, no effect.
• Raising fieldsplit_1_ksp_max_it past 200 — @gthyagi already tried
  this (see earlier Stokes_Constrained is slow in parallel despite low KSP iterations #244 comments), residual stagnates before 200.
• Codim-1 λ sub-DM — retired earlier in this thread (the section
  reduction is already active; the block is already boundary-sized).

Underworld development team with AI support from Claude Code