From 674f83be910fda4ed1ac196561c79509a60a7f6b Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 22 Jul 2026 23:06:13 +0100 Subject: [PATCH 1/2] Automated goal-oriented adaptivity --- firedrake/__init__.py | 1 + firedrake/dmhooks.py | 3 + firedrake/dwr.py | 236 +++++++++++++++++++ firedrake/solving.py | 11 +- firedrake/solving_utils.py | 36 ++- firedrake/variational_solver.py | 7 + tests/firedrake/multigrid/test_snes_adapt.py | 54 +++++ 7 files changed, 343 insertions(+), 5 deletions(-) create mode 100644 firedrake/dwr.py diff --git a/firedrake/__init__.py b/firedrake/__init__.py index 29e9ea9664..369ad7bdd8 100644 --- a/firedrake/__init__.py +++ b/firedrake/__init__.py @@ -160,6 +160,7 @@ def init_petsc(): PointexprOperator, point_expr, MLOperator ) from firedrake.progress_bar import ProgressBar # noqa: F401 +from firedrake.dwr import dwr_marking_callback # noqa: F401 from firedrake.logging import ( # noqa: F401 set_level, set_log_handlers, set_log_level, DEBUG, INFO, diff --git a/firedrake/dmhooks.py b/firedrake/dmhooks.py index 0b82110457..91143f5039 100644 --- a/firedrake/dmhooks.py +++ b/firedrake/dmhooks.py @@ -492,6 +492,9 @@ def _refine_adaptive(dm): coefficient_mapping = {} refined_ctx = refine(ctx, refine, coefficient_mapping=coefficient_mapping) + from firedrake.dwr import DWRMarkingCallback + if isinstance(ctx._marking_callback, DWRMarkingCallback): + refined_ctx._marking_callback = ctx._marking_callback.reconstruct(coefficient_mapping) parent = get_parent(dm) coarsener = get_ctx_coarsener(dm) # Get all DMs from the refined problem diff --git a/firedrake/dwr.py b/firedrake/dwr.py new file mode 100644 index 0000000000..01a9400103 --- /dev/null +++ b/firedrake/dwr.py @@ -0,0 +1,236 @@ +import numpy as np +import ufl +from finat.ufl import BrokenElement, FiniteElement + +from firedrake.assemble import assemble +from firedrake.function import Function +from firedrake.functionspace import FunctionSpace, TensorFunctionSpace +from firedrake.petsc import PETSc +from firedrake.preconditioners.pmg import PMGPC +from firedrake.ufl_expr import TestFunction, TrialFunction, adjoint, derivative +from firedrake.variational_solver import (LinearVariationalProblem, + LinearVariationalSolver, + NonlinearVariationalProblem, + NonlinearVariationalSolver) +from ufl import avg, dS, ds, dx, inner, replace + + +__all__ = ("dwr_marking_callback",) + + +def _replace_arguments(form, *arguments): + return replace(form, dict(zip(form.arguments(), arguments))) + + +def _reconstruct_degree(V, degree): + element = PMGPC.reconstruct_degree(V.ufl_element(), degree) + return V.reconstruct(element=element) + + +def _homogeneous_bcs(bcs, V): + return [bc.reconstruct(V=V, indices=bc._indices, g=0) for bc in bcs] + + +def _both(expr): + return expr("+") + expr("-") + + +def _residual_indicators(F, dual_error, residual_degree, options_prefix): + """Compute the strong residual representation of Rognes and Logg.""" + v, = F.arguments() + V = v.function_space() + mesh = V.mesh().unique() + dim = mesh.topological_dimension + degree = V.ufl_element().degree() + residual_degree + variant = "integral" + + bubble_space = FunctionSpace(mesh, "B", dim + 1, variant=variant) + bubble = Function(bubble_space).assign(1) + if V.value_shape == (): + cell_space = FunctionSpace(mesh, "DG", degree, variant=variant) + else: + cell_space = TensorFunctionSpace(mesh, "DG", degree, + shape=V.value_shape, variant=variant) + cell_trial = TrialFunction(cell_space) + cell_test = TestFunction(cell_space) + cell_residual = Function(cell_space) + cell_problem = LinearVariationalProblem( + inner(cell_trial, bubble * cell_test) * dx, + _replace_arguments(F, bubble * cell_test), cell_residual, + ) + cell_solver = LinearVariationalSolver( + cell_problem, options_prefix=options_prefix + "dwr_cell_" + ) + cell_solver.solve() + + cone_space = FunctionSpace(mesh, "FB", dim, variant=variant) + cone = Function(cone_space).assign(1) + element = BrokenElement(FiniteElement("FB", cell=mesh.ufl_cell(), + degree=degree + dim, variant=variant)) + if V.value_shape == (): + facet_space = FunctionSpace(mesh, element) + else: + facet_space = TensorFunctionSpace(mesh, element, shape=V.value_shape) + facet_trial = TrialFunction(facet_space) + facet_test = TestFunction(facet_space) + facet_residual_hat = Function(facet_space) + facet_rhs = (_replace_arguments(F, facet_test) + - inner(cell_residual, facet_test) * dx) + facet_lhs = (_both(inner(facet_trial / cone, facet_test)) * dS + + inner(facet_trial / cone, facet_test) * ds) + facet_problem = LinearVariationalProblem( + facet_lhs, facet_rhs, facet_residual_hat + ) + facet_solver = LinearVariationalSolver( + facet_problem, options_prefix=options_prefix + "dwr_facet_" + ) + facet_solver.solve() + facet_residual = facet_residual_hat / cone + + indicator_space = FunctionSpace(mesh, "DG", 0) + indicator_test = TestFunction(indicator_space) + indicators = assemble( + inner(inner(cell_residual, dual_error), indicator_test) * dx + + inner(avg(inner(facet_residual, dual_error)), + _both(indicator_test)) * dS + + inner(inner(facet_residual, dual_error), indicator_test) * ds + ) + with indicators.dat.vec as vec: + vec.abs() + return indicators + + +def _dorfler_mark(indicators, fraction): + if not 0 < fraction <= 1: + raise ValueError("marking_fraction must lie in (0, 1]") + local = indicators.dat.data_ro.copy() + if not np.isfinite(local).all(): + raise RuntimeError("DWR error indicators contain non-finite values") + gathered = indicators.comm.allgather(local) + values = np.concatenate(gathered) + total = values.sum() + markers = Function(indicators.function_space()) + if total <= 0: + return markers.assign(1) + ordered = np.sort(values)[::-1] + count = np.searchsorted(np.cumsum(ordered), fraction * total) + 1 + threshold = ordered[min(count - 1, len(ordered) - 1)] + markers.dat.data_wo[:] = local >= threshold + return markers + + +class DWRMarkingCallback: + """Mark cells using an automatically localized dual-weighted residual.""" + + def __init__(self, goal_functional, primal=None, enrichment_degree=None): + if not isinstance(goal_functional, ufl.BaseForm) or goal_functional.arguments(): + raise ValueError("goal_functional must be a 0-form") + self.goal_functional = goal_functional + self._primal = primal + self._enrichment_degree = enrichment_degree + self._high_space = None + if primal is not None: + V = primal.function_space() + self._high_space = _reconstruct_degree( + V, V.ufl_element().degree() + enrichment_degree + ) + + def setup(self, primal, options_prefix): + options = PETSc.Options(options_prefix) + self._primal = primal + self._enrichment_degree = options.getInt("dwr_enrichment_degree", 1) + V = primal.function_space() + self._high_space = _reconstruct_degree( + V, V.ufl_element().degree() + self._enrichment_degree + ) + + def reconstruct(self, coefficient_mapping): + goal = replace(self.goal_functional, coefficient_mapping) + primal = coefficient_mapping[self._primal] + return type(self)(goal, primal, self._enrichment_degree) + + def __call__(self, ctx, current_solution): + return self._mark(ctx, current_solution) + + def _mark(self, ctx, current_solution): + problem = ctx._problem + V = current_solution.function_space() + prefix = ctx.options_prefix or "" + options = PETSc.Options(prefix) + residual_degree = options.getInt("dwr_residual_degree", 1) + marking_fraction = options.getReal("dwr_marking_fraction", 0.5) + high_space = self._high_space + if high_space is None: + raise RuntimeError("DWR marking callback has not been set up") + + dual_low = Function(V, name="dwr_dual_low") + direction = TestFunction(V) + goal_derivative = derivative(self.goal_functional, current_solution, direction) + rhs = assemble(goal_derivative, bcs=_homogeneous_bcs(problem.bcs, V)) + ctx.solve_jacobian_transpose(rhs, dual_low) + if not np.isfinite(dual_low.dat.data_ro).all(): + raise RuntimeError("DWR low-order dual solution contains non-finite values") + + primal_high = Function(high_space, name="dwr_primal_high") + primal_high.interpolate(current_solution) + test, = problem.F.arguments() + test_high = test.reconstruct(function_space=high_space) + F_high = replace(problem.F, {current_solution: primal_high, test: test_high}) + bcs_high = [bc.reconstruct(V=high_space, indices=bc._indices) + for bc in problem.bcs] + high_problem = NonlinearVariationalProblem(F_high, primal_high, bcs_high) + primal_solver = NonlinearVariationalSolver( + high_problem, + options_prefix=prefix + "dwr_primal_", + ) + primal_solver.solve() + if not np.isfinite(primal_high.dat.data_ro).all(): + raise RuntimeError("DWR enriched primal solution contains non-finite values") + + goal_high = replace(self.goal_functional, {current_solution: primal_high}) + dual_high = Function(high_space, name="dwr_dual_high") + direction_high = TrialFunction(high_space) + dual_test_high = TestFunction(high_space) + jacobian_high = derivative(F_high, primal_high, direction_high) + goal_derivative_high = derivative(goal_high, primal_high, dual_test_high) + dual_problem = LinearVariationalProblem( + adjoint(jacobian_high), goal_derivative_high, dual_high, + bcs=_homogeneous_bcs(problem.bcs, high_space), + ) + dual_solver = LinearVariationalSolver( + dual_problem, + options_prefix=prefix + "dwr_dual_", + ) + dual_solver.solve() + if not np.isfinite(dual_high.dat.data_ro).all(): + raise RuntimeError("DWR enriched dual solution contains non-finite values") + + dual_error = dual_high - dual_low + indicators = _residual_indicators( + problem.F, dual_error, residual_degree, prefix + ) + return _dorfler_mark(indicators, marking_fraction) + + +def dwr_marking_callback(goal_functional: ufl.BaseForm): + """Construct a dual-weighted residual cell-marking callback. + + Parameters + ---------- + goal_functional + A scalar UFL 0-form depending on the primal solution. + + Returns + ------- + DWRMarkingCallback + A callback suitable for ``solve(..., marking_callback=...)``. + + Notes + ----- + Options are read from the active solver's PETSc options prefix. The + supported options are ``dwr_enrichment_degree`` (default 1), + ``dwr_residual_degree`` (default 1), and ``dwr_marking_fraction`` + (default 0.5). The auxiliary solvers use the ``dwr_primal_``, + ``dwr_dual_``, ``dwr_cell_``, and ``dwr_facet_`` sub-prefixes. + """ + return DWRMarkingCallback(goal_functional) diff --git a/firedrake/solving.py b/firedrake/solving.py index 4e7e7d2b4f..2225839535 100644 --- a/firedrake/solving.py +++ b/firedrake/solving.py @@ -161,7 +161,7 @@ def _solve_varproblem(*args, **kwargs): eq, u, bcs, J, Jp, M, form_compiler_parameters, \ solver_parameters, nullspace, nullspace_T, \ near_nullspace, \ - options_prefix, restrict, pre_apply_bcs = _extract_args(*args, **kwargs) + options_prefix, restrict, pre_apply_bcs, marking_callback = _extract_args(*args, **kwargs) # Check whether solution is valid if not isinstance(u, Function): @@ -197,7 +197,8 @@ def _solve_varproblem(*args, **kwargs): near_nullspace=near_nullspace, options_prefix=options_prefix, appctx=appctx, - pre_apply_bcs=pre_apply_bcs) + pre_apply_bcs=pre_apply_bcs, + marking_callback=marking_callback) return solver.solve() @@ -291,7 +292,8 @@ def _extract_args(*args, **kwargs): valid_kwargs = ["bcs", "J", "Jp", "M", "form_compiler_parameters", "solver_parameters", "nullspace", "transpose_nullspace", "near_nullspace", - "options_prefix", "appctx", "restrict", "pre_apply_bcs"] + "options_prefix", "appctx", "restrict", "pre_apply_bcs", + "marking_callback"] for kwarg in kwargs.keys(): if kwarg not in valid_kwargs: raise RuntimeError("Illegal keyword argument '%s'; valid keywords \ @@ -336,10 +338,11 @@ def _extract_args(*args, **kwargs): options_prefix = kwargs.get("options_prefix", None) restrict = kwargs.get("restrict", False) pre_apply_bcs = kwargs.get("pre_apply_bcs", True) + marking_callback = kwargs.get("marking_callback", None) return eq, u, bcs, J, Jp, M, form_compiler_parameters, \ solver_parameters, nullspace, nullspace_T, near_nullspace, \ - options_prefix, restrict, pre_apply_bcs + options_prefix, restrict, pre_apply_bcs, marking_callback def _extract_bcs(bcs): diff --git a/firedrake/solving_utils.py b/firedrake/solving_utils.py index 3306e4413a..6896a0d2b5 100644 --- a/firedrake/solving_utils.py +++ b/firedrake/solving_utils.py @@ -218,6 +218,7 @@ def __init__(self, problem, self._post_jacobian_callback = post_jacobian_callback self._post_function_callback = post_function_callback self._marking_callback = marking_callback + self._snes_ref = None self.fcp = problem.form_compiler_parameters # Function to hold current guess @@ -302,7 +303,40 @@ def reconstruct(self, problem=None, mat_type=None, pmat_type=None, **kwargs): for k, v in default_options.items(): if kwargs.get(k) is None: kwargs[k] = v - return _SNESContext(problem, mat_type, pmat_type, **kwargs) + ctx = _SNESContext(problem, mat_type, pmat_type, **kwargs) + ctx._snes_ref = self._snes_ref + return ctx + + def solve_jacobian_transpose(self, rhs, solution): + """Solve the transpose of the current SNES Jacobian. + + Parameters + ---------- + rhs + The dual right-hand side. + solution + The Function in which to store the solution. + """ + if self._snes_ref is None: + raise RuntimeError("The SNES context is not attached to a solver") + snes = self._snes_ref() + if snes is None: + raise RuntimeError("The SNES attached to this context has been destroyed") + with rhs.dat.vec_ro as b, solution.dat.vec_wo as x: + if b.norm() == 0: + x.set(0) + return + ksp = snes.getKSP() + ksp.solveTranspose(b, x) + reason = ksp.getConvergedReason() + if reason < 0: + raise ConvergenceError( + "Transpose Jacobian solve failed to converge after " + f"{ksp.getIterationNumber()} iterations.\nReason:\n " + f"{KSPReasons.get(reason, reason)} " + f"(PC type {ksp.getPC().getType()}, " + f"failure {ksp.getPC().getFailedReason()})" + ) @property def transfer_manager(self): diff --git a/firedrake/variational_solver.py b/firedrake/variational_solver.py index c75dea686f..bb4d4a1d2a 100644 --- a/firedrake/variational_solver.py +++ b/firedrake/variational_solver.py @@ -1,6 +1,7 @@ from __future__ import annotations import ufl +import weakref from itertools import chain from contextlib import ExitStack from types import MappingProxyType @@ -312,6 +313,7 @@ def update_diffusivity(current_solution): self.snes = PETSc.SNES().create(comm=problem.dm.comm) self._ctx = ctx + self._ctx._snes_ref = weakref.ref(self.snes) self._work = problem.u_restrict.dof_dset.layout_vec.duplicate() self.snes.setDM(problem.dm) if marking_callback is not None: @@ -357,6 +359,10 @@ def set_marking_callback(self, callback): """ if not callable(callback): raise TypeError(f"marking callback must be callable, not a {type(callback).__name__}") + from firedrake.dwr import DWRMarkingCallback + if isinstance(callback, DWRMarkingCallback): + with self.inserted_options(): + callback.setup(self._ctx._problem.u, self.options_prefix) self._ctx._marking_callback = callback def get_solution(self): @@ -444,6 +450,7 @@ def solve(self, bounds=None): self.snes.solve(None, work) # The appctx might have been refined self._ctx = dmhooks.get_appctx(self.snes.getDM()) + self._ctx._snes_ref = weakref.ref(self.snes) problem = self._ctx._problem solution = self.snes.getSolution() with problem.u_restrict.dat.vec as u: diff --git a/tests/firedrake/multigrid/test_snes_adapt.py b/tests/firedrake/multigrid/test_snes_adapt.py index 91a04345c6..dd8453eadb 100644 --- a/tests/firedrake/multigrid/test_snes_adapt.py +++ b/tests/firedrake/multigrid/test_snes_adapt.py @@ -20,6 +20,60 @@ def mark_cells(ctx, current_solution): assert solver._ctx._marking_callback is mark_cells +def test_solve_accepts_marking_callback(): + def mark_cells(ctx, current_solution): + M = FunctionSpace(current_solution.mesh(), "DG", 0) + return Function(M).assign(1) + + mesh = UnitSquareMesh(1, 1) + V = FunctionSpace(mesh, "CG", 1) + u = Function(V) + v = TestFunction(V) + + result = solve((u - 1.0)*v*dx == 0, u, marking_callback=mark_cells) + + assert result is u + + +@pytest.mark.skipnetgen +def test_dwr_marking_callback_builds_poisson_markers(): + from netgen.geom2d import SplineGeometry + + geo = SplineGeometry() + geo.AddRectangle((0, 0), (1, 1), bc="boundary") + mesh = Mesh(geo.GenerateMesh(maxh=0.5)) + V = FunctionSpace(mesh, "CG", 1) + u = Function(V) + v = TestFunction(V) + F = inner(grad(u), grad(v))*dx - v*dx + goal = u*dx + direct = {"ksp_type": "preonly", "pc_type": "lu"} + callback = dwr_marking_callback(goal) + dwr_direct = { + f"dwr_{kind}_{key}": value + for kind in ("primal", "dual") + for key, value in direct.items() + } + dwr_local = { + f"dwr_{kind}_{key}": value + for kind in ("cell", "facet") + for key, value in {"ksp_type": "cg", "pc_type": "jacobi"}.items() + } + + problem = NonlinearVariationalProblem( + F, u, bcs=DirichletBC(V, 0, "on_boundary") + ) + solver = NonlinearVariationalSolver( + problem, + solver_parameters={**direct, **dwr_direct, **dwr_local}, + marking_callback=callback, + ) + result = solver.solve() + markers = callback(solver._ctx, result) + + assert markers.function_space().mesh() is mesh + + @pytest.mark.skipnetgen def test_marking_callback_refine_hook_reconstructs_problem(): from netgen.geom2d import SplineGeometry From fb1f77c5244e36ad56a4891b6c6f2340633c103c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 23 Jul 2026 01:42:40 +0100 Subject: [PATCH 2/2] Demo --- .../goal_oriented_adaptivity.py.rst | 136 ++++++++++++++++++ firedrake/dmhooks.py | 9 +- firedrake/dwr.py | 36 +++-- firedrake/solving_utils.py | 3 +- tests/firedrake/demos/test_demos_run.py | 1 + tests/firedrake/multigrid/test_snes_adapt.py | 28 +++- 6 files changed, 189 insertions(+), 24 deletions(-) create mode 100644 demos/goal_oriented_adaptivity/goal_oriented_adaptivity.py.rst diff --git a/demos/goal_oriented_adaptivity/goal_oriented_adaptivity.py.rst b/demos/goal_oriented_adaptivity/goal_oriented_adaptivity.py.rst new file mode 100644 index 0000000000..b2f0ab7c0a --- /dev/null +++ b/demos/goal_oriented_adaptivity/goal_oriented_adaptivity.py.rst @@ -0,0 +1,136 @@ +Goal-oriented mesh adaptivity +============================= + +This demo uses a dual-weighted residual (DWR) estimator to adapt a mesh for +one particular quantity of interest. The example is the weakly symmetric +Hellinger--Reissner elasticity problem studied by Rognes and Logg in +`Automated Goal-Oriented Error Control I +`__. It is a useful test because its mixed +space contains both :math:`H(\mathrm{div})`- and :math:`L^2`-conforming fields. + +The unknowns are the two rows of the stress :math:`\sigma`, the displacement +:math:`u`, and a scalar rotation :math:`\gamma`. We use first-order BDM +elements for each stress row, piecewise constants for the displacement, and +continuous linears for the rotation. The deliberately small initial mesh +keeps the demo quick; increasing its resolution and the number of refinement +steps produces the longer adaptive sequences used in the paper. :: + + from firedrake import * + + mesh = UnitSquareMesh(2, 2) + stress_row = FunctionSpace(mesh, "BDM", 1) + displacement = VectorFunctionSpace(mesh, "DG", 0) + rotation = FunctionSpace(mesh, "CG", 1) + W = stress_row * stress_row * displacement * rotation + + w = Function(W, name="elasticity solution") + sigma0, sigma1, u, gamma = split(w) + tau0, tau1, v, eta = TestFunctions(W) + sigma = as_tensor((sigma0, sigma1)) + tau = as_tensor((tau0, tau1)) + +For shear modulus :math:`\mu` and Lamé parameter :math:`\lambda`, the +compliance tensor is + +.. math:: + + A\sigma = \frac{1}{2\mu}\left(\sigma + - \frac{\lambda}{2(\mu + \lambda)}\operatorname{tr}(\sigma)I\right). + +We manufacture the displacement +:math:`u_0=(xy\sin(\pi y),0)` and use the corresponding body force from the +paper. Prescribed displacement is a natural boundary condition in this +mixed formulation. :: + + x, y = SpatialCoordinate(mesh) + mu = Constant(1.0) + lmbda = Constant(100.0) + compliance = ( + sigma - lmbda/(2*(mu + lmbda))*tr(sigma)*Identity(2) + )/(2*mu) + + body_force = as_vector(( + pi*mu*(2*x*cos(pi*y) - pi*x*y*sin(pi*y)), + (mu + lmbda)*(pi*y*cos(pi*y) + sin(pi*y)), + )) + u0 = as_vector((x*y*sin(pi*y), 0)) + n = FacetNormal(mesh) + + F = ( + inner(compliance, tau) + + dot(div(sigma), v) + + dot(u, div(tau)) + + (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 + +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, +localises the weak residual with bubble and cone functions, and performs +global Dörfler marking. :: + + psi = y*(y - 1) + tangent = as_vector((0, 1)) + goal = psi*dot(dot(n, sigma), tangent)*ds(2) + +The outer solve and all four auxiliary solver families are configured through +PETSc options. Their option prefixes are ``dwr_primal_``, ``dwr_dual_``, +``dwr_cell_``, and ``dwr_facet_``. Here we use direct solves throughout so +that the example focuses on adaptivity. :: + + direct = { + "ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps", + } + solver_parameters = { + **direct, + "snes_adapt_sequence": 1, + "dwr_marking_fraction": 0.5, + **{ + f"dwr_{kind}_{key}": value + for kind in ("primal", "dual") + for key, value in direct.items() + }, + **{ + f"dwr_{kind}_{key}": value + for kind in ("cell", "facet") + for key, value in { + "ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps", + }.items() + }, + } + + initial_dofs = W.dim() + w_adapt = solve( + F == 0, + w, + solver_parameters=solver_parameters, + marking_callback=dwr_marking_callback(goal), + ) + +``solve`` returns the solution on the final adapted mesh. Since the goal's +measure and spatial coordinate belong to the original mesh, we write the same +functional on the final mesh when reporting its value. :: + + adapted_mesh = w_adapt.function_space().mesh().unique() + adapted_sigma0, adapted_sigma1, _, _ = split(w_adapt) + adapted_sigma = as_tensor((adapted_sigma0, adapted_sigma1)) + _, adapted_y = SpatialCoordinate(adapted_mesh) + adapted_normal = FacetNormal(adapted_mesh) + adapted_goal = ( + adapted_y*(adapted_y - 1) + * dot(dot(adapted_normal, adapted_sigma), tangent) + * ds(2, domain=adapted_mesh) + ) + + print(f"degrees of freedom: {initial_dofs} -> {w_adapt.function_space().dim()}") + print(f"weighted shear traction: {assemble(adapted_goal):.8f}") + +The refinement is driven by the error in this boundary traction, rather than +by a global energy norm. More refinement steps can be requested by increasing +``snes_adapt_sequence`` without writing an explicit solve--estimate--refine +loop. diff --git a/firedrake/dmhooks.py b/firedrake/dmhooks.py index 91143f5039..d98dc7485b 100644 --- a/firedrake/dmhooks.py +++ b/firedrake/dmhooks.py @@ -457,15 +457,12 @@ def _refine_adaptive(dm): from firedrake.mg.ufl_utils import refine from firedrake.mg.utils import get_level - # DMAdaptorAdapt() unconditionally destroys its input DM, and each - # adapted input DM remains a level in the AdaptiveMeshHierarchy. - dm.incRef() - ctx = get_appctx(dm) if ctx is None: raise RuntimeError("No _SNESContext found on DM") current_solution = ctx._x - mesh = current_solution.function_space().mesh() + solution_mesh = current_solution.function_space().mesh() + mesh = solution_mesh.unique() hierarchy, level = get_level(mesh) if hierarchy is None: hierarchy = AdaptiveMeshHierarchy(mesh) @@ -489,6 +486,8 @@ def _refine_adaptive(dm): raise ValueError("marking callback must return a DG0 Function or Cofunction") hierarchy.add_mesh(mesh.refine_marked_elements(markers)) + if isinstance(solution_mesh, MeshSequenceGeometry): + solution_mesh.set_hierarchy() coefficient_mapping = {} refined_ctx = refine(ctx, refine, coefficient_mapping=coefficient_mapping) diff --git a/firedrake/dwr.py b/firedrake/dwr.py index 01a9400103..13d7fc7df5 100644 --- a/firedrake/dwr.py +++ b/firedrake/dwr.py @@ -1,8 +1,11 @@ +from __future__ import annotations + import numpy as np import ufl from finat.ufl import BrokenElement, FiniteElement from firedrake.assemble import assemble +from firedrake.exceptions import ConvergenceError from firedrake.function import Function from firedrake.functionspace import FunctionSpace, TensorFunctionSpace from firedrake.petsc import PETSc @@ -35,6 +38,12 @@ def _both(expr): return expr("+") + expr("-") +def _check_finite(function: Function, name: str) -> None: + with function.dat.vec_ro as vec: + if not np.isfinite(vec.norm()): + raise ConvergenceError(f"{name} contains non-finite values") + + def _residual_indicators(F, dual_error, residual_degree, options_prefix): """Compute the strong residual representation of Rognes and Logg.""" v, = F.arguments() @@ -100,12 +109,12 @@ def _residual_indicators(F, dual_error, residual_degree, options_prefix): return indicators -def _dorfler_mark(indicators, fraction): +def _dorfler_mark(indicators: Function, fraction: float) -> Function: if not 0 < fraction <= 1: raise ValueError("marking_fraction must lie in (0, 1]") local = indicators.dat.data_ro.copy() if not np.isfinite(local).all(): - raise RuntimeError("DWR error indicators contain non-finite values") + raise ConvergenceError("DWR error indicators contain non-finite values") gathered = indicators.comm.allgather(local) values = np.concatenate(gathered) total = values.sum() @@ -122,7 +131,9 @@ def _dorfler_mark(indicators, fraction): class DWRMarkingCallback: """Mark cells using an automatically localized dual-weighted residual.""" - def __init__(self, goal_functional, primal=None, enrichment_degree=None): + def __init__(self, goal_functional: ufl.BaseForm, + primal: Function | None = None, + enrichment_degree: int | None = None): if not isinstance(goal_functional, ufl.BaseForm) or goal_functional.arguments(): raise ValueError("goal_functional must be a 0-form") self.goal_functional = goal_functional @@ -135,7 +146,7 @@ def __init__(self, goal_functional, primal=None, enrichment_degree=None): V, V.ufl_element().degree() + enrichment_degree ) - def setup(self, primal, options_prefix): + def setup(self, primal: Function, options_prefix: str) -> None: options = PETSc.Options(options_prefix) self._primal = primal self._enrichment_degree = options.getInt("dwr_enrichment_degree", 1) @@ -144,15 +155,15 @@ def setup(self, primal, options_prefix): V, V.ufl_element().degree() + self._enrichment_degree ) - def reconstruct(self, coefficient_mapping): + def reconstruct(self, coefficient_mapping: dict) -> DWRMarkingCallback: goal = replace(self.goal_functional, coefficient_mapping) primal = coefficient_mapping[self._primal] return type(self)(goal, primal, self._enrichment_degree) - def __call__(self, ctx, current_solution): + def __call__(self, ctx, current_solution: Function) -> Function: return self._mark(ctx, current_solution) - def _mark(self, ctx, current_solution): + def _mark(self, ctx, current_solution: Function) -> Function: problem = ctx._problem V = current_solution.function_space() prefix = ctx.options_prefix or "" @@ -168,8 +179,7 @@ def _mark(self, ctx, current_solution): goal_derivative = derivative(self.goal_functional, current_solution, direction) rhs = assemble(goal_derivative, bcs=_homogeneous_bcs(problem.bcs, V)) ctx.solve_jacobian_transpose(rhs, dual_low) - if not np.isfinite(dual_low.dat.data_ro).all(): - raise RuntimeError("DWR low-order dual solution contains non-finite values") + _check_finite(dual_low, "DWR low-order dual solution") primal_high = Function(high_space, name="dwr_primal_high") primal_high.interpolate(current_solution) @@ -184,8 +194,7 @@ def _mark(self, ctx, current_solution): options_prefix=prefix + "dwr_primal_", ) primal_solver.solve() - if not np.isfinite(primal_high.dat.data_ro).all(): - raise RuntimeError("DWR enriched primal solution contains non-finite values") + _check_finite(primal_high, "DWR enriched primal solution") goal_high = replace(self.goal_functional, {current_solution: primal_high}) dual_high = Function(high_space, name="dwr_dual_high") @@ -202,8 +211,7 @@ def _mark(self, ctx, current_solution): options_prefix=prefix + "dwr_dual_", ) dual_solver.solve() - if not np.isfinite(dual_high.dat.data_ro).all(): - raise RuntimeError("DWR enriched dual solution contains non-finite values") + _check_finite(dual_high, "DWR enriched dual solution") dual_error = dual_high - dual_low indicators = _residual_indicators( @@ -212,7 +220,7 @@ def _mark(self, ctx, current_solution): return _dorfler_mark(indicators, marking_fraction) -def dwr_marking_callback(goal_functional: ufl.BaseForm): +def dwr_marking_callback(goal_functional: ufl.BaseForm) -> DWRMarkingCallback: """Construct a dual-weighted residual cell-marking callback. Parameters diff --git a/firedrake/solving_utils.py b/firedrake/solving_utils.py index 6896a0d2b5..93736c4f12 100644 --- a/firedrake/solving_utils.py +++ b/firedrake/solving_utils.py @@ -307,7 +307,8 @@ def reconstruct(self, problem=None, mat_type=None, pmat_type=None, **kwargs): ctx._snes_ref = self._snes_ref return ctx - def solve_jacobian_transpose(self, rhs, solution): + def solve_jacobian_transpose(self, rhs: Cofunction, + solution: Function) -> None: """Solve the transpose of the current SNES Jacobian. Parameters diff --git a/tests/firedrake/demos/test_demos_run.py b/tests/firedrake/demos/test_demos_run.py index f53a611c2b..69b5ad86ed 100644 --- a/tests/firedrake/demos/test_demos_run.py +++ b/tests/firedrake/demos/test_demos_run.py @@ -51,6 +51,7 @@ Demo(("quasigeostrophy_1layer", "qg_1layer_wave"), ["hypre", "vtk"]), Demo(("saddle_point_pc", "saddle_point_systems"), ["hypre", "mumps"]), Demo(("fast_diagonalisation", "fast_diagonalisation_poisson"), ["mumps"]), + Demo(("goal_oriented_adaptivity", "goal_oriented_adaptivity"), ["mumps"]), Demo(('vlasov_poisson_1d', 'vp1d'), []), Demo(('shape_optimization', 'shape_optimization'), ["adjoint", "vtk"]), Demo(('submesh_reaction_diffusion', 'submesh_reaction_diffusion'), ["netgen", "vtk"]), diff --git a/tests/firedrake/multigrid/test_snes_adapt.py b/tests/firedrake/multigrid/test_snes_adapt.py index dd8453eadb..535b82dce9 100644 --- a/tests/firedrake/multigrid/test_snes_adapt.py +++ b/tests/firedrake/multigrid/test_snes_adapt.py @@ -1,3 +1,5 @@ +import weakref + import pytest from firedrake import * from firedrake import dmhooks @@ -18,6 +20,8 @@ def mark_cells(ctx, current_solution): assert solver.parameters["adaptor_criterion"] == "refine" assert solver._ctx._marking_callback is mark_cells + assert isinstance(solver._ctx._snes_ref, weakref.ReferenceType) + assert solver._ctx._snes_ref() is solver.snes def test_solve_accepts_marking_callback(): @@ -36,13 +40,20 @@ def mark_cells(ctx, current_solution): @pytest.mark.skipnetgen -def test_dwr_marking_callback_builds_poisson_markers(): +@pytest.mark.parallel([1, 2]) +@pytest.mark.parametrize( + ("adapt_option", "criterion"), + (("snes_adapt_sequence", "refine"), + ("snes_adapt_multigrid", "none")), +) +def test_dwr_marking_callback_builds_poisson_markers(adapt_option, criterion): from netgen.geom2d import SplineGeometry geo = SplineGeometry() geo.AddRectangle((0, 0), (1, 1), bc="boundary") mesh = Mesh(geo.GenerateMesh(maxh=0.5)) V = FunctionSpace(mesh, "CG", 1) + old_dim = V.dim() u = Function(V) v = TestFunction(V) F = inner(grad(u), grad(v))*dx - v*dx @@ -65,13 +76,22 @@ def test_dwr_marking_callback_builds_poisson_markers(): ) solver = NonlinearVariationalSolver( problem, - solver_parameters={**direct, **dwr_direct, **dwr_local}, + solver_parameters={ + **direct, + **dwr_direct, + **dwr_local, + adapt_option: 1, + "adaptor_criterion": criterion, + }, marking_callback=callback, ) result = solver.solve() - markers = callback(solver._ctx, result) - assert markers.function_space().mesh() is mesh + assert result.function_space().mesh() is not mesh + assert result.function_space().dim() > old_dim + hierarchy, level = get_level(result.function_space().mesh()) + assert level == 1 + assert hierarchy[0] is mesh @pytest.mark.skipnetgen