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

pbrubeck · 2026-08-03T22:44:56Z

+            if b.norm() == 0:
+                x.set(0)
+                return


Suggested change

if b.norm() == 0:

x.set(0)

return

pbrubeck · 2026-08-03T22:50:18Z

+                if self.snes.getDM() != solution_dm:
+                    # DMAdaptorAdapt() consumed a reference to work, which the
+                    # adapted-away DM keeps as its template global vector.
+                    work.incRef()


isn't better to destroy the work Vec?

pbrubeck · 2026-08-03T22:51:41Z

+    second_goal = assemble(solver.get_goal_functional())
+
+    assert first_dim > old_dim
+    assert isinstance(first_goal, float)


pointless test

pbrubeck · 2026-08-05T22:23:54Z

+      + dot(div(sigma), v)
+      + dot(u, div(tau))


Suggested change

+ dot(div(sigma), v)

+ dot(u, div(tau))

+ inner(div(sigma), v)

+ inner(u, div(tau))

pbrubeck · 2026-08-05T22:24:47Z

+      + (sigma[0, 1] - sigma[1, 0])*eta
+      + gamma*(tau[0, 1] - tau[1, 0])
+  )*dx - dot(body_force, v)*dx - dot(u0, dot(tau, n))*ds


Suggested change

+ (sigma[0, 1] - sigma[1, 0])*eta

+ gamma*(tau[0, 1] - tau[1, 0])

)*dx - dot(body_force, v)*dx - dot(u0, dot(tau, n))*ds

+ inner(sigma[0, 1] - sigma[1, 0], eta)

+ inner(gamma, tau[0, 1] - tau[1, 0])

)*dx - inner(body_force, v)*dx - inner(u0, dot(tau, n))*ds

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().

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