From afad6de46a2f213d7ea758095e9fbd9fce32561d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 16 Aug 2026 19:51:05 +1000 Subject: [PATCH 1/2] Preconditioner defaults: PETSc's GAMG cycle, node aggregation, flexible outer Three defaults that deviated from what the design requires, none carrying a recorded reason. Closes #579 and #576. GAMG cycle. pc_mg_type was "additive" where PETSc's default for GAMG is "multiplicative". Additive applies the levels independently and sums the corrections, so no level sees what another has already removed. git log -S puts its arrival in December 2022, in example scripts, alongside pc_mg_levels=2 -- with one coarse level the two cycles nearly coincide, so it cost nothing where it was written. It was promoted into the solver defaults and the hierarchies got deeper. SolKz at 11 727 unknowns: 243-347 cycles per velocity solve against 74-118. Node aggregation. The velocity operator reached GAMG with block size 1, so it aggregated scalar degrees of freedom rather than nodes. PETSc infers this correctly wherever it can -- with no velocity BCs, or a full vector Dirichlet condition, the field IS and the fieldsplit sub-matrix both come out with block size 2. It collapses to 1 under a COMPONENT-WISE Dirichlet condition, because the field then carries one degree of freedom at a constrained wall node and two inside. That is free slip, so the case PETSc cannot infer is the ordinary one. Setting mat_block_size cannot change the answer -- a preconditioner alters the route and not the fixed point, and the two configurations agree to 2e-7 relative on a 1e-6 solve -- but on free-slip SolKz it is worth 3.48e9 flops against 1.60e9, and with the contrast concentrated in a band rather than spread smoothly, 2.61e11 against 4.44e9. The block size is recorded on the solver as _pc_block_size rather than only at the __init__ call sites, because _apply_preconditioner_options re-applies the bundle on every build and a bundle built without it lists mat_block_size as a stale key and deletes what __init__ set. SNES_MultiComponent keeps block size 1: its unknown is not mesh.dim components in general. Flexible outer Krylov. The Stokes outer ksp_type was never set, so it inherited PETSc's gmres while both fieldsplit sub-blocks are fgmres run to a tolerance. Inside the Schur factorisation the velocity sub-solve computes the search direction, so the operator the outer method applies differs between outer iterations, and GMRES's residual recurrence assumes it does not. Spiegelman notch at refinement 3: 983 velocity iterations per step and DIVERGED_LINEAR_SOLVE, against 58 for fgmres and 23 with a loosened inner tolerance. Raising the velocity iteration cap is inert -- byte-identical residuals -- so the failure is inconsistency, not a shortage of iterations. Set in __init__ and in the strategy setter, which also resets the fieldsplit block. The four call sites that open-coded the GAMG bundle key by key now read it from multigrid_options, which is that module's stated reason for existing and is how one value came to be wrong in four places at once. test_0203 asserted that a fresh solver's first velocity count undercounts the second by a strict margin. The contract is that the first count is a LOWER BOUND and says so; how short it falls depends on where PCSetUp lands relative to the first application of the block, which these defaults move. It is now tight on that problem, so the test asserts the bound rather than the size of the gap. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 83 ++++++++++++------- .../utilities/multigrid_options.py | 65 +++++++++++++-- tests/test_0203_solver_wallclock_guard.py | 20 +++-- 3 files changed, 123 insertions(+), 45 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 3ff11af0..59427964 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -222,6 +222,15 @@ class SolverBaseClass(uw_object): # untouched framework default (eligible for FMG upgrade) apart from an # explicit user override of pc_type, which must be respected. self._pc_managed_value = "gamg" + # Components per node in the field the managed block solves. Subclasses + # whose unknown is a vector set mesh.dim; 1 is right for a scalar. The + # GAMG bundle turns this into `mat_block_size` so that algebraic + # coarsening aggregates nodes rather than scalars — worth a factor of + # two to nine (#579). It has to live here rather than only at the + # __init__ call sites because _apply_preconditioner_options re-applies + # the bundle on EVERY build, and a bundle built without it would list + # `mat_block_size` as a stale key and delete what __init__ had set. + self._pc_block_size = 1 # Latches once the user (or harness) is seen to have set the PC options # themselves. _apply_preconditioner_options() runs on EVERY _build (so # "auto" can re-resolve after a remesh), and without this latch the @@ -915,9 +924,10 @@ class SolverBaseClass(uw_object): f"mesh with refinement >= 1 to enable geometric multigrid.", stacklevel=2, ) - multigrid_options.gamg_bundle().apply( - PETSc.Options(), self.petsc_options_prefix + prefix, - owned=self._managed_pc_options) + multigrid_options.gamg_bundle( + block_size=self._pc_block_size).apply( + PETSc.Options(), self.petsc_options_prefix + prefix, + owned=self._managed_pc_options) self._pc_managed_value = "gamg" def _enforce_galerkin_for_geometric_mg(self): @@ -3254,16 +3264,10 @@ class SNES_Scalar(SolverBaseClass): self.petsc_options["snes_type"] = "newtonls" self._push_managed_option("ksp_type", "gmres") - self._push_managed_option("pc_type", "gamg") - self._push_managed_option("pc_gamg_type", "agg") - self._push_managed_option("pc_gamg_repartition", True) - self._push_managed_option("pc_mg_type", "additive") - self._push_managed_option("pc_gamg_agg_nsmooths", 2) - self._push_managed_option("mg_levels_ksp_max_it", 3) - self._push_managed_option("mg_levels_ksp_converged_maxits", None) + for key, value in multigrid_options.gamg_bundle().settings.items(): + self._push_managed_option(key, value) self.petsc_options["snes_rtol"] = 1.0e-4 - self._push_managed_option("mg_levels_ksp_max_it", 3) if self.verbose == True: self.petsc_options["ksp_monitor"] = None @@ -4167,14 +4171,13 @@ class SNES_Vector(SolverBaseClass): self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_rtol"] = 1.0e-3 self._push_managed_option("ksp_type", "gmres") - self._push_managed_option("pc_type", "gamg") - self._push_managed_option("pc_gamg_type", "agg") - self._push_managed_option("pc_gamg_repartition", True) - self._push_managed_option("pc_mg_type", "additive") - self._push_managed_option("pc_gamg_agg_nsmooths", 2) + # A vector unknown: mesh.dim components per node, which GAMG has to be + # told. See multigrid_options._gamg_settings. + self._pc_block_size = self.mesh.dim + for key, value in multigrid_options.gamg_bundle( + block_size=self._pc_block_size).settings.items(): + self._push_managed_option(key, value) self.petsc_options["snes_rtol"] = 1.0e-3 - self._push_managed_option("mg_levels_ksp_max_it", 3) - self._push_managed_option("mg_levels_ksp_converged_maxits", None) if self.verbose == True: self.petsc_options["ksp_monitor"] = None @@ -5183,14 +5186,13 @@ class SNES_MultiComponent(SolverBaseClass): self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_rtol"] = 1.0e-3 self._push_managed_option("ksp_type", "gmres") - self._push_managed_option("pc_type", "gamg") - self._push_managed_option("pc_gamg_type", "agg") - self._push_managed_option("pc_gamg_repartition", True) - self._push_managed_option("pc_mg_type", "additive") - self._push_managed_option("pc_gamg_agg_nsmooths", 2) + # Block size left at 1. The unknown here is a general multi-component + # field whose components-per-node is not mesh.dim in general, and + # claiming the wrong node size would cost rather than save. Worth + # revisiting per instantiation — see #579. + for key, value in multigrid_options.gamg_bundle().settings.items(): + self._push_managed_option(key, value) self.petsc_options["snes_rtol"] = 1.0e-3 - self._push_managed_option("mg_levels_ksp_max_it", 3) - self._push_managed_option("mg_levels_ksp_converged_maxits", None) if self.verbose == True: self.petsc_options["ksp_monitor"] = None @@ -6025,6 +6027,20 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options["snes_ksp_ew"] = None self.petsc_options["snes_ksp_ew_version"] = 3 + # The outer Krylov must be FLEXIBLE, because both sub-blocks below are + # themselves Krylov solves run to a tolerance. Inside the Schur + # factorisation the velocity sub-solve computes the search direction, so + # the operator the outer method applies differs from one outer iteration + # to the next, and plain GMRES's residual recurrence assumes it does not. + # PETSc's default is `gmres`, and that default is invisible on easy + # problems -- isoviscous SolKz converges in two outer iterations -- and + # expensive on hard ones: on the Spiegelman notch at refinement 3, 983 + # velocity iterations per step and DIVERGED_LINEAR_SOLVE, against 58 for + # fgmres and 23 with a loosened inner tolerance. Raising the velocity + # iteration cap changes nothing (byte-identical residuals), so the failure + # is inconsistency rather than too few iterations. See #576. + self.petsc_options["ksp_type"] = "fgmres" + self.petsc_options["pc_type"] = "fieldsplit" self.petsc_options["pc_fieldsplit_type"] = "schur" self.petsc_options["pc_fieldsplit_schur_fact_type"] = "full" # diag is an alternative (quick/dirty) @@ -6067,13 +6083,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options[f"fieldsplit_{v_name}_ksp_type"] = "fgmres" self.petsc_options[f"fieldsplit_{v_name}_ksp_max_it"] = 200 self.petsc_options[f"fieldsplit_{v_name}_ksp_rtol"] = self._tolerance * 0.1 - self._push_managed_option(f"fieldsplit_{v_name}_pc_type", "gamg") - self._push_managed_option(f"fieldsplit_{v_name}_pc_gamg_type", "agg") - self._push_managed_option(f"fieldsplit_{v_name}_pc_gamg_repartition", True) - self._push_managed_option(f"fieldsplit_{v_name}_pc_mg_type", "additive") - self._push_managed_option(f"fieldsplit_{v_name}_pc_gamg_agg_nsmooths", 2) - self._push_managed_option(f"fieldsplit_{v_name}_mg_levels_ksp_max_it", 3) - self._push_managed_option(f"fieldsplit_{v_name}_mg_levels_ksp_converged_maxits", None) + # The velocity field carries mesh.dim components per node; say so, or + # GAMG aggregates scalars. See multigrid_options._gamg_settings. + self._pc_block_size = self.mesh.dim + for key, value in multigrid_options.gamg_bundle( + block_size=self._pc_block_size).settings.items(): + self._push_managed_option(f"fieldsplit_{v_name}_{key}", value) # Create this dict self.fields = {} @@ -6860,6 +6875,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options["snes_ksp_ew"] = None self.petsc_options["snes_ksp_ew_version"] = 3 + # Flexible for the same reason as in __init__ (#576): the sub-blocks are + # inexact Krylov solves, so the outer operator varies between iterations. + self.petsc_options["ksp_type"] = "fgmres" + self.petsc_options["pc_type"] = "fieldsplit" self.petsc_options["pc_fieldsplit_type"] = "schur" self.petsc_options["pc_fieldsplit_schur_fact_type"] = "full" # diag is an alternative (quick/dirty) diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index 10a5c491..3651cccc 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -281,23 +281,65 @@ def _geometric_mg_settings(coarse, smoother="robust"): return settings -def _gamg_settings(): +def _gamg_settings(block_size=1): """The algebraic-multigrid (GAMG) settings — the fallback whenever no - geometric hierarchy is available.""" - return { + geometric hierarchy is available. + + ``block_size`` is the number of components per node in the field this block + solves. See :func:`gamg_bundle`. + """ + settings = { "pc_type": "gamg", "pc_gamg_type": "agg", "pc_gamg_repartition": True, - "pc_mg_type": "additive", + # PETSc's own default cycle for GAMG, set rather than inherited so that + # a block previously configured for geometric MG (which sets "full") + # cannot leak its cycle through setFromOptions. + # + # This was "additive" from 2022 until #579. Additive multigrid applies + # the levels independently and sums the corrections, so no level sees + # what another has already removed. It entered alongside + # `pc_mg_levels = 2`, where there is one coarse level and the two cycles + # nearly coincide, and stayed as the hierarchy got deeper and they did + # not. Measured on SolKz at 11 727 unknowns, additive cost 243-347 + # cycles per velocity solve against 74-118 multiplicative. + "pc_mg_type": "multiplicative", "pc_gamg_agg_nsmooths": 2, "mg_levels_ksp_max_it": 3, "mg_levels_ksp_converged_maxits": None, } + if block_size > 1: + # Tell GAMG the field has `block_size` components per node, so that it + # aggregates NODES rather than individual scalar degrees of freedom. + # + # PETSc works this out for itself wherever it can, and does: with no + # velocity boundary conditions, or with a full vector Dirichlet + # condition, the field IS and the fieldsplit sub-matrix both come out + # with block size 2. It collapses to 1 under a COMPONENT-WISE Dirichlet + # condition, because the field then genuinely carries one degree of + # freedom at a constrained wall node and two inside, and the matrix is + # not uniformly node-blocked. Block size 1 is the correct description + # there -- and component-wise Dirichlet is free slip, which is the + # standard geodynamics boundary condition, so the case PETSc cannot + # help with is the ordinary one. + # + # Setting it anyway asks GAMG to aggregate on a pairing that does not + # align with nodes at the walls. That cannot change the answer: a + # preconditioner alters the route to the solution and not the solution + # itself, and measured on free-slip SolKz the two configurations agree + # to 2e-7 relative, well inside a 1e-6 solve. What it changes is the + # cost. Free-slip SolKz at 11 727 unknowns: 74-118 cycles per velocity + # solve and 3.48e9 flops against 34-47 and 1.60e9. With the viscosity + # contrast concentrated in a band rather than spread smoothly it is + # larger still -- 2.61e11 flops against 2.97e10, and 32.7 s against + # 4.9 s. + settings["mat_block_size"] = block_size + return settings def _all_keys(): """Every option key any bundle here owns — the basis for the stale lists.""" - keys = set(_gamg_settings()) + keys = set(_gamg_settings()) | set(_gamg_settings(2)) for coarse in GEOMETRIC_MG_COARSE_SOLVERS: keys |= set(_geometric_mg_settings(coarse)) return keys @@ -355,12 +397,21 @@ def geometric_mg_bundle(coarse="redundant", smoother="robust"): return _bundle(_geometric_mg_settings(coarse, smoother)) -def gamg_bundle(): +def gamg_bundle(block_size=1): """Algebraic multigrid — the fallback when no geometric hierarchy exists. + Parameters + ---------- + block_size : int + Components per node in the field this block solves — ``mesh.dim`` for a + velocity or other vector field, ``1`` (the default) for a scalar. On a + vector field this is worth a factor of two to nine; see the commentary + in :func:`_gamg_settings`. Pass ``1`` rather than guessing: telling GAMG + the wrong node size on a scalar block would be a pure loss. + Returns ------- MGSettings Settings and stale keys; see :meth:`MGSettings.apply`. """ - return _bundle(_gamg_settings()) + return _bundle(_gamg_settings(block_size)) diff --git a/tests/test_0203_solver_wallclock_guard.py b/tests/test_0203_solver_wallclock_guard.py index 0aa7e076..dfe82197 100644 --- a/tests/test_0203_solver_wallclock_guard.py +++ b/tests/test_0203_solver_wallclock_guard.py @@ -114,9 +114,17 @@ def test_sub_solve_gauge_sees_work_the_outer_count_cannot(guarded): def test_first_solve_admits_its_count_is_a_lower_bound(): """A fresh solver cannot see the fieldsplit blocks until PETSc has set up the - preconditioner, part-way through its first solve. The count that results is short — - measured, by about half — so the report must say so rather than pass it off as - exact.""" + preconditioner, part-way through its first solve. Whatever the monitor catches is + therefore a LOWER BOUND, and the report must say so rather than pass it off as + exact. + + How short that bound is depends on where PCSetUp falls relative to the first + application of the velocity block, which is a function of the preconditioner + configuration and not something the contract promises. It used to be short by about + half. Since the GAMG and outer-Krylov defaults changed (#579, #576) it is tight on + this problem — same count either way — so this asserts the contract (a bound, and + honestly flagged) rather than the size of the gap, which was incidental. + """ solver, _ = _stokes("lb", lambda v: 1.0) solver.solve() first = solver.solve_report.sub["velocity"] @@ -125,9 +133,9 @@ def test_first_solve_admits_its_count_is_a_lower_bound(): assert not first.complete, "the first solve claimed an exact count it cannot have" assert second.complete - assert second.its > first.its, ( - "the first solve should undercount; if it does not, the 'lower bound' contract " - "is either wrong or no longer needed" + assert second.its >= first.its, ( + "the first solve reported MORE work than the complete count that follows it, so " + "it is not a lower bound at all" ) From 7c386400bdcae18ca3c7a06dc53dd230dc3e9911 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 17 Aug 2026 15:35:24 +1000 Subject: [PATCH 2/2] Declare mat_block_size only where the velocity block is node-blocked MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit MatSetFromOptions hands mat_block_size to PetscLayoutSetBlockSize, which is a hard error rather than a hint it may decline: Arguments are incompatible Local size 67 not compatible with block size 2 Whether it divides is a property of the particular COMBINATION of boundary conditions, not of their kind. Constrained DOFs are absent from the field's global section, so on a 3x3 P2 velocity field: no velocity BCs 98, full vector Dirichlet on every wall 50, component-wise free slip 70 — all even — while (0,0) on Bottom with (0,None) on the other three gives 67. Free slip is safe by an accident of counting, which is why the measurements behind #579 never met this and test_1010_stokesCart did. _withdraw_block_size_if_not_node_blocked drops the option when the block is not node-blocked, immediately before setFromOptions, which is the first point the field decomposition exists. _apply_preconditioner_options re-pushes the bundle on every build, so this runs on every build too. The decision is collective because divisibility is MIXED across ranks: on a 6x6 box the per-rank velocity sizes are [127, 150] at np=2 and [72, 94, 111] at np=3, so a rank-local decision would leave ranks disagreeing about the options DB. test_1022 pins both directions. The second case is the one that matters: a gate that withdrew the option unconditionally would fix the crash and silently undo #579, and the crash test alone cannot tell them apart. test_1021's bundle-ownership test built `owned` from a list that never included gamg_bundle(block_size=2) — the one bundle that sets the key — so every bundle's stale list carried a key `owned` lacked. The sibling is now in the list, which also brings it under the test. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 70 ++++++++++++ tests/test_1021_mg_option_bundle.py | 7 +- tests/test_1022_velocity_block_size_gate.py | 102 ++++++++++++++++++ 3 files changed, 178 insertions(+), 1 deletion(-) create mode 100644 tests/test_1022_velocity_block_size_gate.py diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 59427964..c74b2eb7 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -781,6 +781,67 @@ class SolverBaseClass(uw_object): # Force a full rebuild so the new option bundle is pushed to PETSc. self.is_setup = False + def _withdraw_block_size_if_not_node_blocked(self, field_name=None, prefix=""): + """Remove ``mat_block_size`` when the block it describes is not node-blocked. + + The GAMG bundle declares ``mat_block_size`` so that coarsening aggregates + nodes rather than scalars (#579). ``MatSetFromOptions`` hands that to + ``PetscLayoutSetBlockSize``, which is a HARD ERROR when a rank's local row + count is not divisible by it — PETSc does not treat it as a hint it may + decline: + + Arguments are incompatible + Local size 67 not compatible with block size 2 + + Whether it divides is a property of the particular COMBINATION of + boundary conditions, not of their kind. Constrained degrees of freedom + are absent from the field's global section, so on a 3x3 P2 velocity + field: no velocity BCs gives 98, full vector Dirichlet on every wall + gives 50, and component-wise free slip gives 70 — all even, all fine — + while ``(0, 0)`` on one wall with ``(0, None)`` on the other three gives + 67. There is then no node blocking to declare, and asking for one is not + a lost optimisation but an error. + + The decision is COLLECTIVE. PETSc requires divisibility on every rank, + and a rank-local decision would leave ranks disagreeing about the + contents of the options DB. + + Called immediately before ``setFromOptions`` because that is the first + point at which the field decomposition exists; + ``_apply_preconditioner_options`` re-pushes the bundle on every build, so + this has to run on every build too. + """ + + block_size = getattr(self, "_pc_block_size", 1) + if block_size <= 1: + return + + names, isets, _ = self.dm.createFieldDecomposition() + names = list(names) + if field_name is None: + if len(names) != 1: + return + field_name = names[0] + elif field_name not in names: + return + + local = isets[names.index(field_name)].getLocalSize() + + from mpi4py import MPI + + divides = self.dm.comm.tompi4py().allreduce( + local % block_size == 0, op=MPI.LAND + ) + if divides: + return + + key = f"{self.petsc_options_prefix}{prefix}mat_block_size" + options = PETSc.Options() + if key in options: + options.delValue(key) + + self._pc_block_size_withdrawn = (field_name, block_size, local) + def _apply_preconditioner_options(self): """Push the PETSc option bundle implied by ``self.preconditioner``. @@ -3839,6 +3900,8 @@ class SNES_Scalar(SolverBaseClass): self.dm.setUp() + self._withdraw_block_size_if_not_node_blocked() + self.snes = PETSc.SNES().create(PETSc.COMM_WORLD) self.snes.setDM(self.dm) self.snes.setOptionsPrefix(self.petsc_options_prefix) @@ -4916,6 +4979,8 @@ class SNES_Vector(SolverBaseClass): self.dm.setUp() + self._withdraw_block_size_if_not_node_blocked() + self.snes = PETSc.SNES().create(PETSc.COMM_WORLD) self.snes.setDM(self.dm) self.snes.setOptionsPrefix(self.petsc_options_prefix) @@ -5674,6 +5739,8 @@ class SNES_MultiComponent(SolverBaseClass): self.dm.setUp() + self._withdraw_block_size_if_not_node_blocked() + self.snes = PETSc.SNES().create(PETSc.COMM_WORLD) self.snes.setDM(self.dm) self.snes.setOptionsPrefix(self.petsc_options_prefix) @@ -8766,6 +8833,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): if self._block_constraint_bcs: self._setup_block_fieldsplit_options() + self._withdraw_block_size_if_not_node_blocked( + "velocity", prefix="fieldsplit_velocity_") + self.snes = PETSc.SNES().create(PETSc.COMM_WORLD) self.snes.setDM(self.dm) self.snes.setOptionsPrefix(self.petsc_options_prefix) diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py index 1b60f39e..cea5aee9 100644 --- a/tests/test_1021_mg_option_bundle.py +++ b/tests/test_1021_mg_option_bundle.py @@ -123,7 +123,12 @@ def test_bundles_clear_each_others_keys(): """Bundles share an options prefix, so switching between them must leave no key behind — every key a bundle does not set, a sibling does, and it must appear in that bundle's stale list.""" - bundles = [multigrid_options.gamg_bundle()] + # gamg_bundle(block_size=...) is a sibling like any other: it shares the + # prefix, and it is the one that sets `mat_block_size`. Left out, `owned` + # misses that key while every bundle's stale list carries it, and the + # assertion fails for bundles that are behaving correctly. + bundles = [multigrid_options.gamg_bundle(), + multigrid_options.gamg_bundle(block_size=2)] bundles += [multigrid_options.geometric_mg_bundle(coarse=c) for c in multigrid_options.GEOMETRIC_MG_COARSE_SOLVERS] owned = set().union(*(set(b.settings) for b in bundles)) diff --git a/tests/test_1022_velocity_block_size_gate.py b/tests/test_1022_velocity_block_size_gate.py new file mode 100644 index 00000000..a1c1b735 --- /dev/null +++ b/tests/test_1022_velocity_block_size_gate.py @@ -0,0 +1,102 @@ +"""`mat_block_size` is declared only when the velocity block really is node-blocked. + +The GAMG bundle declares ``mat_block_size`` so coarsening aggregates nodes rather +than scalars (#579). ``MatSetFromOptions`` hands that to +``PetscLayoutSetBlockSize``, which is a hard error — not a hint it may decline — +when a rank's local row count is not divisible: + + Arguments are incompatible + Local size 67 not compatible with block size 2 + +Whether it divides is a property of the particular COMBINATION of boundary +conditions, not of their kind. Constrained degrees of freedom are absent from the +field's global section, so on a 3x3 P2 velocity field: + + no velocity BCs 98 even + full vector Dirichlet, every wall 50 even + component-wise free slip 70 even + (0, 0) Bottom, (0, None) elsewhere 67 ODD + +which is why the free-slip measurements behind #579 never met it and +``test_1010_stokesCart`` did. + +The second test is the one that matters. A gate that withdrew the option +unconditionally would fix the crash and silently discard everything #579 was +for, and the first test alone cannot tell the two apart. +""" + +import pytest +import sympy + +import underworld3 as uw +from petsc4py import PETSc + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _stokes(tag, bcs): + mesh = uw.meshing.StructuredQuadBox(elementRes=(3, 3)) + v = uw.discretisation.MeshVariable(f"Vb{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"Pb{tag}", mesh, 1, degree=1, continuous=True) + + solver = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + solver.constitutive_model = uw.constitutive_models.ViscousFlowModel + solver.constitutive_model.Parameters.shear_viscosity_0 = 1 + solver.bodyforce = sympy.Matrix([0, -1]) + for value, boundary in bcs: + solver.add_dirichlet_bc(value, boundary) + + return solver + + +def _block_size_option(solver): + key = f"{solver.petsc_options_prefix}fieldsplit_velocity_mat_block_size" + options = PETSc.Options() + + return options.getInt(key) if key in options else None + + +def test_an_indivisible_velocity_block_solves(): + """The asymmetric mix that made PCSetUp_FieldSplit refuse the sub-matrix.""" + + solver = _stokes( + "odd", + [ + ((0.0, 0.0), "Bottom"), + ((0.0, None), "Top"), + ((0.0, None), "Left"), + ((0.0, None), "Right"), + ], + ) + + solver.solve() + + assert _block_size_option(solver) is None, ( + "the velocity block is not node-blocked here, so mat_block_size must " + "have been withdrawn" + ) + + +def test_a_divisible_velocity_block_keeps_the_block_size(): + """The control: where the declaration is valid it survives. + + Without this, withdrawing `mat_block_size` on every solve would pass the + test above and quietly undo #579. + """ + + solver = _stokes( + "even", + [ + ((0.0, 0.0), "Bottom"), + ((0.0, 0.0), "Top"), + ((0.0, 0.0), "Left"), + ((0.0, 0.0), "Right"), + ], + ) + + solver.solve() + + assert _block_size_option(solver) == solver.mesh.dim, ( + "a fully Dirichlet velocity block is node-blocked, so GAMG should still " + "be told its block size" + )