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

MULTI-SCALE NUMERICAL SIMULATION OF A TSUNAMI USING MESH ADAPTIVE METHODS Joe Wallwork,1,2 Matthew Piggott,2David Ham,2Hilary Weller3. 1Mathematics of Planet Earth Centre for Doctoral Training, 2Imperial College London, 3University of Reading. INTRODUCTION In 2011 a major earthquake caused the T¯ ohoku tsunami whose leading wave struck the Japanese coast near Fukushima just ten minutes after its genesis. Through efficient numerical simulations, sufficient warning could be provided in future scenarios, allowing evacuation and damage mitigation in coastal areas determined at risk. MAIN OBJECTIVES IInvestigate avenues opened by mesh adaptivity, including guidance by adjoint equations, as in goal-based adaptivity. IImplement a mesh adaptive algorithm using the finite element method (FEM) software provided by Firedrake. Experiment with different a posteriori error estimates. IRun mesh adaptive simulations of a realistic tsunami and make accuracy and efficiency comparisons between approaches. NON-ROTATIONAL SHALLOW WATER EQUATIONS For a region of ocean Ω⊂Ò2, define fluid velocity u:Ω→Ò2, free surface displacement η:Ω→Òand bathymetry b:Ω→Ò. Linearising about the surface at rest ¯ η=0 gives ∂u ∂t +g+η=0,∂η ∂t + + · (bu)=0,(1) for gravitational acceleration g=9.81m s−2. Denote q:=(u,η). MESH ADAPTIVE PROCESS For piecewise linear (Ð1) approximation, the Taylor remainder theorem provides an a priori error estimation result =γvT|H|v,γ=const.,(2) upon which to base an anisotropic adaptive algorithm. We ‘reconstruct’ the Hessian Hof a field frelated to the fluid flow, such as the free surface displacement or fluid speed, using a double L2 projection. Based on (2), a new mesh is generated by modifying the Hessian to ensure symmetric positive definiteness. This provides a metric field, which dictates how the mesh is adapted across the domain. blanIn the isotropic case, metric field values are diagonal matrices whose diagonal entries are the corresponding values of 1 f2. GOAL-ORIENTED ADAPTIVITY For a spatial region A⊂Ω, consider the objective functional J(q)=∫Tend Tstart ∬A η(x,y,t)dxdydt.(3) Using (3) we may consider the free surface displacement near to Fukushima’s Daiichi nuclear power plant, say. Firedrake’s automatic differentiation functionality allows the user to obtain discrete adjoint solutions of (1) associated with (3) in a “discretise-then-differentiate” type approach [Gunzburger, 2002]. blank space blank space blank spa spa Goal-oriented adaptivity seeks to establish a finite element discretisation (in this case, a mesh) which enables us to solve (1) with our error in evaluating (3) being below some tolerance >0. That is, |J(q) − J(qh)| <.(4) We first solve (1) on a relatively coarse mesh, from which adjoint solutions may be extracted and hence a posteriori error estimators constructed. The main error estimator considered in this work is the dual-weighted residual Eh=hRh(qh),λhiL2(Ω),(5) where Rhdenotes the (strong) residual formed by taking the LHS of (1) as a 3-component vector functional and (qh,λh)form our approximation of the primal-dual pair. In the final run, (1) is solved adaptively, with (local) error estimates indicating which regions of mesh should be coarsened or refined. The simplest way of integrating estimators is to create an isotropic metric through scaling the identity matrix by (5). RESULTS The quality of our approximations may be assessed using (4). An ‘exact’ objective functional value for (3) is obtained by iterating over increasingly well-resolved meshes until convergence is attained to three significant figures, as displayed in Figure 1. An objective value J(q)=2.44 ×1013 is achieved on a 196,560 element mesh. 104105 #Elements 2.425 2.430 2.435 2.440 2.445 2.450 2.455 2.460 Objective functional J(u, v, η) = RTend Tstart RAηdxdt ×1013 Figure: Mesh ‘bootstrapping’. We compare (anisotropic) adaption to the free surface Hessian (‘simple adaptivity’) against (isotropic) ‘goal-based’ adaptivity, which incorporates adjoint data to form error indicators using (5). The mixed space Ð1DG −Ð2 is used throughout. Coarse Medium Fine Simple adaptive Goal-based |J(q)−J(qh)| |J(q)| 0.67% 0.34% 0.15% 1.51% 0.20% Mean element count 8,782 20,724 81,902 11,141 11,506 Mean run time (s)11.7 38.1 131.7 405.2 7700.8 The interpolation involved in forming residuals Rh(qh)for the error estimators (5) increases the cost of the goal-based algorithm. Instead, we could solve local boundary value problems on each element, as in the element residual method [Ainsworth & Oden, 1997]. Figure: Initial condition blank space bl[Saito et al., 2011]. Figure: ‘Simple adaptive’ mesh. Figure: ‘Goal-based’ mesh. CONCLUSIONS IThe adjoint problem can be useful for guiding the mesh adaptive process, offering additional information than when only the forward problem (1) is considered. For the T¯ ohoku tsunami, the goal-based approach out-performs the ‘simple adaptive’ method in terms of accuracy and and is competitive regarding element count. IWhilst error and element count are reduced, there is progress to be made in reducing the run time of the mesh adaptive simulations. FURTHER RESEARCH The element residual method should be implemented, to provide a more competitive run time for the goal-based algorithm. Further research will include incorporation of r-adaptivity, allowing ‘mesh movement’. Other ocean modelling applications will be considered for mesh adaptivity, such as storm surges and Gulf Stream separation. REFERENCES Ainsworth & Oden (1997). A posteriori error estimation in finite element analysis. Gunzburger (2002). Perspectives in flow control and optimization. Saito et al. (2011). Tsunami source of the 2011 T¯ ohoku-Oki earthquake, Japan: Inversion analysis based on dispersive tsunami simulations. GitHub github.com/jwallwork23 Centre website mpecdt.org Firedrake firedrakeproject.org Thetis thetisproject.org This work could not have been done without the generous support of