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

Submitted to Imperial College London and the University of Reading in fulfilment of the requirements for the Degree of Master by Research Masters Thesis Multi-scale numerical simulation of a tsunami using mesh adaptive methods Joseph Gregory Wallwork Supervised by Prof. Matthew Piggott & Dr. David Ham. Friday 25th August, 2017 Abstract Mesh adaptive methods are typically categorised as either h-adaptive or r-adaptive. In two dimensions, the former involves 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, and is explored in depth here. This approach 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 aspects of fluid flow, such as a tsunami wave. A hybrid mesh adaptive approach is ideal for tsunami problems, since a large portion of fluid flow, which we would like to accurately resolve, is 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. By applying mesh adaptivity to shallow water problems, this work aims to efficiently generate numerical solutions to tsunami wave propagation problems. The case study of the tsunami which struck Fukushima, Japan, in 2011 is considered, wherein some leading tsunami waves reached the coast in just ten minutes. Numerical results indicate computational cost can be reduced, whilst retaining a sufficiently high accuracy, under mesh adaptivity. Through establishing a highly efficient approach to numerical tsunami simulation, sufficient warning could be provided in future scenarios, allowing for evacuation and damage mitigation in those coastal areas determined to be most at risk. I hereby declare that this Masters thesis is entirely my own work and has not previously been submitted as part of another higher degree at any other university or institution. Wherever external sources have been used, references are provided in the bibliography. Joseph Gregory Wallwork, Imperial College London, Friday 25th August, 2017 ii Contents List of Figures 1 1 Preliminaries 2 1.1 Abbreviations.................................. 2 1.2 Notation..................................... 2 1.3 Notes on computer resources . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Introduction 4 3 Mathematical formulation 6 3.1 Shallowwaterequations ............................ 6 3.2 Forwardproblem ................................ 8 3.3 Boundaryconditions .............................. 9 3.4 Adjointproblem................................. 10 3.5 Finite element problem solving by computer . . . . . . . . . . . . . . . . . 11 4 Methodology and an idealised experiment 12 4.1 Mesh-adaptiveprocess ............................. 12 4.2 Evaluation of hand r-adaptive approaches . . . . . . . . . . . . . . . . . . 14 4.3 Timeintegration ................................ 15 4.4 One-dimensional tsunami test problem . . . . . . . . . . . . . . . . . . . . 16 5 Anisotropic mesh adaptivity 18 5.1 Measuringdistance............................... 18 5.2 Gaugingmeshquality.............................. 21 5.3 Metriccomputation............................... 22 5.4 Hessianrecovery ................................ 24 5.5 Adapting to multiple solution fields . . . . . . . . . . . . . . . . . . . . . . 25 5.6 Metricgradation ................................ 26 5.7 Notes on the adaptive algorithm code . . . . . . . . . . . . . . . . . . . . . 28 5.8 Adaptivitycodetests.............................. 29 5.9 Goal-based mesh adaptivity . . . . . . . . . . . . . . . . . . . . . . . . . . 31 iii 6 Tsunami application 34 6.1 Computational setup for the T¯ohoku tsunami . . . . . . . . . . . . . . . . 34 6.2 Modelverification................................ 36 6.3 Damagequantification ............................. 38 6.4 Modelcomparisons ............................... 39 6.5 Conclusions and future work . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Bibliography 43 Acknowledgements I greatly thank my supervisors Prof. Matthew Piggott and Dr. David Ham for their continued help, support and advice throughout my undergoing of this project. They were always happy to explain concepts I didn’t fully understand and to offer suggestions on how I might proceed with my research. Especial thanks also go to Dr. Alexandros Avdis, Prof. Tatsuhiko Saito and in particular Dr. Nicolas Barral for their extensive knowledge sharing and resource providing. Further thanks to Dr. Stephan Kramer, Dr. Lawrence Mitchell, Dr. Ethan Kubatko, Thomas Gibson and Simon Warder for their sound suggestions relating to computational aspects of the project. Last but certainly not least, thank you to Dr. Anna Radomska for her continued administrative support and always being happy to answer questions and try to resolve issues. List of Figures 2.1 “The Great Wave off Kanagawa” by Hokusai, reinterpreted in the context ofadaptivemeshing. .............................. 4 2.2 Tide and pressure gauge locations near to Fukushima, Japan. . . . . . . . . 5 3.1 Setup of the shallow water equations. . . . . . . . . . . . . . . . . . . . . . 7 4.1 The mesh-adaptive process. . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4.2 Bathymetry and initial conditions for an idealised 1D tsunami simulation. . 16 4.3 Significant regions of ocean in an idealised 1D tsunami simulation. . . . . . 17 5.1 The mapping of an ellipse to the unit circle under a metric. . . . . . . . . . 19 5.2 Meshes adapted with respect to two different scalar functions. . . . . . . . 29 5.3 Mesh adaptive advection and diffusion of a Gaussian source under a constantwindfield.................................. 30 5.4 Mesh adaptive shallow water simulation across a shelf break discontinuity withaGaussiansource. ............................ 31 6.1 Domain geometry used in modelling the T¯ohoku tsunami. . . . . . . . . . . 34 6.2 Initial free surface displacement used in modelling the 2011 T¯ohoku tsunami. 35 6.3 ‘Eyeball norm’ comparisons of the non-rotational and rotational, linear and nonlinear standalone solvers on a fine mesh for the T¯ohoku problem. . . . . 36 6.4 Model verification using time series of free surface measurements and predictions at pressure gauges off the coast of Japan. . . . . . . . . . . . . . . 38 6.5 Damage measure quantification for the T¯ohoku tsunami, alongside the spatial region of importance considered. . . . . . . . . . . . . . . . . . . . . . 39 6.6 A time series comparison of different meshing approaches at pressure gauges off the coast of Japan, alongside error norm values for five approaches to meshing in tsunami modelling. . . . . . . . . . . . . . . . . . . . . . . . . . 40 6.7 Meshes generated after 14 minutes of simulation time, with and without guidance by adjoint information, for the T¯ohoku tsunami problem. . . . . . 41 6.8 Time performance analyses for various meshing approaches to tsunami modelling the T¯ohoku tsunami. . . . . . . . . . . . . . . . . . . . . . . . . 42 1 1. Preliminaries 1.1 Abbreviations AMR: adaptive mesh refinement. DOF: degrees of freedom. FEM: finite element method. GIS: geographical information system. PDE: partial differential equation. SPD: symmetric positive-definite. SW: shallow water. SWEs: shallow water equations. UFL: unified form language. UTM: universal transverse mercator. 1.2 Notation The following items establish all pieces of notation used in this project which are not entirely standard in the literature. •N0:= N∪{0}={0,1,2,3, . . . }denotes the natural numbers with zero included. •For functions f:A→Rand u:A→Rnon A⊆Rn, ∂f ∂u:= h∂f ∂u1. . . ∂f ∂uniT:A→Rn,with the gradient ∂f ∂x=∇f. (1.1) •For functions f:A→Rmand u:A→Rnon A⊆Rn, ∂f ∂u:=     ∂f1 ∂u1. . . ∂f1 ∂un . . ..... . . ∂fm ∂u1. . . ∂fm ∂un    :A→Rm×n,(1.2) provides the Jacobian matrix for the transformation from variables fto variables u. The Jacobian ∂f ∂x=J(f),of fprovides a common example. •If a matrix A∈Rn×nis positive definite then we write A0. 2 •Given a metric space (X, d), a subset U⊆Xthereof has closure U={x∈X| ∀ > 0,∃y∈U|d(x, y)< }, interior U◦=X\(X\U) and boundary ∂U =U\U◦. •We denote by Inthe n×nidentity matrix, for n∈N. •Meshes are denoted by Hand contain three types of entity: (triangular) elements K⊂R2,vertices p∈R2and edges pq, which are often also considered as vectors in R2. •Given some domain, Pkdenotes a degree-kpiecewise polynomial function space defined over the domain, where k∈N. The same notation is used in vector and scalar cases, but the distinction is made clear. 1.3 Notes on computer resources The majority of the computer code used in this project is written in Python, with incorporation of UFL in order to deal with forms. Additional computer packages used are listed below, along with web links for their download and associated documentation. Firedrake: http://firedrakeproject. org GMSH: http://gmsh.info GMT: http://gmt.soest.hawaii.edu Paraview: http://www.paraview.org PlotDigitizer: http://plotdigitizer. sourceforge.net PRAgMaTIc: https://github.com/ meshadaptation/pragmatic QGIS: http://www.qgis.org QMESH: [Avdis et al.,2017] Thetis: http://thetisproject.org In addition, we make use of two geographic databases, referenced below. GSHHG Global Self-consistent, Hierarchical, High-resolution Geography: http://www.soest.hawaii.edu/pwessel/gshhg/ GEBCO General Bathymetric Chart of the Oceans: http://www.gebco.net All Python code for this project can be found on GitHub, at the web address https://github.com/jwallwork23/MResProject. There can be found a number of Python scripts, whose roles are described in Section 5.7, along with directories containing data from the above sources and also plots.1 1Everything needed to run the code is included in the directories of the GitHub page, except for the mesh files. This is for copyrighting purposes, but the meshes can be provided upon request, with the permission of the QMESH developers. 3 2. Introduction Figure 2.1: Interpretation of an adaptively meshes tsunami. Original image “The Great Wave off Kanagawa” by the Japanese artist Katsushika Hokusai (1830-1833) edited by the author. The word tsunami derives from Japanese, where ‘tsu’ refers to the harbour and ‘nami’ refers to a sea wave [Lisitzin,1974]. This language is suggestive of the serious damage that these waves wreak upon their arrival into harbours across Japan. In recent and ancient history, and in particular in Japan, tsunamis have proved to be natural disasters with devastating impacts on human civilisations. This kind of damaging tsunami are common in the seismically active Pacific rim near Japan’s coast, and as a country comprising a collection of small islands, Japan’s coastal cities particularly vulnerable. One of the most memorable tsunamis of recent times, and indeed the most powerful one to strike Japan on record, had its epicentre off the coast of T¯ohoku and struck many parts of the country in 2011, including the Fukushima region. The tsunami caused 15,893 deaths, much devastation to hundreds of thousands of homes and played a large part in causing the level 7 meltdown of the Fukushima Daiichi nuclear power plant.1Although the meltdown event itself caused no deaths, much disruption was caused; an enormous clean-up operation was required, and there was a mass-evacuation of over 160,000 people who lived within a 20 km radius, wherein the yearly dosage of radiation was expected to reach dangerous levels of 20 mSv year−1in the years following the accident.2Only now, six years after the incident, is the radioactivity low enough to begin the clean-up, and still some robots are unable to handle the conditions. Preliminary and later reports of the earthquake, tsunami and following meltdown 1Details obtained from the National Police Agency of Japan: http://www.npa.go.jp/archive/keibi/biki/higaijokyo_e.pdf. 2For further details, see http://fukushimaontheglobe.com. 4 events include those by [Kazama and Noda,2012], [Okada,2011] and [Simons et al.,2011], with further analyses are provided by [Baba et al.,2017] and [Suzuki et al.,2012]. There is an obvious desire to avoid such occurrences in the future, and the Japanese government have even motioned to review its use of nuclear power in its energy grid due to the high risk of these extreme events. Yet in the light of the enormous problem posed by climate change, we should not be too hasty to abandon this relatively clean alternative to carbon dioxide-emitting fossil fuels which has the potential to generate vast amounts of electricity. An alternative approach to the problem is to make tsunami early-warning systems much more efficient. Figure 2.2: Gauge locations near Fukushima, courtesy of [Saito et al., 2011]. The star indicates the earthquake epicentre and small circles indicate aftershock locations. This problem will be considered in this thesis, making use of adaptive meshes in numerical ocean modelling. Mesh adaptivity aims to discretise the domain of ocean of interest in an unstructured way, this not only varies in resolution spatially, but also adapts temporally. In this way, fluid dynamics near to the tsunami wave can be resolved more accurately, whilst minimising the computational effort expended on solving the flow equations in locations distant from both the wave and vulnerable coastal civilisations. The result is an efficient and yet relatively computationally cheap approach. The mathematical model to which we apply the adaptive process is provided by the shallow water equations, which are wellknown to provide a good approximation of large scale oceanic fluid dynamics. These equations may be solved using the finite element method which, unlike the finite difference method, has the advantage of working well on unstructured grids. The central aim of this project is to construct an effective and efficient means of modelling the 2011 T¯ohoku tsunami, which arrives at a solution within the 10 minute timeframe before the waves were first felt, and which sufficiently accurately resolves the coastal bound tsunami wave. Thus, by implementing such a method in the future, vulnerable regions may be given sufficient time to evacuate and damage may be mitigated. Throughout, we refer to the case study tsunami as the T¯ohoku tsunami. The structure of this thesis is as follows. Chapter 3 outlines the mathematical setup, including equations and initial and boundary conditions used. Chapter 4 then goes on to describe the solution approach to solving problems of the form described in Chapter 3 and outlines in brief the general procedure of anisotropic mesh adaptivity, with further detail provided in Chapter 5. Chapters 4 and 5 also include plots and results from the non-adaptive and adaptive algorithms applied to model tsunami problems, respectively. Finally, Chapter 6 uses the content of the preceding chapters to apply anisotropic mesh adaptivity to the T¯ohoku tsunami case study, including some model verification and result comparisons with alternative solution approaches. 5 4. Methodology and an idealised experiment In this chapter we outline and evaluate the general concept of anisotropic mesh adaptivity, which is explained in more detail in Chapter 5, along with the framework within which we solve PDEs using this process. In addition, we consider a model 1D problem to interpret the adjoint problem concept, hinting at how this could prove useful in guiding adaptivity. 4.1 Mesh-adaptive process A fixed mesh defined on a domain Ω has a spatially fixed set of DOFs. If these DOFs are not suitably distributed, or if the problem varies temporally, computational inefficiencies may become apparent. Whilst aiming to avoid such inefficiencies, we would also like to retain a suitably accurate approximation to the true dynamics. As such, it makes sense to consider a finer mesh where our approximation is deemed as poor and a coarser mesh where it is already of a high quality. This is precisely the approach of mesh adaptivity. Of course, we generally do not know the true error made in our approximation. Hence, to perform mesh adaptivity, we must have a means of error estimation. 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 Figure 4.1: The mesh-adaptive process. Moreover, an adaptive mesh varies in time as well as space, meaning the mesh used may change across timesteps, as illustrated in Figure 4.1. The sequence {rk}k∈N⊆N0determines those timesteps in which the mesh is regenerated, which could also be selected adaptively, but which is incremented uniformly in this work. The ability to work with adaptive meshes is particularly attractive for tsunami modelling, where the region of interest (surrounding the coast-bound tsunami wave) moves rather rapidly as time progresses. The task is to ensure both that the path of the tsunami is sufficiently resolved and also the part of the domain containing the wake of the wave becomes more coarsely meshed as it moves away, reducing computational cost. 12 Additionally, we would like to mesh finely near to locations of interest, such as densely populated areas, busy harbours and nuclear power plants. The mesh-adaptive process contains three crucial steps: error estimation, mesh adaption and interpolation, as described in the following subsections. 4.1.1 Error estimation In solving PDEs using FEM we often use P1 approximation, as for the free surface displacement in the SWEs in Chapter 3. For this function space, FEM theory as in [Brenner and Scott,2007] provides the Taylor remainder theorem approximation error result A=γvT|H|v,(4.1) where γ=O(1) is a scalar, v∈Rnis a vector corresponding to a direction and magnitude and His the Hessian of the solution field. As such, it makes sense to consider error estimates for the adaptive process based on the Hessian, described in detail in Chapter 5. 4.1.2 Mesh adaption Having constructed error measures, the next task is to modify the mesh accordingly, providing the mesh adaption. There are various flavours of adaptivity, each with a different approach to this mesh adaption stage. Approaches include p-adaptivity,h-adaptivity and r-adaptivity. We are mainly interested in the latter two approaches. Adaptive mesh refinement (AMR) is a popular type of h-adaptivity. For a summary of the approaches of AMR and p-adaptivity, see [Behrens and Bader,2009] and [Giorgiani et al.,2013], respectively. The cases of hand r-adaptivity are evaluated in Section 4.2. It is also possible to consider combinations of some of these approaches, as discussed in Section 4.2. While the purpose of the previous stage was to flag parts of the mesh for editing, the adaption phase seeks to perform these modifications. 4.1.3 Interpolation of solution Once mesh adaption has been implemented, the final task is to interpolate the solution from the ‘old mesh’ onto the ‘new mesh’, providing another opportunity for accruing error. As such, central to h-adaptive methods is the concept of interpolation error I=kq−Πhqk,in some norm k·k:X→R,(4.2) where again X={f: Ω×(0, T]→R3} 3 q. Interpolation error is made upon approximating the exact solution qby some interpolant Πhqover the mesh, in our case a Lagrange interpolant. Interpolation error does not arise in r-adaptivity, for reasons discussed in Section 4.2. 13 4.2 Evaluation of hand r-adaptive approaches As detailed in Section 5.2, h-adaptivity gauges mesh quality via a functional (5.15) and adapts according to a metric. This SPD tensor field depends on the Hessian of the solution field, in accordance with (4.1). In h-adaptivity, mesh entities may be created or destroyed, thus altering the mesh topology, whilst r-adaptivity fixes the mesh topology, with entities simply moving around the domain, as determined by another functional. The h-adaptive procedure repeatedly evaluates the quality functional following various geometrical and topological operations upon the mesh, aiming towards its minimisation, whereby the mesh has desirable properties. Historically, h-adaptivity is usually considered on Cartesian meshes and the ‘h’ refers to the hwhich usually denotes the size of elements used in FEM. Altering mesh topology can require an excessive amount of interpolation and pose an issue for our problem, since the region surrounding the tsunami, which requires finer resolution, moves with the flow. Excessive interpolation contributes to a heightened computational cost, which is completely against our motive of pursuing adaptive meshing in order to create a more efficient algorithm. This motivates considering the idea of a moving mesh which is able to chase the wave, leaving coarser mesh in its wake. The approach of r-adaptivity is indeed often referred to as mesh movement and can provide a more effective approach for cases in which regions of increased resolution move with the flow. The ‘r’ of r-adaptivity comes from the relocation of the mesh nodes, which can be achieved in a number of ways, discussed in [Piggott et al.,2005]: a variational approach makes use of the Euler-Lagrange equations; it is also possible to apply a socalled mesh smoothing algorithm in order to achieve mesh movement. As noted in [Piggott et al.,2005], the fact r-adaptivity moves the mesh with the flow has the result of reducing large transport velocities, such as those associated with the tsunami. This means there are fewer timestep restrictions than for application of h-adaptivity, providing a clear numerical advantage. On the other hand, unlike with h-adaptivity where new nodes can be created, r-adaptive methods hold the DOF count constant, potentially leading to problems if the dynamics of the problem at hand become rather complex. In these situations, it would be advantageous to be able to ‘inject’ mesh resolution into a certain region – something which is not available in r-adaptivity. In addition, r-adaptivity can produce problematic, so-called ‘tangled meshes’, whereby two or more elements overlap. Clearly, there are disadvantages to the approaches of both hand r-adaptivity, but it appears that many of these are complemented by corresponding advantages for the other approach. As suggested above, it is both desirable and possible to combine the two, yielding a hybrid, hr-adaptive method. One such approach, as described in [Habashi et al.,2000], is to incorporate a moving mesh within a mesh optimisation method such as described earlier. Alternately, h-adaptivity can be incorporated within the structure of moving meshes on a local scale, as in [Lang et al.,2003], to ensure elements are not poorly shaped and neither is the specified error tolerance reached. 14 4.3 Time integration As discussed in Chapter 3, solution of the SWEs requires timestepping, with FEM solves at each iteration. Two candidate timestepping schemes are described in the following. For the purposes of these analyses, consider a simple PDE problem ∂u ∂t =f(u(t), t), u(0) = u0(4.3) defined on some domain Ω ⊆R, for a dependent variable u: Ω →R, forcing term f: Ω →Rand initial value u0∈R. Suppose we have a uniform time discretisation {u(n)}n∈N0of the dependent variable, with (constant) timestep ∆t > 0. One approach is the implicit Euler method, whereby the value of uat the current timestep depends on information available only at the current timestep: u(n+1) −u(n) ∆t=f(u(n+1), t(n+1)).(4.4) Whilst this equation is not as straightforward to solve as with explicit methods, amounting to solving a matrix system, it has the advantage of having a large region of absolute stability. However this approach is only O(∆t) at a particular time tand does not conserve energy. We would prefer to use a higher order, conservative time integrator. In (4.4), the forcing function fis evaluated only at current time values. Another approach is to evaluate at an intermediate time, say halfway between the previous and current timesteps. From this we derive the implicit midpoint rule, given by u(n+1) −u(n) ∆t=fu(n+1) +u(n) 2, t(n)+∆t 2.(4.5) Whilst more computationally expensive than (4.4), this approach has a global error of O(∆t2), meaning the error decays quicker as we consider smaller timesteps. Further, (4.5) is a symplectic integrator in the context of Hamiltonian dynamics, so conserves energy. For explicit time integrators, the spatio-temporal discretisation used must satisfy the Courant-Friedrichs-Levy (CFL) condition (established in [Courant et al.,1928]), c:= ∆t|u| ∆x+|v| ∆y≤1,(4.6) for stability. Since (4.4) and (4.5) are implicit integrators, this condition can be broken without incurring penalties of reduced stability. However, it is still useful to make use of (4.6) when choosing a timestep length. We showed in Section 3.1 that the wave speed of a linear SW wave is given by pg¯ h, where ¯ his the water depth when at rest. As such, for the SWEs, |u|,|v| ≤ √gbmax, where bmax denotes the maximal bathymetry in Ω. In selecting a timestep length, we should ensure it is sufficiently less than the ratio of the minimum tolerated element size and twice this upper bound for the component speeds. 15 4.4 One-dimensional tsunami test problem As a first model SW problem, consider the idealised 1D tsunami problem described in [Davis and LeVeque,2016]. This problem concerns propagation of an initial surface profile caused by a tsunami across a simplified ocean domain of width 400 km. The domain bathymetry bhas a shelf break discontinuity, as illustrated in Subfigure 4.2a: b(x) = (200 m x≤50 km 4,000 m x > 50 km .(4.7) Illustrated in Subfigure 4.2b, the initial condition used is given by η0(x) = (0.4 sin (x−100,000)π 50,000 m 100 km ≤x≤150 km 0 m otherwise .(4.8) As in Section 3.3, (unrealistic) impermeable boundaries are considered, whereby u(4 · 105) = u(0) = 0. This means waves simply reflect off both the coastal and open ocean boundaries. The 1D form of the SWEs used in [Davis and LeVeque,2016] take an alternative linear form µt+g¯ h(x)ηx= 0, ηt+µx= 0,(4.9) where µ=hu represents momentum and we have (again) linearised about a flat surface ¯ηwith zero velocity ¯u= 0. As with (3.3), (4.9) may be written in matrix-vector form qt+A(x)qx=0,where A(x) = "0g¯ h(x) 1 0 #and q="µ η#.(4.10) In 1D tsunami test.py, we also consider the adjoint problem. Given the equation of (4.10), as outlined in [Davis and LeVeque,2016], test by an appropriately sized vector functional λand integrate over the spatio-temporal domain to obtain Z4·105 0Z4200 0 λT(qt+Aqx) dtdx= 0.(4.11) (a) Bathymetry (b) Forward initialisation (c) Adjoint initialisation Figure 4.2: Bathymetry and ‘initial’ conditions used for the 1D tsunami test problem. 16 Integrating by parts in both space and time, Z4·105 0 (λTq)4200 0dx+ZT 0 (λTAq)4·105 0dt−Z4·105 0Z4200 0 qT(λt+ (ATλ)x) dx= 0.(4.12) We consider an objective functional only at the end time, with Tstart =Tend =T, so the integrand of the third term of (4.12) contains the adjoint equation for this problem, λt+ (AT(x)λ)x= 0,where λ=hλµληiT .(4.13) We enforce an ‘initial condition’ of the adjoint problem, as a scaled indicator function λη,0(x)=0.4 1 [10 km,25 km], which approximates a corresponding delta function, as illustrated in Subfigure 4.2c. Again, we impose the boundary condition λu(4·105) = λu(0) = 0. Since the adjoint equation is solved backwards in time, we transform time by t7→ 4200−t, whence ∂t7→ −∂t. Under the same boundary conditions, (4.13) becomes λt−(AT(x)λ)x= 0.(4.14) As in [Davis and LeVeque,2016], we interpret these data by establishing regions of ocean where the forward and adjoint free surface solutions take numerical values with magnitude above some tolerance, say 0.05m, with resulting plots displayed in Subfigures 4.3a and 4.3b, respectively. The plots include a dashed line marking the location, 50 km offshore, of the shelf break discontinuity in bathymetry. Significant regions of the forward problem can be interpreted as areas of ocean to which the initial condition will propagate, whilst the adjoint significance regions, on the other hand, correspond to some approximation of the ‘domain of dependence’ wherein tsunami wave propagative activity will have an effect on the near-coastal area of importance at the end-time of the simulation. Another interpretation may be gleaned by considering regions of ocean wherein the inner product between forward and adjoint problem solutions is significant, as displayed in Subfigure 4.3c. Thus, we gain information concerning which tsunami trajectories are important for our particular problem, giving a first suggestion of where mesh adaption should be performed. Notice the only important regions lie between coast and shelf break. (a) Primal solution η(b) Adjoint solution λη(c) Inner product qTλ Figure 4.3: Regions of ocean where where orward and adjoint shallow water solutions (and the inner product thereof) have magnitude at least 0.05 m for the 1D tsunami test problem. 17 5. Anisotropic mesh adaptivity When a mesh is regenerated using anisotropic mesh adaptivity, the goal is to improve the quality of the worst element of the mesh, in a local sense. In 2D, this is achieved by a performing some combination of the following operations. 1. edge splitting; 2. edge collapsing; 3. edge swapping; 4. node movement. In the former three operations, the approach embodies that of h-adaptivity and in the latter operation it embodies that of r-adaptivity, thereby providing a variant of hr-adaptivity. An optimisation algorithm loops over the nodes of the mesh, and uses the above operations to proposes a new local mesh configuration, based on error measures as mentioned in Subsection 4.1.1. If certain criteria are satisfied, insisting the mesh quality is improved sufficiently across meshes (so it is worth the computational effort associated with regenerating the mesh), the new configuration is accepted. The means of measuring mesh quality is described in Section 5.2, with an alternative approach found in [Pain et al.,2001]. Further intricacies of the mesh adaptive process are also outlined in this chapter. 5.1 Measuring distance Ascalar product is defined as an SPD form, which can be represented by an SPD matrix M, known as a metric. Henceforth in this project, we work solely in 2D and so all meshes Hcan be understood as subsets of R2. As such, for our purposes, a metric provides a map h·,·iM:R2×R2→[0,∞),hx,yiM=xTMy,x,y∈R2.(5.1) Euclidean space E2is obtained by equipping a vector space with such a scalar product. Given the scalar product defined by (5.1), we obtain the corresponding norm by k·kM:R2→[0,∞),kxkM=√xTMx,x∈R2.(5.2) Lengths of mesh edges pq ∈ H can be calculated using (5.2). As commented in [Barral, 2015], angles between vectors can also be computed in the usual way, with the angle θ∈[0,2π) between two vectors x,y∈R2with respect to a metric Mgiven by cos θ=hx,yiM kxkMkykM .(5.3) 18 Further, given an element Kof the mesh has area |K|I2with respect to E2, its corresponding area in the metric space defined by Mcan be found using the determinant: |K|M=pdet(M)|K|I2.(5.4) As Mis symmetric, it is orthogonally diagonalisable, with eigenvalue decomposition M=VTΛV="u1v1 u2v2#"λ10 0λ2#"u1u2 v1v2#(5.5) where λ1, λ2>0, and eigenvectors u= [u1, u2]Tand v= [v1, v2]Tare orthonormal. (1,0) (0,1) M1 2 M−1 2 (h1,0) (0, h2) E2= (R2, I2) (R2, M) Figure 5.1: The mapping of an ellipse to the unit circle under a metric. Combining the above ingredients, we come to an alternative understanding of a metric Min terms of geometry. That is, Mcan be described in terms of ellipses, with its eigenvectors determining the axes in which the ellipse is skewed and the inverse square root of its eigenvalues determining the magnitudes of skew along these axes. In the process of mesh adaption, we seek to obtain a uniform mesh with respect to a metric. Let h1and h2denote the magnitudes of the ellipse axes. Then, for the unit vector e1= [1,0]T∈R2, kh1e1k2 M=hh10iVTΛV"h1 0#=h2 1λ1= 1 ⇐⇒ λ1=1 h2 1 ,(5.6) using (5.5), and similarly for the second standard unit vector, e2= [0,1]T. The metric achieves its purpose of controlling the mesh adaption process by defining the mesh edge lengths we desire at each node of the mesh. In tune with (5.6), for each eigenpair (λi,vi) of M, we take the corresponding edge length in direction vi, as determined by hi=1 √λi .(5.7) Figure 5.1 depicts an ellipse with axes in directions given by the eigenvectors and magnitude given by the inverse square eigenvalues. Application of the square root M1 2of the metric then yields a unit circle (with axes being the standard unit vectors). Here the matrix square root is well-defined since all eigenvalues are positive, meaning we can take the square root through the eigenvalue decomposition, M1 2=VTΛ1 2V=VTdiag(pλ1,pλ2)V. (5.8) To contrast, in isotropic mesh adaptivity, metrics Mare not only SPD, but diagonal, 19 taking the form M= diag(h−2, h−2) for some edge length h > 0. As such, the isotropic approach allows only changes in the size of elements, and not their orientation. The scalar product provided in the Euclidean case by (5.1) is constant across the domain, so distances are measured using the same ellipse, regardless of location. As remarked in [Barral,2015], it is attractive to be able to measure distances in a way dependent on the location within the domain. This is particularly so for tsunami propagation modelling, enabling us to consider mesh elements which have different sizes, using a fine mesh where high resolution is required and a coarse mesh in ‘less important’ regions. For spatially dependent metrics, we consider Riemannian, rather than Euclidean, geometry. We now consider not just one SPD matrix M, but a space M={Mx}x∈Ωthereof, defined on all points of the domain. Locally, each of these metrics defines a scalar product when evaluated on a point x∈Ω. Using these scalar products, we can establish a Riemannian metric space, to which we can extend notions of distance and angle from the Euclidean case. The following extends definitions as in [Pain et al.,2001]. For an element Kof the mesh, consider one of its edges e=pq ∈EK, with endpoints p,q∈R2. The length of ewith respect to metric Mis calculated analytically as `M(e) = Z1 0ppqTM(p+tpq)pq dt≈ k X i=1 ωippqTM(p+αipq)pq,(5.9) using the notation of [Alauzet,2010]. While the strict equality in (5.9) provides a continuous measure of distance along an edge, we need consider a discretised metric Mh={Mp}p∈H, as evaluated at the vertices of the mesh H. Thus we consider a k-point quadrature rule with weights {ωi}k i=1 and Gauss points {αi}k i=1, as in the approximate equality of (5.9). In P1 space, there is linear variation within each element. As such, we consider linear variation of the discretised metric along each edge of the mesh, meaning it is sufficient to calculate length using one-point quadrature, as in [Piggott et al.,2005]: `Me(e) = kpqkMe=ppqTMepq,where Me=1 2(Mp+Mq) (5.10) denotes the edge-centred metric, summed at the endpoints of e. That is, α1=1 2is the single Gauss point. Vertex-wise metrics are calculated using interpolation between meshes. Since M(x)0,∀x∈Ω, in particular Mp, Mq0. Application of the definition of positive-definiteness implies Mp+Mq0, meaning the edge-centred metric is also positive-definite, for each edge of the mesh. Similarly to (5.4), the area of an element Kin Riemannian metric space M={Mp}p∈Ω is given by the integral |K|M=ZKpdet(M(x, y)) dxdy. (5.11) As remarked in [Barral,2015], this may be approximated to first order using (5.4), with 20 Mreplaced by Mapplied at the barycentre of K. With the same modification, (5.3) can be used to calculate the angle between vectors u1and u2in the Riemannian metric space. As already mentioned, the idea of metric-based mesh adaption, as first introduced in [George et al.,1991], is for the Riemannian metric space to be used in computing geometrical quantities which arrange for the new mesh to be a unit mesh with respect to this metric space. This links with the earlier interpretation of moving from an ellipse to a unit circle, in terms of the metric considered. As defined in [Alauzet and Loseille,2016], Kis a unit element if each of its sides have unit length with respect to the governing metric. Hence, K∈ H is a unit element if and only if Kis equilateral with all sides equal to 1. By basic trigonometry, such a triangle has area √3 4with respect to M, and hence Euclidean area √3 4(det(M))−1 2. If we are to consider unit meshes as consisting purely of unit elements, it is clear most domains cannot be filled particularly well. For example, in Euclidean space, where M=I2, a square cannot be filled with equilateral triangles without a number of gaps. It can, however, be filled using right angled triangles. Due to the above, it is necessary for us to relax the constraint that our meshes contain only unit elements. Instead, as in [Barral,2015], we consider a mesh which contains only quasi-unit elements, whose edges each have length in the range [ 1 √2,√2]. In the Riemannian setting, (5.7) carries over to provide a measure of size, hM(v) = kvk2 `M(v),(5.12) with respect to v∈R2, which may be approximated using quadrature as in (5.10). 5.2 Gauging mesh quality Consider a vector J∈RN, where Ndenotes the number of elements on the current mesh, which will likely change value a number of times during a h-adaptive algorithm. We follow [Piggott et al.,2005] in measuring global mesh quality using the ∞-norm, F=kJk∞.(5.13) The use of an ∞-norm means the mesh as a whole is held to have the quality of the worst quality element of the mesh, which we would like to improve, in a local sense. The functional we consider coincides with the one used in PRAgMaTIc (Parallel anisotRopic Adaptive Mesh ToolkIt).1, which provides anisotropic mesh adaptivity in Firedrake for meshes of simplexes and which underpins the code written in this project. In 1Quality functional information described in this section was collected from the associated GitHub page referenced in Section 1.3, in the introduction.tex file of the docs directory. 21 we define the reduced Riemannian space using metric intersection defined in Section 5.5, f M=n∩p∈H f Mp(x)∩f M(x)ox∈Ω.(5.37) As such, we enforce the strongest constraint on element size across metrics in (5.37). Using linear interpolation of the metric, we need only compute (5.37) at each vertex. However, [Alauzet,2010] states that the associated algorithm is of quadratic complexity, so is computationally cumbersome. Instead, we follow the approach of only evaluating intersections over individual edges. For an edge pq ∈ H, we make the approximation f M(p)≈ M(p)∩Mq(p),f M(q)≈ M(q)∩Mp(q).(5.38) 5.7 Notes on the adaptive algorithm code Anisotropic mesh adaptivity functions used in this project are not currently implemented on the master branch of Firedrake. For this reason, we use a different branch3, which will soon be merged into the master. As discussed in Section 4.1, there are three central parts of the adaptive procedure: error estimation, mesh adaption and interpolation. In the following we briefly elaborate on the code required for each of these steps. As discussed in Section 5.3, for computational purposes, error estimation amounts to computing the Hessian of a scalar field. This is achieved by the function construct hessian. The two main modes of Hessian recovery described in Section 5.4 are both implemented in this function, with selection between them available. Having computed a Hessian, the corresponding metric is computed using compute steady metric, making modifications described in Section 5.3. Control of the minimal and maximal element sizes tolerated, as well as the maximal aspect ratio a, is allowed. By allowing larger values of a, we may increase mesh anisotropy. Intersection of metrics by (5.30) is achieved by metric intersection. Functions described here are found in the utility script adaptivity.py The crucial adaptivity construct used in this code is the AnisotropicAdaptation class. Members thereof take as input a mesh and a metric and have two important functions: adapted mesh and transfer solution. The former holds the adapted mesh, under the metric supplied, whilst the latter interpolates fields from the old mesh to the new mesh. As such, the mesh adaption process is encapsulated in adapted mesh. Interpolation using transfer function is only currently supported for Lagrange spaces, and not mixed spaces. Since we consider a mixed Taylor-Hood problem, our purposes require a modification. This is provided in the script interp.py, as the function interp Taylor Hood. As well as extending to the mixed case, this function involves a procedure accounting for when a metric dictates a node is moved outside the domain. An implementation of metric gradation is provided by the function metric gradation. 3Found at https://github.com/taupalosaurus/firedrake/ 28 5.8 Adaptivity code tests Before considering time-dependent PDE problems, which require several remeshing steps over the solution process, we test the adaptive algorithm for some steady problems. The following surfaces considered were introduced in [Olivier,2011]: u1(x, y) = x2+y2, u2(x, y) = tan−10.1 sin(5y)−2x+ tan−10.5 sin(3y)−7x.(5.39) Using sensor tests.py, we are able to generate meshes adapted to (5.39), as displayed in Figure 5.2. In each case double L2projection is used for Hessian reconstruction. The algorithm appears to adapt well to the given sensor functions, with the first example yielding a near uniform mesh in the domain interior, as expected due to ∇∇Tu1= diag(2,2) being constant, although there are some issues on the boundaries. Let us now consider some temporally varying PDE problems. A simple such problem to consider is the advection and diffusion of a pollutant concentration φ, say, under the influence of a constant wind field, u= (1,0) m s−1. Considering a diffusion term with a non-negligible, but still relatively small diffusivity parameter ν > 0, helps to smooth out any small scale spurious structures in the problem. Burgers’ equation is given by ∂φ ∂t +u·∇φ−ν∇2φ= 0.(5.40) (a) Sensor u1(b) Mesh 1 (hmin = 10−3) (c) Sensor u2(d) Mesh 2 (hmin = 10−4) Figure 5.2: Test sensor functions (5.39) and their associated adapted meshes. The initial mesh was uniform with 200 vertices in each direction. Maximal element size and anisotropy were set as hmax = 0.1 and amax = 103, with minimal element size as given. 29 Consider a rectangular domain Ω = [0,4] ×[0,1] m2and a Gaussian initial condition φ(x, y)=0.001 exp(−25((x−α)2+ (y−β)2)),(5.41) with α=β= 0.5. Code for this test case is provided by Burgers test.py. Choosing hmin = 5 mm and hmax = 100 mm as defaults, the algorithm adapts the mesh in such a way that the vast majority of the resolution is focused around the bubble of high concentration, as it is advected along the wind direction, as illustrated in Figure 5.3 and in accordance with expectations. In this experiment we take the diffusivity parameter ν= 10−3m2s−1and perform Hessian reconstruction by double L2projection. In the above, adaption is performed with respect to a single scalar dependent variable, φ, of a scalar PDE (5.40). Incorporating adaptivity into the SWEs is a little more difficult, as these comprise a “2+1” pair of coupled PDEs as in (3.1), with vector and scalar dependent variables.4 For the next test problem, we take a step towards the realistic tsunami case by considering a SW problem with a shelf break discontinuity, as in Section 4.4. Consider the extension of our rectangular domain to a square domain Ω = [0,4] ×[0,4] m2, with a shelf break from 1 cm depth down to 10 cm depth, occurring 50 cm from the left-hand boundary and with the same initial condition as in (5.41), but now with α=β= 2, so that the Gaussian bell is centred within the domain. The script file of simple adaptive SW.py adaptively solves the linear SWEs (3.2) for this setup. As is illustrated in Figure 5.4, 4The SWEs can also be interpreted as a system of three PDEs with three scalar dependent variables or, in the linear case, as a single vector PDE (3.3). Figure 5.3: Advection and diffusion of the initial concentration (5.41) for α=β= 0.5, according to Burger’s equation (5.40), with wind field u= (1,0), purely in the x-direction. 30 Figure 5.4: Linear SW simulation for initial condition (5.41) with α=β= 2 and a shelf break discontinuity. the mesh resolution successfully ‘follows’ the rings propagating outwards from the initial spreading of the initial condition. As we saw in considering the adjoint problem, increased mesh resolution in the latter subplot indicates that the most significant region of the domain for this simulation is the part in the shallow water area ‘near to the shore’. In the SW test case considered here, we adapt only to the fluid speed, not the free surface displacement, or individual components of the fluid velocity. Using the theory as discussed in Section 5.5, and metric intersection, we can adapt to these fields, too. There is, however, an issue with the meshes generated in Figure 5.4. Examining them more closely, it becomes apparent that the denser patches of mesh resolution are often ‘lagging behind’ the crest of the wave itself. This occurs mainly because the mesh is not adapted at every timestep and, when it is, is adapted based on a Hessian reconstructed from the current timestep. In this way, the regions flagged as requiring finer mesh resolution do not necessarily have this requirement for the coming timesteps until the next mesh regeneration. In order to combat this issue, we next consider goal-based adaptivity. 5.9 Goal-based mesh adaptivity Thus far, adaption has been based on a metric calculated from solution data at the current timestep.Goal-based mesh adaptivity seeks to improve this metric, incorporating approximations to the dynamics within regions of the solution field which greatly affect those in the spatial (and temporal) region of interest. For the linear SW case, this involves evaluating an objective functional (3.13) and using the numerical solution of the adjoint equations (3.18). Although the method can be applied more generally, we restrict attention to the linear SW case henceforth, for simplicity. For objective functional (3.13), f(x, t) = 1 A(x)η(x, t) =⇒∂f ∂u=0,∂f ∂η = 1 A(x),(5.42) 31 where Awe have made the (reasonable) assumption that the near-coastal region of importance Ais contained within the solution image Ωs. In this way, the adjoint equations (3.18) have a constant source term in the region A. Given our numerical solution tuple for the SWEs (3.2), we first estimate the error made in the consequent computation of the functional J, providing a measure of what is deemed to be important. In this way, the error estimation stage has the purpose of flagging certain areas of the mesh for refinement, coarsening or movement. To quote pp.1 of [Rannacher,2009], goal-based adaptivity “means the optimization of the mesh and possibly also the discrete model itself on the basis of an a posteriori error estimate”. Such an error estimate takes the general form J(q)−J(qh)≈e(qh) = X K∈H eK(qh),(5.43) where qhdenotes our FEM approximation to the (exact) solution tuple qon a mesh H, with local cell-oriented error indicators eKsumming to give a global error indicator e. Having made an effective choice of error estimator, we seek to minimise the error |J(q)−J(qh)|as in (5.43). To specify the quality of our solution to this problem, we set an upper bound on the error estimate for acceptance. Given an approximation to the solution of the adjoint equations, we make a similar, vectorised version of the argument found in [Power et al.,2006] (to account for the fact we have not just one solution field, but three). In the steady case (where time derivative terms are neglected) we obtain the approximation J(q)−J(qh)≈ hR(qh),λhiΩs,(5.44) which involves both the residual of the forward problem and the approximation to the adjoint problem solution. Thus, in this approach, (5.44) enables us to estimate the associated error in our approximation to the objective functional. There are many other possible choices of error estimate, two of which we discuss in the following. The integral right hand side of (5.44) is approximated using quadrature, in the form of (5.43). We know the exact solution satisfies R(q) = 0. For some linear operator Land what amounts to a source term s, decomposing by R(qh) = Lqh−s, substituting in (5.44) and applying Green’s theorem (ignoring surface terms, as recommended in [Power et al., 2006]) gives an alternative error estimate J(q)−J(qh)≈ hq−qh,L∗λhiΩs.(5.45) For a more accurate error estimate, other approaches are available. Some approaches consider the discrete form of the adjoint, as opposed to the continuous one used above. Another method, described in [Power et al.,2006], is to make ‘coarse-grid’ approximations to the residual. This involves solving the PDE on different mesh resolutions, in order to 32 gain some information concerning how much the residual R(qh) strays from zero. The SWEs (3.2) considered form an unsteady PDE, so it is not so straightforward as to apply steady error estimates at each timestep. An approach to unsteady goalbased anisotropic mesh adaptivity is described in [Belme et al.,2012], which details the extensions required for inner products. For simplicity, however, we follow the approach of [Davis and LeVeque,2016] and consider the simplified error estimate for (5.44), J(q)−J(qh)≈ hqh,λhiΩs.(5.46) The estimate in (5.46) is calculated using vertex-oriented error indicators ep(qh;t) = max τ∈[ e t, t]qh(p, t)Tλh(p, τ),where et= max{t−Tend +Tstart,0}.(5.47) In this approach, at time t, we consider the inner product of forward and adjoint solution approximations over a time range [et, t] within which the objective (3.13) is non-zero. In the isotropic case, the discretised metric field M={Mp}p∈H is simply defined by Mp= diag(ep(qh;t)−2, ep(qh;t)−2). To construct the metric in the anisotropic case, on the other hand, we begin by adapting to a field related to the fluid flow, such as free surface displacement, as usual. We then scale the metric at each vertex by the error indicator. This may be interpreted as replacing the error requirement of (5.17) with e= ep(qh;t)(5.48) Through (5.48), we insist that the error is smaller, and hence that the solution is of better quality in the ‘significant’ regions flagged by the error indicator. In regions of less significance, etakes larger values, meaning that a coarser mesh is allowed. The general goal-based mesh adaptive procedure is summarised by Algorithm 1. Algorithm 1: Goal-based mesh adaptive process /* Fixed mesh adjoint run */ Solve the adjoint problem (backwards in time, from t=T) on a constant, coarse mesh and save the solution data.; /* Adaptive mesh forward run */ for each timestep tn(forward in time, from t= 0)do Estimate the error at tn, using an error estimate (5.43) and the saved data.; Feed the error estimate into the adaptive algorithm by replacing the of (5.17) by the updated error requirement (5.48).; Regenerate the mesh, interpolate variables and solve the forward problem.; end 33 6. Tsunami application As discussed in Chapter 2, this project aims to simulate the tsunami which struck the coast of Fukushima, Japan, in 2011, both accurately and efficiently, using anisotropic mesh adaptivity. The process by which this was achieved is explained in this chapter, with results, analyses and conclusions provided in Sections 6.4 and 6.5. 6.1 Computational setup for the T¯ohoku tsunami Figure 6.1: Domain geometry for the T¯ohoku tsunami. Boundary segments ∂ΩO, ∂ΩFand ∂ΩCare shown in blue, magenta and orange, respectively. White numbers correspond to boundary identity tags. Bathymetry is shown using the QGIS Google Maps plugin. In order to run computations for the example considered, a GIS domain boundary file was constructed using shoreline data from GSHHG, along with bathymetry data from GEBCO. By separating the coastline into parts, we are able to indicate geographical regions where the initial mesh should be finer, such as near to Fukushima. Different GSHHG resolutions are used for the two coastal boundaries near to Fukushima and elsewhere, so small, distant island features are not resolved, otherwise incurring computational expense. Figure 6.1 illustrates the boundary segments used, with bathymetry data visible in the background. Due to the use of different databases for bathymetry and coastlines, with disagreeing sea-level contours, it is necessary to cap the bathymetry so that it is never shallower than 30m, say. The domain boundary ∂Ω comprises of the (disjoint) union of segments ∂ΩO,∂ΩFand ∂ΩC, denoting the open ocean boundary, region of interest surrounding Fukushima and remaining coastal boundary, respectively. A simplification has been made along ∂ΩC, as indicated in Figure 6.1, where a (relatively) small gap between land masses of Hokkaido and Honshu, the Tsugaru Strait, has been closed, as if there exists a sea wall. This is justified as it is sufficiently far from both the earthquake epicentre and the region of interest that any effect the discrepancy may have on the final solution will be negligible. In solving fluids problems on the entire Earth, or when considering timescales of the order of the tides, one should acknowledge both that the Earth is rotating (as manifested in 34 the Coriolis force) and the fact the Earth is (almost) a sphere. We have already discussed the former. With respect to the latter, one approach is to re-define the SWEs in spherical polar co-ordinates. However the Earth is not quite a sphere and there are a number of other effective methods which are regularly used in the computational geosciences. One such approach is transform to UTM coordinates, providing not one single map projection, but a series of projections, taken across the whole globe. The UTM approach divides the Earth into 60 zones, where each band is 6◦in longitude. In addition, zones are divided in the latitudinal direction, given zone letters from the Roman alphabet. UTM coordinates have metric units, taking the horizontal value 500km at zone midpoints. Transformation between latitude-longitude and UTM coordinates is made possible in Python using the functions contained in conversion.py1. Since the domain we consider is contained mainly in zone 54, but also overlaps with zones 53 and 55, we use the force zone number parameter option in our conversion, set to 54, to avoid discontinuities. Figure 6.2: Approximation to the initial surface displacement for the T¯ohoku tsunami used in this project. Using the domain geometry generated in QGIS, we can construct meshes using QMESH online, recently developed at Imperial College [Avdis et al.,2017]. Each boundary segment of the GIS file is given a ‘PhysID’ physical identification tag, so a different level of meshing can be applied within a gradation distance specified in the meshing process. We distinguish between open ocean (PhysID=100) and coastal boundaries (PhysID=200), as displayed in Figure 6.1 - necessary if we were to enforce Dirichlet boundary conditions. We also require a means of interpreting the tsunami initial conditions and bathymetry in Firedrake. Initial condition data was provided by the author of [Saito et al.,2011] as init profile.xyz, containing a set of discretised latitudinal and longitudinal coordinates, along with free surface displacement values. This data was processed using GMT to create a NetCDF surface which could be interpolated in Firedrake, along with the bathymetry field. The resulting initial surface profile is displayed in Figure 6.2. Notice from Figure 6.2 that initial profile of the tsunami is very close to the main Japanese island of Honshu. Witness and measurement reports, such as those used in [Kazama and Noda,2012], [Okada,2011] and [Suzuki et al.,2012], mention that the tsunami wave arrived at many Japanese coastal locations within 30 minutes and at some in just ten minutes. Since we are mainly interested in modelling the wave’s approach to and reaching of the coast, it is reasonable to assume (and fairly clear from experiments in Figure 6.3 overleaf) that the ocean-bound waves will not reach the open ocean boundary within a 25 minute timeframe. Hence, in our application of the SWEs to this case study, a first approximation is to enforce zero flux over the open ocean boundaries. As discussed in Section 3.3, this approach is not particularly physical when waves interact with 1These transformation functions can be downloaded from https://pypi.python.org/pypi/utm. 35 the shoreline and approximates what should be open ocean boundaries as impermeable. Nonetheless, it is sufficient for our purposes, as we are most interested in slightly offshore locations, at gauge locations P02 and P06 as displayed in Figure 2.2. Using code developed in the test problems, an adaptive algorithm for solving the T¯ohoku tsunami SW problem is given by simple adaptive tsunami.py. This script allows for experimentation with different initial QMESH meshes and adaptivity parameters. Adjoint problem solution information is incorporated in goal-based tsunami.py. 6.2 Model verification (a) Non-rot., linear (NL) (b) Non-rot., nonlinear (NN) (c) Rot., linear (RL) (d) Rot., nonlinear (RN) Figure 6.3: Free surface displacement approximations in the T¯ohoku problem after 25 minutes of simulation time under four SW equation sets, where ‘rot.’ denotes ‘rotational’. After establishing the governing SWEs in Section 3.1, we made the simplification of linearity, as in (3.2). We expect the linearised equations to be computationally cheaper to solve than the nonlinear case, hence preferring to use them in our calculations. Before implementing mesh adaptivity in the T¯ohoku tsunami case, let us briefly compare the two approaches to verify the assumption of linearity for ocean-scale SW dynamics. A first comparison of these approaches can be made using an ‘eyeball norm’ consideration of the solutions provided across the entire domain, at a particular point in time. Henceforth an implicit midpoint method is used, with timesteps of length 1 second. Subfigures 6.3a and 6.3b display the shallow water solutions given by the linear and nonlinear, nonrotational approaches after 25 minutes of simulation time, respectively. It is difficult to distinguish any significant differences in the ocean dynamics at this time level, supporting the claim that using the linear standalone solver might be sufficient. Throughout this section, tests are run on a fixed, fine mesh with 97,343 vertices. In the nonlinear case, we take diffusivity parameter ν= 10−3and bottom friction coefficient cb= 0.0025. Additionally, in Section 3.1 the Coriolis effect was assumed to be negligible in tsunami modelling, due to the large difference in timescales and speeds on which the two dynamics act. We can test this assumption numerically by solving the rotational SWEs and comparing the results with the non-rotational case used thus far. As mentioned in Section 36 3.1, in the linear case, the momentum equation corresponding to (3.2) takes the form ∂u ∂t + 2Ω×u+g∇η= 0.(6.1) Equivalent code for the rotational case of (6.1), along with the nonlinear equations, is provided by model verification.py. Evaluating the free surface displacement after 25 minutes, as illustrated in Subfigures 6.3c and 6.3d, yields results indistinguishable by eye to those obtained in the other, non-rotational cases of Figure 6.3. The Coriolis parameter is evaluated vertex-wise using a transformation back into latitude-longitude coordinates. The difficulty in making this analysis is that we do not have the true solution of the fluid dynamics which we seek to approximate. For a more rigorous analysis, we consider a time series of free surface values at two particular locations. A disadvantage of this approach is it only gives a localised, pointwise approximation, rather than a global ones given in Figure 6.3. However, if we choose a location corresponding to an ocean gauge, we can compare model approximations against historical data. Two such gauges are the offshore bottom-pressure gauges P02 and P06, illustrated in Figure 2.2. These pressure gauges, operated by T¯ohoku University, were active between June 2010 and May 2011. We use free surface values from P02 and P06 as calculated in the initial surface inversion analysis of [Saito et al.,2011], as illustrated by the blue lines of Subfigures 6.4a and 6.4b. Given the gauge coordinates (142.5016◦,38.5002◦) and (142.5838◦,38.6340◦), to the nearest tenthousandth of a degree, model outputs may be compared against these data. Using model verification.py, we can run analyses of the various models, yielding the curves in Subfigures 6.4a and 6.4b. The timeseries for the nonlinear solvers are truncated at 25 minutes simulation time, due to extreme computational expense after this point. It is clear from Figure 6.4 there is very little difference between the model results, for this fixed mesh and for these tide gauges. Making use of Python’s clock functionality, we obtain the performance analyses displayed in the table of Subfigure 6.8. Here the averages were taken over two pressure gauge runs considered.2 We observe from Figure 6.8 that the speedup in choosing the linear equations over the nonlinear ones is by a factor of at least 4. In addition, there is a small speedup in considering the linear, non-rotational equations over the linear, rotational equations. It is reasonable to conclude we may exclusively use the linear, non-rotational case henceforth. Our justification is that it provides an similarly accurate and a vastly computationally cheaper approach than the nonlinear cases, and yields a non-negligible saving in computational effort over the linear case wherein the Coriolis effect is considered. It is also clear from Figure 6.4 that the time series for our solution take the right general shape, mimicking the initial rise and fall in free surface rather well, on this high resolution mesh. However, the latter 35 minutes of simulation time is not approximated so well by our model, using any of the equations considered. This should be expected, because, 2Performance tests were run in series on a MacBook Pro with OS X El Capitan. 37 Bibliography F Alauzet. Size gradation control of anisotropic meshes. Finite Elements in Analysis and Design, 46(1):181–202, 2010. Fr´ed´eric Alauzet and Adrien Loseille. A decade of progress on anisotropic mesh adaptation for computational fluid dynamics. Computer-Aided Design, 72:13–39, 2016. Alexandros Avdis, Adam S Candy, Jon Hill, Stephan C Kramer, and Matthew D Piggott. Efficient unstructured mesh generation for marine renewable energy applications. Renewable Energy, submitted, 2017. Baran Aydin. Analytical solutions of shallow-water wave equations. PhD thesis, Citeseer, 2011. Toshitaka Baba, Sebastien Allgeyer, Jakir Hossen, Phil R Cummins, Hiroaki Tsushima, Kentaro Imai, Kei Yamashita, and Toshihiro Kato. Accurate numerical simulation of the far-field tsunami caused by the 2011 T¯ohoku earthquake, including the effects of Boussinesq dispersion, seawater density stratification, elastic loading, and gravitational potential change. Ocean Modelling, 111:46–54, 2017. Nicolas Barral. Time-accurate anisotropic mesh adaptation for three-dimensional moving mesh problems. PhD thesis, Paris 6, 2015. J¨orn Behrens and Michael Bader. Efficiency considerations in triangular adaptive mesh refinement. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 367(1907):4577–4589, 2009. Anca Belme, Alain Dervieux, and Fr´ed´eric Alauzet. Time accurate anisotropic goaloriented mesh adaptation for unsteady flows. Journal of Computational Physics, 231 (19):6323–6348, 2012. S´ebastien Blaise, Amik St-Cyr, Dimitri Mavriplis, and Brian Lockwood. Discontinuous Galerkin unsteady discrete adjoint method for real-time efficient tsunami simulations. Journal of Computational Physics, 232(1):416–430, 2013. S. Brenner and R. Scott. The mathematical theory of finite element methods, volume 15. Springer Science & Business Media, 2007. 44 Richard Courant, Kurt Friedrichs, and Hans Lewy. ¨ Uber die partiellen Differenzengleichungen der mathematischen Physik. Mathematische annalen, 100(1):32–74, 1928. Brisa N Davis and Randall J LeVeque. Adjoint methods for guiding adaptive mesh refinement in tsunami modelling. Pure and Applied Geophysics, 173(12):4055–4074, 2016. Simon W Funke. The automation of PDE-constrained optimisation and its applications. PhD thesis, Imperial College London, UK, 2012. Simon W Funke, Patrick E Farrell, and MD Piggott. Tidal turbine array optimisation using the adjoint approach. Renewable Energy, 63:658–673, 2014. PL George, F Hecht, and MG Vallet. Creation of internal points in Voronoi’s type method. Control adaptation. Advances in engineering software and workstations, 13 (5-6):303–312, 1991. Giorgio Giorgiani, Sonia Fern´andez-M´endez, and Antonio Huerta. Hybridizable discontinuous Galerkin p-adaptivity for wave propagation problems. International Journal for Numerical Methods in Fluids, 72(12):1244–1262, 2013. Wagdi G Habashi, Julien Dompierre, Yves Bourgault, Djaffar Ait-Ali-Yahia, Michel Fortin, and Marie-Gabrielle Vallet. Anisotropic mesh adaptation: towards userindependent, mesh-independent and solver-independent CFD. Part I: general principles. International Journal for Numerical Methods in Fluids, 32(6):725–744, 2000. Motoki Kazama and Toshihiro Noda. Damage statistics (summary of the 2011 off the Pacific coast of T¯ohoku earthquake damage). Soils and Foundations, 52(5):780–792, 2012. Omar Lakkis and Tristan Pryer. A finite element method for nonlinear elliptic problems. SIAM Journal on Scientific Computing, 35(4):A2025–A2045, 2013. Jens Lang, Weiming Cao, Weizhang Huang, and Robert D Russell. A two-dimensional moving finite element method with local refinement based on a posteriori error estimates. Applied Numerical Mathematics, 46(1):75–94, 2003. Eugenie Lisitzin. Sea-level changes, volume 8. Elsevier, 1974. Ramin Moshfegh, Xiangdong Li, and Larsgunnar Nilsson. Gradient-based refinement indicators in adaptive finite element analysis with special reference to sheet metal forming. Engineering computations, 17(8):910–932, 2000. Yoshimitsu Okada. Preliminary report of the 2011 off the Pacific coast of T¯ohoku earthquake. www. bosai. go. jp/e/pdf/Preliminary report110328. pdf, 2011. 45 G´eraldine Olivier. Anisotropic metric-based mesh adaptation for unsteady CFD simulations involving moving geometries. PhD thesis, Universit´e Pierre et Marie Curie-Paris VI, 2011. CC Pain, AP Umpleby, CRE De Oliveira, and AJH Goddard. Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations. Computer Methods in Applied Mechanics and Engineering, 190(29):3771–3796, 2001. MD Piggott, CC Pain, GJ Gorman, PW Power, and AJH Goddard. h, r, and hr adaptivity with applications in numerical ocean modelling. Ocean Modelling, 10(1):95–113, 2005. PW Power, Christopher C Pain, MD Piggott, Fangxin Fang, Gerard J Gorman, AP Umpleby, Anthony JH Goddard, and IM Navon. Adjoint a posteriori error measures for anisotropic mesh optimisation. Computers & Mathematics with Applications, 52(8):1213–1242, 2006. W Pranowo, J Behrens, J Schlicht, and C Ziemer. Adaptive mesh refinement applied to tsunami modelling: Tsunaflash. In The International Conference on Tsunami Warning (ICTW), 2010. Rolf Rannacher. Adaptive finite element discretization of flow problems for goal-oriented model reduction. In Computational Fluid Dynamics 2008, pages 31–45. Springer, 2009. Tatsuhiko Saito, Yoshihiro Ito, Daisuke Inazu, and Ryota Hino. Tsunami source of the 2011 Tohoku-Oki earthquake, Japan: Inversion analysis based on dispersive tsunami simulations. Geophysical Research Letters, 38(7), 2011. Mark Simons, Sarah E Minson, Anthony Sladen, Francisco Ortega, Junle Jiang, Susan E Owen, Lingsen Meng, Jean-Paul Ampuero, Shengji Wei, Risheng Chu, et al. The 2011 magnitude 9.0 T¯ohoku-Oki earthquake: Mosaicking the megathrust from seconds to centuries. Science, 332(6036):1421–1425, 2011. W Suzuki, S Aoi, H Sekiguchi, and T Kunugi. Source rupture process of the 2011 T¯ohokuOki earthquake derived from the strong-motion records. In Proceedings of the 15th Conference on Earthquake Engineering, Lisbon, volume 1650, 2012. Vidar Thom´ee. Galerkin finite element methods for parabolic problems. 2006. M-G Vallet, C-M Manole, Julien Dompierre, S Dufour, and Fran¸cois Guibault. Numerical comparison of some hessian recovery techniques. International Journal for Numerical Methods in Engineering, 72(8):987–1007, 2007. Yu V Vasilevskii and KN Lipnikov. An adaptive algorithm for quasioptimal mesh generation. Computational mathematics and mathematical physics, 39(9):1468–1486, 1999. 46