# E2N: Error Estimation Networks for Goal-Oriented Mesh Adaptation

Joseph G. Wallwork<EMAIL_ADDRESS>Jingyi Lu  Mingrui Zhang
Matthew D. Piggott  Department of Earth Science and Engineering, Imperial
College London, London, SW7 2AZ, UK  Huawei Technologies Co., Ltd, Shenzhen,
China (current address)

###### Abstract

Given a partial differential equation (PDE), goal-oriented error estimation
allows us to understand how errors in a diagnostic quantity of interest (QoI),
or _goal_ , occur and accumulate in a numerical approximation, for example
using the finite element method. By decomposing the error estimates into
contributions from individual elements, it is possible to formulate adaptation
methods, which modify the mesh with the objective of minimising the resulting
QoI error. However, the standard error estimate formulation involves the true
adjoint solution, which is unknown in practice. As such, it is common practice
to approximate it with an ‘enriched’ approximation (e.g. in a higher order
space or on a refined mesh). Doing so generally results in a significant
increase in computational cost, which can be a bottleneck compromising the
competitiveness of (goal-oriented) adaptive simulations. The central idea of
this paper is to develop a “data-driven” goal-oriented mesh adaptation
approach through the selective replacement of the expensive error estimation
step with an appropriately configured and trained neural network. In doing so,
the error estimator may be obtained without even constructing the enriched
spaces. An element-by-element construction is employed here, whereby local
values of various parameters related to the mesh geometry and underlying
problem physics are taken as inputs, and the corresponding contribution to the
error estimator is taken as output. We demonstrate that this approach is able
to obtain the same accuracy with a reduced computational cost, for adaptive
mesh test cases related to flow around tidal turbines, which interact via
their downstream wakes, and where the overall power output of the farm is
taken as the QoI. Moreover, we demonstrate that the element-by-element
approach implies reasonably low training costs.

###### keywords:

error estimation , mesh adaptation , machine learning , neural networks ,
tidal power , Thetis

††journal: Journal of Computational Physics

capbtabboxtable[][0.48]

##  1 Introduction

###  1.1 Goal-Oriented Error Estimation

Numerical simulations across the majority of science and engineering often
have an associated diagnostic _quantity of interest (QoI)_. This could be the
drag or lift on an aeroplane wing, or the level of radiation in a particularly
sensitive region in the vicinity of a nuclear power plant, for example. As
indicated in the abstract, the motivating application utilised in this paper
is in the context of tidal resource assessment for the marine renewable energy
sector, whereby the instantaneous power or total energy output of a tidal farm
provide obvious choices for QoIs [1]. The abstract and idealised nature at
which this is considered makes the results highly relevant to the closely
related wind turbine application.

_Goal-oriented error estimation_ provides a framework for assessing the error
in the chosen QoI due to a particular numerical approximation method for
solving the underlying PDE. This is achieved by means of the _adjoint problem_
associated with the PDE and the QoI. Goal-oriented error estimates are
typically derived from the _dual weighted residual (DWR)_ [2]. The ‘first
order’ DWR result involves the weak residual of the PDE – assuming some finite
element discretisation – but with the test function replaced with the solution
of the adjoint problem. Since this ‘exact’ adjoint solution is not obtainable
in practice, we are unable to evaluate the error estimate without the use of
further approximation methods.

One option is to apply error estimation techniques, such as a posteriori upper
bounds, to construct an alternative formulation, wherein the true adjoint
solution does not appear (see, for example, the ‘difference quotient’ approach
in [2]). Another approach is to approximate the adjoint solution in an
enriched finite element space. This can be done by solving the PDE in a
globally enriched space (see [3], for example), solving PDEs in locally
enriched spaces (see [4], for example), or by applying superconvergent patch
recovery techniques (see [5], for example). In any case, the enrichment
procedure tends to be computationally expensive, although the level of expense
can vary between approaches.

###  1.2 Accelerating Mesh Adaptation using Neural Networks

Goal-oriented error estimation is most commonly deployed for the purposes of
mesh adaptation. By determining local contributions to an estimator (i.e.
_‘error indicators’_), it is possible to establish parts of the mesh that
would benefit from increased resolution, as well as parts of the mesh with
unnecessary degrees of freedom (DoFs). The resulting numerical accuracy versus
efficiency trade-off may be interpreted as a constrained optimisation problem,
whereby we seek to minimise QoI error, subject to having a mesh with
(approximately) some target number of DoFs. However, for the reasons mentioned
above, the goal-oriented error estimation step typically comes at significant
computational expense. In many cases, it can be more computationally demanding
than solving the underlying prognostic PDE, making it a bottleneck within the
goal-oriented mesh adaptation workflow.

Neural networks provide an exciting new way to accelerate advanced
discretisation methods like goal-oriented mesh adaptation. Given the
computational expense of such methods, there is potential for massive gains in
efficiency and widening of access to users that would otherwise view them as
too time intensive. As remarked in [6, p.15], “[the application of neural
networks to emulate mesh adaptation] in the general CFD pipeline is one of the
most risk-averse ways of using machine learning in this context”. That is, the
worst that can happen is that the resulting adapted mesh is of relatively poor
quality and the solver convergence suffers as a consequence. On the other
hand, if the PDE solve step is to be emulated, for example, then an incorrect
application could potentially result in unrealistic simulation data, i.e.
“incorrect physics”.

Whilst this is a relatively new field of research, there are already a number
of existing works. A variety of different approaches have been taken, with
different authors advocating the replacement of different parts of the (goal-
oriented) mesh adaptation workflow by neural networks. These include the test
function choice [7], forward solve [8], adjoint solve [9], derivative recovery
procedure [10], error estimation [11], metric/monitor function/sizing field
construction step [12, 13, 14], and the entire mesh adaptation loop [6, 15,
16].

One of the most relevant references for this work is [13]. In that work, a
neural network is used to ‘solve’ the optimisation problem used to determine a
metric tensor from solution information in the _Mesh Optimization through
Error Sampling and Synthesis (MOESS)_ approach [17]. Each iteration of the
MOESS algorithm requires the consideration of all possible modifications that
can be made to each mesh element and the solution of local PDEs defined on
patches of modified elements and their neighbours. The application of neural
networks means that such sub-meshes and sub-problems do not need to be
explicitly considered, as the mapping from the current state to the optimised
Riemannian metric is ‘learnt’ implicitly during training. This methodology is
tested on steady-state ‘flow past an airfoil’ test cases in [13] and is
extended to use convolutional neural networks in [14].

The authors of [6] point out that, whilst goal-oriented mesh adaptation is a
‘gold-standard’ approach to mesh adaptation, it can be prohibitively
computationally expensive and even infeasible to use, given that many solvers
do not have automated adjoint capability. That work goes further than
emulating the metric construction step of the mesh adaptation pipeline,
covering the forward and adjoint PDE solve steps as well, so that the user of
a trained network does not need to solve adjoint problems at all. Goal-
oriented mesh adaptation is applied to an unsteady compressible flow solver in
a range of rectangular domains with obstacles with randomly generated
geometries, where the QoI is the drag on the object. The domain is interpreted
as a grayscale image comprised of pixels, which may be interpreted as a
rectangular input array so that a deep convolutional neural network can be
applied. The output is a mesh sizing field, which can be likened to a monitor
function or an isotropic metric field.

The works discussed in the previous two paragraphs concern topology-deforming
(i.e. $h$-adaptive) methods. There also exist several papers concerning
topology-preserving (i.e. $r$-adaptive) methods, although none of these
methods use the goal-oriented approach. One example is [8], which uses the
variational mesh movement framework of [18]. In that framework, the metric
tensor now prescribes an anisotropic extension of the ‘mesh density’ concept.
Given the current mesh at a timestep, a metric tensor is constructed based on
the associated PDE solution, often in terms of its gradient or curvature. A
new mesh that is optimal w.r.t. the metric tensor is determined by the
iterative solution of a so-called _mesh movement partial differential equation
(MMPDE)_ [19]. MMPDEs are usually nonlinear and are therefore solved
numerically using iterative methods. At each iteration, the metric must be
recomputed. However, it is undesirable to solve the physical PDE multiple
times during this process from a computational efficiency perspective. As
such, the authors of [8] train a neural network to emulate the PDE solution
itself and use this approximate approach at each iteration of the MMPDE solve.
They then solve the full PDE using conventional methods on the final mesh, so
that physics are accurately represented. The authors of [16] go further and
propose an end-to-end $r$-adaptation method, whereby the solution on the
current mesh is mapped directly to the adapted mesh, removing the need to
solve any MMPDEs. This can be particularly beneficial for highly nonlinear
MMPDEs such as those of Monge-Ampère type, which can be computationally
expensive to solve.

As well as through the use of a remesher, there exist methods for optimising
the input mesh directly. In [20], an adjoint-based differentiable fluid
simulator is embedded into a graph convolution network based model as an
additional layer. With the gradient information computed by the adjoint
method, the input mesh of the simulator is regarded as a trainable parameter,
to be trained jointly with the network. Similarly, [21] and [22] use neural
networks to guide mesh generation. From the user perspective, [6] also
describes a mesh generator, albeit one that emulates a goal-oriented mesh
adaptation loop.

As mentioned above, [11] emulates goal-oriented error estimation using neural
networks. However, that work uses a collocation method based on B-splines to
discretise the strong form of a PDE, as opposed to discretising a weak form
using finite elements. As such, technically it does not do mesh adaptation,
since there is no ‘mesh’ concept. However, the fact that collocation point
locations are optimised means that the discretisation is being adapted based
on the data-driven error estimation process.

###  1.3 Novel Aspects of the Present Work

Previous work by the authors of this paper has applied goal-oriented error
estimation techniques to tidal energy resource assessment with power or energy
output as the QoI (see [3] for the steady-state case and [23] for the time-
dependent case). We build upon this work and seek to accelerate the numerical
experiments presented in [3] using neural networks. For consistency with that
work, we opt to compute the enriched adjoint solution by solving for it on a
globally uniformly refined mesh. For triangular meshes in 2D, this implies an
auxiliary PDE solve on a mesh with four times as many elements as the base
mesh (sometimes referred to as ‘red’ refinement), in which case the error
estimation step becomes a significant portion of the overall computational
cost.

Recall that the authors of [13] train a neural network to emulate a metric
construction strategy. In this work, we instead train a neural network to
emulate the error estimation strategy. This choice is made for a number of
reasons. Firstly, it means that we retain the flexibility to use any metric
formulation – rather than being constrained to the one that the network was
trained on – or indeed to use a different “optimal mesh” concept entirely,
such as one using a monitor function. In a similar vein, we do not need to use
the output of the network for mesh adaptation at all and could perform error
estimation analysis instead. Finally, we seek to apply the neural network in a
focused or ‘surgical’ manner, to a single, computationally expensive component
of the mesh adaptation workflow, as opposed to treating it as a ‘black box’
mesh adaptation algorithm.

A similar ‘focused’ approach is also used in [9], which emulates the adjoint
solve procedure. This is done on the base mesh and the data-driven adjoint
solution is projected into an enriched space, where error indicators are
assembled and thereby used to drive mesh adaptation. The main difference
herein is that we avoid the construction of enriched finite element spaces
altogether, using a neural network to map directly from objects defined on the
base mesh to an error indicator on that same mesh. By avoiding the
construction of enriched spaces, our approach reduces the memory footprint of
the resulting algorithm significantly.

The remainder of the paper is as follows. Section 2 provides background
material on goal-oriented error estimation and mesh adaptation, as well as
introducing the method used for tidal turbine modelling. Section 3 describes
the methodology for emulating the goal-oriented error estimation process,
which is applied to steady-state tidal turbine test cases in Section 4.
Finally, conclusions are drawn and an outlook on future work is made in
Section 5. Details on the derivative recovery procedure and software used to
generate the numerical results presented in this paper are provided in A and
B.

##  2 Background

###  2.1 Forward and Adjoint Problems

Consider a PDE with solution $u$ from a function space $V$ and a variational
formulation,

$\rho(u,v)=0,\quad\forall v\in V,$ (1)

where the residual $\rho:V\times V\rightarrow\mathbb{R}$ is linear in its
second argument. For the purposes of this paper, we consider steady-state
problems, meaning that all derivatives are spatial and there is no time-
dependence.

Suppose there exists a diagnostic QoI, $J:V\rightarrow\mathbb{R}$, which we
would like to accurately evaluate. Given such a QoI, we may formulate the
adjoint equation associated with (1) as

$\mathcal{L}_{u}\>\rho(v,u^{*})=\mathcal{L}_{u}\>J(v),\quad\forall v\in V,$
(2)

where $\mathcal{L}_{w}\>\rho:V\times V\rightarrow\mathbb{R}$ and
$\mathcal{L}_{w}\>J:V\rightarrow\mathbb{R}$ denote the linearisations of
$\rho$ and $J$ about their first arguments, in direction $w\in V$. Here
$u^{*}$ is the so-called _adjoint solution_ , which (2) is linear with respect
to.

In the finite element method, we approximate $u$ in a finite-dimensional
subspace, $V_{h}\subset V$, based upon a mesh $\mathcal{H}_{h}$ of the spatial
domain. This gives rise to a weak formulation for (1),

$\rho(u_{h},v)=0,\quad\forall v\in V_{h},$ (3)

where $u_{h}\in V_{h}\subset V$ is the finite element approximation.
Similarly, we approximate the adjoint solution by solving the discrete adjoint
equation

$\mathcal{L}_{u_{h}}\>\rho(v,u^{*}_{h})=\mathcal{L}_{u_{h}}\>J(v),\quad\forall
v\in V_{h}.$ (4)

Note that the finite element adjoint problem is linearised about $u_{h}$,
rather than $u$.

###  2.2 Error Estimation and Indication

We are primarily interested in goal-oriented error estimators of _dual
weighted residual (DWR)_ type. In particular, the ‘first order’ DWR estimator,
$\mathcal{E}$, which satisfies the relation

$J(u)-J(u_{h})=\underbrace{\rho(u_{h},u^{*}-u^{*}_{h})}_{\mathcal{E}}+R^{(2)},$
(5)

where the remainder term $R^{(2)}$ is quadratic in the forward and adjoint
errors $u-u_{h}$ and $u^{*}-u^{*}_{h}$. Such an error estimator can be
decomposed on an element-wise basis, in order to convey the error contribution
due to each mesh element:

$\mathcal{E}=\sum_{K\in\mathcal{H}_{h}}\mathcal{E}_{K},\qquad\mathcal{E}_{K}:=\rho(u_{h},u^{*}-u^{*}_{h})|_{K}.$
(6)

This decomposition can be used to generate an indicator field, $i_{h}$, which
takes a constant value on each element:

$i_{h}:=\sum_{K\in\mathcal{H}_{h}}i_{K},\qquad
i_{K}(\mathbf{x}):=\left\\{\begin{array}[]{cc}\mathcal{E}_{K}&\mathbf{x}\in
K\\\ 0&\mathbf{x}\not\in K\end{array}\right..$ (7)

The DWR estimator provides a first order approximation to the error in the QoI
by combining residuals of the forward equation and sensitivity information
related to the QoI via the adjoint solution. However, as discussed in Section
1, there is a practical issue: the estimator involves the true adjoint
solution, $u^{*}$, which is in general unknown. As such, it is not possible to
use $\mathcal{E}$ in the form presented in (5) and we instead use an
approximate ‘enrichment method’. The main research question in this paper is:
_can we reliably emulate the (typically) expensive procedure using efficient
machine learning techniques?_

Suppose we choose to evaluate $\mathcal{E}$ by forming and solving the adjoint
problem again on a (globally) uniformly refined mesh and substituting its
solution in place of the true adjoint solution. This accurate – but still
approximate – approach allows us to reliably produce a high quality
representation of the error estimate. However, it is also one of the most
computationally expensive means of doing so (see the discussion in [24,
Section 7.3]). It requires additional data structures, namely a mesh hierarchy
and a means for transferring solution data between its levels. The subscript
$h$ in the notation $\mathcal{H}_{h}$ for the (‘current’) base mesh refers to
the element sizing function, so let $\mathcal{H}_{h/2}$ denote the refined
mesh generated by halving each element’s size in each direction, e.g.
splitting a triangular element into four triangles of equal area in what is
sometimes referred to as ‘red’ refinement. The associated function spaces are
denoted $V_{h}$ and $V_{h/2}$ and satisfy $V_{h}\subset V_{h/2}$, for the
Lagrange elements considered in this paper. Further, denote by $u^{*}_{h/2}\in
V_{h/2}$ the adjoint solution computed on $\mathcal{H}_{h/2}$, which we use to
approximate $u^{*}$. Approximating $\mathcal{E}_{K}$ amounts to the
computation

$\widetilde{\mathcal{E}_{K}}:=\rho\left.\left(P_{h/2}[u_{h}],u^{*}_{h/2}-P_{h/2}[u^{*}_{h}]\right)\right|_{K},\qquad
K\in\mathcal{H}_{h/2},$ (8)

where $P_{h/2}:V_{h}\rightarrow V_{h/2}$ transfers into the enriched space.
The weak residual is computed in this space. The error indicator may then be
approximated in the base space by applying a projection operator
$\Pi_{h}:V_{h/2}\rightarrow V_{h}$:

$\widetilde{i_{h}}:=\Pi_{h}\left[\sum_{K\in\mathcal{H}_{h/2}}\widetilde{i}_{K}\right],\qquad\widetilde{i}_{K}:=\left\\{\begin{array}[]{cc}\widetilde{\mathcal{E}_{K}}&\mathbf{x}\in
K\\\ 0&\mathbf{x}\not\in K\end{array}\right..$ (9)

The reason that we transfer the error indicator information back down to the
base space is that we are primarily motivated by mesh adaptation, which is
driven by such information.

In this work, we choose $P_{h/2}$ to be a prolongation operator (which is
lossless for the finite elements we consider) and $\Pi_{h}$ to be a
conservative projection. As such, successfully emulating the error estimate
evaluation using machine learning is about more than just generating a better
approximation of the adjoint solution; it also involves encoding information
related to the enriched space and the associated interpolation operators.
However, at no point in the data-driven approach are these spaces actually
constructed.

A naïve implementation of $u^{*}_{h/2}$ in (8) involves solving the _forward_
problem in the enriched space and then forming the adjoint problem by
linearising about its solution. The additional cost of solving the enriched
forward problem can be significant, especially for nonlinear problems. In
order to avoid this step during training, we calculate $u^{*}_{h/2}$ in an
approximate sense: form the adjoint problem in the enriched space by
linearising about _the forward solution prolonged from the base space_ ,
$P_{h/2}[u_{h}]$. It was shown in [24, Section 7.5] that this approach yields
similarly effective error estimates to the naïve approach, at a dramatically
reduced computational cost, since prolongation is orders of magnitude cheaper
than global PDE solves.

###  2.3 Tidal Turbine Modelling

In this work, we take as our physical demonstration scenario the depth-
averaged modelling of coastal ocean flows using the _shallow water equations_.
This hydrostatic model is derived from the Navier-Stokes equations and
accounts for nonlinear advection, viscous forces, a free upper surface
inducing pressure gradient like forces, and bed friction (or ‘drag’) due to
sea bed roughness. In this setting, the spatial domain is two-dimensional:
$\Omega\subset\mathbb{R}^{2}$. The prognostic variables are horizontal
velocity $\mathbf{u}:\Omega\rightarrow\mathbb{R}^{2}$ and free surface
elevation $\eta:\Omega\rightarrow\mathbb{R}$. We consider the steady-state
shallow water equations in the non-conservative form,

$\mathbf{u}\cdot\nabla\mathbf{u}+g\nabla\eta+C_{\mathcal{D}}\frac{\|\mathbf{u}\|\mathbf{u}}{\eta+b}=\nabla\cdot(\underline{\boldsymbol{\nu}}\nabla\mathbf{u}),\qquad\nabla\cdot((\eta+b)\mathbf{u})=0,$
(10)

where $g=9.81\,\mathrm{m\,s}^{-2}$ is gravitational acceleration, $b$
[$\mathrm{m}$] is the bathymetry (sea floor topography), $C_{\mathcal{D}}$ is
a dimensionless drag coefficient and $\underline{\boldsymbol{\nu}}$
[$\mathrm{m}^{2}\,\mathrm{s}^{-1}$] is the viscosity tensor. For the purposes
of this paper, we assume an isotropic viscosity coefficient, defined by a
positive scalar field, $\underline{\boldsymbol{\nu}}=\nu\underline{\bf I}$.

Tidal stream turbines are comprised of blades and nacelle positioned on a
vertically oriented shaft, where the plane of blade rotation is perpendicular
to the horizontal component of the fluid flow. As such, they are inherently 3D
objects. However, it is possible to represent them in a 2D setup using a drag
parametrisation, with the presence of a turbine being modelled by a locally
artificially increased drag coefficient. As such, we decompose
$C_{\mathcal{D}}=C_{\mathcal{B}}+C_{\mathcal{T}}$ into a constant background
drag $C_{\mathcal{B}}=0.0025$ due to sea bed friction, and a spatially-varying
drag due to the presence of turbines,

$C_{\mathcal{T}}=\sum_{T\in\mathcal{T}}C_{T}\mathbbm{1}_{T}.$ (11)

Here $\mathcal{T}$ denotes the set of all turbines. In a 2D model, turbine
$T\in\mathcal{T}$ is parametrised over its footprint, which is indicated by
$\mathbbm{1}_{T}$ and has area $A_{T}^{\mathrm{footprint}}$. Clearly, the
turbines extract power based on the area $A_{T}^{\mathrm{swept}}$ that they
sweep in the vertical, not their footprint. As such, we apply the scaling

$C_{T}=\frac{1}{2}\frac{A_{T}^{\mathrm{swept}}}{A_{T}^{\mathrm{footprint}}}c_{T},$
(12)

where the so-called _thrust coefficient_ $c_{T}$ is related to details of the
manufacture and operation of the turbines. We assume a constant thrust
coefficient of 0.8 herein, for simplicity. An additional correction is applied
to account for the fact that the thrust coefficient is defined in terms of the
upstream velocity, rather than the depth-averaged value at the turbine
location which is the value immediately available in the shallow water model
(see [25] for further details).

For the purposes of this paper, we take the QoI to be the power output of the
tidal farm, a proxy of which may be computed as

$J(\mathbf{u},\eta)=\int_{\Omega}\rho\,C_{\mathcal{T}}\|\mathbf{u}\|^{3}\;\mathrm{d}x,$
(13)

where $\rho=1030.0\,\mathrm{kg\,m}^{-3}$ is the assumed density of seawater.
The power output has units of Watts.

The shallow water equations (10) are solved using the coastal ocean model
_Thetis_ [26]. We choose to discretise using the mixed discontinuous-
continuous $\mathbb{P}1_{DG}-\mathbb{P}2$ finite element pair, first
introduced in [27]. The corresponding weak formulation can be found in [3,
Section 3.1]. The resulting nonlinear systems are solved using a Newton method
with line search, with the underlying linear systems solved by application of
a full LU decomposition as a preconditioner [28, 29].

Goal-oriented error estimates for the discretisation described above were
derived in [3, Section 3.2]. We refer to that work for details.

##  3 Methodology

###  3.1 Inputs and Outputs

As discussed above, we seek to emulate DWR error estimation using a machine
learning approach on an element-by-element basis. This means that we seek to
map inputs related to a particular mesh element $K$ of the base mesh to the
associated error contribution, $\mathcal{E}_{K}$. As such, for each element,
there is a single real number to be computed. There are a number of inputs
that we deem to be important, as detailed in the following. All of them are
real-valued, meaning that we effectively seek to obtain a mapping

$\mathrm{E2N}:P\rightarrow\mathbb{R},\qquad P\subseteq\mathbb{R}^{n},$ (14)

where $n$ is the number of scalar input parameters.

##### Coarse DWR approximation

Our objective is to emulate the DWR error estimator (8) without constructing
enriched finite element spaces. A naïve way to do this is to simply replace
the true adjoint solution with the finite element adjoint solution on the base
mesh111Whilst (15) should sum to zero over all mesh elements, since (3) holds
with $v=u^{*}_{h}$, this does not mean that the individual elemental
contributions are zero.:

$\rho(u_{h},u^{*}_{h})|_{K}.$ (15)

This is a low accuracy approximation of a DWR indicator, but is a potentially
useful input feature for the neural network, which is readily calculable on
the base mesh.

##### Physics-based parameters

The next set of input parameters that we consider relate to the flow physics.
Given that the gravitational acceleration and seawater density are held fixed
across test cases, the remaining parameters defining (10) are the viscosity,
drag coefficient and bathymetry. In the case of spatially varying parameter
fields, these are averaged over the element to provide scalar inputs to the
element-wise network.

##### Mesh element parameters

We also consider a number of input parameters related to the mesh element
itself, as detailed in the following. In this work we restrict attention to
meshes comprised of triangular elements.

Mesh element $K$ may be obtained from the reference element $\widehat{K}$
using an affine transformation with Jacobian $\underline{\bf J}_{K}\;$222We
assume $\widehat{K}$ to be the right-angled triangle formed by vertices
$(0,0)$, $(0,1)$ and $(1,0)$. This is in agreement with the convention in
Firedrake.. The matrix $\underline{\bf J}_{K}^{T}\underline{\bf J}_{K}$ is
symmetric positive-definite (SPD), by construction. Symmetry implies the
orthogonal eigendecomposition,

$\underline{\bf J}_{K}^{T}\underline{\bf J}_{K}=\underline{\bf
V}_{K}\underline{\boldsymbol{\Lambda}}_{K}\underline{\bf V}_{K}^{T}.$ (16)

Positive-definiteness implies that the eigenvalues
$\Lambda_{K}=\mathrm{diag}(\lambda_{1,K},\lambda_{2,K})$ are both positive,
meaning we can define $h_{1}:=1/\sqrt{\lambda_{1,K}}$ and
$h_{2}:=1/\sqrt{\lambda_{2,K}}$. Assume w.l.o.g. that the eigenvalues are in
ascending order, so that $h_{1}\geq h_{2}>0$.

A convenient geometric interpretation for (16) is in terms of the minimal
ellipse containing the triangle’s vertices. Such an ellipse is shown in Figure
16, for an arbitrary triangle drawn in blue. The eigenvectors represent the
major semi-axes of the ellipse, with $h_{1}$ being the magnitude of the first
and $h_{2}$ being the magnitude of the second. We define the parameter
$\theta\in[0,2\pi)$ to be the argument of the first eigenvector.

$\widehat{\mathbf{x}}$$h_{1}=\frac{1}{\sqrt{\lambda_{1}}}$$h_{2}=\frac{1}{\sqrt{\lambda_{2}}}$$\theta$

Figure 1:  Ellipse representation of the quantities $h_{1}$, $h_{2}$ and
$\theta$ related to a triangular element. The vector $\widehat{\mathbf{x}}$
points in the direction of the positive $x$-axis.

$\mathbb{P}0$$\mathbb{P}1_{DG}$$\mathbb{P}1$$\mathbb{P}2$

Figure 2:  Diagrams representing the degrees of freedom of various finite
elements. If a node is drawn on an edge then this implies that continuity is
enforced with a neighbouring element.

With the above quantities, we can encapsulate the size of the element with
$d:=h_{1}h_{2}$, the orientation of the element with $\theta$ and the shape of
the element with the aspect ratio $s:=h_{1}/h_{2}$. These three parameters
uniquely define the element. However, as discussed in [13], it is potentially
problematic that orientations $\theta\approx 0$ and orientations
$\theta\approx 2\pi$ are strongly correlated, even though this is not apparent
in polar coordinates. If this parameter were used as input to a neural
network, it could hinder its performance. Instead, we follow the
recommendations of [13] by replacing $\theta$ and $s$ with, respectively,

$\frac{1}{s}\cos^{2}\theta+s\sin^{2}\theta\qquad\text{and}\qquad\left(s-\frac{1}{s}\right)\sin\theta\cos\theta.$
(17)

In addition to the above mesh-based parameters, we also include a parameter
that integrates unity over the intersection of the element boundary and the
domain boundary:

$\int_{\partial K\cap\partial\Omega}\;\mathrm{d}s.$ (18)

Clearly, (18) is zero if the element does not touch the domain boundary.
Otherwise, it has a positive value, whose magnitude is determined by the
length of the shared boundary segment.

##### Solution fields

Given that we seek to perform goal-oriented error estimation, there are many
reasons the forward and adjoint solutions at each DoF make valuable choices of
inputs. We can assume that the adjoint solution from the base space provides a
reasonable first approximation to the enriched version. As such, it may well
be key to guiding the network towards emulating the enrichment process.
Perhaps more importantly, the adjoint solution encodes the QoI, which is
ultimately what we would like to accurately capture. The approximate forward
solution is potentially important because the adjoint solve in the enriched
space involves a linearisation about its prolongation.

An average value of one of the solution fields over an element is typically
insufficient information for the purposes of goal-oriented error estimation,
since intra-element variation often plays an important role. As such, the
input data should reflect this variability. One option is to simply provide
the DoF values. For the (relatively simple) finite elements used in this
paper, all of the DoFs are given by point evaluations. Figure 2 represents the
layouts of these DoFs on an equilateral triangle. The nodal positions are
geometrically accurate for the $\mathbb{P}1$ and $\mathbb{P}2$ diagrams. For
the $\mathbb{P}0$ case, the single DoF is associated with the cell. The
locations of the $\mathbb{P}1_{DG}$ DoFs are identical to the $\mathbb{P}1$
case (at the corners), it is just that inter-element continuity is not
enforced.

A problem with simply supplying the DoFs, however, is that the ordering of the
solution field DoFs on each element is essentially arbitrary in the finite
element package used in this work. This means that the link between the value
of a DoF and its location is lost. To combat this issue, we choose the input
data for a field $f\in\mathbb{P}p_{DG}$ (noting that
$\mathbb{P}p\subseteq\mathbb{P}p_{DG}$) to be its value at the element
centroid, $\mathbf{x}_{K}$ and the corresponding values of up to its $p^{th}$
derivatives. For mixed derivatives, the values are averaged across orderings.
In particular, for a degree 2 field, we take the quantities

$f(\mathbf{x}_{K}),\>\frac{\partial f}{\partial
x}(\mathbf{x}_{K}),\>\frac{\partial f}{\partial
y}(\mathbf{x}_{K}),\>\frac{\partial^{2}f}{\partial
x^{2}}(\mathbf{x}_{K}),\>\frac{1}{2}\left(\frac{\partial^{2}f}{\partial
x\partial y}(\mathbf{x}_{K})+\frac{\partial^{2}f}{\partial y\partial
x}(\mathbf{x}_{K})\right),\>\frac{\partial^{2}f}{\partial
x^{2}}(\mathbf{x}_{K}).$ (19)

Because of the averaging, a degree $p$ field gives rise to
$\frac{1}{2}(p+1)(p+2)$ input data points (the $p^{th}$ triangular number),
which coincides with the DoF count on a triangular element. For a
$\mathbb{P}1_{DG}-\mathbb{P}2$ finite element pair, there are twelve DoFs per
element (six velocity DoFs and six elevation DoFs). These imply twelve input
parameters for the forward solution and twelve for the adjoint.

##### Summary

Overall, we have one input parameter that approximates the DWR estimator
itself on the coarse mesh (see eq. (15)), three input parameters related to
the “physics” of the problem being solved, four related to the mesh element
and 24 input parameters related to the forward and adjoint solution fields.
This gives 32 input parameters in total.

###  3.2 Training

We consider a finite set of test cases, the majority of which are used for
training the neural network (i.e.  _supervised learning_) and the remainder
being reserved as ‘unseen’ problems. For each test case, training data is
generated by applying goal-oriented error _estimation_ (involving the
application of the costly enrichment methods) on an initial mesh constructed
using _gmsh_ [30]. In addition, we apply goal-oriented mesh _adaptation_ based
on the error indicators, until convergence of the QoI or mesh element count
(see [24, Algorithm 3]). A minimum of three mesh adaptation iterations is
used, which we have found to be sufficient for mesh anisotropy to become fully
developed. Training data is extracted from these first three adaptation
iterations, since they are guaranteed to exist for all such test cases.

###  3.3 Network Architecture

Table 4 details the parameter values used for the neural network. We choose a
feed-forward fully connected neural network with a single hidden layer, which
is represented as a graph diagram in Figure 4. A standard sigmoid activation
function is used. All feature data is pre-processed using the arctangent
function, which acts as a normalisation step and avoids biasing towards
extreme values. We opt to use a relatively simple network architecture for
this preliminary work, with the aim of proving that it is possible to
successfully emulate goal-oriented error estimation strategies without
expending significant effort on optimisation.

Network type | Fully connected
---|---
Number of inputs | 32
Number of hidden layers | 1
Number of hidden neurons | 64
Number of outputs | 1
Cost function | Mean square error
Activation function | Sigmoid
Preprocessing function | Arctangent
Learning rate | 0.001
Number of epochs | 2,000
Batch size | 500
Training/validation ratio | 7:3
Optimisation routine | Adam
Software | Pytorch

Figure 3:  A table describing the neural network configuration and parameter
values.

Figure 4:  Graph representation of the network, obtained using [31].

The neural network and its training are implemented in Python using _PyTorch_
[32]. _Adam_ [33] is used as the stochastic optimisation routine.

###  3.4 Metric-based Mesh Adaptation

Whilst the network that we seek to construct relates to error estimation
rather than mesh adaptation, we do use different input meshes to generate the
training data, some of which are obtained following mesh adaptation. Mesh
adaptation is also used in Section 4. As such, we outline the approach in the
following.

##### Background

The metric-based framework was first introduced in [34] and is a mesh
adaptation approach based on _Riemannian metric fields_. A Riemannian metric,
$\mathcal{M}=\\{\underline{\bf M}(\mathbf{x})\\}_{\mathbf{x}\in\Omega}$, is a
tensor field over the spatial domain $\Omega$, which takes the value of an SPD
matrix at each point. The metric provides a concept of mesh optimality which
can in turn be used to drive a mesh adaptation process. That is, at each step
where a mesh modification is to occur, we query which of the available local
operations yields the greatest increase in quality, as measured by a quality
functional related to the metric. (See [24, Equation (5.10)] for example.) For
the 2D metric-based mesh adaptation (or mesh optimisation) tool used in this
work, _Mmg2d_ [35], the available operations are node insertion, node
deletion, edge swapping and node movement. A great advantage of the metric-
based approach is that it enables control of the shape and orientation of
elements, as well as their size, leading to the ability to generate highly
anisotropic elements. Riemannian metrics can also be viewed as continuous
analogues for meshes (see [36, 37] for details) – a duality with many useful
insights for metric-based mesh adaptation. We refer to the above citations and
[38] for further details on the metric-based framework.

##### Goal-Oriented Metrics

It is the metric which drives the mesh adaptation process, so it must be
carefully selected. In this work, we use an anisotropic formulation based on
the goal-oriented error indicators discussed in Subsection 2.2. This is taken
from [39] and is built upon the scaling,

$m(K):=\frac{|K|}{|\widetilde{K}|}\frac{\mathcal{C}\>\mathcal{E}_{K}^{\frac{1}{\alpha+1}}}{\sum_{K\in\mathcal{H}_{h}}\mathcal{E}_{K}^{\frac{1}{\alpha+1}}},$
(20)

where $|K|$ is the area of element $K$ and $\alpha\geq 1$ is a parameter
controlling the extent to which the adapted mesh captures multiple scales. The
parameter value $\alpha=1$ was found to be effective in [3] and so is used
here. Metric complexity is the continuous analogue of the (discrete) mesh
vertex count [36], so the target metric complexity, which we denote by
$\mathcal{C}$, effectively controls the number of vertices in the adapted
mesh. It was proved in [40] that the scaling (20) solves the optimisation
problem of minimising interpolation error, subject to attaining the target
metric complexity.

An element-wise anisotropic metric can be obtained by applying the scaling
(20) to a matrix that encodes the element’s shape and orientation. In this
work, we choose the matrix to be a combination of Hessians. Each Hessian is
recovered from a scalar component of the forward solution using Clément
interpolation [41] (see A for details). The combination is formed by taking
the entry-wise mean. Each Hessian is already symmetric, meaning their mean is,
too. Therefore, like (16), it admits an orthogonal decomposition,

$\underline{\bf H}_{K}=\underline{\bf
V}_{K}\underline{\boldsymbol{\Lambda}}_{K}\underline{\bf V}_{K}^{T}.$ (21)

Assuming w.l.o.g. that the eigenvalues
$\underline{\boldsymbol{\Lambda}}_{K}=\mathrm{diag}(\lambda_{1},\lambda_{2})$
are ordered such that they are ascending in magnitude, we consider a
stretching factor, $s_{K}:=\sqrt{|\lambda_{2}/\lambda_{1}|}$ (much like the
aspect ratio in Subsection 3.1). The anisotropic metric is formed as

$\underline{\bf M}_{K}^{\mathrm{anisotropic}}:=m(K)\>\underline{\bf
V}_{K}\>\mathrm{diag}\left(s_{K},\frac{1}{s_{K}}\right)\>\underline{\bf
V}_{K}^{T}.$ (22)

Similarly as in Subsection 3.1, (22) neatly decomposes the metric into
components affecting size (the scaling $m(K)$), orientation (the eigenvector
matrix $\underline{\bf V}_{K}$) and shape (the stretching factor $s_{K}$).

##### Steady-State Mesh Adaptation Algorithm

The anisotropic metric construction strategy is described above. Building upon
this, we outline the mesh adaptation algorithm in the following. Figure 5
provides a workflow representation.

Solve forward on $\mathcal{H}_{h}^{i}$; evaluate QoI $\mathcal{H}_{h}^{0}$,
$i:=0$ QoI converged? Solve adjoint on $\mathcal{H}_{h}^{i}$ Terminate
Evaluate error indicators on $\mathcal{H}_{h}^{i}$ Increment $i:=i+1$ Mesh
converged? Construct metric of complexity $\mathcal{C}^{i}$ on
$\mathcal{H}_{h}^{i}$ Adapt mesh as $\mathcal{H}_{h}^{i+1}$ yesnoyesno Figure
5:  Diagram illustrating the fixed point iteration loop used for the mesh
adaptation process applied to steady-state PDE problems. The error indication
step is highlighted in a brighter blue than the other computation steps
because it is the one that we seek to accelerate.

We begin with an initial mesh, $\mathcal{H}_{h}^{0}$, where the superscript
indicates the fixed point iteration. Then, in iteration $i\in\mathbb{N}$, we
begin by solving the forward problem in the base finite element space on
$\mathcal{H}_{h}^{i}$ and thereby evaluate the QoI. If it has not converged
since the previous iteration, we progress to solve the adjoint equation in the
base space. The next step is to indicate errors, which is the only place where
the data-driven approach differs from the standard goal-oriented approach.
Given error indicators, we then construct an anisotropic metric. For steady-
state problems, it is often beneficial to ‘spin up’ the target metric
complexity in this step, especially from the perspective of computational
efficiency. As such, $\mathcal{C}^{i}$ denotes the target complexity at step
$i$, and is typically prescribed to plateau at the desired value after three
iterations (see [42] for details). Finally, we adapt (or optimise) the mesh
based on the metric and check for convergence of its element count. If the
overall change in element count is sufficiently small then the algorithm is
terminated. Typically, the process is also terminated (with failure) at this
point if a maximum iteration count is reached. Otherwise, the process
continues. We find that three or more iterations are typically required in
order for anisotropic mesh features to manifest themselves properly, so opt to
turn off the convergence checks for the first three iterations, effectively
imposing a ‘minimum iteration count’.

##  4 Numerical Experiments

###  4.1 Problem Setup

Goal-oriented error estimation and mesh adaptation were previously applied to
steady-state shallow water modelling in [3]. That work focused on two test
cases modelling flow around two turbines positioned in a compass-aligned
channel with uniform inflow velocity from west to east. In one case, the
second turbine was positioned directly behind the first to maximise the
effects of wake interaction. In the other case, the turbines were offset to
the north and south, so that the second turbine was not fully in the wake of
the first. This offsetting gave rise to a heightened power output from the
turbines, in agreement with the literature [43]. By applying isotropic goal-
oriented mesh adaptation, it was found that the power output could be more
accurately captured than with (global) uniform refinement of the initial mesh.
Yet further improvements were found using anisotropic goal-oriented mesh
adaptation.

The aim of the numerical experiments in this section is to replicate the
improved convergence properties of the goal-oriented methods exhibited in [3]
using a data-driven variant, except at significantly reduced computational
cost. In order to achieve this, it is necessary to supply ample training data
for the neural network. Due to the design decision to work on an element-by-
element basis, a significant amount of data is available in a single test
case. The amount of test data is increased by considering a variety of
different test cases. The hope is that this will also improve the
generalisability of the network. Yet more data is obtained by collecting
features on adapted meshes, as well as the initial mesh. For each test case,
we consider the initial mesh and two adapted meshes, implying a total of three
runs per test case.

The domain used in each setup is the same: the rectangle
$\Omega=[1200\,\mathrm{m},500\,\mathrm{m}]$. This is sufficient for the
purposes of modelling steady-state flow in other domains, provided that their
geometries do not introduce flow physics that the network did not “see” during
training. The reason for this is the fact that we work on an element-by-
element basis, meaning that global information (such as domain geometry) is
not available to the network. However, in practice, almost any alteration of
the domain geometry will give rise to different flow physics from those
present in a straight channel and so we are limited to rectangular domains for
now. See Section 5 for further discussion on this.

Whilst the domains used for training are identical, the configuration of
turbines within them are chosen randomly. This reduces the amount of developer
effort required to set up the training data and removes any human bias that
might be introduced by choosing particular configurations. It also means that
the training data can be easily scaled up to an arbitrary number of tidal farm
configurations. For each test case, a random number of turbines between one
and eight (inclusive) is sampled from an integer uniform distribution, i.e.
$U\\{1,8\\}$. For each such turbine, $x$\- and $y$-coordinates are sampled
from the continuous uniform distributions $U(50,1150)$ and $U(50,450)$,
respectively (i.e. at least $50\,\mathrm{m}$ from any of the domain boundaries
in the $\ell_{\infty}$ norm). A proposed turbine location is accepted if it is
at least $50\,\mathrm{m}$ from any other turbines and rejected otherwise. The
selection procedure iterates until all turbines have been positioned in the
farm. Figure 6 shows the first 16 randomly generated setups, out of 100 in
total. Given that we consider the initial mesh and an additional two adapted
meshes per test case, feature data is extracted from a total of 300 meshes.
This gives rise to 711,307 elements across all of our experiments. The feature
and output data associated with these elements are divided into subsets for
training and validation in the ratio 7:3.

Figure 6:  The first 16 setups used for training and validation (i.e. ‘seen’
scenarios). Each plot shows the associated domain boundary and turbine
footprints, along with annotated viscosity coefficient, bathymetry and inflow
velocity.  Figure 7:  The two ‘aligned’ and ‘offset’ setups used for testing
(i.e. ‘unseen’ scenarios). Each plot shows the associated domain boundary and
turbine footprints, along with annotated viscosity coefficient, bathymetry and
inflow velocity.

It should be noted that positioning many turbines in a domain of constrained
size will typically mean that goal-oriented mesh adaptation is less effective
than it would be in the presence of fewer turbines or in the case of modelling
a larger geographical region. This is because there will likely be a higher
proportion of the domain where high mesh resolution is required for an
accurate QoI evaluation. Consequently, the improvement over uniform refinement
may be lessened for the same DoF count. However, these factors do not mean
that training the network on such scenarios is not of value. In fact, the
local approach means that more turbine interaction dynamics can be learnt than
if only small turbine farms were considered.

In order to improve the generalisability of the network, we sample the
viscosity, bathymetry and inflow speed parameters from uniform distributions,
where each parameter is taken to be constant in space. By doing so, we aim to
increase the span of the explored portion of parameter space. Whilst viscosity
is a physical property of fluids, it also provides a parameter which can be
artificially increased in order to stabilise a numerical model. As such, the
viscosity coefficient value is kept as $\mathcal{O}(10^{-1})$ across the test
cases because this helps to maintain stability. In detail, $b\sim U(20,100)$,
$u_{\mathrm{in}}\sim U(0.5,6)$ and $\nu\sim U(0.1,1)$. The values are
annotated in Figure 6. Another means of generalisation could be to choose
spatially varying viscosity, bathymetry and inflow parameters. Note that,
whilst inflow speed is varied across the training data, it is not included as
an input parameter for the network. This is because we would like the network
to be able to generalise to real scenarios where the inflow is so far away
that its value is not of immediate importance, or where different inflow
conditions are used.

The test set is comprised of the two aligned and offset configurations from
[3], which we seek to improve upon. They are shown in Figure 7, along with the
corresponding physical parameters. Using these domain geometries, we generate
quasi-uniform initial meshes using the ‘pack’ algorithm within _gmsh_ , as
shown in Figure 8. In particular, we generate meshes that are as close as
possible to a uniform mesh of right-angled triangles with edge length equal to
the turbine diameter, under the constraint that the turbine footprints are to
be explicitly meshed. On these initial meshes, each turbine is meshed using
two right-angled triangles. The same quasi-uniform meshing strategy is used to
generate initial meshes for the 100 training cases.

(a) Aligned case (b) Offset case

Figure 8:  Quasi-uniform initial meshes used for the ‘aligned’ and ‘offset’
used for testing (i.e. ‘unseen’ scenarios).

###  4.2 Testing

##### Error Indicators

Before making a formal assessment of the data-driven method’s performance, we
first consider the error indicator fields that it generates, in order to
confirm that they have similar structures to those generated by goal-oriented
mesh adaptation.

(a) Aligned, goal-oriented (b) Offset, goal-oriented

(c) Aligned, data-driven (d) Offset, data-driven

Figure 9:  Error indicator fields as computed on the ‘aligned’ and ‘offset’
initial meshes shown in Figure 8. Error indicators computed using standard
goal-oriented mesh adaptation are shown at the top, whilst the corresponding
data-driven indicators are shown at the bottom.

The top row of plots in Figure 9 show the error indicator fields according to
the standard goal-oriented approach. Contributions to the QoI error are
estimated to be largest in magnitude at and upstream of the turbines. This
makes sense because the hydrodynamics in these regions significantly impact
upon the QoI (power output) value. Errors are also indicated along the domain
boundary in the left-hand part of the domain. These occur because the flow
model, Thetis, does not impose boundary conditions exactly; it weakly enforces
them, as is common practice in DG codes. Error contributions are estimated to
be smaller outside of the regions mentioned above, especially downstream of
the turbine farm. This makes sense because we have an advection-dominated
problem, so the downstream hydrodynamics do not impact upon the QoI.

Consider now the bottom row of plots, which are generated using the data-
driven approach. Many of the important features appear to be captured well:
the highest error contributions to the QoI are at and upstream of the
turbines, errors are indicated along the boundary in the left-hand part of the
domain and the lowest errors are indicated downstream of the turbine farm. The
indicator fields are not identical, however. The most notable difference is
that the magnitudes are often not that well captured, especially when the
error estimate values are low. In addition, the error indicator fields tend to
be ‘smoother’ than predicted by the standard goal-oriented approach.

##### Adapted Meshes

The scaling of the data-driven error indicator fields does not need to be
perfect to still produce high quality adaptive meshes on which the QoI can be
accurately approximated. Figure 10 shows the meshes generated on convergence
of the resulting fixed point iteration loops.

(a) Aligned, goal-oriented (82,904 DoFs, 5 iterations) (b) Offset, goal-
oriented (53,617 DoFs, 4 iterations) (c) Aligned, data-driven (74,777 DoFs, 5
iterations) (d) Offset, data-driven (50,741 DoFs, 4 iterations)

Figure 10:  Converged adapted meshes for the ‘aligned’ and ‘offset’ test
cases. Meshes generated using standard goal-oriented mesh adaptation are shown
at the top, whilst the corresponding data-driven meshes are shown at the
bottom. The number of fixed point iterations required for convergence are
provided.

The top row of Figure 10 shows the adapted meshes due to the standard
approach. The influence of the error indicators is clearly visible: the
highest level of mesh resolution is deployed at and upstream of the turbines,
there is also increased resolution along the upstream domain boundaries, and
coarse resolution is used elsewhere, especially downstream. Recall that we use
a metric formulation based on the Hessian of the forward solution. The
influence of the Hessian can also be seen in the anisotropic elements used
between the turbines and downstream of the turbines. These are the regions
where the curvature of the forward solution take the largest values.

The bottom row of Figure 10 shows the corresponding data-driven meshes. The
similarity between the goal-oriented and data-driven adaptive methods is
striking. All of the meshing patterns observed above are also apparent for
these meshes. The main difference is that the data-driven approach leads to
adapted meshes with high anisotropy in the downstream region. This occurs
because the error is not estimated to be as low as with the standard goal-
oriented indicator, meaning that the Hessian’s high curvature in that region
is not ‘cancelled out’ to the same extent.

##### Convergence Analysis

Figure 11 shows QoI error convergence curves associated with various ‘unseen’
test cases. The top row of plots are the ones being considered presently and
the bottom row are discussed later in Subsection 4.3.

In each subfigure, data points are plotted for runs of the different mesh
adaptation methods with increasing overall DoF count. We are able to request
meshes with more DoFs by increasing the target metric complexity in the
scaling (20). This is independent of the error estimation step and so does not
require retraining of the network. Benchmark QoI values are obtained on high
resolution fixed meshes generated by uniformly refining the initial meshes
four times. These values agree with the corresponding QoI values due to three
uniform refinements to at least three significant figures. QoI ‘errors’ may
then be calculated in a relative sense by comparison against these benchmark
values. Lines of best fit are included, to give a sense of the convergence
behaviour.

(a) Aligned case (b) Offset case (c) Reversed (aligned) case (d) Trench
(offset) case

Figure 11:  QoI error convergence curves against DoF count, for ‘unseen’ test
case configurations.

In the case of uniform refinement (blue curve), Subfigures 11(a) and 11(b)
show that higher overall resolution leads to reduced QoI error in general and
that the convergence rate is approximately constant. Applying anisotropic
goal-oriented mesh adaptation (shown in yellow), we observe roughly the same
rate of convergence as with uniform refinement. However, for a given DoF count
in the mixed finite element space, the resulting QoI error values are smaller
by two orders of magnitude. The data-driven version (green) is demonstrated to
be capable of matching the error convergence rate and values in both aligned
and offset cases. It is noteworthy that this is achieved, despite the fact
that the data-driven error indicators do not exactly match those due to the
‘gold standard’ method. We reiterate that the error estimation does not need
to be ‘perfect’ in order to arrive at an effective mesh adaptation routine.

##### Computational Cost at Runtime

So far we have assessed the ability of the network to emulate the goal-
oriented error estimation process. Next we assess its ability to _accelerate_
it, i.e. quantify the associated reduction in computational cost. Figure 12
shows outputs from the same simulations as in Figure 11, except with CPU time
being the independent variable, rather than DoF count.

(a) Aligned case (b) Offset case

Figure 12:  QoI error convergence curves against CPU time, for ‘unseen’ test
case configurations. The numbers of fixed point iterations used to generate
each data point are annotated alongside the data points.

Plotted against CPU time, the rate of decrease is again roughly constant for
the uniform refinement approach. This is not true for the goal-oriented and
data-driven approaches. One reason is that more fixed point iterations are
typically required when the target metric complexity is low, since the
adaptation routine is allowed fewer DoFs than might be preferred. Another
reason is that there are various setup costs associated with mesh adaptation
that are negligible at scale, but often dominate the computational cost of
small problems. As such, the convergence curves are far from monotonic, at
least for the target complexities considered. If errors of 2% or more are
admissible, uniform refinement is the fastest approach. As the admissible
error value is reduced, however, the computational benefit of using adaptive
methods increases. For most error levels or CPU times, the data-driven
approach is shown to come with an additional computational benefit over the
standard goal-oriented approach. That is, the data-driven approach is able to
attain a given error level using fewer DoFs in general.

To investigate the computational cost further, we consider a ‘profiling’
experiment, where we take simulations that give QoI estimates of comparable
accuracy and then contrast the runtimes of the different components of the
associated algorithms. From the aligned test case, we take the data points
generated using a target metric complexity of 3,200, both of which give rise
to errors close to 0.039%. For the offset test case and the same target
complexity, both approaches give rise to errors close to 0.037%. Figure 13
shows the timings for each component. Applying the data-driven approach, the
overall runtime can be approximately halved, compared with standard goal-
oriented mesh adaptation. We do not provide a breakdown of a fixed (uniform)
mesh run that yields a similar QoI accuracy because, as may be expected, such
runs consist solely of a ‘forward solve’.

(a) Aligned (5 iterations, QoI error 0.039%) (b) Offset (4 iterations, QoI
error 0.037%)

Figure 13:  CPU timings for each component of the mesh adaptation algorithms,
for instances of the ‘aligned’ and ‘offset’ test cases. In both cases, a
target complexity of 3,200 is used.

The ‘forward solve’ component corresponds to the solution of the nonlinear
system associated with the discretised shallow water equations and the
subsequent evaluation of the QoI. The two adaptive runs converge in the same
number of fixed point iterations and use the same target complexity, so the
forward solve runtimes are comparable across those bars; similarly for the
‘adjoint solve’ component. The adjoint solves (in the base space/on the base
mesh) take less time than the forward solves because adjoint equations are
always linear, whereas the (forward) shallow water equations are nonlinear and
their numerical solution requires multiple Newton iterations.

For the standard goal-oriented approach, the ‘error estimation’ component
includes the construction of the enriched space, the associated adjoint solve
in that enriched space, assembly of the error indicator in the enriched space
and finally its transfer down to the base space. This accounts for over 56% of
the total runtime, emphasising the motivation for accelerating the error
estimation process. For the data-driven approach, the ‘error estimation’
component includes the extraction of features for the input data and the
subsequent propagation through the network. Together, these account for less
than 11% of the data-driven runtime (less than 5% of the original runtime),
implying a significant cost reduction. In terms of the error estimation step
alone, the data-driven approach yields a 90% saving in computational cost. The
majority of the remaining cost is due to feature extraction, since this
involves some computations related to mesh element geometry, evaluating
coarse-grained approximations to the error estimator and determining gradient
and Hessian values at cell centroids.

The ‘metric construction’ and ‘mesh adaptation’ steps record the relatively
small amount of time spent constructing Riemannian metrics from the error
indicator data and for the mesh optmiser Mmg to use this to generate an
adapted mesh, respectively. Most of the cost associated with the metric
computation is due to Hessian recovery from the forward solution. This is done
in both goal-oriented and data-driven approaches, so the cost is similar in
each case.

##### Computational Cost of Training

In the above, we show that the runtime cost of goal-oriented mesh adaptation
can be reduced significantly by application of neural networks. However,
runtime improvements in efficiency are not the only important factor – the
cost of training should also be considered. The two major costs are generating
the feature and target data and then training the network on these inputs and
outputs. These are shown in Table 1.

Workflow component | Count | Total time taken
---|---|---
Mesh generation | 100 scenarios | 4.7 minutes
Feature generation | 100 scenarios | 95.5 minutes
Network training | 2,000 epochs | 100.8 minutes
| Overall | 201.0 minutes
Table 1:  A table showing the CPU timings for each component of the neural
network training. Temporal values are rounded to the nearest minute.

Given that we train over 100 different scenarios, the cost per scenario is
approximately one minute. For a given amount of input and output data, the
network training cost is then effectively a function of the number of epochs,
taken here to be 2,000. For the experiments presented in this work, the cost
is approximately six seconds per epoch. We argue that the costs are small
prices to pay to double the runtime efficiency.

###  4.3 Generalisation

So far, we have tested the data-driven mesh adaptation approach on two
‘unseen’ cases. For both of these cases, the physics coefficients are well
within the bounds of the parameter spaces that they are randomly sampled from,
i.e. $\nu=0.5\in[0.1,1]$, $b=40\in[10,100]$ and
$u_{\mathrm{inflow}}=5\in[0.5,6]$. In order to test the approach more
thoroughly, we consider some additional test cases whose problem
specifications have key differences from the training cases.

####  4.3.1 Reversed Flow Direction

All cases considered thus far have flow in the positive $x$-direction.
Consider an identical setup to the ‘aligned’ test case, but with the inflow
and outflow boundaries switched and the sign of the inflow velocity reversed:
$u_{\mathrm{inflow}}=-5\,\mathrm{m\,s}^{-1}$.

The top row of Figure 14 shows the converged adaptive meshes for both the
standard goal-oriented and data-driven approaches. These plots exhibit the
same meshing patterns as observed in Subfigures 10(a) and 10(c), but reflected
from left to right. Reflected versions of the plots from those subfigures are
shown in the bottom row, for comparison. The standard goal-oriented meshes
take a very similar form; it is difficult to pick out any notable differences.
The data-driven meshes are also somewhat comparable, although the band of
resolution upstream (i.e. to the right) of the turbines is less apparent in
the ‘reversed’ case. In addition, errors appear to be over-estimated in the
downstream region (i.e. to the left), making the Hessian’s contribution to the
metric even more prominent.

(a) Reversed, goal-oriented (84,281 DoFs) (b) Reversed, data-driven (75,131
DoFs)

(c) Reflected, goal-oriented (82,904 DoFs)

(d) Reflected, data-driven (74,777 DoFs)

Figure 14:  Top row: Converged adapted meshes for a reversed flow version of
the ‘aligned’ test case. The mesh generated using standard goal-oriented mesh
adaptation is shown on the left, whilst the corresponding data-driven mesh is
shown on the right. Both ‘reversed’ cases are run for four fixed point
iterations, for consistency with the original runs. Bottom row: Reflected
versions of the corresponding original results from Subfigures 10(a) and
10(c).

Recall the convergence plot for the ‘aligned’ test case in Subfigure 11(a).
The equivalent plot for the ‘reversed’ test case is shown in Subfigure 11(c).
Whilst the two plots are not identical, they are in good agreement and the
improvement on uniform refinement is clear. As such, whilst the data-driven
approach deploys mesh resolution in a slightly different way in the ‘reversed’
test case, the overall accuracy is still of a high standard.

####  4.3.2 Variable Bathymetry

Next, we consider again the ‘offset’ setup but modify the physical parameters
so that they are outside of the bounds of those used in the training data.
Moreover, we use a spatially varying bathymetry field,

$b\left(\widetilde{y}\right)=160+40\,\widetilde{y}\,\left(1-\widetilde{y}\right),$
(23)

where $\widetilde{y}$ is a linearly rescaled $y$-coordinate so that
$\widetilde{y}=0$ is the bottom boundary and $\widetilde{y}=1$ is the top.
This means $b(\widetilde{y})\not\in[10,100],\>\forall\widetilde{y}$. The
bathymetry forms a parabolic valley, so we refer to this as the ‘trench’ test
case. Note that a consequence of this variable bathymetry is a more complex
flow field. The bathymetry field is represented in $\mathbb{P}0$ space, to
coincide with how features are extracted cell-wise. The other physical
parameters are also set to be well outside of the ranges of values seen
previously: $\nu=2.0\not\in[0.1,1]$; $u_{\mathrm{inflow}}=10\not\in[0.5,6]$.

Figure 15 shows the resulting goal-oriented and data-driven meshes, for which
QoI convergence is achieved after four fixed point iterations. Compared with
the original ‘offset’ test case, the relatively high viscosity value
exacerbates the spread of mesh resolution upstream of the turbines, as
expected. The data-driven adapted mesh is in good agreement with the one
generated using standard goal-oriented adaptation. The only noticeable
difference is that the standard method puts slightly more mesh resolution at
the turbines themselves.

(a) Goal-oriented (54,743 DoFs, 4 iterations) (b) Data-driven (50,792 DoFs, 4
iterations)

Figure 15:  Converged adapted meshes for the ‘offset’ setup with physical
parameters outside the bounds of the training data, including the use of a
variable bathymetry field. The mesh generated using standard goal-oriented
mesh adaptation is shown on the left, whilst the corresponding data-driven
mesh is shown on the right. The number of fixed point iterations required for
convergence are provided.

Convergence analysis for the ‘trench’ test case is shown in Subfigure 11(d),
with the original ‘offset’ results shown in Subfigure 11(b). The trend of the
uniform refinement curve is not as clear for this test case as for the
previous ones. Nonetheless, assuming again that the highest resolution uniform
refinement run gives a good approximation to the true QoI value, we find that
both the standard goal-oriented and the data-driven approaches are able to
achieve comparable accuracy using around an order of magnitude fewer DoFs.

##  5 Conclusion and Outlook

To the best of the authors’ knowledge, this paper presents the first framework
for using neural networks to accelerate specifically the error estimation step
of a goal-oriented mesh adaptation pipeline. The element-based approach taken
in this paper helps to improve the generalisability of the neural network
approach, since the global structure of the problem at hand is not of
significance.

A feed-forward fully connected neural network is trained on 100 randomly
generated tidal turbine modelling test cases, each with varying numbers of
turbines and uniformly distributed physical parameters. The network is then
tested on two established test cases, whose associated parameters fall within
the bounds of the training data. For a given level of accuracy in the QoI
approximation, the data-driven approach is found to match the standard goal-
oriented approach in the sense that it is also able to achieve that accuracy
level using an order of magnitude fewer DoFs than uniform refinement. Further,
the data-driven approach is shown to be capable of halving the total runtime
of standard goal-oriented mesh adaptation (including the forward flow solve
step), by completely avoiding the computationally expensive enrichment step.
The remaining runtime for the data-driven error estimation step (comprised of
feature extraction and network propagation) is 90% lower than for the standard
approach.

Regarding inputs to the neural network, we propose here a feature set that
takes into account the local element geometry, physical parameters and the
variation of the forward and adjoint solution fields over. Derivative and
curvature information are included as well as values, in order to fully
capture the solution fields in each element. It is plausible that the network
would further benefit from having access to information from the surrounding
elements, too, giving a ‘patch-wise’ approach. Perhaps such an approach would
retain generalisability, whilst improving its estimation capability by
accessing more information on error distribution. Note in particular that
patch-wise information is especially relevant to DG discretisations, where
field values tend to differ on each side of inter-element boundaries.

In addition to testing the data-driven mesh adaptation approach on existing
test cases for goal-oriented mesh adaptation that were not used for training,
but which have parameters within the ranges “seen” during training, we perform
generalisation experiments, designed to test the extent to which “unseen”
parameter values and spatially varying bathymetry fields can be handled.
Differences in meshing patterns are observed, but the QoI error convergence
properties are still desirable. Even with these extensions, however, the
applicability of the data-driven approach (as presented) remains limited
within the scope of all scientifically relevant shallow water modelling
problems. In terms of steady-state problems, it would be advantageous to be
able to apply the data-driven approach to scenarios with different boundary
condition formulations, spatially varying viscosity fields and more complex
geometries, including non-convex domains with “holes” or “islands”.
Preliminary investigations undergone as part of this work found that further
training data must be considered if the data-driven approach is to be extended
significantly in any of these directions.

The shallow water equations are often solved in their time-dependent
formulation rather than to steady-state. As such, a major piece of future work
is to train a network to perform goal-oriented error estimation for time-
dependent problems, such as the tidally reversing array example considered in
[23]. Conceptually, this would be very similar to the present work, mapping
feature data to error indicators. Perhaps the training could be done on a
timestep-by-timestep basis, in the same way that we work on an element-by-
element basis in space. The main difference would be that the error indicator
in question is based on the residual of a time-dependent problem, as defined
by the time integrator being used. In particular, it would include an
additional time derivative term.

Building upon the advancements proposed above, the ultimate goal of this work
is to use neural networks to accelerate goal-oriented error estimation and
mesh adaptation in the context of realistic, multi-scale coastal ocean
modelling problems, and indeed a wider range of fluid dynamical and PDE based
problems that can benefit from the use of goal-oriented mesh adaptation. In
such cases, the use of variable resolution is imperative to capture the
multiple scales without extortionately high mesh complexities. Goal-oriented
mesh adaptation provides a way to reduce the mesh complexity, but a data-
driven approach could additionally reduce runtimes.

##  Appendix A Derivative Recovery

The Riemannian metric construction used in this work makes use of the Hessians
of the shallow water solution component fields. Since the velocity components
are approximated in $\mathbb{P}1_{DG}$ space, they are not twice continuously
differentiable. As such, we require the facility to recover derivative
information from a lower order field. We do this using the Clément interpolant
[41], which provides a means of transferring data between $\mathbb{P}0$ spaces
and $\mathbb{P}1$ spaces via weighted averaging over neighbouring cells.

Suppose $u_{h}$ is a (scalar-valued) $\mathbb{P}1_{DG}$ field defined on some
mesh, which approximates a smooth field $u$. Its finite element gradient,
$\nabla u_{h}$, is then a (vector-valued) $\mathbb{P}0$ field. Applying the
Clément interpolant, we obtain a (vector-valued)
$\mathbb{P}1\subseteq\mathbb{P}1_{DG}$ field, $\mathbf{g}_{h}$, which
approximates the gradient $\nabla u$. Applying the same logic again, the
finite element gradient $\nabla\mathbf{g}_{h}$ is a (matrix-valued)
$\mathbb{P}0$ field and its Clément interpolant $\underline{\mathbf{H}_{h}}$
is a (matrix-valued) $\mathbb{P}1$ field, which approximates the Hessian
$\nabla\nabla u$.

Given that velocity is vector-valued, recovering its Hessian essentially means
applying the above procedure to each of its components.

##  Appendix B Software

The shallow water equations are solved using the _Thetis_ coastal ocean model
[26], which is built on top of the _Firedrake_ [44] finite element library.
Firedrake is a Python-based package, which automatically generates C code and
uses _PETSc_ [45, 46] to provide its unstructured mesh concept and to solve
the linear and nonlinear systems which underpin finite element problems.
Metric-based mesh adaptation is achieved using _Mmg_ [35], the coupling of
which with PETSc’s mesh representation was documented in [47].

The software versions used to generate the numerical experiments in this paper
are archived on Zenodo at [48] and [49]. Simulation code is archived at [50].

### Acknowledgements

The authors would like to thank members of the Applied Modelling and
Computation Group (AMCG) at Imperial College London for the many interesting
discussions regarding this work, in particular Stephan C. Kramer and James
Tlhomole.

This work was funded in part under the embedded CSE programme of the ARCHER2
UK National Supercomputing Service (http://www.archer2.ac.uk). Support is also
acknowledged from EPSRC under grant EP/R029423/1.

### CRediT Statement

Joseph G. Wallwork: Conceptualisation, Methodology, Formal Analysis,
Investigation, Writing - Original Draft, Writing – Review & Editing, Software,
Visualisation. Jingyi Lu: Methodology, Investigation, Writing – Review &
Editing. Mingrui Zhang: Methodology, Writing - Review & Editing. Matthew D.
Piggott: Methodology, Writing – Review & Editing, Supervision, Funding
Acquisition.

## References

  * [1] M. D. Piggott, S. C. Kramer, S. W. Funke, D. M. Culley, A. Angeloudis, Optimization of marine renewable energy systems, in: T. M. Letcher (Ed.), Comprehensive Renewable Energy (Second Edition), second edition Edition, Elsevier, Oxford, 2022, pp. 176–220.  doi:10.1016/B978-0-12-819727-1.00179-5.
  * [2] R. Becker, R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta numerica 10 (2001) 1–102.  doi:10.1017/S0962492901000010.
  * [3] J. G. Wallwork, N. Barral, S. C. Kramer, D. A. Ham, M. D. Piggott, Goal-oriented error estimation and mesh adaptation for shallow water modelling, SN Applied Sciences 2 (2020) 1–11.  doi:10.1007/s42452-020-2745-9.
  * [4] V. Dolejší, F. Roskovec, Goal-oriented error estimates including algebraic errors in discontinuous Galerkin discretizations of linear boundary value problems, Applications of Mathematics 62 (6) (2017) 579–605.  doi:10.21136/AM.2017.0173-17.
  * [5] O. C. Zienkiewicz, J. Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, International journal for numerical methods in engineering 24 (2) (1987) 337–357.  doi:10.1002/nme.1620240206.
  * [6] K. Huang, M. Krügener, A. Brown, F. Menhorn, H.-J. Bungartz, D. Hartmann, Machine learning-based optimal mesh generation in computational fluid dynamics, arXiv preprint arXiv:2102.12923 (2021).  doi:=10.48550/arXiv.2102.12923.
  * [7] I. Brevis, I. Muga, K. G. van der Zee, A machine-learning minimal-residual (ML-MRes) framework for goal-oriented finite element discretizations, Computers & Mathematics with Applications 95 (2021) 186–199.  doi:10.1016/j.camwa.2020.08.012.
  * [8] W. Tingfan, L. Xuejun, A. Wei, Z. Huang, L. Hongqiang, A mesh optimization method using machine learning technique and variational mesh adaptation, Chinese Journal of Aeronautics (2021).  doi:10.1016/j.cja.2021.05.018.
  * [9] J. Roth, M. Schröder, T. Wick, Neural network guided adjoint computations in dual weighted residual error estimation, SN Applied Sciences 4 (2) (2022) 1–17.  doi:10.1007/s42452-022-04938-9.
  * [10] A. Khoei, H. Moslemi, M. Seddighian, An efficient stress recovery technique in adaptive finite element method using artificial neural network, Engineering Fracture Mechanics 237 (2020) 107231.  doi:10.1016/j.engfracmech.2020.107231.
  * [11] A. Chakraborty, T. Wick, X. Zhuang, T. Rabczuk, Multigoal-oriented dual-weighted-residual error estimation using deep neural networks, arXiv e-prints (2021) arXiv–2112doi:10.48550/arXiv.2112.11360.
  * [12] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, P. W. Battaglia, Learning mesh-based simulation with graph networks, arXiv preprint arXiv:2010.03409 (2020).  doi:10.48550/arXiv.2010.03409.
  * [13] K. J. Fidkowski, G. Chen, Metric-based, goal-oriented mesh adaptation using machine learning, Journal of Computational Physics 426 (2021) 109957.  doi:10.1016/j.jcp.2020.109957.
  * [14] G. Chen, K. J. Fidkowski, Output-based adaptive aerodynamic simulations using convolutional neural networks, Computers & Fluids 223 (2021) 104947.  doi:10.1016/j.compfluid.2021.104947.
  * [15] J. Yang, T. Dzanic, B. K. Petersen, J. Kudo, K. Mittal, V. Z. Tomov, J. Camier, T. Zhao, H. Zha, T. V. Kolev, R. W. Anderson, D. Faissol, Reinforcement learning for adaptive mesh refinement, arXiv preprint arXiv:2103.01342 (2021).  doi:10.48550/arXiv.2103.01342.
  * [16] W. Song, M. Zhang, J. G. Wallwork, J. Gao, J. Wang, Z. Tian, F. Sun, M. D. Piggott, J. Chen, Z. Shi, X. Chen, M2N: Mesh movement networks for PDE solvers, arXiv preprint (2022).  arXiv:2204.11188, doi:10.48550/arXiv.2204.11188.
  * [17] M. Yano, D. L. Darmofal, An optimization-based framework for anisotropic simplex mesh adaptation, Journal of Computational Physics 231 (22) (2012) 7626–7649.  doi:10.1016/j.jcp.2012.06.040.
  * [18] W. Huang, L. Kamenski, A geometric discretization and a simple implementation for variational mesh generation and adaptation, Journal of Computational Physics 301 (2015) 322–337.  doi:10.1016/j.jcp.2015.08.032.
  * [19] C. J. Budd, W. Huang, R. D. Russell, Adaptivity with moving grids, Acta Numerica 18 (2009) 111–241.  doi:10.1017/S0962492906400015.
  * [20] F. d. A. Belbute-Peres, T. Economon, Z. Kolter, Combining differentiable PDE solvers and graph neural networks for fluid flow prediction, in: International Conference on Machine Learning, PMLR, 2020, pp. 2402–2411.  doi:10.5555/3524938.3525162.
  * [21] Z. Zhang, Y. Wang, P. K. Jimack, H. Wang, MeshingNet: A new mesh generation method based on deep learning (2020).  arXiv:2004.07016, doi:10.1007/978-3-030-50420-5_14.
  * [22] W. Nianhua, L. Peng, C. Xinghua, Z. Laiping, D. Xiaogang, Unstructured mesh size control method based on artificial neural network, Chinese Journal of Theoretical and Applied Mechanics 53 (9) (2021) 1–10.  doi:10.6052/0459-1879-21-334.
  * [23] J. G. Wallwork, A. Angeloudis, N. Barral, L. Mackie, S. C. Kramer, M. D. Piggott, Tidal turbine array modelling using goal-oriented mesh adaptation, Journal of Ocean Engineering and Marine Energy (under review) (2022).  doi:10.31223/X5H06B.
  * [24] J. G. Wallwork, Mesh adaptation and adjoint methods for finite element coastal ocean modelling, Ph.D. thesis, Imperial College London (2021).  doi:10.25560/92820.
  * [25] S. C. Kramer, M. D. Piggott, A correction to the enhanced bottom drag parameterisation of tidal turbines, Renewable Energy 92 (2016) 385–396.  doi:Acorrectiontotheenhancedbottomdragparameterisationoftidalturbines.
  * [26] T. Kärnä, S. C. Kramer, L. Mitchell, D. A. Ham, M. D. Piggott, A. M. Baptista, Thetis coastal ocean model: discontinuous Galerkin discretization for the three-dimensional hydrostatic equations, Geoscientific Model Development 11 (11) (2018) 4359–4382.  doi:10.5194/gmd-11-4359-2018.
  * [27] C. J. Cotter, D. A. Ham, C. C. Pain, A mixed discontinuous/continuous finite element pair for shallow-water ocean modelling, Ocean Modelling 26 (1-2) (2009) 86–90.  doi:10.1016/j.ocemod.2008.09.002.
  * [28] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, J. Koster, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal on Matrix Analysis and Applications 23 (1) (2001) 15–41.  doi:10.1137/S0895479899358194.
  * [29] P. R. Amestoy, A. Guermouche, J.-Y. L’Excellent, S. Pralet, Hybrid scheduling for the parallel solution of linear systems, Parallel Computing 32 (2) (2006) 136–156.  doi:10.1016/j.parco.2005.07.004.
  * [30] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities, International journal for numerical methods in engineering 79 (11) (2009) 1309–1331.  doi:10.1002/nme.2579.
  * [31] A. LeNail, NN-SVG: Publication-ready neural network architecture schematics., Journal of Open Source Software 4 (33) (2019) 747.
URL https://alexlenail.me/NN-SVG/

  * [32] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al., Pytorch: An imperative style, high-performance deep learning library, Advances in neural information processing systems 32 (2019) 8026–8037.  doi:10.5555/3454287.3455008.
  * [33] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).  doi:10.48550/arXiv.1412.6980.
  * [34] P. L. George, F. Hecht, M. G. Vallet, Creation of internal points in Voronoi’s type method. control adaptation, Advances in engineering software and workstations 13 (5-6) (1991) 303–312.  doi:10.1016/0961-3552(91)90034-2.
  * [35] C. Dobrzynski, P. Frey, Anisotropic Delaunay mesh adaptation for unsteady simulations, in: Proceedings of the 17th international Meshing Roundtable, Springer, 2008, pp. 177–194.  doi:10.1007/978-3-540-87921-3_11.
  * [36] A. Loseille, F. Alauzet, Continuous mesh framework part I: well-posed continuous interpolation error, SIAM Journal on Numerical Analysis 49 (1) (2011) 38–60.  doi:10.1137/090754078.
  * [37] A. Loseille, F. Alauzet, Continuous mesh framework part II: validations and applications, SIAM Journal on Numerical Analysis 49 (1) (2011) 61–86.  doi:10.1137/10078654X.
  * [38] C. Pain, A. Umpleby, C. De Oliveira, A. Goddard, Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations, Computer Methods in Applied Mechanics and Engineering 190 (29-30) (2001) 3771–3796.  doi:10.1016/S0045-7825(00)00294-2.
  * [39] J. Carpio, J. L. Prieto, R. Bermejo, Anisotropic “goal-oriented” mesh adaptivity for elliptic problems, SIAM Journal on Scientific Computing 35 (2) (2013) A861–A885.  doi:10.1137/120874606.
  * [40] M. Braack, An adaptive finite element method for reactive-flow problems (Apr 1998).
  * [41] P. Clément, Approximation by finite element functions using local regularization, Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique 9 (R2) (1975) 77–84.
  * [42] J. G. Wallwork, N. Barral, D. A. Ham, M. D. Piggott, Goal-oriented error estimation and mesh adaptation for tracer transport modelling, Computer-Aided Design (2021) 103187doi:10.1016/j.cad.2021.103187.
  * [43] S. Draper, T. Nishino, Centred and staggered arrangements of tidal turbines, Journal of Fluid Mechanics 739 (2014) 72–93.  doi:10.1017/jfm.2013.593.
  * [44] F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. McRae, G.-T. Bercea, G. R. Markall, P. H. Kelly, Firedrake: automating the finite element method by composing abstractions, ACM Transactions on Mathematical Software (TOMS) 43 (3) (2016) 1–27.  doi:10.1145/2998441.
  * [45] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, J. Zhang, PETSc Web page, https://petsc.org/ (2021).
  * [46] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, J. Zhang, PETSc/TAO users manual, Tech. Rep. ANL-21/39 - Revision 3.16, Argonne National Laboratory (2021).
  * [47] J. G. Wallwork, M. G. Knepley, N. Barral, M. D. Piggott, Parallel metric-based mesh adaptation in PETSc using ParMmg (2022).  arXiv:2201.02806, doi:10.48550/arXiv.2201.02806.
  * [48] Thetis coastal ocean model (June 2022).  doi:10.5281/zenodo.6719158.
  * [49] Software used in ‘E2N: Error Estimation Networks for Goal-Oriented Mesh Adaptation’ (June 2022).  doi:10.5281/zenodo.6722155.
  * [50] J. G. Wallwork, Simulation code for ‘E2N: Error Estimation Networks for Goal-Oriented Mesh Adaptation’ (July 2022).  doi:10.5281/zenodo.6811761.