Goal-Oriented Metric-Based Mesh Adaptive Tidal Farm Modelling MARINE 2021

The modelling of a tidal array farm is an inherently multi-scale endeavour. It requires the simultaneous resolution of tidal processes across tens or hundreds of kilometres of coastal ocean (including estuaries, or even entire seas), the hydrodynamics in the neighbourhood of the farm (hundreds of metres), the wakes of individual turbines (metres, or tens of metres) and device hydrodynamics (sub-metre). As such, the construction of an accurate, computationally eﬃcient numerical model requires careful consideration of the underlying discretisation. In this paper, we apply time-dependent mesh adaptation techniques based on the Riemannian metric framework to an idealised tidal array and assess the quality of the resulting approximations. Whilst classical hierarchical mesh adaptation methods modify mesh element/cell size in order to improve resolution locally, the metric-based approach also allows for control of element shape and orientation, which can be especially advantageous for advection-dominated problems. Metrics are normalised in such a way that the resulting discretisation is multi-scale in both space and time as per Alauzet and Olivier (2010). Typically, metrics are constructed from recovered derivatives of solution ﬁelds, such as ﬂuid vorticity. Alternatively, metrics may be derived from goal-oriented error estimates, enabling accurate estimation of a diagnostic quantity of interest (QoI). In the context of tidal farm modelling, one clear QoI is the power output. Building upon the idealised steady-state test case considered in Wallwork et al. (2020), which represents turbines using a drag parametrisation in a depth-averaged shallow water model, we demonstrate here that goal-oriented mesh adaptation can be used to obtain an accurate approximation of tidal farm power output using relatively few overall degrees of freedom.


INTRODUCTION
Goal-oriented error estimation and mesh adaptation techniques have the primary objective of minimising the error accrued when evaluating a diagnostic quantity of interest (QoI). In the context of modelling a tidal farm -a collection of tidal stream turbines designed to extract energy from high velocity tidal currents -this could be the total energy output of the farm, the associated profit or a measure of its environmental impact. Given a particular choice, error estimators such as the dual weighted residual (Becker & Rannacher, 2001) can be used to deduce error indicator fields conveying how much error may be attributed to different components of the discretisation (such as particular finite elements and timesteps). These error indicators 'localise' the error contributions and can be used to drive mesh adaptation routines which construct new discretisations better representing the QoI.
In this paper, error indicator data is interpreted by a hybrid hr-adaptive mesh adaptation toolkit, Pragmatic (Barral et al., 2016), in terms of the Riemannian metric framework. This 'optimal mesh model', first introduced in George et al. (1991), defines a Riemannian metric space based on error indicator data and updates the mesh so that its elements are 'unit' when viewed in that space, i.e. they have edges of (near) uniform length, and consequently are appropriately multi-scale in physical space. For details on metric-based mesh adaptation, see Alauzet (2011a, 2011b), Pain et al. (2001), and Piggott et al. (2009).
The tidal farm modelling framework used in this paper parametrises turbines as patches of increased bottom friction in a depth-averaged 2D shallow water model. This model was chosen because it allows large, multi-scale coastal ocean modelling problems to be represented at a greatly reduced computational cost, compared with a 3D model with the same resolution in the horizontal direction. This is applicable subject to some key assumptions, such as limited vertical accelerations, well-mixed water column and relatively negligible vertical dimensions compared to the horizontal. Goal-oriented metric-based mesh adaptation was previously applied to such a model in the steady-state case in Wallwork et al. (2020). It was found to lead to a reduction in error accrued when estimating array power output, for the same number of degrees of freedom. This paper extends the previous work to the time-dependent, tidally varying case, illustrating how goal-oriented adaptation can be used to resolve the varying multi-scale hydrodynamics within and around an idealised tidal array, assuming that the 2D hydrodynamics are most important.
The organisation of this paper is as follows. Section 2 describes the depth-averaged tidal turbine modelling framework, including details on the model parameters and discretisation used. Section 3 describes a goal-oriented metric-based mesh adaptation method, which is framed around minimising the error accrued when evaluating the energy output of a tidal farm. The mesh adaptation method of Section 3 is applied to an idealised tidal turbine array test case in Section 4. To the best of the authors' knowledge, this presents the first application of goal-oriented mesh adaptation methods to a time-dependent tidal farm modelling problem. Finally, conclusions are drawn and future work is proposed in Section 5.

Shallow Water Equations
In this paper, coastal hydrodynamics are modelled using the time-dependent, non-rotational, nonlinear shallow water equations in non-conservative form: The prognostic variable tuple q := (u, η) is comprised of (depth averaged) fluid velocity, u, and free surface elevation, η, and is defined on a two-dimensional spatial domain, Ω ⊂ R 2 , and some time interval. The total water depth, H = η + b, arises from a bathymetry field, b, which is fixed in time.
Bottom friction is represented in the model using a quadratic drag parametrisation. The dimensionless drag coefficients c b and c t convey the background drag and drag due to the turbines, respectively, the latter of which is described in Subsection 2.5. Numerical experiments in this paper assume a constant background drag, c b = 0.0025.

Discretisation
Given a mesh of the domain, spatial discretisation of (1) is achieved using the P1 DG − P1 DG equal order DG discretisation due to the Thetis coastal ocean model (Kärnä et al., 2017). Thetis uses the SIPG method described in Hillewaert (2013) for the purposes of representing viscosity. Whilst Thetis can be run on triangular or quadrilateral meshes, we opt to use triangular elements in this paper, since the mesh adaptation toolkit we use assumes simplicial meshes. Thetis is based on the Firedrake Python finite element package (Rathgeber et al., 2016), which automatically generates C code and uses PETSc (Balay et al., 2019) to solve linear and nonlinear systems.
Meshes are denoted by the symbol H in this work and are interpreted as finite collections of elements, K ⊂ Ω. Time integration is performed using the Crank-Nicolson method with implicitness θ = 0.5.

Viscosity
Whilst viscosity is a physical property of a fluid, it also provides a tool to be deployed by the numerical modeller. Treating it consistently with the eddy viscosity approximation, artificially increasing viscosity smooths out turbulent effects and -when used appropriately -can improve the stability of a numerical scheme. By smoothing out turbulent effects, the flow becomes increasingly laminar. This is undesirable in applications which aim to capture vortices in the wake of tidal turbines, since such features are inherently turbulent and therefore cannot be captured without an appropriately small viscosity.
We address the above concerns by adopting a mesh-dependent viscosity -one which is sufficiently small to capture vortices in the wake of turbines, yet sufficiently large elsewhere to ensure stability of the spatial discretisation. Our specification is made via the mesh Reynolds number, Re K := h K U/ν, where U is a characteristic fluid speed and h K is the circumradius of an isotropic mesh element K ∈ H. Given a small viscosity, ν target , which we would like to impose near the turbines and a maximum tolerated mesh Reynolds number, Re max , we set This element-wise quantity is projected into (vertex-wise) P1 space, which has an additional smoothing effect that can be beneficial for the shallow water solver. For simplicity, the characteristic speed is chosen to be constant in both time and space in this work, providing an upper bound. Consequently, the mesh Reynolds number is overestimated. This is preferable to underestimation because it puts more emphasis on stability, rather than vortex capture.

Tides
Tidal forcings appear as boundary conditions for the free surface elevation on forced boundary segments Γ F ⊂ ∂Ω. (Impermeable) free-slip conditions are applied on the unforced boundaries, Γ freeslip = ∂Ω\Γ F . The use of free-slip conditions is a rather simplistic approach to modelling the dynamics on any coast other than tall cliffs. However, they are sufficient for the purposes of the numerical experiments presented in Section 4. Taken together, we have boundary conditions, where (u, η) now denote the fluid velocity and elevation values on the domain boundary. Since we are using a DG discretisation, these values are weakly enforced.

Tidal Turbine Parametrisation
Suppose we have a tidal farm, F, interpreted as a finite set of turbines, with each turbine T ∈ F interpreted as a subset of the spatial domain. Each turbine is meshed explicitly, so that T ⊂ H. The initial (un-adapted) mesh for the experimental setup presented in Section 4 includes 6 m fine resolution in order to capture the turbine footprints.
With the representation above, the turbines may be modelled using a drag parametrisation which acts within its 'footprint'. The drag coefficient c t in (1) depends on the thrust coefficient C T = C T (u) in each turbine footprint, T . The thrust coefficient relates the drag force exerted on the flow by turbine T to the area it sweeps in the vertical, A swept , via In general, the thrust coefficient is a function of fluid velocity, which may depend on a cut-in speed, for example. However, it is assumed to be constant across all turbines within a farm for the purposes of this paper. The turbine drag coefficient, c t , is calculated by scaling the thrust according to the ratio of area swept in the vertical and footprint area, A footprint , since it is the area swept by the turbine blades which determines the energy extraction: where 1 T is an indicator function which is unity within the footprint of turbine T and zero elsewhere.
It should be remarked that the thrust coefficient as described is based on an upstream velocity. In order to account for the fact that the velocity at the turbine is depth-averaged velocity, we use the thrust coefficient correction recommended in Kramer and Piggott (2016).
Given that the turbines are deployed in water of density ρ > 0 (assumed to be the constant ρ = 1030 kg m −3 herein), a proxy for the energy output of the tidal farm is given by where [T start , T end ] is the time interval we seek to accurately capture energy output over. The integrands of the time integrals on the RHS of (6) (i.e. the space integrals) provide proxies for the power output of the individual turbines.

Error Estimation
The application of goal-oriented mesh adaptation presented in this paper seeks to accurately approximate the energy output of the farm. That is, we seek to minimise δE : is the exact solution of the shallow water problem and q h = (u h , η h ) is a finite element approximation thereof.
Given the definition of energy output in (6), it is likely that mesh adaptation based on a metric derived from the fluid speed would go some way to provide an accurate energy output approximation. For example, an anisotropic metric may readily be derived from a recovered Hessian of this field. However, it is possible that such a strategy will involve applying increased refinement in regions of the domain far from the tidal farm, where the (curvature of) fluid velocity happens to be large, but that flow and its structure has little or no impact upon the energy output. For example, dynamics downstream of the farm and/or individual turbines if that flow does not impinge on further downstream turbines. This implies a computational inefficiency which we would prefer to avoid.
Goal-oriented error estimation provides a means of generating metrics which admit accurate energy output estimates, without wasting resolution unnecessarily on regions of the domain (and intervals in time) where the hydrodynamics are less impactful upon it. It does this using the adjoint problem associated with the prognostic equation set -in this case the shallow water equations (1). The adjoint problem depends on the QoI -in this case the energy output -and conveys how its sensitivities propagate across the space-time domain. Let q * denote the solution of the associated adjoint problem and q * h denote solution of the adjoint of the discrete shallow water model. Note that q * h lives in the same DG function space as q h .
For the purposes of the following presentation, we restrict attention to the finite element problem solved in a single timestep. The shallow water problem may then be expressed weakly in the form where r(·, ·) is referred to as the weak residual on the single timestep. Due to the pioneering work of Becker and Rannacher (2001) the energy output error may be approximated to first order by the error estimate which is known as the dual weighted residual (DWR). That is, the error accrued in evaluating the energy output may be approximated by a known expression involving the 'forward' finite element approximation, q h , and the error in the corresponding 'adjoint' approximation, q * h .

Error Indication
In order to apply mesh adaptation, we need to first deduce local contributions to the global quantity (8), so that a mesh adaptation method can make local modifications to the discretisation, as appropriate. The simplest means of doing so is to consider element-wise contributions. We refer to Wallwork et al. (2020) for details on how this can be performed in the case of shallow water tidal farm model.
Without further treatment, the adjoint error term e * := q * − q * h cannot be evaluated directly, since it contains the exact adjoint solution tuple. One approach is to substitute q * with a higher order approximation, such as may be obtained by solving the adjoint problem in a (globally or locally) enriched finite element space, or through application of a technique such as superconvergent patch recovery. An alternative approach is proposed in Becker and Rannacher (2001), which has a much lower computational cost than such methods in general. Its application to time-dependent tracer transport problems was considered in Wallwork et al. (2021). The so-called difference quotient method reformulates the error indicator on each element, rather than attempting to approximate the true adjoint solution. The resulting error indicators take the form where (u h , v h ) stand for the x-and y-components of u h , respectively, and ∂K denotes the edge set of an element K ∈ H with circumradius h K . Here (Ψ u , Ψ v , Ψ η ) denote components of the strong residual and (ψ u , ψ v , ψ η ) are flux terms due to the integration by parts and the DG discretisation. In practice, the Laplacians of the adjoint solution components are constructed using a recovery method. In this work, the gradient is obtained using an L 2 projection, the finite element derivative of which is an element-wise quantity. Note that the weighting term evaluates flux terms in the forward problem at the finite element adjoint solution.
There are a number of caveats associated with error indicator (9). Firstly, its sum over all elements consistently overestimates the true QoI error, with the overestimation increasing with mesh size. As such, it is not useful as an error estimator. Nevertheless, its localisation to individual elements provides useful information for mesh adaptation. This is because the metrics constructed from such error estimators (as detailed in the following subsection) are normalised in space, meaning that global scale factors are unimportant. Secondly, two terms which act to symmetrise the viscosity operator in the SIPG method must be dropped when constructing an error indicator of the form (9) from the P1 DG − P1 DG formulation. Finally, the estimator uses results which have as yet only been proved for elliptic problems. Despite these drawbacks, the results in Section 4 show that difference quotients can still give promising results when applied to the hyperbolic problem investigated here.

The Metric-Based Framework
Given an error indicator field, the metric-based mesh adaptation framework uses the concept of a Riemannian metric to dictate the transformations which should be applied to a mesh. In the twodimensional case considered here a Riemannian metric M = {M(x)} x∈Ω is a tensor field corresponding to a SPD 2 × 2 matrix at each point of the domain. Since the metric is defined point-wise, it gives rise to a Riemannian metric space, (R 2 , M).
The core idea of the metric-based mesh adaptation method is to use this Riemannian metric space within the mesher. Adaptation operations are applied to the mesh until it becomes a unit mesh when viewed in the Riemannian metric space. That is, each of its elements have edges of unit length when mapped into the metric space. In general, it is not possible to tessellate a 2D domain with equilateral triangles, so this definition is relaxed slightly to ensure a quasi-unit mesh. Note that, whilst the transformed mesh is quite regular, its representation in the (Euclidean) physical space can be rather distorted.
One of the great advantages of metric-based mesh adaptation is that it allows for control of element shape and orientation, as well as size. It is possible to define anisotropic goal-oriented metrics and perform anisotropic goal-oriented mesh adaptation based on scalar error indicators such as (9). For example, the a posteriori anisotropic metric due to Carpio et al. (2013) was applied to a steady-state tidal farm modelling problem in Wallwork et al. (2020). However, anisotropic metrics are not applied to time-dependent problems in this paper, for brevity; they will be used in future work. Instead, the isotropic case is assumed, whereby the metric is the identity matrix scaled by E K in each element. The metric-based mesh adaptation toolkit used herein expects metrics defined in P1 space (i.e. vertexwise), so an additional projection step is applied.
For further details on metric-based mesh adaptation applied to steady-state problems, we refer to Alauzet (2011a, 2011b) and Pain et al. (2001).

Time-Dependent Case
An additional implementation detail in the time-dependent goal-oriented case is the requirement to solve adjoint equations over a sequence of meshes. We take a fixed-point iteration approach, which involves partitioning the simulated time period into N subintervals, A different mesh is associated with each subinterval, and mesh-to-mesh interpolation is applied to transfer data from one subinterval to the next. The fixed-point iteration approach taken in this paper is based on that outlined in Belme et al. (2012); we refer to that work for details on the algorithm. The main difference in the paper presented here is that it uses goal-oriented metrics derived from an a posteriori error result (8) rather than an a priori one. A notable feature of both implementations is that they use space-time normalisation methods, which allow for the discretisation to be multi-scale in both space and time. In particular, DOFs are distributed across the temporal domain so as to improve the QoI accuracy. This means that the meshes associated with some subintervals may be much finer or coarser than others.
One of the key tunable parameters for metric-based mesh adaptation is the metric complexity. For steady-state problems, it is the analogue of the mesh vertex count in (continuous) metric space. For time-dependent problems, it is the continuous analogue of the sum of all mesh vertex counts over all timesteps. By increasing the target metric complexity, we allow for heightened overall mesh resolution. Whilst this typically implies a heightened computational cost, it should also imply a reduction in QoI error (provided that the metric is chosen appropriately).

Software for Mesh Adaptation
In this work, adjoint equations associated with (1) are formulated and solved using the discrete adjoint approach. This approach amounts to differentiating through the discretised model, as opposed to the 'continuous' model (as in the continuous adjoint approach). The discrete adjoint method is made available in Firedrake and Thetis by the dolfin-adjoint package (Farrell et al., 2013).
The Pragmatic (Barral et al., 2016) anisotropic mesh adaptation toolkit is used to take an input mesh and a metric defined upon it and generate a new mesh. In 2D mode, Pragmatic applies four different types of mesh transformation: vertex insertion, vertex removal, edge swapping and local Laplacian smoothing. The first three modify the mesh topology, whilst the fourth does not. Mesh-to-mesh solution transfer is achieved using the conservative interpolation operator provided by libsupermesh (Farrell & Maddison, 2011;Maddison et al., 2016).

Idealised Tidal Array
The numerical experiments presented in Wallwork et al. (2020) involve a simple idealised array of two turbines. In those experiments, the problem was solved to steady state, meaning there was no tidal forcing and the flow was actually laminar. However, the setup is useful for illustrating the impact that the locations of the turbines within the array can have upon the total power output. In that paper, the power output was found to be 15% lower when the turbines were aligned in the flow compared with when they were were offset by one turbine diameter in opposite directions orthogonal to the background flow. In addition, it was illustrated that the application of goal-oriented metric-based mesh adaptation can lead to more accurate approximations of array power output for similar numbers of DOFs.
In this paper, we consider the extension to the time-dependent case and a larger tidal farm. The test case was originally proposed in Divett et al. (2013) and consists of an array of fifteen turbines arranged in three rows and five columns. That work also demonstrated the importance of array configuration in tidal farm design; a staggered turbine layout was found to extract 54% more energy from the flow than a centred, aligned configuration. 1 This is because a staggered array gives the turbine wakes more opportunity to recover before interacting with another downstream turbine, as well as turbines being able to exploit accelerated bypass flow. Tidal turbines act to extract energy from the flow, so their wakes are effectively momentum deficits. For the version of the tidal array test case underpinning the numerical experiments in this paper, we consider the centred, aligned configuration alone, which implies strongly nonlinear interactions between turbines and wakes.
An initial mesh of the rectangular domain Ω = [−1500, 1500] × [−500, 500] with increased mesh resolution in the array region was generated using gmsh (Geuzaine & Remacle, 2009). Figure 1 shows the spatial domain, including the forced and free-slip boundary segments, the five columns of turbines and a zoom region used in subsequent plots.
The resulting boundary conditions are given by substitution of (11) in (3). Note that the tidal period is 10% of the M2 tidal constituent; it was reduced in (Divett et al., 2013) in order to emphasise vorticity. Figure 2 shows snapshots of the fluid velocity field over a range of time levels over the single tidal cycle [0, T tide ], as computed on a high resolution mesh with 184,896 elements and 554,688 DOFs in the P1 DG − P1 DG function space. Starting from a near zero initial velocity and a linear elevation field which satisfies the boundary conditions, the flow speed increases Westward for the first half of the interval, before reversing to flow Eastward. At time levels 0.2 T tide and 0.4 T tide , the West-most column of turbines experiences relatively laminar flow, whilst the flow around the other turbines is much more turbulent. At time levels 0.6 T tide and 0.8 T tide , it is the East-most column which experiences relatively laminar flow. The velocity field at t = T tide is stored to disk, along with the corresponding elevation. These 'spun-up' hydrodynamics are used as initial conditions for subsequent mesh adaptive simulations on the interval [T tide , 1.5 T tide ].

Mesh Adaptation Based on Energy Output
In this subsection, isotropic goal-oriented metrics are applied to the tidal array test case, with the aim of minimising the error in energy output. The simulated interval [T tide , 1.5 T tide ] is divided into 40 subintervals of equal length. A target metric complexity of 2,000,000 is imposed. Figure 3 shows At t = T tide , the spun-up hydrodynamics are interpolated as an 'initial condition'. Accordingly, the mesh associated with the first subinterval is adapted so that these hydrodynamics may be accurately captured. The goal-oriented metric does not attempt to capture the entire field, however: only the upstream hydrodynamics are held to be of importance; little mesh resolution is used surrounding the fifth (i.e. East-most) columns of turbines. The nonlinear interactions between turbines and the fact that turbines act to remove momentum from the flow means that the first (i.e. West-most) column of turbines extracts the most energy when the flow is Eastward. This can be seen in Subfigure 4a, which shows the power output of each turbine column under an extension of the high resolution fixed mesh simulation of Subsection 4.1 for another half tidal cycle. Turbine columns to the East do occasionally extract similar amounts of energy, but they generally extract much less and the levels fluctuate due to the turbulence. Moreover, it is not the instantaneous power output values that are of interest, but the areas under the curves in Figure 4. That the darkest (West-most) curve has by far the largest area indicates why the mesh adaptation algorithm gives so much precedence to the first column. Given that the first column is by far the most important and that we have an advection-dominated problem, dynamics downstream of this column are held to be far less important than the dynamics surrounding and upstream of those turbines. Turbines in the second column generate the least power because the flow conditions in the wake of the first turbine column are close to laminar, meaning that they experience a significant momentum deficit for much of the simulation. That the first column generates the most power and the second column generates the least power is consistent with results reported in Divett et al. (2013). Similar mesh structures as at t = T tide are apparent in the snapshots at t = 1.125 T tide and t = 1.25 T tide . Interestingly, little mesh resolution is deployed for the purposes of capturing vortices. In fact, most of the mesh resolution is used to capture the relatively laminar flow conditions found West of the second turbine column.
The arrow of time implies that dynamics at a particular time are only important from then onwards.
At the start of the simulation, the dynamics can potentially impact the entire solution trajectory, whereas towards the end of the simulation the impact is much more limited. This explains why so much more resolution is deployed at t = T tide than at t = 1.5 T tide , despite the velocity (and hence the contribution to the energy output) being near-zero in both cases. Figure 4 suggests that the power curve is slightly out of phase with the tidal forcing, with the outputs of all columns becoming near-zero before the tide turns. This is likely due to the significant drag forces inherent in the tidal farm. As a consequence, the power outputs of all columns are significantly larger at t = 1.125 T tide than at t = 1.375 T tide . Combined with the above argument about the arrow of time, this explains why the mesh used at time t = 1.375 T tide is so much coarser than that at t = 1.125 T tide . At t = 1.375 T tide and t = 1.5 T tide , the pre-specified maximum metric magnitude is being imposed within the turbine footprints.
Consider now the columnar power output curves shown in Subfigure 4b, which are due to the goaloriented mesh adaptation simulation. Despite the fact that even the adapted mesh with the highest element count stated (59,628 in Subfigure 3b) has less than a third of the high resolution fixed mesh element count, the power output curve for the first column agrees well with that shown in Subfigure 4a. The general trend of the second column's power output curve is also well represented, even though it contributes the least. This is partly because its upstream conditions are well resolved and partly because its variability is fairly small. The power curves for the remaining turbine columns are not well captured. This is for a number of reasons. Firstly, the power curves of these columns are much more difficult to capture than the first two columns because of the turbulent conditions there, of which the high variability in those curves is a symptom. Secondly, the target complexity is relatively small for a problem of this size; if it were increased then it is likely that the mesh adaptation algorithm would deploy more DOFs for the purpose of capturing columns 3-5. Note that the variation in columns 3-5 is reduced in Subfigure 4b, compared with Subfigure 4a. This is likely due to numerical diffusion introduced by the coarse mesh resolution in that region.
To conclude this subsection, we propose further avenues for investigation based on the numerical experiments presented above. Firstly, goal-oriented mesh adaptation is only applied over a single flood tide. In future work, it would be interesting to see how the adaptation algorithm acts to deploy mesh resolution over a sequence of flood and ebb tides. It would also be interesting to see the different ways in which mesh resolution is deployed for alternative configurations of the same farm, such as those considered in Divett et al. (2013). It was revealed in that paper that, unlike the centred and aligned configuration considered here, the second column of a staggered array generates the most energy, due to accelerated bypass flow. As such, it seems reasonable to expect that the goal-oriented meshes would look quite different. Given that the x-component of the flow conditions is much more significant than the y-component, it is plausible that incorporating anisotropy into the metric would help to more effectively resolve the dynamics. Finally, the above results only consider one fixed mesh resolution level and one target metric complexity. Consideration of multiple values for each of these parameters would allow for more detailed convergence analysis.

CONCLUSIONS
In this paper, we have applied a goal-oriented metric-based mesh adaptation method to the problem of simulating time-dependent hydrodynamics in an idealised tidal farm test case. This is the first time goal-oriented mesh adaptation has been used in the context of time-dependent tidal farm modelling. In Section 3, an accurate energy output assessment of the farm is sought via goal-oriented error estimation techniques. Section 4 demonstrates how the goal-oriented adaptation algorithm deploys mesh resolution over the spatio-temporal domain in order for this aim to be achieved, whilst retaining relatively few degrees of freedom in the adapted meshes.
In future work, we will apply goal-oriented metric-based mesh adaptation framework to more complex problems such as proposed tidal power infrastructure projects, with spatially varying bathymetry data and realistic tidal forcings. In turn, the mesh adaptation capability can be used to accelerate design optimisation calculations (Culley et al., 2016;Funke et al., 2014).