diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 3c0f92cd..7e442b5f 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -210,6 +210,13 @@ class SolverBaseClass(uw_object): # helper below is a no-op. self._preconditioner = "auto" self._pc_option_prefix = None + # An explicit `preconditioner="fmg"` on a single-field solver cannot + # take the native route (#276), so it is honoured via custom-P + # transfers over the mesh's own dm_hierarchy instead (#478). This flag + # carries the request from _apply_preconditioner_options (build time) + # to custom_mg.build_transfers (first solve); it is re-derived on + # every resolution, same staleness rule as _pc_resolved. + self._pc_single_field_geo_requested = False # The pc_type value this helper last managed. Subclasses that opt in set # their __init__ default ("gamg"); used in "auto" mode to tell an # untouched framework default (eligible for FMG upgrade) apart from an @@ -238,6 +245,22 @@ class SolverBaseClass(uw_object): # resolved would be exactly the stale-but-authoritative-looking summary this # reporting exists to prevent. self._pc_resolved = False + # Readable record of every preconditioner fallback / degrade / guard-skip + # decision taken for this solver, keyed by site name — see the public + # `pc_fallbacks` property. Written ONLY through _record_pc_fallback (the + # record is WRITTEN, never inferred — same doctrine as + # _push_managed_option). Reset rule: cleared each time + # _apply_preconditioner_options re-resolves (the same staleness rule as + # _pc_resolved); solve-time sites (custom_mg, rotated_bc) record after + # that, so the record always describes the CURRENT resolution. + # Reason vocabulary (fixed; tests assert on it): + # "unavailable" — the requested configuration could not be built here + # "declined" — available in principle, but a policy chose otherwise + # "build_failed" — an attempted build raised and a fallback was used + # "check_skipped" — a correctness guard did not run (its failure mode + # is sanctioned, but the skip is now on the record) + # "forced" — a required key overrode a user/unset value + self._pc_fallbacks = {} # Custom multigrid prolongation hierarchy (see set_custom_mg / # utilities.custom_mg). None => standard FMG/GAMG path, unchanged. @@ -270,6 +293,48 @@ class SolverBaseClass(uw_object): out[key] = self.petsc_options.getString(name) return out + def _record_pc_fallback(self, site, *, requested, installed, reason, detail=""): + """Record one preconditioner fallback / degrade / guard-skip decision. + + The single write path into ``pc_fallbacks`` (records are WRITTEN, never + inferred). ``reason`` must come from the fixed vocabulary documented at + ``_pc_fallbacks`` in ``__init__``. Runs on every rank — the record is + state, not output, so it must not be rank-gated the way warnings are. + """ + self._pc_fallbacks[site] = dict(requested=requested, installed=installed, + reason=reason, detail=detail) + + @property + def pc_fallbacks(self): + """Every preconditioner fallback the current resolution took, by site. + + A dict keyed by site name (e.g. ``"single_field_gate"``, + ``"no_hierarchy"``, ``"custom_mg.build"``); each value is a dict with + ``requested`` (what was asked for), ``installed`` (what actually runs), + ``reason`` (one of ``"unavailable"``, ``"declined"``, ``"build_failed"``, + ``"check_skipped"``, ``"forced"``) and ``detail``. Empty means the + resolved preconditioner is exactly what was requested and every guard + ran — the clean-solve state a test can assert on. + + The record is reset whenever the preconditioner options re-resolve (a + rebuild/remesh), and solve-time sites (``utilities.custom_mg``, + ``utilities.rotated_bc``) re-record each solve, so it always describes + the current resolution. Solvers that manage their own PC options + (``_pc_option_prefix is None``) never clear at rebuild; their sites are + custom_mg-only, which re-record per solve. + + Returns + ------- + dict + A copy — mutating it does not affect the solver. + + See Also + -------- + preconditioner_settings : the option values the managed block resolved to. + strategy : a readable summary of the same resolution. + """ + return {site: dict(rec) for site, rec in self._pc_fallbacks.items()} + @property def _user_overridden_pc_options(self): """The managed-block keys the USER set, as (key, value) pairs. @@ -660,9 +725,19 @@ class SolverBaseClass(uw_object): - ``"auto"`` (default) — use geometric Full Multigrid (FMG) when the mesh carries a genuine refinement hierarchy (``len(mesh.dm_hierarchy) > 1``, i.e. built with ``refinement >= 1``), - otherwise fall back to algebraic multigrid (GAMG). - - ``"fmg"`` (alias ``"mg"``) — force geometric multigrid. Requires a + otherwise fall back to algebraic multigrid (GAMG). On a single-field + (scalar/vector) solver ``"auto"`` keeps GAMG even with a hierarchy — + the decline is recorded in :attr:`pc_fallbacks`. + - ``"fmg"`` (alias ``"mg"``) — geometric multigrid. Requires a refinement hierarchy; warns and falls back to GAMG if none exists. + On a single-field solver the native FMG path is unreliable (#276), + so the request is honoured via custom-P transfers built over the + same hierarchy (``utilities.custom_mg``) — installed on the live PC + at the first solve. If that build fails (e.g. a deformed mesh whose + coarse levels kept reference coordinates), the solve degrades to + GAMG with a readable :attr:`pc_fallbacks` record. This is a + *preference*; ``custom_mg.set_custom_fmg`` is the *demand* form and + raises on failure instead. - ``"gamg"`` — force algebraic multigrid (the historical default). Geometric multigrid is inherently robust to mesh anisotropy (it is built @@ -687,6 +762,13 @@ class SolverBaseClass(uw_object): f"preconditioner must be 'auto', 'fmg', or 'gamg' (got {value!r})" ) self._preconditioner = choice + # A hierarchy cached by an earlier RESOLUTION (the auto/"fmg" install) + # must not outlive a new explicit choice: auto_inject_custom_mg + # re-installs solver._custom_mg unconditionally, which would leave a + # later preconditioner="gamg" unreachable with a clean pc_fallbacks + # record. A user registration (set_custom_fmg) is a demand and is kept. + if isinstance(self._custom_mg, dict) and self._custom_mg.get("auto_cached"): + self._custom_mg = None # Force a full rebuild so the new option bundle is pushed to PETSc. self.is_setup = False @@ -704,6 +786,12 @@ class SolverBaseClass(uw_object): # Every path from here is a resolution decision, including "the user owns # these options, leave them alone". self._pc_resolved = True + # Fresh resolution => fresh fallback record (same staleness rule as + # _pc_resolved). Solve-time sites re-record after this. The custom-P + # reroute request is re-derived below for the same reason (a remesh + # can collapse the hierarchy it depends on). + self._pc_fallbacks.clear() + self._pc_single_field_geo_requested = False opts = self.petsc_options @@ -750,22 +838,44 @@ class SolverBaseClass(uw_object): # on the common curved-shell cases (issue #276) as well as some flat # high-degree ones. The Stokes velocity sub-block (prefix # "fieldsplit_velocity_") is the validated, robust native-FMG path and is - # unaffected. So never auto-route a single-field solver to native FMG — - # fall back to GAMG. Geometric MG on a scalar/vector solver is available, - # robustly, via ``underworld3.utilities.custom_mg.set_custom_fmg`` (own - # barycentric/RBF prolongation + Galerkin coarse operators; no injection). + # unaffected. So a single-field solver never routes to NATIVE FMG. Two + # routes instead (#478): + # * explicit `preconditioner="fmg"` is honoured via custom-P transfers + # over the mesh's own dm_hierarchy (no DMCreateInjection anywhere, so + # the err62 failure mode cannot arise). The options DB deliberately + # keeps GAMG as the safe base configuration — the custom-P PCMG is + # installed on the LIVE PC at the first solve (auto_inject_custom_mg + # -> custom_mg.build_transfers, requested-native source), the same + # shape the adapt-child pickup uses. If the transfer build fails, the + # solve degrades to that GAMG base, recorded in `pc_fallbacks`. + # * "auto" keeps GAMG (a default change needs its own validation + # campaign, per #478) — the decline is recorded, never warned. + # `set_custom_fmg` remains the DEMAND form (raises on failure); + # `preconditioner="fmg"` is a PREFERENCE (degrades, loudly and readably). if want_fmg and prefix == "": - if self._preconditioner == "fmg" and uw.mpi.rank == 0: - import warnings - warnings.warn( - f"[{self.name}] preconditioner='fmg' is not supported on a " - f"single-field (scalar/vector) solver: native geometric FMG " - f"needs DMCreateInjection, which PETSc cannot reliably build " - f"on a refined DMPlex (issue #276). Falling back to GAMG. For " - f"geometric MG on this solver use " - f"underworld3.utilities.custom_mg.set_custom_fmg().", - stacklevel=2, - ) + if self._preconditioner == "fmg": + self._pc_single_field_geo_requested = True + self._record_pc_fallback( + "single_field_gate", + requested="native geometric FMG (preconditioner='fmg')", + installed="custom-P geometric MG (resolved at first solve)", + reason="declined", + detail="native single-field FMG needs DMCreateInjection, " + "which PETSc cannot reliably build on a refined " + "DMPlex (#276); the request is honoured via custom-P " + "transfers over the mesh hierarchy instead, degrading " + "to GAMG (recorded) if the transfer build fails") + else: + # "auto" declines silently by design (never a new warning) — + # but the decline is the direction that matters (#484), so it + # is on the record. + self._record_pc_fallback( + "single_field_gate", + requested="geometric FMG (preconditioner='auto', hierarchy present)", + installed="gamg", + reason="declined", + detail="single-field native FMG is fragile (#276); " + "auto never routes there") want_fmg = False # The option VALUES live in utilities.multigrid_options, which is the @@ -789,14 +899,22 @@ class SolverBaseClass(uw_object): owned=self._managed_pc_options) self._pc_managed_value = "mg" else: - if self._preconditioner == "fmg" and n_levels <= 1 and uw.mpi.rank == 0: - import warnings - warnings.warn( - f"[{self.name}] preconditioner='fmg' requested but the mesh " - f"has no refinement hierarchy; falling back to GAMG. Build the " - f"mesh with refinement >= 1 to enable geometric multigrid.", - stacklevel=2, - ) + if self._preconditioner == "fmg" and n_levels <= 1: + self._record_pc_fallback( + "no_hierarchy", + requested="geometric FMG (preconditioner='fmg')", + installed="gamg", + reason="unavailable", + detail="the mesh has no refinement hierarchy; build it with " + "refinement >= 1 to enable geometric multigrid") + if uw.mpi.rank == 0: + import warnings + warnings.warn( + f"[{self.name}] preconditioner='fmg' requested but the mesh " + f"has no refinement hierarchy; falling back to GAMG. Build the " + 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) @@ -829,6 +947,13 @@ class SolverBaseClass(uw_object): return gkey = f"{prefix}pc_mg_galerkin" if (not opts.hasName(gkey)) or opts.getString(gkey) == "none": + self._record_pc_fallback( + "galerkin_forced", + requested=f"{gkey} unset (or 'none')", + installed="both", + reason="forced", + detail="UW3 installs no coarse-DM operator callbacks, so " + "geometric MG requires Galerkin RAP coarse operators") if uw.mpi.rank == 0: import warnings warnings.warn( @@ -5844,6 +5969,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self._tolerance = 1.0e-4 self._strategy = "default" + # Owned-option latch state (see _resolve_owned_option): the value THIS + # solver last pushed per key, and any user value latched per key. + # Ownership is RECORDED, never inferred — the same doctrine as + # _managed_pc_options. + self._owned_option_pushes = {} + self._owned_option_user = {} + # Participate in the auto FMG/GAMG switch on the velocity fieldsplit # block (see the `preconditioner` property). The velocity pc/mg keys # set below are the GAMG default; _apply_preconditioner_options() @@ -6456,14 +6588,27 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): """ Solver convergence tolerance for the Stokes saddle-point system. - Setting this value automatically configures PETSc tolerances for the - coupled velocity-pressure solve using Schur complement fieldsplit: - - ``snes_rtol``: Set to ``tolerance`` - - ``ksp_atol``: Set to ``tolerance * 1e-6`` - - ``fieldsplit_pressure_ksp_rtol``: Set to ``tolerance * 0.1`` - - ``fieldsplit_velocity_ksp_rtol``: Set to ``tolerance * 0.033`` + Setting it configures the PETSc tolerances of the coupled + velocity-pressure Schur-fieldsplit solve. The keys fall into two + ownership classes (#483): + + **OWNED** — re-asserted before every solve, *unless you set the key + explicitly, after which your value is honoured* (the same latch that + makes ``snes_max_it`` reachable): + + - ``snes_rtol`` = ``tolerance`` + - ``ksp_atol`` = ``tolerance * 1e-6`` + + **DERIVED at set time** — written once when you assign ``tolerance`` + (the class table ``_TOLERANCE_DERIVED_KEYS``), then yours to override: - Also enables Eisenstat-Walker adaptive tolerance (``snes_ksp_ew``). + - ``fieldsplit_pressure_ksp_rtol`` = ``tolerance * 0.1`` + - ``fieldsplit_velocity_ksp_rtol`` = ``tolerance * 0.033`` + + Also enables Eisenstat-Walker adaptive tolerance (``snes_ksp_ew``), + which re-picks the outer ``ksp_rtol`` every Newton step — so to steer + the linear solve via ``ksp_rtol`` you must first switch + ``snes_ksp_ew`` off. Returns ------- @@ -6483,8 +6628,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): #: BELOW the tolerance demanded of the outer solve. These factors are that margin. #: Their existence is principled; their size is inherited convention, so they are #: DEFAULTS a user may override rather than values this property owns outright. - _INNER_RTOL_MARGIN = {"fieldsplit_pressure_ksp_rtol": 0.1, - "fieldsplit_velocity_ksp_rtol": 0.033} + #: Subclasses declare their own table (Stokes_Constrained derives the outer + #: ksp_rtol and the Eisenstat-Walker pins instead — a real design difference: + #: EW pinning owns its outer accuracy). `_INNER_RTOL_MARGIN` is the + #: historical name for this class's table, kept as an alias. + _TOLERANCE_DERIVED_KEYS = {"fieldsplit_pressure_ksp_rtol": 0.1, + "fieldsplit_velocity_ksp_rtol": 0.033} + _INNER_RTOL_MARGIN = _TOLERANCE_DERIVED_KEYS @tolerance.setter def tolerance(self, value): @@ -6495,57 +6645,82 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options["ksp_atol"] = self._tolerance * 1.0e-6 - # Setting the tolerance re-derives the inner margins from it — that is this + # Setting the tolerance re-derives the margins from it — that is this # property's job, and a user who changes the tolerance expects it. What must NOT # happen is `solve()` re-deriving them on every call: it did (pyx `solve()` # round-trips `self.tolerance` immediately before `setFromOptions()`), which # overwrote any user value between it being set and PETSc reading it and made # both documented options silently unreachable (#477). `solve()` now re-asserts # only the outer keys, via `_reassert_outer_tolerances`. - for key, margin in self._INNER_RTOL_MARGIN.items(): + self._derive_tolerance_margins() + + def _derive_tolerance_margins(self): + """Write this class's DERIVED tolerance keys from the current tolerance. + + One mechanism, two tables: each saddle-point class declares + ``_TOLERANCE_DERIVED_KEYS`` ({option key: margin factor}) and this + method applies it. Derivation happens at SET time only — a user who + overrides a derived key afterwards keeps it (#477/#483); `solve()` + never re-derives these. + """ + for key, margin in self._TOLERANCE_DERIVED_KEYS.items(): self.petsc_options[key] = self._tolerance * margin - def _resolve_snes_max_it(self, default): - """The nonlinear iteration cap for this solve, honouring a user-set - ``snes_max_it``. + def _resolve_owned_option(self, key, default): + """The value this solve should push for an option the solver OWNS, + honouring a user-set value. - ``solve()`` pushes this option before every solve, so it has to be able to tell - its OWN previous push from a value the user set — otherwise the first solve makes - the option permanently unreachable. Call once, before this solve pushes anything. + ``solve()`` re-pushes the owned keys before every solve, so the + resolver has to tell its OWN previous push from a value the user set — + otherwise the first solve makes the option permanently unreachable + (the #477 failure shape; ruling D18, generalised for #483). Call once + per key, before this solve pushes anything. """ - pushed = getattr(self, "_snes_max_it_pushed", None) - user = getattr(self, "_snes_max_it_user", None) - if self.petsc_options.hasName("snes_max_it"): + pushed = self._owned_option_pushes.get(key) + user = self._owned_option_user.get(key) + if self.petsc_options.hasName(key): try: - current = int(self.petsc_options.getInt("snes_max_it")) + current = type(default)(self.petsc_options.getString(key)) except Exception: return default if pushed is None or current != pushed: # Never pushed by us, or the user has moved it since. Latch: from here on # the option is theirs, because the next solve will read back OUR push of # THEIR value and would otherwise mistake it for our own default. - self._snes_max_it_user = current + self._owned_option_user[key] = current return current if user is not None: return user return default + def _push_owned_option(self, key, value): + """Push an owned option and remember what was pushed, so a later solve + can tell this solver's own value apart from a user override.""" + self.petsc_options.setValue(key, value) + self._owned_option_pushes[key] = value + + def _resolve_snes_max_it(self, default): + """The nonlinear iteration cap for this solve, honouring a user-set + ``snes_max_it`` (see ``_resolve_owned_option`` for the mechanism).""" + return self._resolve_owned_option("snes_max_it", int(default)) + def _push_snes_max_it(self, value): - """Push ``snes_max_it`` and remember what was pushed, so a later solve can tell - this solver's own value apart from a user override.""" - value = int(value) - self.petsc_options.setValue("snes_max_it", value) - self._snes_max_it_pushed = value + self._push_owned_option("snes_max_it", int(value)) def _reassert_outer_tolerances(self): """Re-push the OUTER tolerance keys before a solve, leaving the inner margins be. `solve()` may have changed `snes_max_it` and the SNES type for a Picard warm-up, - so the outer settings are re-asserted before the real solve. The sub-block rtols - are deliberately excluded: they belong to whoever set them last, which may be the + so the outer settings are re-asserted before the real solve. The keys are OWNED + (re-pushed each solve) but ownership is polite: a user who explicitly sets + `snes_rtol` or `ksp_atol` is honoured from then on — before #483 both were + silently discarded here every solve, the worst of the reachability middle + grounds (documented as settable, actually owned). The sub-block rtols are + deliberately excluded: they belong to whoever set them last, which may be the user (#477). Overwriting them here is what made them unsettable.""" - self.petsc_options["snes_rtol"] = self._tolerance - self.petsc_options["ksp_atol"] = self._tolerance * 1.0e-6 + for key, derived in (("snes_rtol", float(self._tolerance)), + ("ksp_atol", float(self._tolerance) * 1.0e-6)): + self._push_owned_option(key, self._resolve_owned_option(key, derived)) # The Eisenstat-Walker flags are NOT re-asserted here. solve() never changes # them, so they stay in the options DB from the `tolerance` setter and # setFromOptions picks them up regardless — while re-asserting would make @@ -6851,10 +7026,20 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): the block solver's default ``newtonls`` defect-corrects a linear system in many steps when the Schur approximation is stiff. - Default ``False`` (opt-in). On the fieldsplit/iterative path it is - bit-identical on uniform ``mu`` and cracks the moderate-contrast wall; - but a monolithic ``lu`` solve factorizes the Pmat (``pc_use_amat`` is a - no-op there), so this term is not inert for direct solves — hence opt-in. + Default ``False`` (opt-in). **Reachability** (#486, from the PETSc + fieldsplit source): the flag swaps only the *Pmat* (h,h) block, so it + is live in exactly two regimes — + + - ``pc_fieldsplit_schur_precondition = "a11"``: the Schur + preconditioner is the grouped ``[p,h]`` Pmat block, which carries + the swap; + - a monolithic direct factorisation (``pc_type = lu``/``cholesky``) + of the Pmat. + + Under ``Stokes_Constrained``'s own defaults (``selfp`` + + ``diag_use_amat``) the Pmat (h,h) block is never read by the Schur + preconditioner and the flag is INERT — setting it there records a + ``multiplier_schur_pc`` entry in :attr:`pc_fallbacks` and warns. """ return self._multiplier_schur_pc @@ -8225,6 +8410,55 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # Lagrange-multiplier rows (block-constrained Stokes). Guarded: no-op # for ordinary Stokes. Register the interior screening residual and the # diagonal mass Jacobian/preconditioner on each multiplier's field. + # + # multiplier_schur_pc reachability (#486, resolved by tracing PETSc + # fieldsplit.c): the flag swaps only the Pmat (h,h) block below. With + # schur_precondition=selfp the Schur preconditioner Sp is assembled by + # MatSchurComplementGetPmat from sub-matrices split out of the AMAT + # whenever pc_fieldsplit_diag_use_amat is set (jac->mat[1], not + # jac->pmat[1]) — so under this class's defaults (selfp + + # diag_use_amat) the swapped block is provably never read. It IS read + # under schur_precondition=a11 (Sp = the grouped [p,h] Pmat block) and + # under a monolithic direct factorisation of the Pmat. An explicit + # opt-in silently doing nothing is exactly the #477 class -> record + # AND warn (unlike the auto declines, this one earns the warning). + if self._multipliers and self._multiplier_schur_pc: + _opts = self.petsc_options + _schur_pre = (_opts.getString("pc_fieldsplit_schur_precondition") + if _opts.hasName("pc_fieldsplit_schur_precondition") + else "") + _pc_type = (_opts.getString("pc_type") + if _opts.hasName("pc_type") else "") + # Read the VALUE, not the key's presence: diag_use_amat set to + # "false" means the Pmat block IS read and the opt-in is live + # (measured: Schur pre differs by rel-Frobenius 0.30 flag-on/off). + _diag_amat = _opts.getBool("pc_fieldsplit_diag_use_amat", False) + if (_schur_pre != "a11" and _diag_amat + and _pc_type not in ("lu", "cholesky")): + self._record_pc_fallback( + "multiplier_schur_pc", + requested="1/mu multiplier Schur mass (Pmat h,h block)", + installed="unread — selfp builds Sp from the Amat A11 block", + reason="declined", + detail=f"pc_fieldsplit_schur_precondition=" + f"'{_schur_pre or 'selfp'}' with diag_use_amat: the " + f"Pmat (h,h) block never reaches the Schur " + f"preconditioner; set schur_precondition='a11' (or " + f"factorise the Pmat directly) to make the opt-in " + f"live") + if uw.mpi.rank == 0: + import warnings + warnings.warn( + f"[{self.name}] multiplier_schur_pc=True has no effect " + f"under pc_fieldsplit_schur_precondition=" + f"'{_schur_pre or 'selfp'}' with diag_use_amat: the " + f"1/mu multiplier Schur mass is written into the Pmat " + f"(h,h) block, which selfp never reads (Sp is built " + f"from the Amat). Use " + f"petsc_options['pc_fieldsplit_schur_precondition'] = " + f"'a11' to make it live. See solver.pc_fallbacks.", + stacklevel=2, + ) for k, mvar in enumerate(self._multipliers): fid = mvar._solver_field_id # Operator (Amat) (lambda,lambda) block is ALWAYS the true screening eps so diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index d03c27f2..2f129b8e 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -2401,21 +2401,40 @@ def __init__( # preconditioner partition-dependent, which only cancels once the TRUE # residual is driven down). Pinning EW's initial = max rtol to the solver # tolerance makes the outer fgmres iterate until genuinely converged, so - # the velocity is partition-independent to round-off. Kept in sync by the - # `tolerance` setter below. - self.petsc_options["snes_ksp_ew_rtol0"] = self._tolerance * 1.0e-1 - self.petsc_options["snes_ksp_ew_rtolmax"] = self._tolerance * 1.0e-1 + # the velocity is partition-independent to round-off. Applied through the + # class derived-key table (kept in sync by the `tolerance` setter below); + # the table also writes ksp_rtol = tolerance * 0.1, which at the class + # default tolerance equals PETSc's own ksp_rtol default (and EW re-picks + # ksp_rtol per step regardless). + self._derive_tolerance_margins() return + #: Constrained derives DIFFERENT keys from the tolerance than the base + #: Stokes table: the outer ``ksp_rtol`` and the Eisenstat-Walker pins, + #: NOT the inner fieldsplit margins. That is a real design difference — + #: EW pinning owns this class's outer accuracy (see ``__init__``) — made + #: explicit here rather than hand-rolled in a second code path (#483). + _TOLERANCE_DERIVED_KEYS = {"ksp_rtol": 0.1, + "snes_ksp_ew_rtol0": 0.1, + "snes_ksp_ew_rtolmax": 0.1} + @property def tolerance(self): """Solver tolerance (see :class:`SNES_Stokes_SaddlePt.tolerance`). - Overridden so that, in addition to ``snes_rtol`` / ``ksp_rtol`` / - ``ksp_atol``, the Eisenstat-Walker initial and max relative tolerances are - pinned to ``tolerance * 0.1`` — otherwise EW's default (0.3) under-solves - the ill-conditioned augmented constrained system on a linear solve and the - velocity becomes partition-dependent (see ``__init__``). + Same two ownership classes as the base property (#483): + + **OWNED** (re-asserted each solve unless you set the key explicitly, + after which your value is honoured): ``snes_rtol`` = ``tolerance``, + ``ksp_atol`` = ``tolerance * 1e-6``. + + **DERIVED at set time** (the class table ``_TOLERANCE_DERIVED_KEYS``; + yours to override afterwards): ``ksp_rtol``, ``snes_ksp_ew_rtol0`` + and ``snes_ksp_ew_rtolmax``, all ``tolerance * 0.1`` — the EW pins + replace the base class's inner fieldsplit margins because EW's + default (0.3) under-solves the ill-conditioned augmented constrained + system on a linear solve and the velocity becomes + partition-dependent (see ``__init__``). """ return self._tolerance @@ -2423,10 +2442,8 @@ def tolerance(self): def tolerance(self, value): self._tolerance = value self.petsc_options["snes_rtol"] = value - self.petsc_options["ksp_rtol"] = value * 1.0e-1 self.petsc_options["ksp_atol"] = value * 1.0e-6 - self.petsc_options["snes_ksp_ew_rtol0"] = value * 1.0e-1 - self.petsc_options["snes_ksp_ew_rtolmax"] = value * 1.0e-1 + self._derive_tolerance_margins() def solve(self, *args, **kwargs): """Solve the constrained Stokes system (see :meth:`SNES_Stokes.solve`). @@ -2517,20 +2534,31 @@ def saddle_preconditioner(self): automatically. The grouped :math:`[p,\\lambda]` Schur preconditioner is formed by - ``selfp`` from the operator blocks, and the pressure mass it needs is the - ``1/viscosity`` (``1/constitutive_model.K``) term supplied automatically. - There is nothing for the user to set; this property is inert and assigning - to it raises. (The base :class:`SNES_Stokes` keeps a settable + ``selfp`` **from the operator (Amat) blocks alone**: with + ``diag_use_amat`` set (this class's default), selfp assembles + :math:`S_p \\approx A_{11} - A_{10}\\,\\mathrm{diag}(A_{00})^{-1}A_{01}` + from Amat sub-blocks and never reads the Pmat — so the automatic + ``1/viscosity`` pressure mass participates only if you override + ``pc_fieldsplit_schur_precondition = "a11"`` (an earlier version of + this docstring claimed selfp used it; that was drifted, see #486). + There is nothing for the user to set; this property is inert and + assigning to it raises. (The base :class:`SNES_Stokes` keeps a settable ``saddle_preconditioner`` as an advanced override.) """ + # TODO(BUG): should Constrained's selfp see the Pmat blocks at all? + # Under selfp + diag_use_amat the 1/mu pressure mass (_pp_G0) and the + # multiplier Schur mass (multiplier_schur_pc) are both provably unread + # by the Schur preconditioner (PETSc fieldsplit.c trace, #486). Whether + # Sp should instead be built with the Pmat A11 block is a solver-design + # question — out of scope for the #486 instrumentation, not changed here. return None @saddle_preconditioner.setter def saddle_preconditioner(self, value): raise AttributeError( "Stokes_Constrained does not use `saddle_preconditioner`: the Schur " - "preconditioner is built automatically (selfp + the 1/viscosity mass " - "from constitutive_model.K). Remove this assignment." + "preconditioner is built automatically (selfp, assembled from the " + "operator's Amat blocks). Remove this assignment." ) def _viscosity_scale(self): diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index f9d02b81..a933c010 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -705,6 +705,13 @@ def _install_velocity_block_transfers(solver, Ps, verbose=False): snes.computeJacobian(x0, J, Pmat) except PETSc.Error: # fallback: throwaway max_it=0 solve assembles + splits the operator + solver._record_pc_fallback( + "custom_mg.velocity_block_assembly", + requested="direct Jacobian assembly (computeFunction/computeJacobian)", + installed="throwaway max_it=0 assembly route; same PC installed", + reason="build_failed", + detail="snes.computeJacobian raised; the operator is assembled by " + "a zero-iteration solve instead") saved = (solver.petsc_options.getString("snes_max_it") if solver.petsc_options.hasName("snes_max_it") else None) solver.petsc_options["snes_max_it"] = 0 @@ -855,8 +862,14 @@ def build(self, solver): except Exception: # Sanctioned swallow: setUp can fail on a not-yet-fully-configured # SNES (pre-solve injection). The install paths call setUp again; - # the finest map then reads the DM's current global section. - pass + # the finest map then reads the DM's current global section. The + # skip is recorded so "the section was finalized" is checkable. + solver._record_pc_fallback( + "custom_mg.presolve_setup", + requested="pre-build snes.setUp() (finalize the DM section)", + installed="deferred to install-time setUp", + reason="check_skipped", + detail="setUp raised on the not-yet-fully-configured SNES") coords, maps, ncomp = [], [], [] for k, mesh in enumerate(self.level_meshes): @@ -935,7 +948,15 @@ def _assert_finest_matches_operator(solver, finest_map, parallel): try: op_n = int(solver.snes.getJacobian()[0].getSize()[0]) except Exception: - return # can't read operator -> skip + # Sanctioned: no readable operator to check against — record the + # skipped guard rather than silently waiving it (#484). + solver._record_pc_fallback( + "custom_mg.finest_operator_check", + requested="finest reduced-map vs operator span check", + installed="unchecked", + reason="check_skipped", + detail="could not read the assembled operator") + return if op_n <= 0: return if parallel: @@ -958,6 +979,53 @@ def install(self, solver, verbose=False): _install_transfers(solver, self.transfers, verbose=verbose) +class _DMLevelView: + """A mesh-shaped view over one DM of a native ``refine()`` hierarchy. + + ``CustomMGHierarchy`` consumes *meshes*, but the requested-native source in + :func:`build_transfers` (#478) has only the raw ``mesh.dm_hierarchy`` DMs. + This adapter provides exactly what a coarse level is asked for: ``.dm`` + (used by ``_clone_dm_with_solver_discretisation`` — the hierarchy DMs carry + the boundary labels ``refine()`` propagates, which is all copyDS needs) + plus ``_get_coords_for_basis``, delegated UNBOUND to + ``discretisation.Mesh`` so there is one implementation. That method reads + only ``self.dm`` and the four scalars copied here (``dim``, ``cdim``, + ``isSimplex``, ``qdegree`` — shared by every level of a uniform + refinement), which is what makes the unbound call safe. No full ``Mesh`` + is built: a coarse MG level needs no variables, caches, or registries. + """ + + def __init__(self, dm, fine_mesh): + # A CLONE, never the hierarchy DM itself: the base (gmsh-imported) + # level may need its coordinate field repaired below, and the shared + # hierarchy DMs also feed the native Stokes FMG route. + self.dm = dm.clone() + self.dim = fine_mesh.dim + self.cdim = fine_mesh.cdim + self.isSimplex = fine_mesh.isSimplex + self.qdegree = fine_mesh.qdegree + + # The gmsh-imported BASE level carries section-only coordinates (its + # coordinate field is a PetscContainer, no PetscFE), and + # DMCreateInterpolation from such a source silently returns a ZERO + # matrix (measured) — _get_coords_for_basis would then hand every + # level-0 node the coordinate (0,0) and the transfer build collapses. + # refine() gives the child levels an FE coordinate space; give the + # clone's base the same: P1 Lagrange on the identical vertex layout. + cdm = self.dm.getCoordinateDM() + field = cdm.getField(0) + fobj = field[0] if isinstance(field, tuple) else field + if not isinstance(fobj, PETSc.FE): + fe = PETSc.FE().createLagrange(self.dim, self.cdim, self.isSimplex, + 1, self.qdegree, comm=PETSc.COMM_SELF) + cdm.setField(0, fe) + cdm.createDS() + + def _get_coords_for_basis(self, degree, continuous): + from underworld3.discretisation import Mesh + return Mesh._get_coords_for_basis(self, degree, continuous) + + # --------------------------------------------------------------------------- # # Entry points # --------------------------------------------------------------------------- # @@ -988,8 +1056,9 @@ def set_custom_fmg(solver, coarse_meshes, *, builder="barycentric", def build_transfers(solver, field_id=None): """The custom-P prolongations this solver should drive, built and ready to - install — from either a solver-set hierarchy (``set_custom_fmg``) or a - **mesh-owned** one (a ``mesh.adapt`` refinement child). + install — from a solver-set hierarchy (``set_custom_fmg``), a **mesh-owned** + one (a ``mesh.adapt`` refinement child), or a **requested-native** one + (explicit ``preconditioner="fmg"`` on a single-field solver, #478). This is the shared resolution rule for every route that can drive custom-P multigrid: the standard solve path via :func:`auto_inject_custom_mg`, and the @@ -997,12 +1066,20 @@ def build_transfers(solver, field_id=None): ``utilities.rotated_bc`` and so never reaches the standard injection hook (#467). Both must answer "which hierarchy does this solver get?" the same way. + Resolution order: **solver-set > mesh-owned > requested-native.** The first + is a DEMAND (build errors raise — the user registered it explicitly); the + other two are PREFERENCES, built through the same opportunistic arm + (barycentric, RBF retry, degrade to the solver's default preconditioner — + every step recorded in ``solver.pc_fallbacks``). + A refinement child carries ``mesh._custom_mg_coarse_meshes`` (the static coarse tail), so a :class:`CustomMGHierarchy` ``[*coarse, solver.mesh]`` targeting ``field_id`` (0 for the Stokes velocity block, None for scalar/vector) is built lazily on first solve — every solver on an adapted - mesh drives geometric MG with no per-solver call. A solver-set hierarchy (if - present) always wins. + mesh drives geometric MG with no per-solver call. The requested-native + source instead wraps the mesh's own ``dm_hierarchy`` tail in + :class:`_DMLevelView` adapters — the same coarse levels native FMG would + use, driven through injection-free custom-P transfers. Parameters ---------- @@ -1037,28 +1114,46 @@ def build_transfers(solver, field_id=None): # so it is faithful to the operator on adapt children too — including scalar # semi-Lagrangian advection-diffusion (which earlier had to be skipped). coarse = getattr(solver.mesh, "_custom_mg_coarse_meshes", None) - if coarse is None: + if coarse is not None: + # An EXPLICIT preconditioner choice beats the opportunistic pickup. Before + # this guard, `solver.preconditioner = "gamg"` on an adapt child was + # silently clobbered back to the custom-P PCMG at solve time (measured: + # both arms of test_0842's fmg-vs-gamg comparison ran pc_type=mg), so a + # user could not opt out and any FMG-vs-GAMG comparison was vacuous. + # `_pc_user_override` is the same statement in the other spelling: the + # solver's option manager has latched "the user owns this block's pc_type" + # (they wrote a pc_type of their own into petsc_options), and an + # opportunistic pickup must stand down for exactly the same reason. + # "auto" (the default) still picks up the mesh-owned hierarchy. + # NOTE the arity: this function returns a 2-tuple, never bare None — a bare + # `return` here is what turned the gate into a TypeError at the call site + # when this hunk migrated from auto_inject_custom_mg (which returns nothing) + # during the #488 x #471 merge. + if (getattr(solver, "_preconditioner", "auto") == "gamg" + or getattr(solver, "_pc_user_override", False)): + return None, None + level_tail = list(coarse) + builder = getattr(solver.mesh, "_custom_mg_builder", "barycentric") + elif getattr(solver, "_pc_single_field_geo_requested", False): + # Requested-native source (#478): an explicit `preconditioner="fmg"` + # on a single-field solver. The gate in _apply_preconditioner_options + # set the flag only when the mesh reported a hierarchy, but re-check + # here — a remesh between build and solve can collapse it, and this + # arm must degrade readably, never crash a solve. + hierarchy_dms = list(getattr(solver.mesh, "dm_hierarchy", []) or []) + if len(hierarchy_dms) < 2: + solver._record_pc_fallback( + "custom_mg.requested_native", + requested="custom-P geometric MG over mesh.dm_hierarchy", + installed="default preconditioner", + reason="unavailable", + detail="the refinement hierarchy is gone (collapsed by a " + "remesh between build and solve)") + return None, None + level_tail = [_DMLevelView(dm, solver.mesh) for dm in hierarchy_dms[:-1]] + builder = "barycentric" + else: return None, None # nothing to inject - - # An EXPLICIT preconditioner choice beats the opportunistic pickup. Before - # this guard, `solver.preconditioner = "gamg"` on an adapt child was - # silently clobbered back to the custom-P PCMG at solve time (measured: - # both arms of test_0842's fmg-vs-gamg comparison ran pc_type=mg), so a - # user could not opt out and any FMG-vs-GAMG comparison was vacuous. - # `_pc_user_override` is the same statement in the other spelling: the - # solver's option manager has latched "the user owns this block's pc_type" - # (they wrote a pc_type of their own into petsc_options), and an - # opportunistic pickup must stand down for exactly the same reason. - # "auto" (the default) still picks up the mesh-owned hierarchy. - # NOTE the arity: this function returns a 2-tuple, never bare None — a bare - # `return` here is what turned the gate into a TypeError at the call site - # when this hunk migrated from auto_inject_custom_mg (which returns nothing) - # during the #488 x #471 merge. - if (getattr(solver, "_preconditioner", "auto") == "gamg" - or getattr(solver, "_pc_user_override", False)): - return None, None - - builder = getattr(solver.mesh, "_custom_mg_builder", "barycentric") # Retry with the RBF builder before abandoning geometric MG. The # barycentric builder has LOCAL support: it re-triangulates the coarse # DOF cloud and locates each fine DOF in one simplex, so a coarse DOF @@ -1074,7 +1169,7 @@ def build_transfers(solver, field_id=None): _attempts = [builder] + (["rbf"] if builder != "rbf" else []) h = Ps = None for _i, _b in enumerate(_attempts): - h = CustomMGHierarchy(list(coarse) + [solver.mesh], builder=_b, + h = CustomMGHierarchy(level_tail + [solver.mesh], builder=_b, field_id=field_id) try: Ps = h.build(solver) @@ -1082,6 +1177,13 @@ def build_transfers(solver, field_id=None): except Exception as exc: # pragma: no cover - defensive import warnings if _i + 1 < len(_attempts): + solver._record_pc_fallback( + "custom_mg.transfer_builder", + requested=_b, + installed=f"{_attempts[_i + 1]} (DENSE transfer)", + reason="build_failed", + detail=f"{exc}; the RBF rescue is a performance cliff — " + f"its transfer is dense (nnz/row == n_coarse), see #424") warnings.warn( f"custom_mg: {_b} transfer build failed ({exc}); " f"retrying with the '{_attempts[_i + 1]}' builder, which " @@ -1093,9 +1195,15 @@ def build_transfers(solver, field_id=None): f"it as a performance cliff and fix the cause, not the " f"symptom (#424).") continue + solver._record_pc_fallback( + "custom_mg.build", + requested=f"custom-P geometric MG ({' -> '.join(_attempts)})", + installed="default preconditioner", + reason="build_failed", + detail=str(exc)) warnings.warn( - f"custom_mg: mesh-owned FMG build failed ({exc}); using the " - "solver's default preconditioner.") + f"custom_mg: opportunistic custom-P FMG build failed ({exc}); " + "using the solver's default preconditioner.") return None, None return h, Ps @@ -1142,6 +1250,14 @@ def auto_inject_custom_mg(solver, field_id=None): # default preconditioner (round-3b annulus finding, 2026-07). if op_n > 0 and pr != op_n: import warnings + solver._record_pc_fallback( + "custom_mg.dimensional_guard", + requested="custom-P geometric MG hierarchy", + installed="default preconditioner", + reason="unavailable", + detail=f"finest transfer {pr}x{pc} is incompatible with the " + f"operator (size {op_n}); set_custom_fmg() an explicit " + f"hierarchy to override") warnings.warn( "custom_mg: mesh-owned adapt-mesh FMG transfer is incompatible " f"with this solver's operator (transfer {pr}x{pc}, operator {op_n}); " @@ -1149,10 +1265,22 @@ def auto_inject_custom_mg(solver, field_id=None): "set_custom_fmg() an explicit hierarchy to override.") return except Exception: - pass # can't check -> don't block working cases + # Sanctioned: an unreadable operator must not block working cases — + # but the skipped guard is on the record. + solver._record_pc_fallback( + "custom_mg.dimensional_guard", + requested="finest-transfer vs operator size check", + installed="unchecked (hierarchy installed anyway)", + reason="check_skipped", + detail="could not read the assembled operator to check the " + "finest transfer against it") h.install(solver, verbose=False) - solver._custom_mg = {"mode": "hierarchy", "hierarchy": h, "verbose": False} + # auto_cached marks this as a RESOLUTION product (auto/fmg install), not a + # user registration: the preconditioner setter drops it so a later explicit + # choice re-resolves instead of re-injecting this hierarchy unconditionally. + solver._custom_mg = {"mode": "hierarchy", "hierarchy": h, "verbose": False, + "auto_cached": True} def inject_custom_mg(solver): diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 0e2c9b17..adb502a8 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -1294,6 +1294,27 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal # (or a test) can assert on it instead of inferring it from timings "velocity_pc": "custom-FMG" if custom_Pl is not None else "GAMG", "schur_pre": "1/mu-mass" if Mp is not None else "selfp"} + # Mirror the degraded arms into the solver-wide fallback record (#484): + # the ctx keys above stay authoritative for this path's own tests, but + # "was anything substituted?" must be answerable in ONE place for any + # solver. Recorded only when a substitution actually happened. + if custom_Pl is None: + solver._record_pc_fallback( + "rotated.velocity_pc", + requested="custom-P geometric MG on the rotated velocity block", + installed="GAMG", + reason="unavailable", + detail="no multigrid hierarchy (set_custom_fmg or a mesh-owned " + "adapt tail) is available to the rotated path") + if Mp is None: + solver._record_pc_fallback( + "rotated.schur_pre", + requested="1/mu pressure-mass Schur preconditioner", + installed="selfp + jacobi", + reason="unavailable", + detail="the native pressure-mass Pmat block could not be built; " + "selfp degrades on curved/deformed boundaries and " + "variable viscosity") else: ksp = ctx["ksp"] nsp = ctx["nsp"] diff --git a/tests/parallel/ptest_1020_fmg_single_field_parallel.py b/tests/parallel/ptest_1020_fmg_single_field_parallel.py new file mode 100644 index 00000000..4244469e --- /dev/null +++ b/tests/parallel/ptest_1020_fmg_single_field_parallel.py @@ -0,0 +1,56 @@ +"""Parallel (MPI) test: the single-field explicit-fmg unlock (#478) at np>1. + +A distribute-then-refine hierarchy is co-partitioned (refine() never moves +points across the decomposition), so the requested-native custom-P transfers +are rank-local and the route must work unchanged in parallel: + + * explicit ``preconditioner="fmg"`` on a scalar solver -> LIVE pc_type "mg" + over every hierarchy level, converged; + * ``"auto"`` keeps GAMG with the decline recorded (route parity with serial). + +Run: + + cd tests/parallel + mpirun -np 2 python ./ptest_1020_fmg_single_field_parallel.py +""" + +import underworld3 as uw + +rank = uw.mpi.rank + +mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=0.2, refinement=2, qdegree=2, +) +assert len(mesh.dm_hierarchy) == 3 + + +def poisson_on(name): + T = uw.discretisation.MeshVariable(name, mesh, 1, degree=2) + p = uw.systems.Poisson(mesh, T) + p.constitutive_model = uw.constitutive_models.DiffusionModel + p.constitutive_model.Parameters.diffusivity = 1.0 + p.add_dirichlet_bc(0.0, "Bottom") + p.add_dirichlet_bc(1.0, "Top") + return p + + +# --- explicit fmg: custom-P geometric MG on the live PC --------------------- +p_fmg = poisson_on("Tfmg") +p_fmg.preconditioner = "fmg" +p_fmg.solve() +assert p_fmg.snes.getConvergedReason() > 0, "fmg-routed solve did not converge" +pc = p_fmg.snes.getKSP().getPC() +assert pc.getType() == "mg", f"expected live pc 'mg', got {pc.getType()!r}" +assert pc.getMGLevels() == len(mesh.dm_hierarchy) +assert "custom_mg.build" not in p_fmg.pc_fallbacks, "the transfer build degraded" + +# --- auto: GAMG + the recorded decline (route parity with serial) ----------- +p_auto = poisson_on("Tauto") +p_auto.solve() +assert p_auto.snes.getConvergedReason() > 0 +assert p_auto.snes.getKSP().getPC().getType() == "gamg" +rec = p_auto.pc_fallbacks["single_field_gate"] +assert rec["reason"] == "declined" and rec["installed"] == "gamg" + +uw.pprint("ptest_1020_fmg_single_field_parallel: PASS") diff --git a/tests/test_1020_fmg_single_field_lockout.py b/tests/test_1020_fmg_single_field_lockout.py index 3e19500e..ba1d7a24 100644 --- a/tests/test_1020_fmg_single_field_lockout.py +++ b/tests/test_1020_fmg_single_field_lockout.py @@ -1,34 +1,41 @@ #!/usr/bin/env python3 -"""Regression: native geometric FMG is locked out for single-field solvers (#276). +"""Explicit geometric MG on a single-field solver is honoured, not locked out (#478). Native geometric FMG relies on ``DMCreateInjection`` between refined DMPlex levels. PETSc can build that for the Stokes velocity sub-block but NOT for a single-field (scalar/vector) discretisation on a refined DMPlex — it fails at -solve time with err62 ("Could not locate matching functional for injection"). - -So ``preconditioner="fmg"``/``"auto"`` on a scalar/vector solver must fall back -to GAMG (and solve), never route to native FMG and crash. Geometric MG on such -a solver remains available via ``utilities.custom_mg.set_custom_fmg`` (covered by -test_1016). +solve time with err62 ("Could not locate matching functional for injection", +#276). The old gate therefore declined ``preconditioner="fmg"`` to GAMG with a +warning — locking geometric MG out of every scalar/vector solver even though +the robust custom-P route (no injection anywhere) was sitting in +``utilities.custom_mg``. + +Now: explicit ``"fmg"`` routes to custom-P transfers over the mesh's own +``dm_hierarchy`` (installed on the LIVE PC at first solve; GAMG stays in the +options DB as the degrade base). ``"auto"`` is deliberately unchanged (a +default flip needs its own validation campaign) — it keeps GAMG and records +the decline in ``pc_fallbacks``. """ import sympy import pytest import underworld3 as uw +from underworld3.utilities import custom_mg pytestmark = [pytest.mark.level_2, pytest.mark.tier_a] @pytest.fixture def annulus_hierarchy(): - # refinement>=1 => a multi-level dm_hierarchy (the FMG trigger) + # refinement>=1 => a multi-level dm_hierarchy; annulus + qdegree=3 is the + # exact geometry x degree combination that err62'd under native FMG (#276) return uw.meshing.Annulus( radiusOuter=1.0, radiusInner=0.5, cellSize=1.0 / 6.0, refinement=2, qdegree=3 ) -def _poisson(mesh): - T = uw.discretisation.MeshVariable("T", mesh, 1, degree=3) +def _poisson(mesh, name="T"): + T = uw.discretisation.MeshVariable(name, mesh, 1, degree=3) p = uw.systems.Poisson(mesh, T) p.constitutive_model = uw.constitutive_models.DiffusionModel p.constitutive_model.Parameters.diffusivity = 1.0 @@ -37,33 +44,46 @@ def _poisson(mesh): return p -def test_scalar_fmg_falls_back_and_solves(annulus_hierarchy): - """preconditioner='fmg' on a scalar solver falls back to GAMG and solves - (previously err62 DMCreateInjection).""" - assert len(annulus_hierarchy.dm_hierarchy) > 1 # hierarchy present +def test_scalar_fmg_routes_to_custom_geometric_mg(annulus_hierarchy): + """preconditioner='fmg' on a scalar solver now RUNS geometric MG (custom-P; + previously silently declined to GAMG) — read off the LIVE PC, with every + level of the native hierarchy driven.""" + assert len(annulus_hierarchy.dm_hierarchy) > 1 p = _poisson(annulus_hierarchy) p.preconditioner = "fmg" - p.solve() # must not raise err62 - assert p.snes.getConvergedReason() > 0 - # resolved to GAMG, NOT native geometric mg - assert p.petsc_options.getString("pc_type") == "gamg" - - -def test_scalar_auto_uses_gamg_not_native_fmg(annulus_hierarchy): - """preconditioner='auto' on a scalar solver with a hierarchy must NOT select - native FMG (it would crash) — it silently uses GAMG.""" - p = _poisson(annulus_hierarchy) - # 'auto' is the default; solve and confirm no native-mg pc_type - p.solve() + p.solve() # native routing would err62 here; custom-P must not assert p.snes.getConvergedReason() > 0 - assert p.petsc_options.getString("pc_type") != "mg" + pc = p.snes.getKSP().getPC() + assert pc.getType() == "mg", "explicit fmg is still being declined" + assert pc.getMGLevels() == len(annulus_hierarchy.dm_hierarchy) + # the reroute is recorded, and the build did NOT degrade + assert p.pc_fallbacks["single_field_gate"]["reason"] == "declined" + assert "custom-P" in p.pc_fallbacks["single_field_gate"]["installed"] + assert "custom_mg.build" not in p.pc_fallbacks + assert "custom_mg.transfer_builder" not in p.pc_fallbacks + + +def test_scalar_fmg_degrades_readably_when_the_build_fails( + annulus_hierarchy, monkeypatch): + """Negative control for the unlock: if the transfer build fails, the solve + must still converge on the GAMG base left in the options DB, and say so.""" + def broken(coarse_coords, fine_coords): + raise RuntimeError("synthetic builder failure (test probe)") + + monkeypatch.setitem(custom_mg._BUILDERS, "barycentric", broken) + monkeypatch.setitem(custom_mg._BUILDERS, "rbf", broken) + p = _poisson(annulus_hierarchy, name="Tdeg") + p.preconditioner = "fmg" + with pytest.warns(UserWarning, match="custom-P FMG build failed"): + p.solve() + assert p.snes.getConvergedReason() > 0 # stability oracle + assert p.snes.getKSP().getPC().getType() == "gamg" + assert p.pc_fallbacks["custom_mg.build"]["reason"] == "build_failed" -def test_cartesian_box_scalar_also_falls_back(): - """The lockout is single-field-wide, not geometry-specific: native single- - field FMG is fragile across geometry×degree×refinement (it errors 62 even on - some flat high-degree boxes), so a scalar solver on a flat box also falls - back to GAMG and solves rather than risk the crash.""" +def test_cartesian_box_scalar_fmg_also_unlocked(): + """The unlock is single-field-wide, not geometry-specific: the flat + high-degree box that also err62'd under native FMG runs custom-P too.""" box = uw.meshing.UnstructuredSimplexBox( minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.4, refinement=1, qdegree=3 ) @@ -75,14 +95,28 @@ def test_cartesian_box_scalar_also_falls_back(): p.add_dirichlet_bc(0.0, "Bottom") p.add_dirichlet_bc(1.0, "Top") p.preconditioner = "fmg" - p.solve() # must not raise err62 + p.solve() assert p.snes.getConvergedReason() > 0 - assert p.petsc_options.getString("pc_type") == "gamg" + pc = p.snes.getKSP().getPC() + assert pc.getType() == "mg" + assert pc.getMGLevels() == len(box.dm_hierarchy) + + +def test_scalar_auto_still_uses_gamg_and_records_the_decline(annulus_hierarchy): + """'auto' is deliberately NOT flipped by #478 (that default change needs its + own np2/np4 + adaptivity validation): it keeps GAMG, and the decline is the + recorded migration probe.""" + p = _poisson(annulus_hierarchy, name="Tauto") + p.solve() + assert p.snes.getConvergedReason() > 0 + assert p.snes.getKSP().getPC().getType() == "gamg" + assert p.pc_fallbacks["single_field_gate"]["reason"] == "declined" + assert p.pc_fallbacks["single_field_gate"]["installed"] == "gamg" def test_stokes_velocity_fmg_still_selected(annulus_hierarchy): - """The lockout is single-field only: the Stokes velocity sub-block still gets - native geometric FMG on the same hierarchy.""" + """The single-field reroute does not touch Stokes: the velocity sub-block + still gets native geometric FMG on the same hierarchy.""" v = uw.discretisation.MeshVariable("v", annulus_hierarchy, annulus_hierarchy.dim, degree=2) pp = uw.discretisation.MeshVariable("p", annulus_hierarchy, 1, degree=1) stokes = uw.systems.Stokes(annulus_hierarchy, velocityField=v, pressureField=pp) @@ -92,3 +126,39 @@ def test_stokes_velocity_fmg_still_selected(annulus_hierarchy): stokes._build(False, False, None) # the velocity sub-block keeps native geometric multigrid assert stokes.petsc_options.getString("fieldsplit_velocity_pc_type") == "mg" + + +def test_options_db_and_record_stay_honest_between_build_and_solve(annulus_hierarchy): + """Between _build and the first solve the options DB deliberately says gamg + (the safe degrade base) while the plan is custom-P mg. The #471 gating rule + plus the record's 'resolved at first solve' phrasing keep the report honest; + after the solve the DB carries the installed bundle.""" + p = _poisson(annulus_hierarchy, name="Thon") + p.preconditioner = "fmg" + p._build() + assert p.petsc_options.getString("pc_type") == "gamg" # degrade base + rec = p.pc_fallbacks["single_field_gate"] + assert "resolved at first solve" in rec["installed"] + p.solve() + assert p.preconditioner_settings["pc_type"] == "mg" # bundle written + assert p.snes.getKSP().getPC().getType() == "mg" + + +def test_switch_to_gamg_after_fmg_is_honoured(annulus_hierarchy): + """Setting preconditioner='gamg' AFTER an fmg solve must install a real + GAMG. The #534 review measured the failure this pins: the fmg install + cached solver._custom_mg, auto_inject re-installed it unconditionally, and + the explicit gamg choice was silently unreachable — with pc_fallbacks + EMPTY, a false-clean record (the recorder's own contract violated).""" + p = _poisson(annulus_hierarchy, name="Tsw") + p.preconditioner = "fmg" + p.solve() + assert p.snes.getKSP().getPC().getType() == "mg" + p.preconditioner = "gamg" + p.solve() + pc = p.snes.getKSP().getPC() + assert pc.getType() == "gamg", ( + "explicit gamg after an fmg solve was overridden by the cached " + "custom-P hierarchy") + # honoured explicitly => no fallback story to tell + assert "single_field_gate" not in p.pc_fallbacks diff --git a/tests/test_1022_pc_fallback_observability.py b/tests/test_1022_pc_fallback_observability.py new file mode 100644 index 00000000..12c5b8de --- /dev/null +++ b/tests/test_1022_pc_fallback_observability.py @@ -0,0 +1,213 @@ +#!/usr/bin/env python3 +"""Every preconditioner fallback leaves a readable record (#484). + +The solver ecosystem degrades gracefully in many places (single-field FMG +gate, missing hierarchy, transfer-build failures, guard skips, forced +Galerkin). Before #484, 10 of the 12 fallback sites left no queryable state — +a warning at best, silence at worst — so neither a user nor a test could ask +"did I get what I asked for?". Now every site writes into +``solver.pc_fallbacks`` through one recorder, with a fixed reason vocabulary: +``unavailable``, ``declined``, ``build_failed``, ``check_skipped``, ``forced``. + +House rule: every probe is proven to fire AND proven silent on healthy state. +""" +import numpy as np +import pytest + +import underworld3 as uw +from underworld3.utilities import custom_mg + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +REASONS = {"unavailable", "declined", "build_failed", "check_skipped", "forced"} + + +def _poisson_on(mesh, degree=1, name="T"): + T = uw.discretisation.MeshVariable(name, mesh, 1, degree=degree) + p = uw.systems.Poisson(mesh, T) + p.constitutive_model = uw.constitutive_models.DiffusionModel + p.constitutive_model.Parameters.diffusivity = 1.0 + p.add_dirichlet_bc(0.0, "Bottom") + p.add_dirichlet_bc(1.0, "Top") + return p + + +def _box(cellSize=0.4, refinement=0, qdegree=2): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=cellSize, refinement=refinement, qdegree=qdegree) + + +@pytest.fixture +def adapt_child(): + mesh = _box(cellSize=0.2, refinement=1) + + def metric(points): + r = np.linalg.norm(np.asarray(points) - 0.5, axis=1) + return 1.0 / np.where(r < 0.25, 0.05, 0.2) ** 2 + + return mesh.adapt(metric, max_levels=1) + + +# --------------------------------------------------------------------------- # +# The global negative control: a clean solve records NOTHING +# --------------------------------------------------------------------------- # +def test_clean_solve_has_empty_fallback_record(): + """No hierarchy, default preconditioner: the resolution is exactly what was + asked for, so the record must be empty — this is the state every other + test's 'record present' assertion is measured against.""" + p = _poisson_on(_box(refinement=0)) + p.solve() + assert p.snes.getConvergedReason() > 0 + assert p.pc_fallbacks == {} + + +def test_explicit_gamg_records_nothing(): + """Explicit 'gamg' on a hierarchy mesh is honoured verbatim — no decline, + no record (the negative control for the auto-decline probe).""" + p = _poisson_on(_box(refinement=1)) + p.preconditioner = "gamg" + p.solve() + assert p.snes.getConvergedReason() > 0 + assert "single_field_gate" not in p.pc_fallbacks + assert "no_hierarchy" not in p.pc_fallbacks + + +# --------------------------------------------------------------------------- # +# Gate + no-hierarchy sites (the pyx _apply_preconditioner_options sites) +# --------------------------------------------------------------------------- # +def test_fmg_without_hierarchy_records_unavailable_and_still_warns(): + p = _poisson_on(_box(refinement=0)) + p.preconditioner = "fmg" + with pytest.warns(UserWarning, match="no refinement hierarchy"): + p.solve() + rec = p.pc_fallbacks["no_hierarchy"] + assert rec["reason"] == "unavailable" + assert rec["installed"] == "gamg" + assert p.petsc_options.getString("pc_type") == "gamg" + # negative control: a hierarchy mesh must NOT carry this record + p2 = _poisson_on(_box(refinement=1), name="T2") + p2.preconditioner = "fmg" + p2.solve() + assert "no_hierarchy" not in p2.pc_fallbacks + + +def test_auto_decline_is_recorded_but_stays_silent(recwarn): + """'auto' on a single-field solver with a hierarchy declines geometric FMG + by design. The decline is now on the record — but it must NOT gain a + warning (auto's silence is a stability contract).""" + p = _poisson_on(_box(refinement=2)) + p.solve() + assert p.snes.getConvergedReason() > 0 + rec = p.pc_fallbacks["single_field_gate"] + assert rec["reason"] == "declined" + assert rec["installed"] == "gamg" + fmg_warnings = [w for w in recwarn.list + if "FMG" in str(w.message) or "GAMG" in str(w.message)] + assert fmg_warnings == [], "the auto decline must not warn" + + +def test_reason_vocabulary_is_closed(): + """Whatever a solve records, the reasons come from the documented set.""" + p = _poisson_on(_box(refinement=1)) + p.preconditioner = "fmg" + p.solve() + for site, rec in p.pc_fallbacks.items(): + assert rec["reason"] in REASONS, (site, rec) + assert set(rec) == {"requested", "installed", "reason", "detail"} + + +# --------------------------------------------------------------------------- # +# custom_mg sites (transfer builder ladder, guards) +# --------------------------------------------------------------------------- # +def test_transfer_builder_failure_records_the_rbf_rescue(adapt_child, monkeypatch): + """Barycentric build fails -> RBF rescue: recorded as build_failed, the + solve still converges on geometric MG. Degree 2 so the exact recorded + (vertex-level) transfers cannot bypass the builder.""" + def broken(coarse_coords, fine_coords): + raise RuntimeError("synthetic barycentric failure (test probe)") + + monkeypatch.setitem(custom_mg._BUILDERS, "barycentric", broken) + p = _poisson_on(adapt_child, degree=2, name="Trbf") + with pytest.warns(UserWarning, match="retrying with the 'rbf' builder"): + p.solve() + assert p.snes.getConvergedReason() > 0 + rec = p.pc_fallbacks["custom_mg.transfer_builder"] + assert rec["reason"] == "build_failed" + assert "rbf" in rec["installed"] + assert "custom_mg.build" not in p.pc_fallbacks # the rescue succeeded + + +def test_total_transfer_failure_records_and_solves_on_default_pc( + adapt_child, monkeypatch): + def broken(coarse_coords, fine_coords): + raise RuntimeError("synthetic builder failure (test probe)") + + monkeypatch.setitem(custom_mg._BUILDERS, "barycentric", broken) + monkeypatch.setitem(custom_mg._BUILDERS, "rbf", broken) + p = _poisson_on(adapt_child, degree=2, name="Tfail") + with pytest.warns(UserWarning, match="custom-P FMG build failed"): + p.solve() + # stability oracle: the degrade path must still converge + assert p.snes.getConvergedReason() > 0 + rec = p.pc_fallbacks["custom_mg.build"] + assert rec["reason"] == "build_failed" + assert rec["installed"] == "default preconditioner" + assert p.snes.getKSP().getPC().getType() != "mg" + + +def test_unpatched_adapt_child_solve_records_no_build_failure(adapt_child): + """Negative control for both probes above: the healthy adapt-child pickup + installs geometric MG with no build-failure record.""" + p = _poisson_on(adapt_child, degree=2, name="Tok") + p.solve() + assert p.snes.getConvergedReason() > 0 + assert p.snes.getKSP().getPC().getType() == "mg" + assert "custom_mg.transfer_builder" not in p.pc_fallbacks + assert "custom_mg.build" not in p.pc_fallbacks + assert "custom_mg.finest_operator_check" not in p.pc_fallbacks + + +def test_finest_operator_guard_skip_is_recorded(): + """The finest-operator guard's early return (operator unreadable) was a + silent swallow; it now records check_skipped. Driven directly: a solver + with no SNES yet is exactly the unreadable-operator state.""" + p = _poisson_on(_box(refinement=0), name="Tguard") + assert not hasattr(p, "snes") or p.snes is None or True # un-built solver + custom_mg.CustomMGHierarchy._assert_finest_matches_operator( + p, finest_map=np.arange(4), parallel=False) + rec = p.pc_fallbacks["custom_mg.finest_operator_check"] + assert rec["reason"] == "check_skipped" + + +# --------------------------------------------------------------------------- # +# Forced Galerkin site +# --------------------------------------------------------------------------- # +def test_user_mg_without_galerkin_is_forced_and_recorded(): + p = _poisson_on(_box(refinement=1), name="Tgal") + p.petsc_options["pc_type"] = "mg" + with pytest.warns(UserWarning, match="forcing pc_mg_galerkin=both"): + p._build() + rec = p.pc_fallbacks["galerkin_forced"] + assert rec["reason"] == "forced" + assert p.petsc_options.getString("pc_mg_galerkin") == "both" + # negative control: setting the key yourself satisfies the requirement + p2 = _poisson_on(_box(refinement=1), name="Tgal2") + p2.petsc_options["pc_type"] = "mg" + p2.petsc_options["pc_mg_galerkin"] = "both" + p2._build() + assert "galerkin_forced" not in p2.pc_fallbacks + + +# --------------------------------------------------------------------------- # +# Reset rule +# --------------------------------------------------------------------------- # +def test_record_resets_when_the_preconditioner_re_resolves(): + p = _poisson_on(_box(refinement=0), name="Trst") + p.preconditioner = "fmg" + with pytest.warns(UserWarning, match="no refinement hierarchy"): + p.solve() + assert "no_hierarchy" in p.pc_fallbacks + p.preconditioner = "gamg" + p.solve() + assert p.pc_fallbacks == {}, "stale records must not survive re-resolution" diff --git a/tests/test_1023_saddle_tolerance_ownership.py b/tests/test_1023_saddle_tolerance_ownership.py new file mode 100644 index 00000000..c96bdee7 --- /dev/null +++ b/tests/test_1023_saddle_tolerance_ownership.py @@ -0,0 +1,167 @@ +#!/usr/bin/env python3 +"""Tolerance-key ownership on the saddle-point solvers is honest (#483). + +Two defects, both in the middle ground between "owned" and "settable": + +1. ``snes_rtol`` / ``ksp_atol`` were re-pushed by every solve + (``_reassert_outer_tolerances``), silently discarding a user's explicit + value — documented as settable, actually owned. They are now routed + through the same latch that made ``snes_max_it`` reachable (ruling D18): + the framework keeps re-asserting them until the user sets one, after which + the user's value is honoured. + +2. ``Stokes`` and ``Stokes_Constrained`` derived different option keys from + ``tolerance`` through two unrelated hand-rolled code paths. Both now run + one mechanism (``_derive_tolerance_margins``) over a per-class table + (``_TOLERANCE_DERIVED_KEYS``), and the difference is documented: the base + class derives the inner fieldsplit margins, Constrained derives the outer + ``ksp_rtol`` and the Eisenstat-Walker pins. + +Every arm reads LIVE PETSc objects (or the options DB where EW makes the live +value per-iteration) after TWO solves — the second solve is the one that used +to clobber. +""" +import sympy +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _mesh(): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25, qdegree=3) + + +def _stokes(mesh, constrained=False, suffix=""): + v = uw.discretisation.MeshVariable(f"v{suffix}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"p{suffix}", mesh, 1, degree=1) + cls = uw.systems.Stokes_Constrained if constrained else uw.systems.Stokes + stokes = cls(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((1.0, 0.0), "Top") + stokes.bodyforce = sympy.Matrix([0, -1]) + return stokes + + +def _solve_twice(stokes): + stokes.solve() + stokes.solve() # the second solve is the one that used to clobber + + +# --------------------------------------------------------------------------- # +# Reachability: user overrides survive the solve, on BOTH classes +# --------------------------------------------------------------------------- # +def test_stokes_outer_and_inner_overrides_are_honoured(): + mesh = _mesh() + stokes = _stokes(mesh, suffix="a") + stokes.tolerance = 1e-6 + stokes.petsc_options["snes_rtol"] = 3.0e-3 + stokes.petsc_options["ksp_atol"] = 7.0e-9 + stokes.petsc_options["fieldsplit_velocity_ksp_rtol"] = 0.02 + stokes.petsc_options["fieldsplit_pressure_ksp_rtol"] = 0.005 + _solve_twice(stokes) + # outer keys, live off the SNES/KSP + assert stokes.snes.getTolerances()[0] == pytest.approx(3.0e-3) + assert stokes.snes.getKSP().getTolerances()[1] == pytest.approx(7.0e-9) + # inner fieldsplit rtols, live off the sub-KSPs (field 0 = velocity) + vel_ksp, pres_ksp = stokes.snes.getKSP().getPC().getFieldSplitSubKSP() + assert vel_ksp.getTolerances()[0] == pytest.approx(0.02) + assert pres_ksp.getTolerances()[0] == pytest.approx(0.005) + + +def test_stokes_defaults_still_owned_without_overrides(): + """Negative control: with NO user overrides the ownership must still work — + the latch must not simply have stopped the framework pushing.""" + mesh = _mesh() + stokes = _stokes(mesh, suffix="b") + stokes.tolerance = 1e-6 + _solve_twice(stokes) + assert stokes.snes.getTolerances()[0] == pytest.approx(1e-6) + assert stokes.snes.getKSP().getTolerances()[1] == pytest.approx(1e-12) + vel_ksp, pres_ksp = stokes.snes.getKSP().getPC().getFieldSplitSubKSP() + assert vel_ksp.getTolerances()[0] == pytest.approx(1e-6 * 0.033) + assert pres_ksp.getTolerances()[0] == pytest.approx(1e-6 * 0.1) + + +def test_constrained_outer_overrides_and_ew_pins_are_honoured(): + mesh = _mesh() + stokes = _stokes(mesh, constrained=True, suffix="c") + stokes.tolerance = 1e-6 + stokes.petsc_options["snes_rtol"] = 2.0e-3 + stokes.petsc_options["ksp_atol"] = 5.0e-9 + stokes.petsc_options["snes_ksp_ew_rtol0"] = 4.0e-5 # derived-key override + _solve_twice(stokes) + assert stokes.snes.getTolerances()[0] == pytest.approx(2.0e-3) + assert stokes.snes.getKSP().getTolerances()[1] == pytest.approx(5.0e-9) + # EW re-picks the live ksp rtol per iteration, so the DB is the contract + # for the pins: the user's value must survive both solves. + assert float(stokes.petsc_options.getString("snes_ksp_ew_rtol0")) \ + == pytest.approx(4.0e-5) + # the pin the user left alone keeps its tolerance-derived value + assert float(stokes.petsc_options.getString("snes_ksp_ew_rtolmax")) \ + == pytest.approx(1e-6 * 0.1) + + +def test_constrained_defaults_still_owned_without_overrides(): + mesh = _mesh() + stokes = _stokes(mesh, constrained=True, suffix="d") + stokes.tolerance = 1e-6 + _solve_twice(stokes) + assert stokes.snes.getTolerances()[0] == pytest.approx(1e-6) + assert stokes.snes.getKSP().getTolerances()[1] == pytest.approx(1e-12) + for key in ("ksp_rtol", "snes_ksp_ew_rtol0", "snes_ksp_ew_rtolmax"): + assert float(stokes.petsc_options.getString(key)) \ + == pytest.approx(1e-6 * 0.1), key + + +# --------------------------------------------------------------------------- # +# Order rule: set-time derivation vs later user override +# --------------------------------------------------------------------------- # +def test_tolerance_wins_when_set_after_the_user_key_and_vice_versa(): + mesh = _mesh() + # arm 1: user key BEFORE tolerance= -> the set-time derivation wins + s1 = _stokes(mesh, suffix="e") + s1.petsc_options["snes_rtol"] = 3.0e-3 + s1.tolerance = 1e-6 + _solve_twice(s1) + assert s1.snes.getTolerances()[0] == pytest.approx(1e-6) + # arm 2: user key AFTER tolerance= -> the user wins across two solves + s2 = _stokes(mesh, suffix="f") + s2.tolerance = 1e-6 + s2.petsc_options["snes_rtol"] = 3.0e-3 + _solve_twice(s2) + assert s2.snes.getTolerances()[0] == pytest.approx(3.0e-3) + + +# --------------------------------------------------------------------------- # +# snes_max_it: the original D18 latch, now running on the shared mechanism +# --------------------------------------------------------------------------- # +def test_snes_max_it_latch_regression(): + mesh = _mesh() + s1 = _stokes(mesh, suffix="g") + s1.petsc_options["snes_max_it"] = 7 + _solve_twice(s1) + assert s1.snes.getTolerances()[3] == 7 + # negative control: untouched -> the framework default (50) is live + s2 = _stokes(mesh, suffix="h") + _solve_twice(s2) + assert s2.snes.getTolerances()[3] == 50 + + +# --------------------------------------------------------------------------- # +# Docs/code lockstep: the tables are exactly what the docstrings promise +# --------------------------------------------------------------------------- # +def test_derived_key_tables_match_the_documentation(): + assert uw.systems.Stokes._TOLERANCE_DERIVED_KEYS == { + "fieldsplit_pressure_ksp_rtol": 0.1, + "fieldsplit_velocity_ksp_rtol": 0.033, + } + assert uw.systems.Stokes_Constrained._TOLERANCE_DERIVED_KEYS == { + "ksp_rtol": 0.1, + "snes_ksp_ew_rtol0": 0.1, + "snes_ksp_ew_rtolmax": 0.1, + } diff --git a/tests/test_1024_multiplier_schur_pc.py b/tests/test_1024_multiplier_schur_pc.py new file mode 100644 index 00000000..324ca305 --- /dev/null +++ b/tests/test_1024_multiplier_schur_pc.py @@ -0,0 +1,162 @@ +#!/usr/bin/env python3 +"""multiplier_schur_pc: verified live where it can be, instrumented where it is inert (#486). + +The flag swaps only the Pmat (h,h) block of each Lagrange-multiplier field +(``_setup_solver``: the DS JacobianPreconditioner term becomes pressure's +``1/mu`` mass instead of the screening ``eps`` mass). Tracing PETSc's +fieldsplit.c resolves where that block is actually read: + +* ``pc_fieldsplit_schur_precondition = "a11"`` — the Schur preconditioner is + the grouped ``[p,h]`` **Pmat** block: the flag is LIVE (matrix-level oracle + below, immune to "both converge in 4 iterations"). +* class defaults (``selfp`` + ``diag_use_amat``) — Sp is assembled from + **Amat** sub-blocks; the Pmat (h,h) block is provably never read: the flag + is INERT. That inertness (the issue's original measurement) is codified as + a negative control pinning the PETSc semantics, and the solver now records + + warns when the opt-in cannot reach the PC. +""" +import numpy as np +import pytest +import scipy.sparse as sp +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _schur_pre_matrix(solver): + """The assembled Schur preconditioner matrix Sp, as scipy CSR.""" + schur_ksp = solver.snes.getKSP().getPC().getFieldSplitSubKSP()[1] + Sp = schur_ksp.getOperators()[1] + i, j, v = Sp.getValuesCSR() + return sp.csr_matrix((v, j, i), shape=Sp.getSize()) + + +def _rel_diff(A, B): + denom = sp.linalg.norm(A) + 1e-300 + return sp.linalg.norm(A - B) / denom + + +def _constrained(mesh, name, schur_pre, flag): + """Open-top box, TWO constraint arcs (Left/Right), viscosity contrast, + linear (ksponly) so the extracted PC is free of nonlinear noise.""" + xx, yy = mesh.X + v = uw.discretisation.MeshVariable(f"v{name}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"p{name}", mesh, 1, degree=1) + s = uw.systems.Stokes_Constrained(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = sympy.exp( + sympy.log(1000.0) * xx) # 3 decades across the box + s.bodyforce = sympy.Matrix( + [0.0, sympy.sin(sympy.pi * xx) * sympy.cos(sympy.pi * yy)]) + s.add_dirichlet_bc((0.0, 0.0), "Bottom") + s.add_constraint_bc(0.0, "Left", normal=sympy.Matrix([[-1.0, 0.0]])) + s.add_constraint_bc(0.0, "Right", normal=sympy.Matrix([[1.0, 0.0]])) + s.petsc_options["snes_type"] = "ksponly" + if schur_pre is not None: + s.petsc_options["pc_fieldsplit_schur_precondition"] = schur_pre + s.multiplier_schur_pc = flag + return s + + +@pytest.fixture(scope="module") +def solves(): + """All four (regime x flag) arms on one mesh, warnings captured per arm.""" + import warnings as _w + + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0, 0), maxCoords=(1, 1), cellSize=0.2, qdegree=3) + out = {} + for key, schur_pre, flag in ( + ("selfp_off", None, False), + ("selfp_on", None, True), + ("a11_off", "a11", False), + ("a11_on", "a11", True)): + s = _constrained(mesh, key, schur_pre, flag) + with _w.catch_warnings(record=True) as caught: + _w.simplefilter("always") + s.solve() + out[key] = (s, [str(w.message) for w in caught]) + return out + + +# --------------------------------------------------------------------------- # +# Live regime: under a11 the flag provably changes the Schur preconditioner +# --------------------------------------------------------------------------- # +def test_flag_is_live_under_a11(solves): + """Matrix-level oracle: the assembled Schur pre differs flag-on vs + flag-off under a11 — the DS JacobianPreconditioner path for multiplier + fields reaches the PC. (If this ever finds NO difference, that path is + broken upstream of the flag — file that as its own bug.)""" + Sp_off = _schur_pre_matrix(solves["a11_off"][0]) + Sp_on = _schur_pre_matrix(solves["a11_on"][0]) + assert _rel_diff(Sp_on, Sp_off) > 1e-6, ( + "multiplier_schur_pc changed nothing under a11 — the multiplier " + "JacobianPreconditioner path is broken upstream of the flag") + + +def test_probe_is_regime_aware_under_a11(solves): + """Where the flag IS live, there must be no inertness record or warning.""" + s_on, warned = solves["a11_on"] + assert "multiplier_schur_pc" not in s_on.pc_fallbacks + assert not [w for w in warned if "multiplier_schur_pc" in w] + + +# --------------------------------------------------------------------------- # +# Inert regime codified (the issue's measurement, now a pinned control) +# --------------------------------------------------------------------------- # +def test_flag_is_inert_under_class_defaults_and_says_so(solves): + """selfp + diag_use_amat: Sp identical flag-on/off (pins the fieldsplit.c + semantics we traced — if PETSc changes FieldSplitSchurPre, this fails), + and the solver records + warns that the opt-in cannot reach the PC.""" + Sp_off = _schur_pre_matrix(solves["selfp_off"][0]) + s_on, warned = solves["selfp_on"] + Sp_on = _schur_pre_matrix(s_on) + assert _rel_diff(Sp_on, Sp_off) == 0.0, ( + "selfp read the Pmat (h,h) block — the PETSc fieldsplit semantics " + "this class's defaults rely on have changed") + rec = s_on.pc_fallbacks["multiplier_schur_pc"] + assert rec["reason"] == "declined" + assert "selfp" in rec["installed"] + assert [w for w in warned if "multiplier_schur_pc" in w], ( + "an explicit opt-in doing nothing must warn") + + +def test_no_flag_means_no_record_and_no_warning(solves): + """Negative control: with the flag off, the probe stays silent.""" + s_off, warned = solves["selfp_off"] + assert "multiplier_schur_pc" not in s_off.pc_fallbacks + assert not [w for w in warned if "multiplier_schur_pc" in w] + + +# --------------------------------------------------------------------------- # +# Setter contract +# --------------------------------------------------------------------------- # +def test_toggling_after_a_solve_forces_ds_reregistration(solves): + s_off, _ = solves["selfp_off"] + assert s_off.is_setup + s_off.multiplier_schur_pc = True + assert not s_off.is_setup, ( + "toggling multiplier_schur_pc must force the DS term to re-register") + s_off.multiplier_schur_pc = False # restore for any later use + + +def test_probe_reads_the_flag_value_not_its_presence(): + """diag_use_amat set to FALSE under selfp means the Pmat (h,h) block IS + read — the opt-in is live and the inertness warning must stay silent. + The #534 review measured the Schur pre differing by rel-Frobenius 0.30 + flag-on/off in this regime while the probe (reading hasName, not the + value) still recorded 'declined'.""" + import warnings as _warnings + + mesh = uw.meshing.StructuredQuadBox(elementRes=(6, 6)) + s = _constrained(mesh, "valrd", None, True) + s.petsc_options["pc_fieldsplit_diag_use_amat"] = False + with _warnings.catch_warnings(record=True) as caught: + _warnings.simplefilter("always") + s.solve() + assert "multiplier_schur_pc" not in s.pc_fallbacks, ( + "probe read the key's presence, not its value: flag=False makes the " + "opt-in live, not inert") + assert not [w for w in caught if "multiplier_schur_pc" in str(w.message)]