Automated goal-oriented adaptivity by pbrubeck · Pull Request #5273 · firedrakeproject/firedrake

pbrubeck · 2026-07-22T22:08:00Z

Description

Automates dual-weighted residual goal-oriented adaptive mesh refinement.

dwr_marking_callback(goal_functional) API.
Configuration through prefixed PETSc.Options() keys.
Enriched primal/dual solves via _SNESContext.solve_jacobian_transpose(), residual localization, and Dörfler marking.
Top-level solve(..., marking_callback=...) forwarding.
Demo from Rognes and Logg 2012

pbrubeck · 2026-07-31T16:42:00Z

+        f"dwr_{kind}_{key}": value
+        for kind in ("cell", "facet")
+        for key, value in {"ksp_type": "cg", "pc_type": "jacobi"}.items()


don't use dict comprehensions or f-strings, these need to be human readable.

pbrubeck · 2026-07-31T16:42:21Z

+    geo = SplineGeometry()
+    geo.AddRectangle((0, 0), (1, 1), bc="boundary")
+    mesh = Mesh(geo.GenerateMesh(maxh=0.5))


use UnitSquareMesh

pbrubeck · 2026-07-31T16:43:31Z

+@pytest.mark.parametrize(
+    ("adapt_option", "criterion"),
+    (("snes_adapt_sequence", "refine"),
+     ("snes_adapt_multigrid", "none")),


This requires https://gitlab.com/petsc/petsc/-/merge_requests/9447

Suggested change

("snes_adapt_multigrid", "none")),

# ("snes_adapt_multigrid", "none")),

pbrubeck · 2026-08-03T11:30:34Z

+    def solve_jacobian_transpose(self, rhs: Cofunction,
+                                 solution: Function) -> None:


Where should this function live?

NonlinearVariationalSolver: exposes it to the user, needs adjoint

_SNESContext: too obscure to be useful to a user

DWRMarkingCallback: prevents code reusability

DWRMarkingCallback estimates the dual-weighted-residual error indicator for a user-supplied goal functional and turns it into the DG0 marker Function that NonlinearVariationalSolver's existing marking-callback adaptive-refinement machinery expects, so goal-oriented mesh adaptivity is driven the same way as any other marking callback. Estimating the DWR indicator needs a dual (adjoint) solve on an enriched space, and a primal solve's Jacobian is symmetric only in specific cases, so _SNESContext grows solve_jacobian_transpose() to solve the dual problem against the transposed primal Jacobian. Finding "the" ksp behind solve_jacobian_transpose() and PMG's/adaptive-refinement's reconstructed _SNESContexts is unreliable via a weakref carried across every reconstruct() call (PETSc returns a fresh Python wrapper on every getKSP()/getDM(), and the weakref has to be explicitly re-attached each time); composing it directly on the DM instead (dm.setAttr("_ksp", ksp)) survives reconstruction for free since PETSc's attribute compose/query operates on the underlying PetscObject. _refine_adaptive() carries the composed ksp forward onto each newly-refined DM, alongside the parent/ctx-coarsener/appctx propagation it already does. NonlinearVariationalSolver.set_marking_callback() now recognises a DWRMarkingCallback and runs its setup() against the primal solution and options prefix; get_goal_functional() exposes the (possibly-adapted) goal functional for inspection after solve(). Also along the way: - add_hooks/SetupHooks now record the appctx a saved hook stack was built for and only replay it when the current appctx matches. Adaptive refinement replaces the appctx on every adapted solve, so replaying a stale hook stack anchored the DM chain on an already-torn-down root DM and surfaced as PETSc error 101 out of DMRefine() on a solver's second solve(). - NonlinearVariationalSolver.solve() no longer caches its PETSc work vector across calls (the problem may have been reconstructed onto an adapted mesh since the last solve) and incRef()s it whenever DMAdaptorAdapt() has swapped the solution DM out from under the solve, since DMAdaptorAdapt() steals a reference to it as the adapted-away DM's template global vector; without the extra ref a second collection of the same vector double-frees it and segfaults the interpreter. - _refine_adaptive() snaps to solution_mesh.unique() and, for a MeshSequenceGeometry, calls set_hierarchy() after adding the refined mesh, so goal-oriented adaptivity works on mesh sequences too.

Walks through DWRMarkingCallback on a Poisson problem: setting a goal functional, attaching the callback to a solver via set_marking_callback(), and inspecting the adapted mesh and goal value after solve().

Add a global DWR estimate of the error in the goal functional, split into the error committed by discretising and the error committed by not solving the algebraic system exactly, and stop adapting once it meets a tolerance. * -dwr_atol/-dwr_rtol stop the loop once |eta| < max(atol, rtol*|J(u_h)|). A marking callback that returns None now means "stop adapting": the refine hook hands the same DM back, and a SNES convergence test reports immediate convergence, so the rest of -snes_adapt_sequence costs nothing. PETSc's DMAdaptor has no tolerance of its own to do this with. * -dwr_monitor reports the estimate once per cycle. An exact_solution kwarg adds the true error and an effectivity index to that output. * Warn when the solver error estimate exceeds the discretisation error estimate, since refining then cannot help. * A marking callback is no longer required: without one, adaptive refinement degenerates to uniform refinement, as grid sequencing asks. Three fixes this uncovered: * dwr_ options were read from the reconstructed context, whose prefix is renamed after the multigrid level it becomes, so every one of them silently reverted to its default after the first refinement. * Function's second positional argument is val, not name, so the enriched-order primal solve raised on every mark. * Interpolation between a space and itself is the identity. PETSc asks for it when the adaptor hands back the DM it was given. Port the p-Laplacian goal-based adaptivity demo from #4893, with netgen replaced by a plain mesh, and give it and the elasticity demo names that say which is linear and which is not. The elasticity demo needs MUMPS null-pivot detection for its indefinite system. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>

pbrubeck · 2026-08-11T12:02:45Z

+  u.interpolate(0.99*u_exact)
+
+  v = TestFunction(V)
+  quadrature = {"degree": degree + 4}


Suggested change

quadrature = {"degree": degree + 4}

quadrature = {"degree": 2*degree + 2}

pbrubeck · 2026-08-11T12:04:12Z

+The DWR machinery localises the residual by projecting it onto cell-bubble and
+facet-bubble spaces. Those two auxiliary solves take their own options under
+the ``dwr_cell_`` and ``dwr_facet_`` prefixes. Both operators are well
+conditioned, so a diagonally preconditioned solve is enough. ::
+
+  solver_parameters.update({
+      "dwr_cell_ksp_type": "cg",
+      "dwr_cell_pc_type": "jacobi",
+      "dwr_facet_ksp_type": "cg",
+      "dwr_facet_pc_type": "jacobi",
+  })
+


These are the defaults, I don't see why one would ever want to change them

Suggested change

The DWR machinery localises the residual by projecting it onto cell-bubble and

facet-bubble spaces. Those two auxiliary solves take their own options under

the ``dwr_cell_`` and ``dwr_facet_`` prefixes. Both operators are well

conditioned, so a diagonally preconditioned solve is enough. ::

solver_parameters.update({

"dwr_cell_ksp_type": "cg",

"dwr_cell_pc_type": "jacobi",

"dwr_facet_ksp_type": "cg",

"dwr_facet_pc_type": "jacobi",

})

pbrubeck · 2026-08-11T12:05:25Z

+
+There is no loop to write here, unlike in the ad hoc implementations such a
+method usually needs. PETSc composes the solver, the estimator, the marker and
+the refiner, exactly as ``-snes_grid_sequence`` composes a solver with uniform


-snes_grid_sequence is not used in the repo, we should not refer to it.

pbrubeck · 2026-08-11T12:11:36Z

+    PETSc's `DMAdaptor` runs a fixed number of ``-snes_adapt_sequence`` steps
+    and has no error tolerance of its own. Once the marking callback declines
+    to mark anything, the mesh stops changing, so each remaining step would
+    re-solve a problem that is already solved. Reporting convergence
+    immediately makes those steps cost nothing.


we should add the relevant petsc changes to support early termination

_SNESContext.solve_jacobian_transpose found its KSP by reading a "_ksp" attribute stamped onto a DM. That was wrong twice over. The attribute was a strong reference, and KSPSetDM() makes the KSP hold the DM, so the pair formed a cycle that the garbage collector cannot break, with nothing to tear it down. A DM is also shared by every solver built on the same FunctionSpace, since it comes from the cached dof_dset. A second solver on that space overwrote the first solver's entry, and the first solver then solved against the second's operator with no error. Give the context a weak reference to its SNES instead, and resolve the KSP from it per call. A context is one-to-one with a solver, so nothing aliases; the reference is weak because the SNES owns the DM that owns the context during a solve. Contexts rebuilt by reconstruct() for field splits and coarse levels deliberately do not inherit it: the outer SNES's Jacobian describes a different problem. Rename the method to solve_jacobian(b, x, transpose=False) and report failure through a new check_ksp_convergence(), so a failed Jacobian solve reads like a failed solve(). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>

pefarrell · 2026-08-11T13:52:01Z

+
+Our goal is the weighted average shear traction on the right boundary.  Its
+exact value is approximately :math:`-0.06029761071`.  The DWR callback
+linearises this functional, solves the low- and enriched-order dual problems,


I was expecting this discussion of the DWR callback above a code block that uses the DWR callback?

pefarrell · 2026-08-11T13:52:30Z

+
+Our goal is the weighted average shear traction on the right boundary.  Its
+exact value is approximately :math:`-0.06029761071`.  The DWR callback
+linearises this functional, solves the low- and enriched-order dual problems,


Is it necessary to emphasise that the DWR callback 'linearises this functional'? Doesn't seem important to mention.

pefarrell · 2026-08-11T13:52:58Z

+  goal = psi*dot(dot(n, sigma), tangent)*ds(2)
+
+The Hellinger--Reissner system is indefinite.  We therefore ask MUMPS to detect
+null pivots (``icntl_24``) rather than give up on the first one it meets.  We


Does it work with the default icntl_24 setting?

pefarrell · 2026-08-11T13:53:36Z

+Each auxiliary solve the estimator performs has its own options prefix, so
+none of them is silently configured by the outer solver's options.
+``dwr_cell_`` and ``dwr_facet_`` localise the residual onto cells and facets;
+both mass-like operators are well conditioned, so a diagonally preconditioned


Aren't these solves block-diagonal (localised to each cell and facet)? That seems more important to mention than 'well conditioned'

They are diagonal. They only have one dof per cell/facet

pefarrell · 2026-08-11T13:54:09Z

+``snes_adapt_sequence`` bounds the number of SOLVE--ESTIMATE--MARK--REFINE
+cycles.  ``dwr_rtol`` stops the loop early once the estimated error in the goal
+falls below that fraction of :math:`|J(w_h)|`, and ``dwr_atol`` sets an absolute
+tolerance instead.  ``dwr_monitor`` reports the estimate once per cycle, split


Be explicit about what happens if both are set

pefarrell · 2026-08-11T13:54:48Z

+
+``snes_adapt_sequence`` bounds the number of SOLVE--ESTIMATE--MARK--REFINE
+cycles.  ``dwr_rtol`` stops the loop early once the estimated error in the goal
+falls below that fraction of :math:`|J(w_h)|`, and ``dwr_atol`` sets an absolute


Be explicit: it is below the specified fraction of the current uncorrected functional value with the not-enriched discretisation?

pbrubeck · 2026-08-12T07:38:37Z

+
+  solver_parameters.update({
+      "dwr_cell_ksp_type": "cg",
+      "dwr_cell_pc_type": "jacobi",
+      "dwr_facet_ksp_type": "cg",
+      "dwr_facet_pc_type": "jacobi",
+  })
+


These are the defaults

pbrubeck · 2026-08-12T08:28:30Z

+    subprocess.run([sys.executable, "-c", _COLLECT_AFTER_ADAPT], check=True)
+
+
+def _jacobian_solver(V, u):


move to different test file

An OptionsManager deletes every option under its prefix whose name matches one of its defaults, including the options that the database already held. The DWR callback rebuilds its enriched, cell and facet solvers on every adapted mesh, so only the first mesh saw the requested dwr_cell_, dwr_facet_ and dwr_enriched_ options; the rest silently used preonly and lu. The callback now captures the options under its prefix when a solver attaches it, passes them to those solvers as solver_parameters, and carries that copy through the reconstruction onto each adapted mesh. _adapt_marked_cells leaves dm_plex_transform_type out of its parameters when the user has already set it, so the same deletion cannot take that choice away either. Also rewrite the docstrings and comments of this branch to the AGENTS.md prose rules: active verbs, explicit relative pronouns, one idea per sentence. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

pbrubeck commented Jul 22, 2026

View reviewed changes

Comment thread firedrake/dwr.py Outdated

pbrubeck commented Jul 24, 2026

View reviewed changes

Comment thread firedrake/dmhooks.py

pbrubeck commented Jul 24, 2026

View reviewed changes

Comment thread firedrake/dmhooks.py Outdated

pbrubeck commented Jul 28, 2026

View reviewed changes

Comment thread demos/goal_oriented_adaptivity/goal_oriented_adaptivity.py.rst Outdated

Base automatically changed from pbrubeck/mg-redist to main July 31, 2026 10:37