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
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,
where
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¶
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.
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.