3. Goal-oriented mesh adaptation

3.1. Error indicators

Goal-oriented mesh adaptation is the term used to refer to mesh adaptation methods that are driven by goal-oriented error estimates. Recall from the error estimation section of the manual that the fundamental result is the dual weighted residual (DWR) (6), which can be written in approximate form

(1)\[J(u)-J(u_h)\approx\rho(u_h,(u_h^*)^+),\]

where \((u_h^*)^+\) denotes an approximation of the adjoint solution from a (finite-dimensional) superspace of the base finite element space, i.e. \(V_h^+\supset V_h\).

Mesh adaptation is all about wisely using varying resolution to equidistribute errors. This means increasing resolution where errors are deemed to be high and reducing it where errors are deemed to be low. We cannot immediately deduce spatially varying information from (1) as it is currently written. A simple, but effective way of extracting such information is to compute the element-wise contributions,

(2)\[\rho(u_h,(u_h^*)^+)=\sum_{K\in\mathcal{H}}i_K,\]

where

(3)\[i_K:=\rho(u_h,(u_h^*)^+)|_K.\]

is the so-called error indicator for element \(K\). The second output of indicate_errors() contains error indicator fields – element-wise constant fields, whose value on \(K\) is the indicator \(i_K\).

Note that so far we have only discussed how to estimate and indicate errors for steady-state problems. The extension to the time-dependent case is similar, in the sense that we take a timestep-by-timestep approach in time like how we take an element-by-element approach in space. Avoiding some minor details, we can apply all of the subsequently discussed methodology to the weak residual associated with a single timestep. For simple time integration schemes like Implicit Euler, the main difference will be that the weak residual also includes a term that discretises the time derivative.

3.2. Adapting based on error indicators

Once error indicator fields have been computed, there are many different ways to perform mesh adaptation. One common approach is to use a hierarchical type approach, where the mesh resolution is locally increased in elements where the indicator value falls below a pre-specified tolerance and is locally decreased if the indicator values are already lower than required. This is sometimes called “adaptive mesh refinement (AMR)”, although the literature is not entirely consistent on this. The terms “quad-tree” and “oct-tree” are used when it is applied to quadrilateral and hexahedral meshes, respectively. Sadly, this form of mesh adaptation is not supported in Firedrake.

As described in the Animate documentation, metric-based mesh adaptation is another approach which is currently being integrated into Firedrake and is supported on developer branches. In that case, we need to construct a Riemannian metric from an error indicator field. Goalie provides a number of different driver functions for doing this. The simplest is isotropic_metric(), which takes an error indicator field and constructs a metric that specifies how a mesh should be adapted purely in terms of element size (not orientation or shape). Alternative anisotropic formulations, which combine error indicator information with curvature information from the Hessians of solution fields are provided by anisotropic_metric(). See the API documentation (and references therein) for details.

The metric-based framework has been the main backend used for goal-oriented mesh adaptation during the development of Goalie. However, this does not mean other approaches would not work. Mesh movement (or \(r\)-adaptation) has been implemented in Firedrake on a number of occasions (such as [CWK+22, MCB18]) and could plausibly be driven by goal-oriented error estimates.

3.3. Fixed point iteration loops

In some mesh adaptation approaches, it is common practice to adapt the mesh multiple times until convergence is attained, in some sense. This is often the case under metric based mesh adaptation, for example. Goalie includes two methods to facilitate such iterative adaptation approaches, as described in the following.

In the non-goal-oriented case, there is the fixed_point_iteration() method, which accepts a Python function that describes how to adapt the mesh as an argument. This provides the flexibility to use different adaptation routines. The function should take as input the MeshSeq instance and the dictionary of solution fields from each timestep. For example, it could compute the Hessian of the solution field at each timestep using compute_hessian(), ensure that it is a metric using enforce_spd() and then combine this information using averaging or intersection. In each iteration, the PDE will be solved over the sequence of meshes (with data transferred inbetween) using solve_forward() and a Hessian-metric will be constructed on each mesh. All of the meshes are then adapted.

The iteration is terminated according to AdaptParameters when either the pre-specified maximum iteration count maxiter is reached, or the meshes no longer change substantially. This is the case when none of the corresponding element counts changes more than the relative tolerance element_rtol.

Note

The convergence criteria are not actually checked until the minimum iteration count miniter has been met.

In the goal-oriented case, the fixed_point_iteration() method takes a similar form, except that indicate_errors() is called, rather than solve_forward(). That is, the forward problem is solved over all meshes in the sequence, then the adjoint problem is solved over all meshes in reverse and finally goal-oriented error indicators are computed. As such, the adaptor function depends on these error indicators, as well as the MeshSeq and solution field dictionary (which contains solutions of both the forward and adjoint problems).

Like with the simpler case, the goal-oriented fixed point iteration loop is terminated if the maximum iteration count or relative element count convergence conditions are met. However, there are two additional convergence criteria defined in GoalOrientedParameters. Convergence is deduced if the QoI value changes between iterations by less than qoi_rtol. It is similarly deduced if the error estimate value changes by less than estimator_rtol.

3.4. References

[CWK+22]

Mariana CA Clare, Joseph G Wallwork, Stephan C Kramer, Hilary Weller, Colin J Cotter, and Matthew D Piggott. Multi-scale hydro-morphodynamic modelling using mesh movement methods. GEM-International Journal on Geomathematics, 13(1):1–39, 2022. doi:10.1007/s13137-021-00191-1.

[MCB18]

Andrew TT McRae, Colin J Cotter, and Chris J Budd. Optimal-transport–based mesh adaptivity on the plane and sphere using finite elements. SIAM Journal on Scientific Computing, 40(2):A1121–A1148, 2018. doi:10.1137/16M1109515.