Quantifying Computational Efficiency of Adaptive Mesh Refinement for Shallow Water Solvers

Non-uniform, dynamically adaptive meshes are a useful tool for reducing computational complexities for geophysical simulations that exhibit strongly localised features such as is the case for example for tsunami, hurricane or typhoon prediction. Theoretical insight for meshbased numerical methods, however, is largely restricted to uniform meshes as they allow for a traditional definition of convergence and computational complexities linearly multiply with mesh resolution. To gain more insight into adaptive meshes and build a theoretical framework for their assessment, we present and discuss a number of mesh metrics that we apply to simulations on an adaptive triangular mesh. The latter is driven by physics-based refinement indicators that capture relevant physical processes and determine the areas of mesh refinement/coarsening. The mesh metrics take into account a number of characteristics of numerical simulations such as numerical errors, spatial resolution, as well as computing time and allow for a novel definition of convergence which can be used to demonstrate a convergence speed-up of adaptive simulations. Furthermore, we obtain numerical evidence for computational overhead as well as the distribution of the numerical errors across the various mesh resolutions.


Introduction
A large part of research in the geosciences depends on computer simulations. The latter are oftentimes used to predict or simulate experiments that would be too difficult, costly, dangerous or, straight out, impossible to carry out in the real world. The complexity of these experiments requires solving mathematical equations to a very high level of detail. Running these detailed simulations with a spatially uniform resolution, however, leads to long computing times. This, in turn, means that we require a lot of computational resources to carry out this type of research which ultimately limits the number and type of experiments that we are capable of simulating.
In recent years element-based models have gained popularity for solving equations describing geophysical phenomena. These models include continuous and discontinuous finite elements, see for example Brenner and Scott (2000); Hesthaven and Warburton (2008) as well as finite volumes, see for example Castro et al. (2012); LeVeque (2002) . The common feature of all these models is that the domain of interest is decomposed into polygonal cells. In these cells the computer model then determines values for the quantities of interest. As an example, for a simulation of flooding of a coastline, these quantities could be the water height D and the horizontal momentum of the water Du = (Du 1 , Du 2 ) . More detailed simulations require smaller cells. Hence, very detailed simulations are so costly because they require computations on a lot of (small) cells. As described in Behrens (2006) this does not mean that smaller scale processes are necessarily resolved. The quality of the resolution of underlying physics heavily depends on the numerical model that is used.
One approach to reduce computational cost is that of non-uniform resolution, i.e. the use of cells that have variable size within the domain. These non-uniform meshes can be static such as the ones used in Giraldo (2000) with a 2D discontinuous Galerkin (DG) model to solve advection problems or in Luettich and Westerink (1995) for the simulation of storm surge and have been continuously improved over the past years. In contrast to these locally refined but static meshes, adaptive meshes can capture evolving features as presented in Behrens (2006). Especially suited for the simulation of spatially localised and moving phenomena, these non-uniform meshes can take the form of nested meshes as in Mandli and Dawson (2014) for the simulation of hurricane storm surge, or of unstructured meshes, see for example Piggott et al. (2008). The general idea of all the above is identical: Areas of interest are refined while areas that are considered not as crucial for the simulation output are kept relatively coarse. This has the advantage of reducing the overall computational cost of simulations while still retaining a high resolution in relevant areas. In this manuscript we will consider adaptive triangular meshes as described in Behrens et al. (2005) and use bisection along a marked edge as a refinement strategy which has advantages in terms of computational efficiency at the cost of not allowing for fully unstructured mesh refinement.
Note, that for dynamical phenomena, the area of interest of a simulation may move with time, hence complicating the determination of the area to be refined. Therefore it is important to identify a good measure that captures the movement of the area of interest. We will refer to this measure as a refinement indicator η τi with the understanding that τ i denotes the ith element of the mesh and it will be used to automate the mesh manipulation process as we will describe in more detail in subsection 3.1. Automatic mesh refinement has been shown to improve simulation accuracy Dietachmayer and Droegemeier (1992) for idealised meteorological applications even when disjoint areas of high resolution are used. Furthermore it was shown in Iselin et al. (2002) that for advection problems, they can reduce diffusive and phase errors. The needed refinement indicators require insight into the physical problem and ideally represent the sensitivities of the numerical simulation. The quality of a non-uniform simulation will depend on how well the refinement indicator captures the numerical error. Two different types of refinement indicators have been discussed in the literature: (a) heuristic/ physics-based error indicators, see for example Müller et al. (2013); and (b) adjoint-based refinement indicators as in Beckers et al. (2019); Farrell et al. (2013).
The goal of using adaptive mesh refinement (AMR) for numerical simulations is to reduce the number of degrees of freedom. Methods for AMR, however, cause additional computational overhead for grid management and manipulation. A priori estimates that quantify if the reduction of degrees of freedom outweighs this overhead are generally not possible to obtain. Furthermore direct comparisons of uniform and adaptive simulations are difficult because the traditional definition of convergence does not apply. Hence, this study attempts a structured quantification of the efficiency of AMR using mesh metrics. This goes beyond the first attempt to quantify mesh quality which has been made in Gu et al. (2001), with a sole focus on posteriori error estimates. The set of mesh metrics that we propose will correlate many parameters of adaptive simulations such as resolution, numerical errors, and computing time to name a few. This has the advantage that it allows us to directly quantify how large the overhead of AMR is and how large the area of refinement may be in order for AMR to be more efficient than uniform mesh refinement. A definition of numerical convergence taking into account the numbers of degrees of freedom will furthermore allow to compare convergence properties of uniform and adaptive simulations and reveals that AMR can lead to improved convergence rates.
The remainder of this manuscript is organised as follows. In section 2 we define the set of metrics for measuring various aspects of computational efficiency of AMR, section 3 then describes the numerical 2D shallow water model that we used for this study and provides detail on the adaptive mesh refinement and the required mesh parameters. In section 4 we then investigate three test cases: Subsection 4.1 shows how the metrics can be applied to a localised travelling vortex. Throughout the simulation the percentage area that is refined until the smallest mesh width is almost constant. We discuss percental error distribution on the mesh levels and how AMR can influence the relationship of resolution to degree of freedoms depending on mesh parameters. In subsection 4.2 we then study the mesh metrics for a test case that involves wetting and drying and a very sensitive refinement indicator. In this case the percentage area changes during evolution, hence making it more challenging for the automated AMR. In section 4.3 we finally show a test that one could expect to not be suited for AMR. A centred impulse at time t = 0 that evolves and produces waves throughout the entire domain. Even here, we can show that the adaptive mesh captures the initial impulse and the subsequent propagating waves and due to the temporal localisation of the impulse, we achieve a higher computation efficiency for the entire simulation. Overall, we find that the use of the adaptive mesh improves computational efficiency with respect to all mesh metrics. Finally section 5 summarises the findings of the paper and discusses how the presented metrics will be of use for numerical modellers who seek to assess the computational efficiency of their models.
-Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.

Metrics of Adaptivity
In order to set the scope for our evaluation of properties of adaptive numerical methods, we introduce a number of evaluation criteria -metrics -especially selected for adaptive mesh refinement methods. While in principle all the accuracy metrics applied to uniform grid methods are available for adaptive methods as well, we want to focus on the evaluation of adaptive methods in this study. By introducing a number of different viewpoints on how to measure performance of an algorithm or method, we also counteract the common problem of comparability of different implementations. The same method implemented in two different ways can expose very different performance characteristics. Performance of adaptive algorithms can be assessed from different viewpoints. In general, it is possible to focus on accuracy, on computational requirements such as run time or memory, or on resolution. It is also interesting to assess the overhead imposed by the additional capability to adaptively control the mesh. So, in the remainder of this section, we will introduce three groups of metrics and relate them to common evaluation hypotheses.

Metrics focusing on accuracy
Accuracy can only be measured for certain test cases for which analytical solutions exist. In case of unknown analytical solutions, convergence behaviour can be used as a replacement, assuming that the algorithm is correct and converges to the exact solution. In that case a (costly) high resolution uniformly computed solution can serve as a reference. In both cases, the error can be measured in an appropriate norm. Commonly, the 1 , 2 and ∞ error which we denote with e 1 , e 2 , and e ∞ respectively is used, defined as follows: where u and u h are the exact/reference and computed solution respectively, and Ω represents the computational domain. Of course u needs to be projected appropriately to the discrete space in order to be able to compare it with u h . Furthermore, we introduce the symbol n Ω to denote the number of degrees of freedom in a domain Ω. In order to derive a comparable metric for different methods, one can now ask how many degrees of freedom (DOF) a given algorithm needs for a fixed error. This is a common task in engineering, where a certain accuracy is required and the cost for achieving it is to be minimised. One can also state this question inversely: given a fixed n Ω , what is the error a method generates. Since in most cases the number of DOF directly relates to the computational cost, both questions are closely related to corresponding questions measuring the run time or memory requirements with respect to the error.
In all the above cases we consider ratios of certain values: • The ratio of error over DOF: r d fixed = eρ nΩ for ρ ∈ {1, 2, ∞}, • The ratio of DOF over error: r e fixed = nΩ eρ = r −1 d fixed for ρ ∈ {1, 2, ∞}, • The ratio of time over error r t−to−sol = t eρ , which corresponds to a time-to-solution with given accuracy and ρ ∈ {1, 2, ∞}.
In the naming of the metrics d fixed stands for number of DOF fixed, and e fixed for error fixed. This is to emphasise the different perspectives that a user of the metrics might be interested in. The metric r d fixed stresses that numerical errors can differ for equal (or fixed) numbers of DOF depending on the location of DOFs, where as r e fixed stresses that the same numerical error can be achieved with different numbers of DOFs.

Metrics focusing on resolution
In many geoscientific or physics-based multi-constituent applications, accuracy is hard to determine, so resolution plays an important role. Usually it is measured in the (local) mesh size or shortest wave length resolvable. In Gu et al. (2001), an attempt was made to quantify how well flow was resolved. An important question for this study is, how many DOF are needed in order to resolve certain features.
-Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.
Another question is, given a certain number of DOF (limited e.g. by the given computing infrastructure), how fine can the (local) resolution be and therefore the size of the physical model. An interesting question for method inter-comparison is what error can be achieved by which resolution. Of course, this question can only be answered if a converging reference solution is known. Now, in order to rigorously define metrics related to resolution, a definition of resolution needs to be given. Since in this manuscript triangular meshes are considered, we use the equivalent measures shortest edge length and radius of inscribed circle. In both cases we will use the expression h = min(h τ ), where h τ denotes either the length of edge τ or the radius of the inscribed circle in cell τ . The minimum is taken globally over all edges or cells, respectively. Note that in higher order methods, h τ needs to consider the distance between DOF and not the shortest edge length of cells, since waves can be resolved within cells.
Then we can define the following metrics, again as ratios: • The ratio of resolution over error r eff−res = h eρ for ρ ∈ {1, 2, ∞}, • The ratio of resolution over DOF r res = h nΩ .

Metrics focusing on computational resources
Trying to assess computational demands is a highly challenging task, since inefficient implementations can hamper an objective assessment of such characteristics for different methods. However, it is one of the few really relevant metrics, since in practical applications one is interested in obtaining a most accurate solution with minimal computational effort. Hence, of importance here is the metric already mentioned above as r t−to−sol . Furthermore, one can assess the time an implementation needs per DOF. By this, an efficient implementation can be characterised. It is also interesting to look at this ratio for changing problem sizes. Optimally, the ratio should keep constant for increasing problem sizes, but often cache effects or non-optimal algorithmic complexity inhibits this. Interesting questions arise also from memory requirements. The amount of memory necessary to achieve a certain error or a certain resolution is an important measure for adaptively refined algorithms.
The following ratios may be helpful to assess computational efficiency of adaptive algorithms: • The ratio of time over error r t−to−sol (see above), • The ratio of time over DOF r t−per−DOF = t nΩ , • The ratio of memory required for resolution r mem−res = M B h , where M B is the memory in bytes.

Further remarks
One very important information for assessing the quality of adaptive mesh methods is the amount of overhead the adaptive control imposes on the algorithm, compared to a non-adaptive method. This is often only assessable by rigorous profiling of a running program. A second important prerequisite for a rigorous assessment of adaptive methods is an accurate and robust error estimator or refinement criterion. A description of our refinement strategy can be found in section 3.1. If the local adaptation strategy selects the wrong (or insufficient) area for better resolution, then a comparison is often corrupted. For example, none of the error-based metrics defined in subsections 2.1 and 2.2 could be utilised; in this case the error will always be bad, since the refinement strategy does not capture the relevant areas. It is therefore paramount to compare a uniform grid solution to an equally resolved adaptive solution first and to make sure the difference with respect to solution error is not dominating.

The Numerical Model and Adaptive Mesh Refinement
As a representative model for this study we chose a discontinuous Galerkin (DG) model that solves the depth-integrated shallow water equations in two dimensions which can be written in flux form where the prognostic variables are U = (D, Du) : the water depth D and the 2D momentum Du defined on a spatial domain of interest Ω ⊂ R 2 . Spatial coordinates are denoted as x = (x, y) ∈ Ω. ∂U ∂t =: U t Version: November 28, 2019 -Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.
denotes the temporal derivative and ∇· = ∂ ∂x · + ∂ ∂y · is the divergence operator. F is a flux function and S a source term defined as where g = 9.81m s −2 is the acceleration due to gravity, ⊗ a vector product in R 2 , and I 2 is the 2 × 2 identity matrix. For the reader familiar with shallow water models, we note that we chose to denote the water height with D instead of the more commonly used h, which in this manuscript, we have reserved for the resolution as we think it improves readability. The source term models a temporally constant bathymetry b = b(x). Throughout this paper vector valued quantities are indicated by a bold print while all other quantities are assumed to be scalar.
The discretisation that we are using is a discontinuous Galerkin (DG) discretisation obtained through three steps: (a) decomposing the computational domain Ω = i τ i into conforming triangles with i the index of the triangles; (b) approximating the prognostic variables U = k U k (t)φ k (x) by linear Lagrange polynomials with k the corresponding index for the basis functions; and (c) integrating the resulting equations in space against test functions. In the present model the test functions coincide with the basis functions φ k used for the approximation, so that the resulting semi-discrete system reads for all triangles τ i and j = 1, 2, 3 leading to a second order accurate spatial discretisation. Here, F * is a numerical approximation of the flux at the cell interfaces. In our model, we used Rusanov's flux (see for example Toro (2009)). Note that we integrated the flux term by parts twice to obtain the strong form of the equations as in Hesthaven and Warburton (2008). In this form, we integrate the jump over the edges which has desirable properties with respect to wellbalancing as shown in Beisiegel (2014) and more recently in Vater et al. (2019). After numerical integration and re-organisation, the system (1) can be written as where H denotes the discretised version of the right hand side which includes the fluxes as well as the source terms. As this model is used for flood simulations, slope limiting is required to prevent spurious velocities at wet/dry interfaces. The slope limiter we are using is velocity-based and described in detail in Vater et al. (2019). In the latter study we show that it can be successfully applied to tsunami benchmarks and to model flood scenarios. It is almost parameter-free and robust when applied to unstructured meshes. The system (2) can be solved using a strong stability preserving (SSP) multi-stage Runge Kutta method. We used Heun's method (RK22) for this study.

Computational Efficiency and Adaptive Meshes
The main focus of this study is on characteristics of the underlying adaptive mesh and not on the numerical model presented in the previous subsection. The computational model described in the previous section and in more detail in Beisiegel et al. (2019) uses the grid generator amatos (see Behrens et al. (2005)) to create dynamically adaptive and conforming triangular meshes. Finer triangles are obtained by bisecting coarse triangles along a marked edge (usually the longest edge) as suggested in Rivara (1984). The dynamic mesh adaptation process involves problem-dependent refinement indicators η τi , such as the absolute value of the gradient of fluid height at time t: for each element τ i as in Behrens (2006), to control the element-wise refinement and coarsening. Using problem dependent and user-defined tolerances 0 ≤ θ crs < θ ref ≤ 1 the mesh manipulation is then carried out as follows: with η max = η max (t) := max τi⊂Ω η τi (t) the maximum value of the refinement indicator over all elements. The values 0 ≤ θ crs < θ ref ≤ 1 mark the fraction of the maximum error below/above which an element Version: November 28, 2019 -Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.
is considered for coarsening/refinement. The size of the smallest and largest element of the mesh is determined through mesh levels λ crs ≤ λ ref . Starting from an initial (coarse) triangular mesh, the mesh generator will uniformly refine λ crs times to establish the coarsest mesh level and then, using the refinement indicator, refine areas of interest until the desired finest mesh level is reached using the tolerances as described in (3). The mesh node values are then interpolated or restricted after modification using the known Lagrange basis functions φ k in each element with an additional mass correction that has been implemented for the coarsening to ensure mass conservation.
The amatos mesh, created as described above, has the further advantage that it uses a cache-efficient space-filling curve-ordering of elements (see Behrens and Bader (2009)), which allows fast access of neighbouring elements: A feature that is particularly beneficial for localised numerical methods such as the presented discontinuous Galerkin since elements only communicate over edges.
For convenience, the meshes are kept conforming, i.e. free of hanging nodes, throughout the simulation. We stress that this is not required by the method itself. Hanging nodes would require to combine two or more Riemann solutions over one (coarse) edge as is for example recently done in Kopera and Giraldo (2014) and Hermann et al. (2011).

Application of Metrics and Numerical Tests
The mesh metrics defined in section 2 are a computational tool for quantifying the efficiency of adaptive simulations. In the following, we will compare the computational efficiency of adaptive simulations using the mesh described in subsection 3.1 with a corresponding uniform simulation. We will focus on three dynamically different problems: a quasi-stationary travelling vortex in subsection 4.1, for which the percental area of elements on the finest mesh level stays constant, long wave resonance in a paraboloid basin in subsection 4.2, for which the dynamics of the problem determine the percental area of fine resolution elements and finally a centred impulse in subsection 4.3 which only exhibits a localised feature of interest initially. We will apply the metrics to all mentioned test cases and demonstrate a gain in computational efficiency for the adaptive simulations -even for the centred impulse whose temporal behaviour would seem unsuited for AMR due to the expansion of a large number of waves during evolution.

A Strongly Localised Feature: A Quasi-Stationary Vortex. Mesh Metrics and Error Distribution
In where r = x − c 2 is the radial distance from c, the parameters v max = 0.5 and r m = 0.45 are maximum tangential velocity and the vortex radius respectively and the scaling factor s is defined as , where r vm = 1 2 −2 + 2 1 + 4r 4 m is the radius of maximum winds. With a background height, h bg = 1, and a background velocity, where x and y are the spatial coordinates x = (x, y) . The final time of the simulation is T end = 4s, i.e. until the vortex reaches its starting position again. The analytical solution for this test can be found in Vater (2013). We ran a number of uniform and adaptive simulations with mesh levels λ ref ranging from 6 to 14.
(B) A fixed coarse mesh level λ crs = 4. Higher resolved simulations were obtained by increasing λ ref ; The numerical water depth at t = 1s with a corresponding adaptive mesh is depicted in Figure 1. We observe that after an initial calibration period, the mesh follows the vortex and refines only the area of interest. The calibration period is needed because the initial conditions are not exactly balanced on a discrete level, so that gravity waves are emitted which vanish as soon as the model reaches a balanced state. Over time, the percentage of elements on a given mesh level do not show a lot of fluctuation as shown in Figure 2. Since the vortex is not changing in size over time, this is an indication that it is well captured by the adaptive mesh and internal iteration to flag elements for refinement and coarsening are kept minimal which adds to the computational efficiency of the mesh manipulation.
The proposed mesh metrics have been plotted in Figure 3. The solid line depicts results obtained on a uniform mesh while the two non-solid lines refer to adaptively refined meshes: The dashed line represents strategy (A) and the dashed-dotted line represents strategy (B). The adaptive simulations both used the same refinement indicator η τi as stated above.
We observe that both sets of adaptive simulations, using strategy (A) and (B), overall, give better results than the set of uniform simulations. In particular, the top left panel shows r d fixed for all three sets of simulations and indicates an accelerated convergence with respect to the number of degrees of freedom n Ω . This modified definition of numerical convergence uses n Ω instead of the traditional mesh width h. For uniform simulations both definitions lead to an equivalent result. The advantage of this modified definition is that it improves the comparability of adaptive and uniform simulations. Using linear least squares regression analysis, we find that the mean slope of the solid line for r d fixed is 1.04 while the mean slope of the dashed dotted is 1.41 -an increase of over 40%.
Moreover, we note, that using an adaptive mesh we consistently achieve the same L 2 -error with less CPU time as can be seen from the metric r t−to−sol (top middle panel of Figure 3) indicating an increase in computational efficiency. Per degree of freedom, however, the adaptive simulations are computationally more expensive. This can be seen from the metric r t−per−DOF (bottom middle panel in Figure 3). There, we show that for a given number of degrees of freedom n Ω , the uniform simulation (solid line) requires less CPU time per degree of freedom. The reason for this is the computational overhead caused by mesh manipulation and management when dynamically adaptive meshes are used. We remark, though, that the overhead produced by adaptive simulations is constant (the plotted lines are merely shifted by a constant) and approximately does not increase with increasing n Ω . The metric r eff−res (top right panel of Figure 3) furthermore shows that for a decreasing mesh width h, i.e. increasing λ ref , the adaptive simulations yield a higher effective resolution. This means that given a fixed number of degrees of freedom n Ω , adaptive simulations resolve finer features if the underlying equations model them. As an example Behrens (2006) mentions that you cannot expect to resolve turbulence with a fine resolution if your equation does not model turbulence. Finally metric r res gives insight about how the local resolution changes by increasing the number of degrees of freedom (bottom left panel of Figure 3). Since this depends on the systematic ways in which λ crs and λ ref are chosen, we elaborate on this in the following paragraphs.
We investigate the two refinement strategies (A) and (B) with more analytical rigour here, looking at r res as our primary metric. When considering an adaptive mesh, we observe that elements of a number of mesh levels between λ crs and λ ref are present, whereas in uniform meshes only one mesh level occurs.
Let us look at the number of elements as an indicator for the number of degrees of freedom (or unknowns), determining the computational work load. The exact number of unknowns can be determined by Euler's formula for plane graphs that relates the number of vertices, edges and cells by #vertices − #edges + #cells = 2, considering the location of unknowns. In our DG method, all unknowns are located on cells (elements), which means our n Ω is directly related to the number of cells.
The refinement strategy implemented by the grid generator amatos bisects elements, thus a refinement step generates two elements out of the elements marked for refinement. With this knowledge we can now -Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al. study the effect of the different refinement strategies. Uniform refinement will increase the refinement level by one (i.e. λ ref → λ ref +1, note that for uniform meshes λ crs = λ ref ) and will double the number of elements (i.e. n Ω → 2 · n Ω ). Refinement strategy A will increase both, λ crs and λ ref by one, but will leave d λ unchanged. This means all elements of levels λ crs through λ ref will be refined, which is all elements altogether. Therefore, n Ω → 2 · n Ω as in the uniform refinement case. This is confirmed by the metric r res as it shows the dashed line (adaptive strategy A) being parallel to the solid line (uniform simulation) -the number of degrees of freedom n Ω evolve in the same way. In contrast to this, refinement strategy B will only increase λ ref , in other words λ crs → λ crs , λ ref → λ ref + 1, and d λ → d λ + 1. Note that only those elements found by the refinement indicator will be refined. Most likely, these will be only those that are already on the (previous) finest level. Therefore, the additional number of elements depends on the fraction of the domain found relevant by the error indicator.
Let us assume that we start with a uniform mesh in which about 10% of the elements are flagged for refinement. Then only those 10% elements are refined which means the number of unknowns is doubled only for those 10% (i.e. n Ω → n Ω + 0.1 · n Ω = 1.1 · n Ω ). Using this formula recursively for increasing numbers of refinement levels shows how dramatic the reduction of complexity can be.
To conclude our investigation of this test case, we studied the spatial distribution of the numerical error. The goal is to minimise the numerical error through capturing it within the fine resolution area of the simulation. It is expected that adaptive meshes lead to a more uniform distribution of the numerical error compared to uniform meshes if the refinement indicator sufficiently captures the areas where numerical error is large. We show examples of the distribution of the numerical error in Figure 4. Overall, we observe that both simulations yield a localised error, i.e. the majority of the error is captured in the high resolution area of the simulation. The adaptive simulation (left panel of Figure 4) shows a few elements with a higher error close to the transition to coarser elements as expected while the remaining error is equally distributed throughout the fine mesh area. We note that the uniform simulation used 131072 elements as opposed to the adaptive simulation which used on average 9476 elements, or less than 10% the number of elements compared to the uniform one.
To further illustrate the different behaviours of the two adaptive strategies (A) and (B) and to support the hypothesis that keeping λ crs as small as possible leads to efficiency benefits, we consider error distributions for two representative examples. Tables 2 and 1 show the distribution (per cent) of elements on a certain mesh level (left) as well as the distribution of numerical L 2 errors per mesh level (right). The columns of the tables correspond to one simulation with mesh levels as stated, e.g. the first column in table 2 corresponds to an adaptive simulation with mesh parameters λ crs = 4 and λ ref = 6. Using θ crs = 0.1, θ ref = 0.2 for all simulations, we observe that the majority of elements resulting from refinement strategy (B), are on the finest mesh level in comparison to configuration (A), where the majority of elements are both on the finest and coarsest mesh level. We furthermore see that for configuration (B), the largest part of the error is also to be found on the finest mesh level while the error on level 4 is decreasing with increasing λ ref and all intermediate mesh levels λ crs ≤ λ ≤ λ ref only contain small errors. For configuration (A) we obtain the majority of errors on the finest mesh level as well. The error on the coarsest level, however, is increasing with increasing the overall resolution.

Refinement Indicators For Highly Sensitive Problems: Wetting and Drying and a Dynamically Changing Area of Maximum Refinement.
In this section we explore the behaviour of our mesh metrics to a challenging problem that includes the modelling of wetting and drying as well as a high resolution area that is dynamically changing in size according to model dynamics. Present high-order information, when not controlled properly through the refinement indicator η τi , will propagate into coarser mesh regions where it can cause large numerical error. This is partially due to the non-optimality of the marking strategy that we employ. As shown in Dörfler (1996) for Poisson's equation, using the L 2 -norm instead of the maximum norm in the definition of η τi , we obtain the least amount of refined triangles at the cost of a possibly larger number of iterations for the mesh manipulation. The presented computational model described in section 3, see also Beisiegel et al. (2019) is capable of simulating wetting and drying using slope limiters during the evolution. This will prevent spurious velocities from developing close to the wet/dry interface. The fluid depth D will not be set to zero in these areas unless below a small cut off tolerance = 10 −K with K ∈ N which may lead to small water films remaining in dry areas -an effect that makes the wetting and drying very sensitive to adaptive mesh refinement of the wet area of the domain Ω -an aspect, that we will describe in the following subsection.
-Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.

Long Wave Resonating in a Parabolic Basin
This nonlinear problem can be found in Lynett et al. (2002) in a scaled version and its original analytical solution was first determined in Thacker (1981). In a domain Ω = [0, 8000] 2 , a parabolic basin of shape b(x) = h 0 (1 − r 2 a 2 ) with h 0 = 1 the centre of the basin, a = 2500 the distance from the centre to the shoreline and r = x 2 the radius, the initial water depth is described as where the constants are determined as Similar to the experiment in subsection 4.1, we ran sets of simulations with uniform and adaptive mesh refinement with mesh parameters ranging from 8 ≤ λ crs , λ ref ≤ 17, θ ref = 0.05, and θ crs = 0.001 with a time step of ∆t = 16 s for h = 207.11 m which corresponds to λ ref = 8 and smaller ∆t for larger mesh levels in a way that the CFL condition is still fulfilled. For the adaptive simulations, we chose the refinement strategy labelled strategy (A) in the previous subsection 4.1 with d λ = 3 to ensure that the wet area of the domain is completely captured within the high resolution area of the adaptive mesh. For a more detailed description on how this is achieved, we would like to refer the reader to subsection 4.2.2.
As a stable refinement indicator we chose η τi = D τi , the mean value of the fluid depths over each element τ i ⊂ Ω. Figure 5 shows an example of the simulated fluid depth D and a corresponding adaptive mesh at time t = 2000s. We observe that the fluid is well captured by the adaptive mesh and the dry area is kept coarse throughout the simulation, hence minimising computations in areas which have no influence on the dynamics of the problem. Applying the metrics defined in section 2 to all simulations of this test case, we obtain Figure 6. Depicted are results obtained on a uniform mesh (solid line) and on an adaptive mesh (dashed line). We observe that the adaptive simulation consistently achieves a lower numerical error for a given number of degrees of freedom n Ω as demonstrated by the metric r d fixed in the top left panel of Figure 6. Using the modified definition of convergence that uses n Ω instead of the more commonly used h, we see that numerical convergence is accelerated by about 10% in this case which is significantly lower than for the -Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.  -Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.
test case presented in subsection 4.1. We attribute this to the lacking smoothness of the solution at the wet/dry interface In line with the findings in subsection 4.1, the mesh refinement using strategy (A) makes it possible to achieve a finer spatial resolution for a given number of degrees of freedom as indicated by r res in the bottom left panel of Figure 6. This plot furthermore shows that the corresponding plotted line of the adaptive simulation (dashed line) is a linear translate of the uniform simulation as would be expected from the use of strategy (A) (see also final paragraphs of subsection 4.1). Furthermore, as can be seen by r eff−res in the top right panel of Figure 6, the effective resolution is consistently higher for a given numerical error.
The rapid speed of change of the wet/dry interface and with that the high resolution part of the mesh are challenging for this physics-based and automated type of adaptive mesh refinement. Changing the percental area of the high resolution area of the mesh requires a high number of mesh manipulations, i.e. changing the size of the domain that is refined until the finest mesh level λ ref . This is in contrast to the previous test case in subsection 4.1 that comprised a propagating local feature that did not change in size over time. Although there is a significant amount of time spent in mesh manipulation and management, the metric r t−to−sol (top middle display of Figure 6) shows that the time to solution is smaller for the adaptive simulations while the time per degree of freedom (see r t−per−DOF in bottom middle display of Figure 6) is slightly larger for the adaptive simulation. Confirming the findings from the previous subsection.
To conclude, Figure 7 furthermore shows that the adaptive mesh captures the majority of the numerical error in the high resolution part of the mesh (right display) and that compared to the uniform simulation, the point wise error is smaller.

A Note on Grid Parameters
The adaptive mesh refinement that is used in this study heavily depends on a good choice of refinement indicator η τi as described in subsection 3.1 as the automated mesh manipulation process is driven entirely by it. Hence, it is important, that the refinement indicator is a good proxy for numerical error. For example, in the presence of wetting and drying (see previous subsection 4.2), it is important to only finely resolve the wet part of the domain. The reason being that artificial waves -in our case present close to the wet/dry interface -might not be completely filtered out by the slope limiter in the sense of D ≡ 0, and, hence, might be artificially amplified if energy from high-order modes is propagated into lower-order modes without additional filtering and with that causing numerical error.
Moreover, the adaptive mesh refinement described in subsection 3.1 shows the dependence of the automatic refinement and with that of the simulation quality on four key parameters λ crs , λ ref , θ crs , θ ref . Our observations confirm that the mesh refinement is sensitive towards the choice of the λs which determine the range of spatial resolution of the simulation, and of the θs as shown in Figure 8, which determine the fraction of the maximum error below/above which an element is considered for coarsening/refinement. We stress that the mesh manipulation is based on physics-based, heuristic refinement indicators, i.e. regions are considered for refinement or coarsening where these indicators are large or small. Therefore, good knowledge of the underlying problem is required to fulfill the above requirement that refinement indicators (proxies) correlate well with numerical error. As said before, the choice of the λs determines the coarsest and finest mesh width or in other words the length of the waves that are resolved in the simulation. In general, a higher resolution can be achieved by increasing the number of degrees of freedom n Ω which can be initiated by mesh refinement (h-refinement) as in this study or through increasing the order of the used polynomials (p-refinement). We note that it was found in Jameson (2000) that the number of points needed for a fixed numerical error is not constant with respect to the order of the numerical scheme -in our case the order of the polynomials φ k -but rather declines with increasing order. In theory, d λ can be arbitrarily large. However, sufficient physical space surrounding the feature of interest is required such that elements on all mesh levels from λ crs to -Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al. λ ref may envelop it. This is because mesh levels continuously increase from coarse to fine and in the absence of the required space, the feature of interest either will not be refined with the desired high resolution or the coarse area will be refined more than desired, hence the maximum benefit from AMR might not be obtained. For that reason, we solely chose strategy (B) in subsection 4.2.
The choice of the θs can be particularly crucial in cases where the refinement is determined through a refinement indicator that is very sensitive towards small perturbations as the one described in section 4.2. In this case an increase of the refinement area (a decrease of θ ref ) can lead to an increase of numerical error as small errors at the wet/dry are amplified when they are higher resolved by the finer mesh. This supports that the adaptive mesh refinement can only be as good as the underlying refinement indicator.

An Only Initially Localised Feature: A Centred Impulse
This subsection investigates the viability and benefit of using AMR for phenomena that are only initially local and subsequently spread throughout the entire domain in a way that capturing it is impossible while keeping the number of elements on the finest grid level λ ref constant. To illustrate this further, let us consider the following: Let Ω = [2, 2] 2 be a square domain with periodic boundary conditions. The initial water height is defined by D(x, 0) = 1.5 · e − x−xc 2 2 α with a shape parameter α = 0.001, and x c = (1, 1) the centre of the domain. The bathymetry is flat b(x) = 0 and velocities are assumed to be zero, i.e. u = 0. Figure 9 shows two snapshots of the numerical solution of the fluid height D on a uniformly refined grid with λ ref = 13. It is visible that the initially localised elevated fluid height propagates through the entire domain. We have run simulations with a number of varying spatial resolutions with corresponding mesh levels λ ref between 6 and 13. There is no analytical solution available for this test case, so in order to compute the numerical error we chose the numerical solution on a uniform mesh with λ ref = 13 as a reference solution. Figure 10 shows an example of such a comparison with an adaptive mesh with λ crs = 7 and λ ref = 10. The adaptive mesh is plotted in the third row of Figure 10 and shows that it captures the waves even with a growing area of maximum refinement. Analysing the numerical error over time, we observe that the point wise ∞ error is largest at t = 0 which is why we chose to compare numerical results at this time step t = 0. The difference between the adaptive simulation and the numerical reference is shown in the bottom row in Figure 10. We have chosen different colour scales for each plot in this row to increase readability. Using the numerical error as described above, we can compute the metrics defined in Section 2. The result can be found in Figure 11 and allows for many of the same conclusions as the previous sections. Using r d fixed , we can see that even in the more conservative ∞ norm, the same error can be achieved with a lower number of degrees of freedom. Furthermore, r res allows for the conclusion -Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.  -Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.  -Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.
that a smaller mesh width is achieved with less DOF (n Ω ). Both metrics, r t−to−sol and r eff−res however indicate that the adaptive simulation does not lead to any improvements in computational efficiency. Times to solution are identical and so is the effective resolution. The reason for this can be seen from r t−per−DOF : Since the refined area (area of the mesh, where the mesh level is λ ref ) increases over time, additional computational overhead is created. This overhead increases with increased number of DOF because additional elements are flagged for refinement and manipulated. However, this does not lead to the conclusion that adaptive simulations are not useful for this type of application. In fact, the cost per DOF increases over time with increasing area of maximum refinement as can be seen from Figure 12. It shows that towards the end of the simulation the cost per DOF increases by a factor of 4. We conclude that especially for highly time sensitive applications such as tsunami propagation adaptive simulations do lead to an increased efficiency especially during the important first moments of the simulation. This is furthermore supported by Figure 13 which shows (left display) that the number of elements is increasing over time in order to capture the ocurring waves. Moreover, as can be seen from the right display, we observe that the cpu time is not linearly increasing over time. This is an indication that AMR is still useful even if only temporarily the feature of interest is spatially localised.

Discussion and Concluding Remarks
In this study we have introduced a set of metrics that can be used to measure the computational efficiency of adaptive mesh refinement (AMR). Our adaptive mesh is structured and refined by bisection along longest edges for efficiency reasons as elaborated on in Behrens and Bader (2009). In a comparison with uniform simulations using the same numerical model, we showed that adaptive simulations are computationally more efficient with respect to the metrics. Using the Discontinuous Galerkin model described in Beisiegel et al. (2019) as an example, we studied a set of three test cases: (1) A spatially localised quasi-stationary travelling vortex; (2) a challenging long wave resonance in a paraboloid basin that comprises wetting and drying; and (3) a centred impulse that is only initially spatially localised.
While the test case 1 (see also subsection 4.1) comprised a roughly constant sized area of high resolution the test cases 2 (see also subsection 4.2) and 3 (see also subsection 4.3 comprised a highly refined area that dynamically changed in size. In Test case (3) even beyond the point where it could be labelled as localised.
Our focus then was on dynamically adaptive triangular meshes that were refined by bisection using physics-based, or heuristic, refinement indicators η τi as described in subsection 3.1. Insight into the dynamics of the numerical test problems allowed us to find indicators that were good proxies for numerical errors. The latter is important as the quality of the AMR strongly depends on the quality of the refinement indicator.
The metrics correlate numerical error to numbers of degrees of freedom n Ω (DOF) and CPU time as well as local (minimal) resolution h and made it possible to gain new insights into AMR. Using the metric r d f ixed = e 2 /n Ω , we consistently find that simulations using AMR achieve the same numerical 2 error with significantly fewer DOF. Re-defining convergence by using n Ω instead of mesh width h, the metric r d f ixed can furthermore be used to compare convergence properties of uniform and adaptive simulations. For all test cases we find that the dynamically adaptive mesh leads to an accelerated convergence and in case of the stationary vortex in subsection 4.1 even by up to 40%. This has an impact on the total run time of the solution as well. The metric r t−to−sol = t/e 2 shows that the same 2 accuracy of simulation result can be achieved with less computational effort when using AMR.
One of the major arguments against AMR is that mesh manipulation and management add additional computational overhead to the overall cost of the simulation. Using the metric r t−per−DOF = t/n Ω we studied the computational cost per degree of freedom and find that although adaptive simulations are more expensive per DOF, the additional cost is independent of n Ω . In fact, we find that the overhead appears to be constant per DOF.
Since AMR is known as a tool that allows for the consistent simulation of multi-scale phenomena where small scale features interact with larger scales, we studied local resolution h using the two mesh metrics r res = h/n Ω and r ef f −res = e 2 /h. These show that using AMR achieves the same 2 error with a smaller spatial resolution, i.e. at a given error, you will resolve finer scale phenomena. Furthermore, as r res shows, AMR achieves a higher spatial resolution given a fixed number of DOF. In a comparison of two different -Non-peer reviewed Earth-Arxiv preprint --Submitted to Mathematics of Computation -Beisiegel et al.
mesh refinement strategies for the quasi-stationary vortex in subsection 4.1, we furthermore found that AMR is most efficient if the coarse mesh level λ crs is kept as small (i.e. as coarse) as possible. This is because for a set of adaptive simulations with constant d λ the metric r res evolves in a same way as for the uniform simulations. Finally, we demonstrated in subsection 4.3 that for even only initially localised phenomena, we can achieve a computationally more efficient simulation using AMR. In practice we expect this to be especially important for highly time critical applications such as tsunami propagation. The presented metrics may in future help to investigate computational and numerical properties of AMR methods more rigorously. It was our intention to demonstrate their usefulness along one example code. It would now be interesting to see other AMR code developers adopt these metrics and see how different numerical schemes and different implementations expose their characteristic properties.