Multi-scale numerical simulation of a tsunami using mesh adaptive methods

Multi-scale numerical simulation of a tsunami using mesh adaptive methods An exploration of anisotropic mesh adaptivity Joe Wallwork Mathematics of Planet Earth CDT Supervised by Matthew Piggott and David Ham [email protected] Abstract Mesh adaptive methods are typically categorised as either h-adaptive or r-adaptive. In two dimensions, the former involve operations altering the number of mesh degrees of freedom by the insertion or deletion of element edges, while the latter hold both the number of degrees of freedom and the mesh topology fixed and redistribute mesh entities (vertices, edges and elements) geometrically. Anisotropic mesh adaptivity seeks to incorporate aspects of both hand r-adaptivity, providing a hybrid (hr) approach. As such, anisotropic mesh adaptivity benefits from the h-adaptive ability to completely regenerate a mesh before it gains tangled nodes, as well as the r-adaptive ability to allow degrees of freedom to follow features of fluid flow, such as a tsunami wave. A hybrid mesh adaptive approach is ideal for tsunami problems, since the important features we would like to accurately resolve are clustered in a relatively small region of ocean, which itself moves as time progresses. A standalone finite element shallow water solver is constructed for solving tsunami modelling problems, along with an anisotropic mesh adaptivity library capable of adaption both to fields related to the flow (such as fluid speed) and as guided by adjoint solution data. Introduction In 2011 a tsunami caused by a major earthquake lead to enormous damage wreaked on the coast of Fukushima, Japan, with leading tsunami waves reaching the Japanese coast after just ten minutes. Through highly efficient numerical tsunami simulations, sufficient warning could be provided in future scenarios, allowing evacuation and damage mitigation in coastal areas determined most at risk. Main Objectives 1. Build a standalone shallow water solver, using the automated finite element method code writing software provided by Firedrake. 2. Investigate avenues opened by mesh adaptivity, including guidance of the mesh adaptive process using adjoint equations, as in goal-based adaptivity. 3. Implement a mesh adaptive algorithm within the standalone shallow water solver. Experiment with different methods of Hessian reconstruction, normalisation and metric gradation. 4. Run mesh adaptive simulations for a realistic tsunami case study and make numerical comparisons between approaches, in terms of both accuracy and efficiency. Non-rotational shallow water equations We linearise about a flat surface ¯ η, typically taken as zero. For an ocean domain Ω⊂R2, define fluid velocity u: Ω →R2, free surface displacement η: Ω →Rand bathymetry b: Ω →R. Then ∂u ∂t +g∇η=0,∂η ∂t +∇·((¯ η+b)u)=0,(1) where g= 9.81m s−2denotes gravitational acceleration. Model verification experiments were carried out on a realistic ocean domain displayed in Figure 1, using the initial condition shown in Figure 2. These experiments showed very little difference in the approximations made by (1) and the corresponding nonlinear or rotational counterparts. Figure 1: Problem domain. Figure 2: Initial condition [2]. Figure 3: Gauge locations [2]. Adjoint shallow water equations Correspondingly, for adjoint variables (λu, λη):Ω→R3, we have −∂λu ∂t −b∇λη=∂J∗ ∂u,−∂λη ∂t −g∇·λu=∂J∗ ∂η ,(2) for an objective functional, J. For a spatial region of importance A⊂Ω, this is given by J(u, η) = ZTend Tstart ZA η(x, y, t) dxdydt. (3) The objective (3) allows us to consider the free surface displacement near to Daiichi nuclear power plant, say, thereby evaluating the extent to which damage is caused on the Japanese coast. Mesh adaptive process Solver H0,q0 Error model Remesher Interpolant At timestep trk, we start with mesh Hkand solution qk. Hi,qi Hi+1,qi Hi+1,qi+1 For piecewise linear (P1) approximation, the Taylor remainder theorem error result =γvT|H|v,(4) provides an error estimate upon which to base an adaptive algorithm. We compute the Hessian of a field related to the fluid flow, such as free surface displacement or fluid speed. Special reconstruction is required, such as double L2projection. Based on (4), a new mesh is generated by modifying the Hessian to obtain an error metric, which dictates the way a mesh is adapted across the domain. Results Five meshing approaches were considered: increasingly refined fixed meshes (i), (ii) and (iii); a ‘simple adaptive’ approach (iv), adapting to the free surface Hessian; and a goal-based adaptive approach (v), incorporating adjoint solution information. Time series at two offshore pressure gauges shown in Figure 3 are displayed in Figure 4 and 5, along with errors, run times and mean vertex counts. 0 5 10 15 20 25 Time elapsed (mins) −2 −1 0 1 2 3 4 5 Free surface (m) Gauge measurement Mesh approach (i) Mesh approach (ii) Mesh approach (iii) Mesh approach (iv) Mesh approach (v) Figure 4: Timeseries at gauge P02. 0 5 10 15 20 25 Time elapsed (mins) −1 0 1 2 3 4 5 Free surface (m) Gauge measurement Mesh approach (i) Mesh approach (ii) Mesh approach (iii) Mesh approach (iv) Mesh approach (v) Figure 5: Timeseries at gauge P06. Norm type Coarse Medium Fine Simple adaptive Goal-based L1norm 0.398 0.291 0.321 0.315 0.255 P02 L2norm 0.480 0.331 0.355 0.351 0.326 L∞norm 1.190 0.788 0.780 0.659 1.035 L1norm 0.324 0.128 0.111 0.163 0.106 P06 L2norm 0.460 0.172 0.136 0.225 0.152 L∞norm 1.368 0.547 0.390 0.785 0.680 Mean vertex count 3,126 17,086 97,343 5,896 3,756 Mean run time 39.9 s 239.3 s ≈4 m 1701.4 s ≈28 m 1602.8 s ≈27 m 2383.4 s ≈40 m The goal-based approach first solves (2) on a relatively low resolution, fixed mesh. Then (1) is solved adaptively, using adjoint solution data to establish pointwise a posteriori error indicators. This establishes which regions are significant and can be fed into adaptive procedure. ‘Simple adaptive’ and goal-based meshes are displayed in Figures 6 and 7, respectively. Figure 6: ‘Simple adaptive’ approach. Figure 7: ‘Goal-based’ approach. Conclusions 1. Anisotropic mesh adaptivity provides a means of maintaining both a low vertex count during tsunami modelling and a sufficiently well-resolved region surrounding important dynamics. 2. The adjoint problem can be useful for guiding the mesh adaptive process, offering additional information than in the case where only the forward problem is considered. For the T¯ ohoku tsunami, the goal-based approach out-performs the ‘simple adaptive’ method in accuracy and vertex count. Further Research I will further research anisotropic mesh adaptivity as a PhD project, including a more thorough incorporation of r-adaptivity into the approach described here. I will consider a number of other applications for mesh adaptivity in the geosciences, such as storm surges and Gulf Stream separation. Whilst my MRes computations were performed on a standalone solver of my own design, the PhD work will be part of the Thetis project, created for solving coastal, near-estuarine and ocean FEM problems. References [1] Power et al. Adjoint a posteriori error measures for anisotropic mesh optimisation. 2006. [2] Saito et al. Tsunami source of the 2011 Tohoku-Oki earthquake, Japan: Inversion analysis based on dispersive tsunami simulations. 2011. [3] Davis & LeVeque. Adjoint methods for guiding adaptive mesh refinement in tsunami modelling. 2016. GitHub github.com/jwallwork23 Centre website mpecdt.org Firedrake firedrakeproject.org Thetis thetisproject.org