Adaptive discontinuous Galerkin method for tsunami modeling and - PowerPoint PPT Presentation

S Abarbanel Memorial 2018 Adaptive discontinuous Galerkin method for tsunami modeling and prediction on a global scale Jan S Hesthaven EPFL-SB-MATH-MCSS Jan.Hesthaven@epfl.ch B. Bonev (EPFL, CH) F. Giraldo (NPS, US) M. Hajihassanpour (Sharif, Iran) M. A. Kopera (UC Santa Cruz, US)

Motivation Cost ‣ Fukushima (Japan, 2011) 
 Created by deep sea earth 500’000 people evacuated ~ $650B quakes or large meteors ‣ Indian Ocean (2004) 
 ~225’000 killed, ~ $20B

A Tsunami warning system Once an earth quake and a possible tsunami is detected, a model is needed to predict propagation and issue a warning and evacuation The goal is this work is develop such a model on a global scale

Challenges ‣ What is the right physical model ? ‣ Tsunami’s are very long wave lengths - depends on generation 
 ‣ How do we solve the equations ? ‣ We need flexibility and accuracy over long time ‣ How do we deal with land-ocean interfaces ‣ Embedding is the only realistic approach ‣ We must avoid artificial waves ‣ Scheme must well balanced ‣ What about boundaries ? ‣ Local or global model ? ‣ Efficiency ? ‣ A global model must be adaptive

The physical model We consider the shallow water model - long waves Draft Dra Surface Elevation z ϕ + τ Geoid x ϕ − τ u Topography The Shallow Water Equations in one dimension: ∂ t ϕ + ∂ ∂ ∂ x ( ϕ u ) = 0 ∂ t ( ϕ u ) + ∂ ∂ ϕ u 2 + 1 = − ∂ 1 2 ϕ 2 2 ∂ x τ ∂ x

The physical model Draft We express it on a more general form the SWE in conservation form We write the shalloe water equations as a system of m conservation laws on Ω – x : ; ∂ q ∂ t + Ò · F ( q ) = S ( x , q ) in Ω ◊ (0 , Œ ) (3) q = q 0 on Ω ◊ { t = 0 } , F ( q ) and S ( x , q ) are the flux and source terms and q = q ( x , t ) : Ω ◊ (0 , Œ ) æ R m is the unknown solution with suitable initial condition q 0 ( x ) Note that it is inhomogeneous, i.e., in steady state we must have raft ∂ q ∂ t = 0 ⇔ Ò · F ( q ) = S ( x , q ) . in the so-called lake at rest solution, It is important to respect this at the discrete level

The physical model Draft To avoid problems with artificial boundaries, we consider the full Draft global model. We use Cartesian coordinates with a 4D state vector ϕ w $ T q = # ϕ ϕ u ϕ v r and the associated flux S T S T S T ϕ u ϕ v ϕ w D ϕ u 2 + 1 2 ϕ 2 ϕ uv ϕ uw V ˆ V ˆ V ˆ F ( q ) = i + j + k . W X W X W X ϕ v 2 + 1 2 ϕ 2 ϕ uv ϕ vw U U U ϕ w 2 + 1 2 ϕ 2 ϕ uw ϕ vw $ # $ The source term is C ( x , u ) = − 2 ω z ϕ Coriolis force x × u , i.e., x · u = 0, R 2 µ ( x , q ) = 1 ˜ S ( x , q ) = C ( x , u ) − ϕ Ò τ + µ x . � ϕ ∇ τ + ∇ · ˜ � R 2 x · F , Bottom topography Restriction to sphere

The computational model Divide Ω into non-overlapping quadrilateral, curvilinear elements D k ∈ G s.t. Draft Draft Ω = t K k =1 D k . Approximate solution lies in the set of piecewise polynomials: v ∈ L 2 ( Ω ) | ∀ D k ∈ G , v | D k ∈ P ( D k ) * V ( Ω , G ) := ) so that K n q k n q ( x , t ) ≈ q N ( x , t ) = N ( x , t ) . k =1 2 , and respective Lagrange ba With GLL nodes { ξ i } constructed on I = [ − 1 , 1] 2 , and respective Lagrange basis With GLL nodes constructed on functions L i ( ξ ) : ( N +1) 2 ÿ q k q k N ( x ) = N ( x i ) L j ( ξ ( x )) . j =1 ( − 1 , +1) (+1 , +1) D k ξ = Ψ ( x ) I x = Ψ − 1 ( ξ ) η z ( − 1 , − 1) (+1 , − 1) y ξ x

The computational model Draft Draft We use the discontinuous Galerkin formulation Weak Form/Green’s Form 1 ∂ q N ⁄ j 2 n · F ∗ L i d x = − N L i d x . − F N · Ò − S N ˆ ∂ t D ∂ D Dra An additional integration by parts yields the strong form: An additional integration by parts yields the strong form: Strong Form/Divergence Form 1 ∂ q N ⁄ j 2 + Ò · F N − S N L i d x = n · ( F N − F ∗ N ) L i d x . ˆ ∂ t ∂ D D Suitable numerical representations of the flux and source terms F N , S N are: ( N +1) 2 ( N +1) 2 ÿ ÿ F N ( q N ) = F ( q ( x j )) L j ( x ) , S N ( x , q N ) = S ( x j , q ( x j )) L j ( x ) , j =1 j =1 We use the local Lax-Friedrichs flux = 1 − c N , q + N ) + F N ( q + q + ! " ! N ) " ! " F ∗ F N ( q − q − N − q − . N N N 2 2

Nonconforming grids We consider a general setup F ∗ , 1 q 1 N , q 0 q 1 N , P 1 0 q 0 = F ∗ � � � � , Fluxes evaluated at children N N N N F ∗ , 2 q 2 N , q 0 q 2 N , P 2 0 q 0 � � = F ∗ � � , N N N N Maximum rule = 1 F ∗ , 0 q 0 N , q 1 N ⊕ q 2 P 0 P 1 0 q 0 N , q 1 + P 0 P 2 0 q 0 N , q 2 1 F ∗ 2 F ∗ � � � � � � �� . N N N N N N 2

Well balanced property Draft Geophysical Systems often have steady-state solutions q such that Draft ∂ q ∂ t = 0 ⇔ Ò · F ( q ) = S ( x , q ) . We are mostly interested in the so-called lake at rest solution, which is given as ϕ ( x , t ) = ϕ 0 − τ ( x ) , u ( x , t ) = 0 . Most practical situations involve perturbation of this solution. Definition $ | denote the numerical representation of the ‘lake at rest’ solution, Let q N = # ϕ N , ϕ u N such that ∀ x i : ϕ N ( x i ) = max { ϕ 0 − τ N ( x i ) , 0 } (20) ϕ u N ( x i ) = 0 . (21) We call a scheme with right-hand side R N ( q N ) well-balanced, i ff it exactly preserves the ‘lake at rest’ solution, i.e. R N ( q N ) = 0 . (22)

Well balanced property t f Replacing the integrals with the quadrature rules yields the discrete version of the weak form a ⁄ j n · F ∗ R ( q N ) = F N · Ò L i − S N L i d x + N d x ˆ ∂ D D r n · F ∗ N L i ] =: R N ( q N ) . ≈ Q D [ F N · Ò N L i − S N L i ] + Q ∂ D [ ˆ (23) D Remark In general, for the weak form , we can only guarantee the well-balanced property if numerical integration is exact! Due to the rational Jacobian, this is not the case here.

Well balanced property Draft The strong form reads ⁄ j n · ( F N − F ∗ R ( q N ) = ( Ò · F N − S N ) L i d x − N ) L i d x ˆ ∂ D D n · ( F N − F ∗ N ) L i ] =: R N ( q N ) . ≈ Q D [( Ò · F N − S N ) L i ] − Q ∂ D [ ˆ Proposition Let G be a conforming mesh of Ω and let the approximation of the bathymatry τ N ∈ V ( Ω , G ) be continuous on Ω . Moreover let ϕ 0 > τ N , i.e. there are no dry areas in Ω . Under these circumstances, the strong form of the scheme is well-balanced. Individual conditions for well-balancedness As a consequence, the strong form decouples the well-balancing condition into Ò · F N = S N , which has to be satisfied on all nodes and the flux consistency condition n · ( F N − F ∗ N ) = 0 ˆ on the interface nodes. This can also be extended to non-conforming grids

Wetting/Drying interface Draft Some of the issues at the wet/dry interface: • How do we maintain positivity on nearly dry nodes? • How do we evaluate ( ϕ u ) 2 for small ϕ ? ϕ • How do we guarantee stability near dry regions? • How can we keep the well-balanced property at the wet/dry interface? • Some possibile strategies • Grid conforming to the wet/dry interface. + Accurate treatment of the interface. − Expensive re-meshing and treatment of boundary conditions required. • Fixed mesh but dry cells are turned o ff . + Simple to handle. − Sudden inclusion/exclusion of the dry elements breaks conservation. • Keep a thin layer on drying nodes. + Not very expensive and avoids the sudden inclusion/exclusion. − Treatment of artificial pressure gradients due to the dry nodes.

Positivity limiter Draft Idea: Ensure positivity of the average ⁄ Draft ϕ = ϕ ( x )d x D k in each cell, then rescale the nodal values in order to be positive 7 . CFL-like condition Assuming exact integration, we can show that under the condition J e J ∆ tc ≤ ω 1 (35) 2 where c is the signal velocity, ω 1 the first weight of the numerical integration and J , J e are the Jacobians of the volume and edge parametrizations respectively. If we use a convex combination of Euler forward steps, this property is retained. Strong-stability preserving Runge-Kutta (SSPRK) schemes are such schemes 8 .

Dr Positivity limiter ϕ N ϕ ∗ N ϕ N D k Figure: Application of the positivity limiter We use the linear rescaling around the positive average ; < ϕ N θ = min 1 , m = min { ϕ N ( x i ) } , ϕ N − m i to ensure positivity of the nodal values 9 . The rescaled solution is then: ϕ ∗ N = θ ( ϕ N − ϕ N ) + ϕ N u ∗ N = θ ( u N − u N ) + u N

D Wetting/Drying interface ϕ + τ ϕ N + τ N τ N τ D k ¯ D k ˜ D k Dr (a) exact solution (b) numerical approximation Remark The positivity limiter ensures positivity, but if we do not allow negative waterheights, we can not ensure ϕ N Ò ϕ N = − ϕ N Ò τ N . Therefore well-balancedness is lost. Solution We replace Ò with Ò N , where Ò N is a finite-di ff erence approximation of the gradient operator.