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

▶
adaptive discontinuous galerkin method for tsunami
SMART_READER_LITE
LIVE PREVIEW

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

Dec 16, 2022 468 likes •880 views

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

slide-1
SLIDE 1

Jan S Hesthaven EPFL-SB-MATH-MCSS Jan.Hesthaven@epfl.ch

Adaptive discontinuous Galerkin method for tsunami modeling and prediction on a global scale

S Abarbanel Memorial 2018
  • B. Bonev (EPFL, CH)
  • F. Giraldo (NPS, US)
  • M. Hajihassanpour (Sharif, Iran)
  • M. A. Kopera (UC Santa Cruz, US)
slide-2
SLIDE 2

Motivation

Cost

  • Fukushima (Japan, 2011)


500’000 people evacuated ~ $650B

  • Indian Ocean (2004)


~225’000 killed, ~ $20B Created by deep sea earth quakes or large meteors

slide-3
SLIDE 3

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

slide-4
SLIDE 4

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
slide-5
SLIDE 5

The physical model

We consider the shallow water model - long waves

Draft

x z −τ ϕ ϕ + τ Geoid Topography Surface Elevation u

Dra

The Shallow Water Equations in one dimension: ∂ ∂tϕ + ∂ ∂x(ϕu) = 0 ∂ ∂t(ϕu) + ∂ ∂x 1 ϕu2 + 1 2ϕ22 = − ∂ ∂xτ
slide-6
SLIDE 6

The physical model

Draft

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, Œ) q = q0
  • n Ω ◊ {t = 0},
(3) F (q) and S(x, q) are the flux and source terms and q = q(x, t) : Ω ◊ (0, Œ) æ Rm is the unknown solution with suitable initial condition q0(x)

We express it on a more general form 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

slide-7
SLIDE 7

The physical model

To avoid problems with artificial boundaries, we consider the full global model. We use Cartesian coordinates with a 4D state vector

Draft

q = # ϕ ϕu ϕv ϕw$T

and the associated flux

Draft

F (q) = S W U ϕu ϕu2 + 1 2ϕ2 ϕuv ϕuw T X Vˆ i + S W U ϕv ϕuv ϕv2 + 1 2ϕ2 ϕvw T X Vˆ j + S W U ϕw ϕuw ϕvw ϕw2 + 1 2ϕ2 T X Vˆ k. # $

The source term is

D r

$

˜ S(x, q) = C(x, u) − ϕÒτ + µx. Coriolis force Bottom topography Restriction to sphere

i.e., x · u = 0, C(x, u) = −2ωzϕ R2 x × u, µ(x,q) = 1 R2 x · ϕ∇τ + ∇ · ˜ F
  • ,
slide-8
SLIDE 8

The computational model

Draft

Divide Ω into non-overlapping quadrilateral, curvilinear elements Dk ∈ G s.t. Ω = tK k=1 Dk. Approximate solution lies in the set of piecewise polynomials: V(Ω, G) := ) v ∈ L2(Ω) | ∀Dk ∈ G, v|Dk ∈ P(Dk)* so that q(x, t) ≈ qN(x, t) = K

n

k=1 qk N(x, t). With GLL nodes constructed on 2, and respective Lagrange ba

Draft

n With GLL nodes {ξi} constructed on I = [−1, 1]2, and respective Lagrange basis functions Li(ξ): qk N(x) = (N+1)2 ÿ j=1 qk N(xi)Lj(ξ(x)). x z y Dk ξ η I ξ = Ψ(x) x = Ψ−1(ξ) (−1, −1) (−1, +1) (+1, +1) (+1, −1)
slide-9
SLIDE 9

The computational model

Draft

We use the discontinuous Galerkin formulation Weak Form/Green’s Form ⁄ D 1 ∂qN ∂t − F N · Ò − SN 2 Li dx = − j ∂D ˆ n · F ∗ N Li dx. An additional integration by parts yields the strong form:

Draft

An additional integration by parts yields the strong form: Strong Form/Divergence Form ⁄ D 1 ∂qN ∂t + Ò · F N − SN 2 Li dx = j ∂D ˆ n · (F N − F ∗ N) Li dx.

Dra

Suitable numerical representations of the flux and source terms F N, SN are: F N(qN) = (N+1)2 ÿ j=1 F (q(xj))Lj(x), SN(x, qN) = (N+1)2 ÿ j=1 S(xj, q(xj))Lj(x), We use the local Lax-Friedrichs flux F ∗ N ! q− N, q+ N " = 1 2 ! F N(q− N) + F N(q+ N)" − c 2 ! q+ N − q− N " .
slide-10
SLIDE 10

Nonconforming grids

We consider a general setup

F ∗,1 N
  • q1
N,q0 N
  • = F ∗
N
  • q1
N,P1 0q0 N
  • ,
F ∗,2 N
  • q2
N,q0 N
  • = F ∗
N
  • q2
N,P2 0q0 N
  • ,

Fluxes evaluated at children

F ∗,0 N
  • q0
N,q1 N ⊕ q2 N
  • = 1
2
  • P0
1 F ∗ N
  • P1
0q0 N,q1 N
  • + P0
2 F ∗ N
  • P2
0q0 N,q2 N
  • .

Maximum rule

slide-11
SLIDE 11

Well balanced property

Draft

Geophysical Systems often have steady-state solutions q such that ∂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.

Draft

Definition Let qN = # ϕN, ϕuN $| denote the numerical representation of the ‘lake at rest’ solution, such that ∀xi: ϕN(xi) = max {ϕ0 − τN(xi), 0} (20) ϕuN(xi) =0. (21) We call a scheme with right-hand side RN(qN) well-balanced, iff it exactly preserves the ‘lake at rest’ solution, i.e. RN(qN) = 0. (22)
slide-12
SLIDE 12

Well balanced property

D r a f t

Replacing the integrals with the quadrature rules yields the discrete version of the weak form R(qN) = ⁄ D F N · ÒLi − SNLi dx + j ∂D ˆ n · F ∗ N dx ≈ QD [F N · ÒNLi − SNLi] + Q∂D [ˆ n · F ∗ N Li] =: RN(qN). (23) 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.
slide-13
SLIDE 13

Well balanced property

Draft

The strong form reads R(qN) = ⁄ D (Ò · F N − SN)Lidx − j ∂D ˆ n · (F N − F ∗ N) Li dx ≈ QD [(Ò · F N − SN)Li] − Q∂D [ˆ n · (F N − F ∗ N) Li] =: RN(qN). 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 = SN, which has to be satisfied on all nodes and the flux consistency condition ˆ n · (F N − F ∗ N) = 0
  • n the interface nodes.

This can also be extended to non-conforming grids

slide-14
SLIDE 14

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 off.
+ 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.
slide-15
SLIDE 15

Positivity limiter

Draft

Idea: Ensure positivity of the average ϕ =

⁄

Dk ϕ(x)dx in each cell, then rescale the nodal values in order to be positive7.

Draft

CFL-like condition Assuming exact integration, we can show that under the condition Je J ∆tc ≤ ω1 2 (35) where c is the signal velocity, ω1 the first weight of the numerical integration and J,Je 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 schemes8.
slide-16
SLIDE 16

Positivity limiter

Dr

Dk ϕ∗ N ϕN ϕN Figure: Application of the positivity limiter We use the linear rescaling around the positive average θ = min ; 1, ϕN ϕN − m < , m = min i {ϕN(xi)} to ensure positivity of the nodal values9. The rescaled solution is then: ϕ∗ N = θ(ϕN − ϕN) + ϕN u∗ N = θ(uN − uN) + uN
slide-17
SLIDE 17

Wetting/Drying interface

D

¯ Dk ˜ Dk ϕ + τ τ (a) exact solution Dk ϕN + τN τN (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.

Dr

Solution We replace Ò with ÒN, where ÒN is a finite-difference approximation of the gradient
  • perator.
slide-18
SLIDE 18

Wetting/Drying interface

Note: This is locally 1st order only - but in very few cells

D

We call the node xi a wet node if ϕ(xi) > ϕmin and consequently, a dry node if ϕ(xi) ≤ ϕmin. Then, we define the finitie difference operator Dξf(ξi) :=

Y _ _ _ ] _ _ _ [

f(ξi+1)−f(ξi−1) ξi+1−ξi−1 if ξi+1, ξi, ξi−1 are wet nodes f(ξi+1)−f(ξi) ξi+1−ξi if ξi+1, ξi are wet nodes, but ξi−1 is dry f(ξi)−f(ξi−1) ξi−ξi−1 if ξi, ξi−1 are wet nodes, but ξi+1 is dry
  • therwise
and similiarly Dη. In the physical domain, these operators translate to Dx := ∂ξ ∂xDξ + ∂η ∂xDη. Using these operators we replace the exact gradient operator with ÒN := # Dx Dx Dx

$|.

slide-19
SLIDE 19

1D tests

: × → F(q) =
  • ϕu
ϕu2 + 1 2ϕ2
  • S(x,q) =
  • −ϕ∇τ
  • .

We consider the 1D shallow water equation

D

−2 −1 1 2 0.1 0.2 0.3 0.4 x h t = 0.75 −2 −1 1 2 −0.5 0.5 x hu t = 0.75

D

100 101 102 103 10−8 10−5 10−2 1
  • 2
|Ω|/h L2 Error 100 101 102 10−13 10−7 10−1 1
  • 2
1
  • 3
1
  • 4
1
  • 5
|Ω|/h L2 Error
slide-20
SLIDE 20

1D tests

A sloping beach test

D

−1 1 −0.2 0.2 0.4 x h + b t = 1.0 −1 1 −2 −1 1 2 ·10−3 x hu t = 1.0

D

0.2 0.4 0.6 0.8 1 10−15 10−10 t
  • rel. mass error
0.2 0.4 0.6 0.8 1 10−17 10−10 10−3 t
  • rel. energy error
slide-21
SLIDE 21

1D tests

Draft

101 102 103 10−4 10−2 100 1
  • 1
1
  • 1
|Ω|/h L2 Error 101 102 103 10−5 100 1
  • 2
1
  • 4
|Ω|/h L2 Error
slide-22
SLIDE 22

1D tests

slide-23
SLIDE 23

Putting it all together

We have presented the elements of a tool that

  • solves the shallow water equation on the sphere
  • is a high-order accurate, non-conforming DG methods with 


support for adaptive solution.

  • is well balanced
  • deals robustly with wetting/drying
  • initial simple tests suggest that it performs well

Does it work ?

slide-24
SLIDE 24

A fantasy problem

An event east of Africa Source is a sudden displacement of the bottom

slide-25
SLIDE 25

Validation

So it looks good — but is it right ? Two different models for earth quake excitation

  • Static excitation 


a sudden displacement of the sea floor 


  • Dynamic excitation


a dynamic displacement of sea floor

ied as 휁(풙, 푡) = ∑ 푓 푖(푡, 푡0 푖 , 푡0 푖+1) 푂푖(풙푖,푙푖, 푤푖,푑푖,. . . ), 푁 푖=1 unctions 푓 푖(푡, 푡0 푖 , 푡0 푖+1) can be used [28, 29, 38]. Here, a l 푓 푖(푡,푡0 푖 , 푡0 푖+1) = ⎩ ⎪ ⎪ ⎨ ⎪ ⎪ ⎧ if 푡 < 푡0 푖 푡 − 푡0 푖 푡0 푖+1 − 푡0 푖 if 푡0 푖 ≤ 푡 ≤ 푡0 푖+1 1 if 푡0 푖+1 < 푡 .

The parameters are recovered from measurements 
 and earth quake models

slide-26
SLIDE 26

First Validation case - Tohoku event

Japan, 2011. Static model for earth quake

slide-27
SLIDE 27

First Validation case - Tohoku event

What does the grid look like ? The benefits of the non-conforming grid is clear

slide-28
SLIDE 28

First Validation case - Tohoku event

D

1 2 3 4 5 6 7 8 9 10 t = 0 d t = 0.3 d

The DART system was put in place after the Thailand event. 
 We use it for validation now

slide-29
SLIDE 29

First Validation case - Tohoku event

Draft

−1 1 2 DART 21418 −0.5 0.5 1 DART 21413 −0.2 0.2 0.4 DART 21416 −0.2 0.2 Wave amplitude DART 52402 −0.2 0.2 DART 52405 −0.2 0.2 DART 46408 0:00 6:00 12:00 −0.1 0.1 DART 52403 0:00 6:00 12:00 −0.1 0.1 0.2 Hours after earthquake DART 52406 0:00 6:00 12:00 −0.2 0.2 DART 51407

Blue - data Red - N=2 Green - N=4

slide-30
SLIDE 30

First Validation case - Tohoku event

Static or dynamic excitation ? Differences are minimal except in near field

slide-31
SLIDE 31

First Validation case - Tohoku event

Static or dynamic excitation ?

time (hr) height (m) 2 4 6 8 10 12
  • 2.5
  • 2
  • 1.5
  • 1
  • 0.5
0.5 1 1.5 2 Buoy data N=4 Dynamic N=2 Dynamic 0.2 0.4 0.6 0.8 1 1.2
  • 1
1 2 time (hr) height (m) 2 4 6 8 10 12
  • 2.5
  • 2
  • 1.5
  • 1
  • 0.5
0.5 1 1.5 2 Buoy data N=4 Static N=2 Static 0.2 0.4 0.6 0.8 1 1.2
  • 1
1 2

Buoy 1 Buoy 4

slide-32
SLIDE 32

Second Validation case - Thailand event

Very complex and long exitation Involves 14 faults and a 10 min excitation
 Longest earth quake ever recorded and the 3rd strongest (9.1-9.3)

slide-33
SLIDE 33

Second Validation case - Thailand event

Excitation

shallow waters. Note that in both the Tohoku and Sumatra-Andaman tsunamis this ratio is relatively small.
slide-34
SLIDE 34

Second Validation case - Thailand event

Grid used Temporal adaptivity not used

comparison to the Tohoku tsunami simulation is required. 325 326
slide-35
SLIDE 35

Second Validation case - Thailand event

slide-36
SLIDE 36

Second Validation case - Thailand event

Static or dynamic excitation ? Same scale as before - differences are substantial in near field

slide-37
SLIDE 37

Second Validation case - Thailand event

Costal validation Delay caused by bottom friction

have a satisfactory agreement with the measurements obtained by the satellites. time (hr) height (m) 1 2 3 4 5 6
  • 3
  • 2
  • 1
1 2 3 Measurements (Hanimaadhoo) N = 4 Dynamic N = 4 Static time (hr) height (m) 1 2 3 4 5 6
  • 3
  • 2
  • 1
1 2 3 Measurements (Gan) N = 4 Dynamic N = 4 Static Lat SSH (m) 10
  • 1
  • 0.5
0.5 1 Measurements N=4 Dynamic Fujii and Satake Gopinathan et al.

Cross validation against past work

slide-38
SLIDE 38

Second Validation case - Thailand event

Satellite validation

Lat SSH (m) 10
  • 1
  • 0.5
0.5 1 Measurements N=4 Dynamic N=2 Dynamic Lat SSH (m) 10
  • 1
  • 0.5
0.5 1 Measurements N=4 Static N=2 Static Lat SSH (m)
  • 20
  • 10
  • 1
  • 0.5
0.5 1 Measurements N=4 Dynamic N=2 Dynamic Lat SSH (m)
  • 20
  • 10
  • 1
  • 0.5
0.5 1 Measurements N=4 Static N=2 Static

Jason-1

1.9h after quake

TOPEX

2h after quake
slide-39
SLIDE 39

Second Validation case - Thailand event

Static or dynamic excitation ? Jason-1 TOPEX/Poseidon

Buoy L2-error static model L2-error dynamic model Hanimadhoo 0.6234 0.5959 Gan 0.4353 0.3995 Jason-1 0.2047 0.2048 TOPEX/Poseidon 0.2775 0.1527
slide-40
SLIDE 40

Thank you !