Goal-oriented error estimation and mesh adaptation for shallow water modelling

This study presents a novel goal-oriented error estimate for the nonlinear shallow water equations solved using a mixed discontinuous/continuous Galerkin approach. This error estimator takes account of the discontinuities in the discrete solution and is used to drive two metric-based mesh adaptation algorithms: one which yields isotropic meshes and another which yields anisotropic meshes. An implementation of these goal-oriented mesh adaptation algorithms is described, including a method for approximating the adjoint error term which arises in the error estimate. Results are presented for simulations of two model tidal farm configurations computed using the Thetis coastal ocean model (Kärnä et al. in Geosci Model Dev 11(11):4359–4382, 2018). Convergence analysis indicates that meshes resulting from the goal-oriented adaptation strategies permit accurate QoI estimation using fewer computational resources than uniform refinement.


Introduction
Coastal modelling problems are typically multi-scale, often with a strongly direction-dependent flow. As such, anisotropic mesh adaptation is an attractive prospect for providing both accuracy and reduced computational cost; by adapting the mesh such that its anisotropy is aligned with the flow, the number of degrees of freedom (DOFs) required to yield an accurate solution is reduced. The metric based approach to anisotropic mesh adaptation was first introduced in [20] and uses Riemannian metric fields to control not only the size of mesh elements, but also their shape and orientation. This approach was shown to be particularly suited to multi-scale ocean modelling in [31]. For a review of recent progress in the field of anisotropic mesh adaptation, see [1].
Further, coastal modelling problems often involve a diagnostic quantity of interest (QoI) of greater importance than the solution itself. Goal-oriented error estimation represents the error accrued in computing the QoI in terms of PDE residuals and solutions of associated adjoint equations. Used within a mesh adaptation algorithm, such estimators can yield meshes permitting accurate QoI approximation. The majority of goal-oriented error estimators are based upon the pioneering work of [9,10]. More recently, researchers have developed goal-oriented estimators which take account of discontinuous Galerkin (DG) discretisations (see [15,22]). Integration into the metric based mesh adaptation framework has also been developed (see [12,27,32]). The depth-averaged shallow water equations are often used in coastal ocean models, providing approximations 2 Goal-oriented mesh adaptation

Metric-based mesh adaptation
In this paper, Riemannian metric fields are used to drive the mesh adaptation process. These tensor fields, often referred to simply as metrics, are symmetric positive-definite (SPD) linear forms defined pointwise, which give rise to all of the geometrical quantities necessary to perform mesh adaptation. In this work, metrics are derived using the error estimates described in Sect. 2.2. For details on the metric-based approach, see [1,7,20,29,34].
The spatial domain, denoted by ⊂ ℝ n , is assumed to have piecewise smooth boundary ∶= . For a mesh H of the domain, we denote mesh elements by K ∈ H and the edge set of element K by K . In this work we restrict attention to triangular mesh elements, for simplicity. Denote the set of all edges which are not on the domain boundary (internal edges) by int .

Error estimation
The shallow water equations may be written in the 'residual form' (1) ( ) = 0, with solution living in a space of functions denoted V. Throughout this paper, we shall refer to (1) as the forward equation and to as the forward solution. The diagnostic QoI, J, is a functional which maps members of V onto the real number line. Associated with (1) and the QoI is an adjoint equation, The adjoint solution * also lives in V and conveys the propagation of sensitivities of the QoI to perturbations in the forward solution. For a finite dimensional subspace V h ⊂ V , (1)-(2) have Galerkin approximations given by where ⟨⋅, ⋅⟩ is the L 2 inner product and h and * h are finite element approximations to the forward and adjoint solutions. Typically, integration by parts is applied in constructing the weak residuals ( h , ⋅) and * ( * h , ⋅) from the inner products.
We consider the classical a posteriori goal-oriented error estimate known as the dual weighted residual, due to [9,10]. Therein, the error result is presented, where the remainder term R (2) is quadratic in the forward and adjoint errors ∶= − h and * ∶= * − * h . Element-wise dual weighted residual error indicators and a global error estimator may be derived as

Goal-oriented metrics
There are many potential ways to construct metric tensor fields using the scalar error indicator (6), one being to appropriately scale an identity matrix; this is referred to as an isotropic approach, since resulting meshes are relatively isotropic. For problems with strong directional dependence, anisotropic metrics can be beneficial, allowing control of the shape and orientation of mesh elements, as well as size.
In [37], an isotropic metric was compared with two approaches to constructing anisotropic metrics from goaloriented error estimates (based on the work of [27,32]) for advection-diffusion problems discretised using continuous finite elements. In this subsection, we utilise a different approach, developed in [12], which straightforwardly permits discontinuous discretisations.

Isotropic metric
First, consider the isotropic case. Let |K | denote the reference element volume and |K| denote the volume of an arbitrary element K ∈ H . Constructing an element-based isotropic metric in n-dimensions amounts to the scaling [12] where n is the element-wise identity metric. The desired element volume |K | is chosen to minimise interpolation error, such that the metric complexity N > 0 is smaller than some desired complexity, N max . Metric complexity can be viewed as the continuous analogue of the mesh vertex count [1,7]. The proof that |K | solves this optimisation problem is given in [11]. The parameter ≥ 1 which arises in its solution is not known a priori. However, [12] states that its influence on (7) is negligible, provided that we are sufficiently close to the optimal volume. We have found = 1 to be an effective choice in practice. It is worth remarking that |K | is mesh-dependent. This means that meshes generated using metric K are heavily influenced by the mesh upon which the metric was constructed. However, as detailed in Sect. 3.4, the mesh adaptation algorithm used in this work iteratively solves the PDE, evaluates (7) and adapts the mesh until convergence criteria are met. As such, the dependence on the initial mesh diminishes as the algorithm progresses. This is in agreement with what we have found in numerical experiments. For further discussion of the coupled mesh adaptation-PDE solution process, see [1].

Anisotropic metric
As with the approaches considered in [37], the anisotropic metric construction of [12] uses a recovered Hessian of the prognostic variables. In this instance, we compute the element-averaged Hessian K on an element K by solving an auxiliary finite element problem (for details, see [29]).
As a symmetric matrix, the Hessian has an orthogonal eigen-decomposition K = K K K T , with eigenvalue The stretching factor associated with the element-averaged Hessian is defined by We construct an element-wise metric by modifying the eigenvalues appropriately. In the two dimensional case, [12] A vertex-based metric is obtained from (9) by projection from ℙ0 to ℙ1 . In this work, the projection is applied using a Galerkin projection, which amounts to averaging the element-wise values surrounding a vertex.

Shallow water equations
In this work we consider the nonlinear shallow water equations for velocity Suppose the exact solution = ( , ) lives in a function space V with finite dimensional subspace V h . For all test functions = ( , ) ∈ V h , we have a Galerkin formulation of (10) given by where h ∈ V h is a finite element approximation to . The weak form is decomposed into advection, gravity, viscosity, drag and continuity terms. Taking all test functions, (11) forms a nonlinear system of equations which is solved using a Newton iteration. A direct solver is used in this work to solve the linear system of equations that arises in the update step of Newton's method based on the Jacobian of (11) [3,4].
For V h , we select the mixed finite element space ℙ1 DG − ℙ2 . That is, velocity is piecewise linear and discontinuous across elemental boundaries, whilst elevation is piecewise quadratic and continuous. This element pair was shown to be suitable for shallow water modelling in [13].
Clearly, functions from the velocity space are discontinuous across internal edges. Thus, internal edge integrals whose integrand involves such functions are not well defined. This necessitates the introduction of the following restriction operators, which assign a single value on each edge. For an edge ∈ int , arbitrarily label each side with + and −, giving rise to normal vectors ̂ ± . Define the restriction operators given by the average {{v}}| ∶= 1 Ignoring boundary condition implementation (for brevity), Thetis uses the discretisation, ∑ n j=1 S ij T ij for n-dimensional tensor functions and . If the elevation was chosen from a discontinuous space (which is also a valid choice in Thetis) then the gravity term would be integrated by parts and we would need Riemann solutions for and .
Lax-Friedrichs stabilisation [26] is applied in the advection term with parameter = 1 The interior penalty parameter used in the viscous term is chosen in line with [17], depending upon the polynomial degree, variation of and the minimal angle in each mesh element. For simplicity, we set as the largest element-wise value. It is important to update the penalty parameter whenever the mesh is adapted. Boundary conditions are imposed weakly, therefore contributing additional terms to (12).

Goal-oriented error estimate
As discussed in Sect. 2.2, the construction of goal-oriented error estimators from (12) typically involves integrating by parts on each element. Doing so enables us to derive from (5)-(6) error indicators of the form where ( h ) is the strong PDE residual (3), ( h ) concatenates the residuals associated with the boundary conditions and E DG K ( h , * ) contains flux terms arising from the DG discretisation.
As in [10], (13) is written in this compact form so that the structure of the error indicator is transparent. This formulation has the advantage of separating out different components of the error estimator as regards the part of the problem they relate to. The first term assesses how well the PDE is solved on element interiors, the second term assesses the extent to which the boundary conditions have been weakly imposed and the final term describes the magnitude of the flux terms contributed by the DG discretisation of the velocity space. Since boundary conditions are neglected in (12), we are primarily interested in the nature of the flux term, , which conveys the smoothness of the discrete solution. Analysis of these contributions is useful in informing model development and assessing whether mesh features arise due to approximation error, boundary conditions or discretisation choice.
In order to derive an error indicator of the form (13), we need to appropriately integrate (12) by parts and substitute test functions for the adjoint error. No integration by parts is required for the gravity or drag terms, since it is not used in their derivation from corresponding terms in (12). The application of integration by parts on a particular element K results in a term involving an integral over its edge set, K . In the DG method itself, these edge integrals must be restricted, so as to yield a unique value on each edge. However, since we seek element-based error indicators, there is no need to restrict quantities which are discontinuous across elemental boundaries.
Integrating by parts in the advection, viscosity and continuity terms yields where ̂ ∈ {̂ + ,̂ − } is the outward-pointing normal to K . The notation used in (14) can be simplified by application of three restriction identities, as follows. Consider a scalar function f, vector function and an edge shared by two adjacent elements K + , K − ∈ H whose outward-pointing normals are ̂ + and ̂ − , respectively. Then Replacing test functions with the adjoint error in (14) and applying (15), we obtain error indicators of the form (13), where the flux term is given by Here * = ( * , * ) is the exact adjoint solution, * h = ( * h , * h ) is a finite element approximation thereof and * = ( * , e * ) is the corresponding error.

Approximation of the adjoint error
Note that (13) and (16) contain the (unknown) exact adjoint solution, * . In practice, it suffices to approximate it in an enriched space, V + h ⊃ V h . We again choose V + h as ℙ1 DG − ℙ2 , but defined on a uniformly refined mesh. That is, a single, global iso-ℙ2 refinement is made, which amounts to inserting vertices wherever a quadrature node would exist in a quadratic element. The adjoint equation is then derived using a linearisation about the projected forward solution, + h h ∈ V + h . Whilst this approach is shown to be effective in Sect. 4.3, it implies an additional computational cost associated with solving the adjoint equation on a mesh with four times as many elements. In future work, a more efficient approach will be implemented, such as the one described on pp.590-593 of [15], which solves local PDEs to approximate * .

Implementation details
The anisotropic metric formulation described in Sect. 2.3 relies on the provision of an approximate Hessian. Hessian recovery techniques typically seek second derivatives for scalar fields. In this case, there are a number of options for fields to recover a Hessian from-the free surface elevation, velocity components and the fluid speed being the most obvious candidates. Individual Hessians may be combined, using a strategy such as metric averaging or metric superposition (see pp. 131-138 of [7] for details and an investigation of the differences). In this work, we superpose the free surface elevation and velocity component (16) Hessians, so that the anisotropy of each field is accounted for. The eigenvalues of the resulting metric are then modified according to (9) in order to account of the error indicator (13), giving rise to a goal-oriented metric. Mesh adaptation is performed using Pragmatic [8,34], which is a mesh optimisation library primarily performing h-like adaptive operations, with some local Laplacian smoothing. The mesh adaptation workflow used is identical to that described in Algorithm 1 of [37]. The adaptation loop is terminated if the QoI value, number of mesh elements or error estimator (13) change by less than 0.1% from one iteration to the next.
Previously published studies, such as [19], use the dolfin-adjoint package [18], which automates the solution of the discrete adjoint of a finite element problem expressed in Unified Form Language (UFL) [2]. In this work, we wish to be able to discretise the adjoint problem in a different way than the forward problem. Thus, we exploit the automatic differentiation capabilities of UFL in order to avoid an error-prone manual calculation.
Given a weak form PDE 'F == 0' with finite element solution 'q' and QoI 'J' , the adjoint solution, 'q_star' , may be computed in only a few lines of code: Writing the adjoint equation using Firedrake solve calls enables us to solve the adjoint equation in V + h , as opposed to V h . To do this, we replace q by q_∈ V + h and F by F_ in the above, where F_ is defined by prolonging the variables used in F from V h to V + h . The Firedrake installation used for all simulations documented in Sect. 4 is archived at [23,36], with all simulation code archived at [38].

Tidal turbine modelling
Marine renewable energy is an active area of research in coastal ocean modelling (for example, see [14,16,19,28,35]). In particular, tidal power presents an opportunity to generate large amounts of low-carbon electricity in coastal countries such as the UK. One major advantage of tidal power over other renewable energy sources is that tides are highly predictable, meaning that power is generated reliably. The shallow water depth in which tidal turbines are deployed has a significant impact on wake recovery and hydrodynamic blockage effects, meaning that the positioning of turbines can be very important. By modifying the shallow water equations to account for tidal turbines positioned within the domain, recent research formulated tidal array positioning as a PDE-constrained optimisation problem [19]. Solving this problem gives the configuration with maximum power, with the potential to incorporate penalties in the optimisation functional to account for financial [14] and environmental impact factors [16]. A parametrisation based approach is used, meaning that the turbines are modelled using a density function d = d( ).
Modifying (10) to account for a set of tidal turbines T amounts to choosing an appropriate drag coefficient C d . Suppose turbine T has thrust coefficient c T , area A T and footprint indicated by T . For a background drag C b = 0.0025, The binary footprint function T ∶ → {0, 1} is unity in the region where turbine T is deployed and zero elsewhere.
As such, it specifies the spatial region where the turbine drag is active. The thrust coefficent used in (17) is based on an upstream velocity, whereas is the depth-averaged velocity at the turbine. Hence, we correct the thrust coefficient using the rescaling recommended in [25]. We use the following proxy for the power output of the tidal array: This provides a QoI for goal-oriented error estimation and has units of Watts.

Problem setup
Whilst realistic tidal turbine applications are inherently time-dependent, we consider a steady-state test case to highlight the impact of the turbine position on power output. For a simple tidal farm with two turbines, T 1 and T 2 , we consider two configurations: one where T 1 is directly upstream of T 2 and another where the turbines are offset by one turbine diameter to south and north, respectively. Our numerical experiments involve applying the goal-oriented mesh adaptation approaches described in Sect. 2 to provide accurate approximations to the total tidal farm power output in each case. Coarse initial meshes which take account of the tidal turbines are generated using gmsh [21] and shown in Fig. 1a, b. Figure 1c, d show the magnitude of the velocity (i.e. the speed) given by solving the test case on meshes which have been refined three times using iso-ℙ2 refinement.
Observe that, in the aligned case, the momentum deficit at T 2 is more significant than at T 1 . This is because the wake of T 1 has not fully recovered and so the fluid speed meeting T 2 is lower than that which meets the first. That is, T 2 is fully in the wake of T 1 and thus overall experiences slower flow and consequently generates less power [cf. (18)]. In the offset case, the momentum deficit is similar at each turbine, suggesting that the presence of T 1 has less impact upon the fluid speed meeting T 2 . The result is that the total power output of the aligned array configuration is lower than that of the offset configuration, as shown in the QoI values displayed on the fourth row of Table 1. Table 1 illustrates the convergence of QoI values under iso-ℙ2 (uniform) refinement to four significant figures. Final values, J 0 = 19.7174 kW and J 1 = 23.1905 kW , present benchmark values to approximate using adaptive meshes. In agreement with expectations, the converged power output in the aligned configuration is 15% lower than in the offset case. This illustrates how the  positioning of turbines within an array can significantly impact power output.

Convergence analysis under goal-oriented mesh adaptation
Having obtained benchmark values for the power output under uniform refinement in each turbine configuration, we apply mesh adaptation in an attempt to achieve convergence to these values using fewer DOFs. Figure 2 shows meshes generated by the two adaptation strategies, each with low resolution downstream of the second turbine. This should be expected since we have an advection-dominated problem, meaning the power output of the array is independent of the downstream dynamics. We observe high mesh resolution surrounding and upstream of the turbines.
The moderate advection-dominance manifests in the anisotropy of the meshes shown in Fig. 2c, d. The maximum aspect ratio is reported to be higher in the aligned configuration, in which case almost all of the dynamics to which the QoI value is sensitive lie in a narrow horizontal band through the middle of the domain. The moderate mesh resolution near the domain boundaries in the tidal farm and upstream regions is due to residual contributions from the weakly enforced boundary conditions. Figure 3 is useful in understanding how each component of the error indicator (13) contributes towards the mesh adaptation. We observe that the PDE residual is most significant surrounding the turbines, whilst the error indicator contributions due to flux terms are generally in the upstream region. As has already been noted, the fact that we weakly impose boundary conditions means that there are significant flux term contributions near to the boundary. Figure 4 illustrates the convergence of evaluated QoI values to the benchmark values established in Table 1 under each adaptation strategy with increasing DOF count. This validates the adaptive solution strategies. The DOF count is increased by increasing the target mesh complexity in (7).
For a given number of DOFs, isotropic goal-oriented adaptation generally offers an improvement in QoI approximation accuracy over uniform refinement. Anisotropic goal-oriented adaptation is shown to yield yet more accurate power output estimates, with the improvement over uniform refinement of the initial mesh being significant. Fig. 2 Example meshes generated using goal-oriented mesh adaptation. In the aligned case, adaptation is applied to the mesh shown in Fig. 1a in order to generate a, c. In the offset case, adaptation is applied to the mesh shown in Fig. 1b

Discussion
The numerical experiments conducted in Sect. 4 validate the two goal-oriented mesh adaptation strategies considered. The experimental results show that the anisotropic approach in particular is able to achieve notably higher QoI estimation accuracy than is obtained using the same number of DOFs under uniform refinement. The isotropic approach also offers some improvement over uniform refinement for these test cases.
Consider again the goal-oriented meshes shown in Fig. 2. Observe that, in each case, the tidal farm and other well-resolved regions correspond to a small proportion of the domain, with low resolution used away from these regions. One can imagine that extending the problem size to something akin to a realistic tidal turbine modelling application would imply high resolution deployed in an even smaller proportion of the domain. Goal-oriented mesh adaptation is most effective when the QoI is insensitive to the dynamics in the majority of the domain and this becomes increasingly true as the domain size increases relative to the tidal farm. As remarked above, goal-oriented adaptation strategies deploy high resolution mostly in upstream regions for advection-dominated problems. Any combination of widening the channel, extending the downstream region and increasing the advection-dominance of the problem (as per the advection/viscosity relationship determined by the Reynolds number) would imply a leftwards shift of the goal-oriented convergence curves shown in Fig. 4, relative to the uniform refinement curve. It is likely that all three of these would be present in a realistic application, meaning that goal-oriented mesh adaptation would require relatively few DOFs to provide an accurate power output estimate.
Mixed finite element methods are becoming increasingly popular discretisation choices for shallow water problems [13]. Through the derivation shown in Sect. 3.2, we illustrate how to perform goal-oriented error estimation for mixed discretisations with discontinuous components. The derivation of this estimate follows the procedure used in [22], which considers DG discretisations for advection-diffusion-reaction problems.
The anisotropic mesh adaptation algorithm used in this work is based on that described in [12], which also focuses on advection-diffusion-reaction problems. This elementbased approach is advantageous because it straightforwardly permits the incorporation of discontinuous goal-oriented error indicators, such as arise from DG discretisations. As in the numerical experiments considered in [12], we observe convergence of the QoI to reference values computed on a high resolution mesh. In agreement with the numerical experiments considered in [22], we found that, for a given number of DOFs, the relative error in evaluating the QoI was smaller under anisotropic adaptation than under isotropic adaptation.

Conclusion
The main achievement of this paper is the formulation of a goal-oriented error estimate for the nonlinear shallow water equations solved using a mixed discontinuous/ continuous finite element method, along with the implementation of isotropic and anisotropic mesh adaptation algorithms using this estimate.

Compliance with ethical standards
Conflicts of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.