Make the fall off the multigrid loud, and move the penalty into the total stress (#625) - #632
Merged
Merged
Conversation
Implements the FMG-conditional penalty: penalty=10 when the custom-P multigrid will run, 0 otherwise, with the bet confirmed after setup and the harmful GAMG+penalty pairing warning and recording in pc_fallbacks. All five pairing cases behave as specified (tests/test_0206). It should NOT be merged as it stands. Three existing tests fail, and they are right: test_1017, test_0835 and test_0836 each assert that the FMG solution and the GAMG solution of the same problem agree to 1e-4 -- that a preconditioner changes the path, not the answer. Making an operator term depend on the preconditioner breaks that invariant, and the measured disagreement is 5.1e-4 to 1.4e-3. Kept as the record of what the conditional design costs. Underworld development team with AI support from Claude Code
NOT READY TO MERGE — see the blocker below. The penalty is now a discretisation choice applied unconditionally (`Stokes.DEFAULT_PENALTY = 10`), not one selected by the preconditioner. Selecting it from the solver was tried first and rejected: test_1017, test_0835 and test_0836 each assert that the FMG and GAMG solutions of the same problem agree to 1e-4 — a preconditioner changes the path, not the answer — and a conditional penalty broke all three by 5.1e-4 to 1.4e-3. Those three pass again here. The real defect is addressed instead. A velocity block that drops off the custom-P multigrid because the mesh has no hierarchy now warns and records `velocity.fell_back_from_fmg` in pc_fallbacks, naming the remedy. That fallback is the start of the recurring failure: GAMG degrades under refinement until it hits its iteration cap, a capped velocity solve corrupts S = -B A^-1 B^T, and the pressure block stalls behind it (976 s against 25.6 s at h=1/30, #625). An explicit preconditioner="gamg" is a decision, not a fallback, and is left alone. Measured basis (SolCx, eta 1e6, P2-P0disc, 2592 cells): under FMG the penalty is 21% faster, cuts the Schur count per application 59 -> 18 and total velocity iterations 546 -> 270. The accuracy cost is a consistent perturbation rather than a changed answer — same convergence rate, gap shrinking 1.102 -> 1.087 -> 1.066 over three refinements. BLOCKER, found by the suite: test_1018's spherical dynamic topography moves 0.4192 -> 0.3021, 28%, on the vertex-sampled surface value (the facet-integrated surface and CMB values are unaffected). With the penalty active the recovered p is the Lagrange multiplier, so the constraint reaction sigma_nn needs the documented correction p_mech = p - lambda*mu*(div u), and the rotated free-slip reaction path does not apply it. Defaulting the penalty before that is fixed would silently change everyone's dynamic topography. 516 of 517 solver tests pass; the failure above is the one. Underworld development team with AI support from Claude Code
The penalty was added at the weak-form assembly point, F1 = stress + penalty term, OUTSIDE the stress definition. So the operator being solved carried it while `stress`, `stress_1d` and everything downstream did not, and each consumer had to apply `p_mech = p - lambda*mu*(div u)` by hand -- which the penalty docstring duly tells users to do for pressure-dependent laws. The term is ISOTROPIC, a modification of the pressure, so it belongs in the total stress. `stress` is now `stress_deviator - (p - lambda*mu*(div u))*I` and F1 is just `self.stress`. The weak form is unchanged; every consumer of the total stress now inherits the correction through the same route. `stress_deviator` keeps excluding it, correctly -- it is not deviatoric. That also keeps it out of the viscoelastic history, which tracks the deviator via `constitutive_model.flux`, where it does not belong. This does NOT fix the test_1018 topography shift, and the byte-identical result says why: `boundary_normal_traction` recovers sigma_nn from the Cartesian nodal reaction r_c = A.u - b, not from the symbolic stress. That residual already carried the penalty, so that path was never inconsistent. See the next commit message for what the 28% actually is. 24 related tests pass. With DEFAULT_PENALTY at its shipped value the refactor is a no-op; it only bites when a penalty is set. Underworld development team with AI support from Claude Code
Everything for the change is in place; the value is the one thing held back. At 10 the spherical dynamic topography recovered from the rotated free-slip reaction drops 0.4192 -> 0.3021 -- 28% -- on the VERTEX-SAMPLED value, while the facet-integrated value stays correct (test_1018, whose docstring is "reaction loads must be divided by boundary mass to recover POINTWISE stress"). So augmentation corrupts the de-smearing from reaction loads to pointwise stress, presumably because lambda*mu*(div u) is non-zero cell-by-cell for P2-P0 and averages out over a facet integral but not at a vertex. Not a bookkeeping error that the p_mech move fixed: that recovery reads the Cartesian nodal reaction r_c = A.u - b, which already carried the penalty, and the measured value was byte-identical before and after. Dynamic topography is the main product of that machinery, so the default stays 0 until it is understood. test_0206's invariant test patches DEFAULT_PENALTY to 7.0 for its duration, so it still means something while the shipped value is 0 -- otherwise it would pass by comparing 0 to 0 however the value were chosen. 727 solver tests pass. Underworld development team with AI support from Claude Code
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What lands
Two things that stop the recurring "the Schur solve wandered again" session, plus
the machinery for the third, held back deliberately.
1. A velocity block that falls off the multigrid says so. FMG needs a mesh
hierarchy. Without one the block drops to GAMG with nothing said — measured,
refinement=0gives one hierarchy level andgamg,refinement=2givesmg.That fallback is the start of the failure: GAMG degrades under refinement until
it hits its iteration cap, a capped velocity solve corrupts
S = -B A^-1 B^T,and the pressure block stalls behind it. Measured at h=1/30, changing only
fieldsplit_velocity_ksp_max_it200 -> 5000: 976 s -> 25.6 s, outer 44 -> 2,pressure 200/application -> 30, identical answer.
Now warns and records
velocity.fell_back_from_fmginpc_fallbacks, naming theremedy. An explicit
preconditioner="gamg"is a decision, not a fallback, and isleft alone — a warning that fires on a deliberate choice trains people to ignore
it.
2. The grad-div penalty moves into the total stress, out of
F1. It wasadded at the assembly point,
F1 = stress + penalty term, outside the stressdefinition — so the operator carried it while
stress,stress_1dand everydownstream consumer did not, and each had to apply
p_mech = p - lambda*mu*(div u)by hand. The term is isotropic, amodification of the pressure, so
stressis nowstress_deviator - (p - lambda*mu*(div u))*IandF1is justself.stress.The weak form is unchanged; consumers inherit the correction through the one
route.
stress_deviatorkeeps excluding it, correctly — it is not deviatoric, whichalso keeps it out of the viscoelastic history (which tracks the deviator through
constitutive_model.flux).3.
Stokes.DEFAULT_PENALTY, held at 0. The measured case for turning it onis strong: under FMG,
penalty = 10is 21% faster, cuts the Schur count perapplication 59 -> 18 and total velocity iterations 546 -> 270, and the accuracy
cost is a consistent perturbation — same convergence rate, gap shrinking
1.102 -> 1.087 -> 1.066 over three refinements.
It is held anyway, because at 10 the spherical dynamic topography from the
rotated free-slip reaction drops 0.4192 -> 0.3021, 28%, on the
vertex-sampled value while the facet-integrated value stays correct
(
test_1018, whose docstring is "reaction loads must be divided by boundarymass to recover pointwise stress"). Augmentation corrupts the de-smearing from
reaction loads to pointwise stress. Dynamic topography is the main product of
that machinery, so the constant stays 0 until that is understood; everything
else is in place and the change is this one value.
Rejected design, recorded
Selecting the penalty from the preconditioner — on for FMG, off for GAMG — was
implemented first and rejected.
test_1017,test_0835andtest_0836eachassert that the FMG and GAMG solutions of the same problem agree to 1e-4: a
preconditioner changes the path, not the answer. A conditional operator term
broke all three, by 5.1e-4 to 1.4e-3. Kept on
experiment/penalty-paired-to-preconditioner, marked do-not-merge.Testing
727 solver tests pass (
test_100*,test_101*,test_102*,test_020*,test_083*,test_084*), includingtest_1018and the three invariant tests.tests/test_0206_automatic_penalty_pairing.py— 5 tests: the default ispreconditioner-independent (patching
DEFAULT_PENALTYto 7.0 for the duration,so it still means something while the shipped value is 0), an explicit penalty is
honoured, falling off the multigrid is loud and recorded, a hierarchy means no
warning (the negative control), and an explicit GAMG choice is not nagged about.
Not covered
The penalty appears with the same
F1bolt-on shape inStokes_Constrained(line 3408) and
NavierStokes(line 4930). The refactor is scoped toSNES_Stokes; those carry the same split and are untouched.Underworld development team with AI support from Claude Code