diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index f3df7cf6..c42bbb47 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 @@ -772,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``. @@ -915,9 +985,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 +3325,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 @@ -3830,6 +3895,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) @@ -4162,14 +4229,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 @@ -4906,6 +4972,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) @@ -5176,14 +5244,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 @@ -5663,6 +5730,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) @@ -6016,6 +6085,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) @@ -6058,13 +6141,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 = {} @@ -6935,6 +7017,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) @@ -8822,6 +8908,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/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" ) 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" + )