Demos for goal-oriented adaptivity

demos/demo_references.bib

@@ -601,3 +601,35 @@ @article{Hu:2014
   doi={10.1007/s10915-014-9821-5}
 }
 
+@article{Becker:2001,
+  doi = {10.1017/S0962492901000010},
+  year = 2001,
+  volume = {10},
+  pages = {1--102},
+  author = {Roland Becker and Rolf Rannacher},
+  title = {An optimal control approach to a posteriori error estimation in finite element methods},
+  journal = {Acta Numerica}
+}
+
+@article{Rognes:2013,
+  doi = {10.1137/10081962X},
+  year = 2013,
+  volume = {35},
+  number = {3},
+  pages = {C173--C193},
+  author = {Marie E. Rognes and Anders Logg},
+  title = {Automated goal-oriented error control {I}: stationary variational problems},
+  journal = {{SIAM} Journal on Scientific Computing}
+}
+
+@incollection{Endtmayer:2024,
+  doi = {10.1016/bs.aams.2024.08.003},
+  year = 2024,
+  volume = {59},
+  pages = {19--108},
+  author = {Bernhard Endtmayer and Ulrich Langer and Thomas Richter and Andreas Schafelner and Thomas Wick},
+  title = {A posteriori single- and multi-goal error control and adaptivity for partial differential equations},
+  booktitle = {Error Control, Adaptive Discretizations, and Applications, Part 2},
+  series = {Advances in Applied Mechanics},
+  publisher = {Elsevier}
+}

demos/goal_oriented_nonlinear_p_laplacian/goal_oriented_nonlinear_p_laplacian.py.rst

@@ -0,0 +1,196 @@
+Goal-based adaptivity for stationary boundary value problems
+============================================================
+
+.. rst-class:: emphasis
+
+    The dual-weighted residual (DWR) method designs global and local error
+    estimators for a goal functional :math:`J(u)`, where :math:`u` solves a
+    partial differential equation. Deriving DWR for a
+    specific problem usually takes substantial expertise: the appropriate
+    adjoint equation, the residual formulation, and so on. This demo shows
+    how Firedrake implements DWR automatically for stationary boundary value
+    problems.
+
+    The demo was contributed by `Patrick Farrell
+    <mailto:patrick.farrell@maths.ox.ac.uk>`__, based on the MSc project of
+    `Joseph Flood <mailto:josephdflood01@gmail.com>`__.
+
+The dual-weighted residual method :cite:`Becker:2001` designs global and local
+error estimators for the error in a goal functional :math:`J(u)`. While
+implementing DWR by hand involves substantial expertise, the high-level
+symbolic UFL representation of the problem permits the *automation* of DWR
+:cite:`Rognes:2013`.
+
+Here we apply DWR to a nonlinear stationary boundary value problem, the
+:math:`p`-Laplacian:
+
+.. math::
+    -\nabla \cdot \left( |\nabla u|^{p-2} \nabla u \right) = f \text{ in } \Omega, \quad u = 0 \text{ on } \partial \Omega.
+
+We solve it on the unit square with a known analytical solution, so that we can
+compute effectivity indices for our error estimates. Any Firedrake mesh can be
+refined in place, so we start from a uniform mesh:
+
+.. code-block:: python
+
+  from firedrake import *
+
+  mesh = UnitSquareMesh(8, 8)
+
+  degree = 2
+  V = FunctionSpace(mesh, "CG", degree)
+  x, y = SpatialCoordinate(mesh)
+
+  p = Constant(5)
+  u_exact = x*(1-x)*y*(1-y)*exp(2*pi*x)*cos(pi*y)
+  f = -div(inner(grad(u_exact), grad(u_exact))**((p-2)/2) * grad(u_exact))
+
+The manufactured solution vanishes on all four edges, so the boundary condition
+is homogeneous. The problem is strongly nonlinear, so for the purposes of this
+demo we cheat and start from an initial guess close to the exact solution. The
+integrand is not polynomial, so we ask for a richer quadrature rule than the
+degree of the space alone would suggest:
+
+.. code-block:: python
+
+  u = Function(V, name="Solution")
+  u.interpolate(0.99*u_exact)
+
+  v = TestFunction(V)
+  quadrature = {"degree": 2*degree + 2}
+  F = (inner(inner(grad(u), grad(u))**((p-2)/2) * grad(u), grad(v))*dx(**quadrature)
+       - inner(f, v)*dx(**quadrature))
+  bcs = DirichletBC(V, 0, "on_boundary")
+
+To apply goal-based adaptivity we need a goal functional. Here we use the
+integral of the normal derivative of the solution over the top boundary, which
+`UnitSquareMesh` numbers ``4``:
+
+.. code-block:: python
+
+  n = FacetNormal(mesh)
+  J = inner(grad(u), n)*ds(4)
+
+The estimate splits in two: the error committed by discretising, and the error
+committed by not solving the algebraic system exactly. Because DWR measures
+that algebraic error too, a very coarse tolerance on the nonlinear solver
+costs us nothing:
+
+.. code-block:: python
+
+  solver_parameters = {
+      "snes_atol": 1.0e-6,
+      "snes_rtol": 1.0e-1,
+      "snes_linesearch_type": "l2",
+      "ksp_type": "preonly",
+      "pc_type": "lu",
+  }
+
+The DWR machinery localises the residual by projecting it onto cell-bubble and
+facet-bubble spaces.  The resulting local problems are diagonal because each
+has one degree of freedom per cell or facet.
+
+Now the adaptive loop itself. ``snes_adapt_sequence`` is the *maximum* number
+of SOLVE--ESTIMATE--MARK--REFINE cycles. ``goal_atol`` stops it early once the
+estimated error in the goal is small enough, and ``goal_rtol`` gives a tolerance
+relative to the absolute value of the current goal functional on the
+uncorrected, non-enriched solution.  If both tolerances are set, the loop stops
+when the estimate falls below the larger threshold.  ``goal_marking_fraction``
+is the fraction of the total estimated error that Dörfler marking must
+capture, and ``goal_monitor`` reports the estimate at each cycle:
+
+.. code-block:: python
+
+  solver_parameters.update({
+      "adaptor_criterion": "refine",
+      "snes_adapt_sequence": 10,
+      "goal_atol": 5.0e-3,
+      "goal_marking_fraction": 0.5,
+      "goal_monitor": None,
+  })
+
+The adaptive loop is handled inside the solver.  It runs the configured
+SOLVE--ESTIMATE--MARK--REFINE cycle, and the marking callback supplies the
+cells to refine.  We also pass the exact solution, which lets the monitor
+report the true error and an effectivity index.  Of course, this is not needed
+in general:
+
+.. code-block:: python
+
+  initial_dofs = V.dim()
+  problem = NonlinearVariationalProblem(F, u, bcs)
+  solver = NonlinearVariationalSolver(
+      problem,
+      solver_parameters=solver_parameters,
+      marking_callback=GoalOrientedMarker(J, exact_solution=u_exact),
+  )
+  uh = solver.solve()
+
+``solver.solve()`` returns the solution on the final adapted mesh, in the
+original degree-2 space.  The enriched solution that the estimator computes is
+not returned.  ``solver.get_marking_callback().error_estimate`` is the
+estimate :math:`\eta` of the error in the goal for the returned solution, and
+it is the value that stopped the loop:
+
+.. code-block:: python
+
+  print(f"degrees of freedom: {initial_dofs} -> {uh.function_space().dim()}")
+  print(f"error estimate: {solver.get_marking_callback().error_estimate:.6e}")
+
+The loop terminates after eight refinements, short of the ten it was allowed,
+because the estimate fell below ``goal_atol``. The effectivity indices
+
+.. math::
+
+    I = \frac{\eta}{J(u) - J(u_h)}
+
+printed by ``goal_monitor`` stay close to one throughout, between roughly
+:math:`0.83` and :math:`1.30`. The estimate is therefore a faithful measure of
+the error it steers by. Tightening ``goal_atol`` refines further and drives the
+estimate down accordingly.
+
+Because the two parts of the estimate are reported separately, it is easy to see
+which one is limiting. On the last cycle this run prints
+
+.. code-block:: text
+
+    DWR: the solver error exceeds the discretisation error. Tighten the solver tolerances.
+
+By then the mesh is fine enough that the discretisation error falls below the
+algebraic error left by a ``snes_rtol`` of :math:`10^{-1}`. Further refinement
+is of no use: the algebraic error limits the accuracy of the goal, not the
+mesh. Tightening ``snes_rtol`` and ``snes_atol`` is what helps. Here the
+estimate already sits within ``goal_atol``, so the warning is safe to ignore.
+
+Finally we plot the adapted mesh. The refinement concentrates near the top
+boundary, where the goal functional lives. Refinement driven by a global energy
+norm would instead spread evenly over the domain:
+
+.. code-block:: python
+
+  import matplotlib.pyplot as plt
+  from firedrake.pyplot import triplot
+
+  fig, axes = plt.subplots()
+  triplot(uh.function_space().mesh().unique(), axes=axes)
+  axes.set_aspect("equal")
+  axes.legend(loc="center left", bbox_to_anchor=(1.05, 0.5))
+  fig.savefig("goal_oriented_nonlinear_p_laplacian_mesh.png", bbox_inches="tight")
+
+.. image:: goal_oriented_nonlinear_p_laplacian_mesh.png
+   :width: 60%
+   :alt: The adapted mesh, refined near the top boundary.
+   :align: center
+
+The estimator used here weights the primal residual by the dual error. Richer
+variants exist. Symmetrising the estimate with the adjoint residual leaves a
+remainder cubic in the errors rather than quadratic :cite:`Endtmayer:2024`.
+That costs two further solves per cycle.
+
+A python script version of this demo can be found :demo:`here
+<goal_oriented_nonlinear_p_laplacian.py>`.
+
+.. rubric:: References
+
+.. bibliography:: demo_references.bib
+   :filter: docname in docnames

docs/source/advanced_tut.rst

@@ -27,6 +27,7 @@ element systems.
    1D Vlasov-Poisson equation using vertical independent function spaces.<demos/vp1d.py>
    Degree-independent multigrid convergence using patch relaxation.<demos/poisson_mg_patches.py>
    Multigrid on adaptively-refined mesh hierarchies.<demos/adaptive_multigrid.py>
+   Goal-oriented adaptivity for the nonlinear p-Laplacian.<demos/goal_oriented_nonlinear_p_laplacian.py>
    Monolithic multigrid with Vanka relaxation for Stokes.<demos/stokes_vanka_patches.py>
    Vertex/edge star multigrid relaxation for H(div).<demos/hdiv_riesz_star.py>
    Auxiliary space patch relaxation multigrid for H(curl).<demos/hcurl_riesz_star.py>

tests/firedrake/demos/test_demos_run.py

@@ -52,6 +52,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_nonlinear_p_laplacian", "goal_oriented_nonlinear_p_laplacian"), ["matplotlib"]),
     Demo(('vlasov_poisson_1d', 'vp1d'), []),
     Demo(('shape_optimization', 'shape_optimization'), ["adjoint", "vtk"]),
     Demo(('submesh_reaction_diffusion', 'submesh_reaction_diffusion'), ["netgen", "vtk"]),

tests/firedrake/multigrid/test_snes_adapt.py

@@ -404,6 +404,34 @@ def test_goal_oriented_marker_reconstructs_exact_solution():
     assert solver.get_marking_callback().error_estimate != 0.0
 
 
+@pytest.mark.parallel([1, 2])
+def test_goal_oriented_marker_mixed_space():
+    mesh = UnitSquareMesh(2, 2)
+    W = FunctionSpace(mesh, "RT", 1) * FunctionSpace(mesh, "DG", 0)
+    w = Function(W)
+    sigma, u = split(w)
+    tau, v = TestFunctions(W)
+    F = (inner(sigma, tau) + inner(u, div(tau)) + inner(div(sigma), v) + inner(1, v))*dx
+    n = FacetNormal(mesh)
+    goal = inner(sigma, n)*ds(2)
+
+    num_refinements = 2
+    parameters = _goal_poisson_solver_parameters("snes_adapt_sequence", "refine", num_refinements)
+    parameters["pc_factor_mat_solver_type"] = "mumps"
+    problem = NonlinearVariationalProblem(F, w)
+    solver = NonlinearVariationalSolver(
+        problem, solver_parameters=parameters, marking_callback=GoalOrientedMarker(goal),
+    )
+    result = solver.solve()
+
+    hierarchy, level = get_level(result.function_space().mesh())
+    assert level == num_refinements
+    assert result.function_space().ufl_element() == W.ufl_element()
+    assert result.function_space().dim() > W.dim()
+    assert np.isfinite(solver.get_marking_callback().error_estimate)
+    assert solver.get_marking_callback().error_estimate != 0.0
+
+
 @pytest.mark.parallel([1, 2])
 def test_adaptive_refine_without_marking_callback_is_uniform():
     mesh, V, problem, _ = _goal_poisson_problem(n=2)