Well-Balanced and Positivity Preserving DG Schemes for Shallow Water - PowerPoint PPT Presentation

Well-Balanced and Positivity Preserving DG Schemes for Shallow Water Flows with Shock Capturing by Adaptive Filtering Procedures Sigrun Ortleb 1 , Andreas Meister 1 , Thomas Sonar 2 1 Mathematics and Natural Sciences, University of Kassel, Germany 2 Computational Mathematics, TU Braunschweig, Germany HYP 2012

The 2D Shallow Water Equations H water height b bottom elevation ∂ u ∂ t + ∇ · F ( u ( x , t )) = s ( u , x , t ) u = ( H , Hv 1 , Hv 2 ) T , ∂ x f 1 + ∂ ∂ Where ∇ · F = ∂ y f 2 ,     Hv 1 Hv 2 Hv 2 1 + 1 2 gH 2  ,  , f 1 = f 2 = Hv 1 v 2   Hv 2 2 + 1 2 gH 2 Hv 1 v 2 (0 , − gHb x , − gHb y ) T s =

SWE – Challenges for a numerical scheme Shock capturing, Well-balancedness, Positivity preservation Xing, Zhang, Shu ’10: 3rd order DG scheme on cartesian grids, WB + PP, Limiters Bryson, Epshteyn, Kurganov, Petrova ’10: 2nd order central-upwind scheme on triangular grids, WB + PP, Limiters In this talk: 3rd order WB and PP DG scheme on triangular grids Shock capturing by (low cost) filtering procedures ∗ , + Positivity preservation for implicit time integration ∗ Joint work with: ∗ Andreas Meister (Univ. Kassel) + Thomas Sonar (TU Braunschweig)

The DG Scheme for SWE Ω u h ( · , t ) T h = { τ 1 , τ 2 , . . . , τ # T h } Triangulation Triangulation of computational domain Piecewise polynomial approximation u h ( x , t ) of degree N Variational formulation in space DG advantages: conservative method, highly local, unstructured grids easy to implement

The DG Scheme for SWE Variational formulation � � � � d u Φ d x + F ( u ) · n Φ d σ − F ( u ) · ∇ Φ d x = s Φ d x dt ↑ ↑ ↑ ↑ τ i ∂τ i τ i τ i F num ( u − h , u + h , n ) F ( u h ) u h s h

The DG Scheme for SWE Variational formulation � � � � d F num ( u − h , u + u h Φ d x + h , n ) Φ d σ − F ( u h ) ·∇ Φ d x = s h Φ d x dt τ i ∂τ i τ i τ i Use orthogonal polynomials q ( N ) � u i ˆ u h | τ i ( x , t ) = k ( t )Φ k ( x ) , q ( N ) = ( N + 1)( N + 2) / 2 k =1 Semidiscrete equation for coefficients → explicit RK method � � � � d u i F num ( u − h , u + / � Φ k � 2 dt ˆ = − h , n ) Φ k d σ + F ( u h ) · ∇ Φ k d x L 2 k ∂τ i τ i � s h Φ k d x / � Φ k � 2 + L 2 τ i Quadrature rules

The DG Scheme for SWE Variational formulation � � � � d F num ( u − h , u + u h Φ d x + h , n ) Φ d σ − F ( u h ) ·∇ Φ d x = s h Φ d x dt τ i ∂τ i τ i τ i Use orthogonal polynomials q ( N ) � u i ˆ u h | τ i ( x , t ) = k ( t )Φ k ( x ) , q ( N ) = ( N + 1)( N + 2) / 2 k =1 Semidiscrete equation → Time integration: RK method � � � � d u i F num ( u − h , u + / � Φ k � 2 dt ˆ = − h , n ) Φ k d σ + F ( u h ) · ∇ Φ k d x L 2 k ∂τ i τ i � s h Φ k d x / � Φ k � 2 + L 2 τ i Quadrature rules

Orthogonal Proriol-Koornwinder-Dubiner Basis � � � 1 − s � l 21 + r Φ lm ( r , s ) = P 0 , 0 P 2 l +1 , 0 1 − s − 1 ( s ) , l + m ≤ N m l 2 P α,β Jacobi polynomials n Defined on reference element T T = { ( r , s ) ∈ R 2 | − 1 ≤ r , s ; r + s ≤ 0 } Proriol (1957), Koornwinder (1975), Dubiner (1991)

New Modal Filtering Approach Damping step: Carried out after each time step � � 2 p � � l + m − const · N ∆ t u i , mod u i ˆ = exp ˆ lm lm h i N � − α i · η 2 p � � Multiplication by filter function σ ( η ) = exp Resulting from SV modified equation on T Meister, Ortleb, Sonar: - On Spectral Filtering for Discontinuous Galerkin Methods on Unstructured Triangular Grids, NMPDE, ’11 - New Adaptive Modal and DTV Filtering Routines for the DG Method on Triangular Grids applied to the Euler Equations, GEM, ’12

Shock Indication and Adaptivity u i , mod u i ˆ = σ ( η ) · ˆ lm lm � − s i α i η 2 p � σ ( η ) = exp Modified filter function Resolution indicator ω res � u i lm ) 2 (ˆ � ˜ if ˜ s i > tol s i l + m = N ω res = s i = � i lm ) 2 + ǫ 0 otherwise u i (ˆ l + m < N tol = 0 . 01 � 1000(5 N 4 + 1) ω res � ˜ s i = min , 1 i Persson, Peraire ’06; Barter, Darmofal ’07

Circular Dam-Break Problem Computational domain: [0 , 40] × [0 , 40] Initial conditions: � 2 . 5 if ( x − 20) 2 + ( y − 20) 2 ≤ 2 . 5 2 H ( x , y , 0) = 0 . 5 else Boundaries: outflow conditions

Circular Dam-Break Problem Modally filtered DG solution: Water height N = 6 , K = 2200. α i = 2 N ∆ t h i , p = 1. t = 0 . 4 t = 0 . 7 t = 1 . 4

Circular Dam-Break Problem Modally filtered DG solution: Water height N = 6 , K = 2200. α i = 2 N ∆ t h i , p = 1. t = 3 . 5 t = 4 . 7

DTV Filtering on Non-Cartesian Graphs U (0) = [ u h ( x j )] j ∈ DTV nodes for k = 0,1,2,. . . U ( k +1) = DTV ( U ( k ) ) Meister, Ortleb, Sonar ’11 t = 0 . 4: before DTV filtering after DTV filtering

Circular Dam-Break Problem DTV filtered solutions: Water height N = 6 , K = 2200. 49 subtris, λ = 1. t = 0 . 4 t = 0 . 7 t = 1 . 4

Circular Dam-Break Problem DTV filtered solutions: Water height N = 6 , K = 2200. 49 subtris, λ = 1. t = 3 . 5 t = 4 . 7

Well-balanced Scheme � � � � d F WB ( u − ∗ , u + ∗ , H − , n ) Φ d σ − u h Φ d x + F ( u h ) ·∇ Φ d x = s h Φ d x dt τ i ∂τ i τ i τ i   0 � ( H − ) 2 − ( H − ∗ ) 2 � ∗ , u + ∗ , u + n 1 g F WB ( u − ∗ , H − , n ) = F num ( u − ∗ , n ) +   2 � ( H − ) 2 − ( H − ∗ ) 2 � n 2 g 2 Hydrostatic reconstruction Audusse et al. ’04 u ± ( H ± ∗ , H ± ∗ · ( v 1 ) ± , H ± ∗ · ( v 2 ) ± ) T , = ∗ � 0 , H ± + b ± − max � b − , b + �� H ± = max ∗

Well-balanced Scheme � � � � d ∗ , u + F WB ( u − u h Φ d x + ∗ , H − , n ) Φ d σ − F ( u h ) ·∇ Φ d x = s h Φ d x dt τ i ∂τ i τ i τ i   0 � ( H − ) 2 − ( H − ∗ ) 2 � n 1 g F WB ( u − ∗ , u + ∗ , H − , n ) = F num ( u − ∗ , u + ∗ , n ) +   2 � ∗ ) 2 � ( H − ) 2 − ( H − n 2 g 2 Hydrostatic reconstruction Audusse et al. ’04 Higher order quadrature necessary � 1 � 2 gH 2 � � � x Φ d x = τ gHH x Φ d x = − τ gHb x Φ d x τ � � � − 1 � lm ) 2 · lm ) 2 + ǫ WB filter: Indicator ω res u i u i = (ˆ (ˆ i l + m = N l + m < N lm = ˆ lm + ˆ u i H i b i now depending on ˆ lm

Small Perturbation Test Smooth bottom elevation on Ω = [0 , 2] × [0 , 1] : b ( x ) = 0 . 8 e − 5( x − 0 . 9) 2 − 50( y − 0 . 5) 2 Initial water height � 1 − b ( x ) + 0 . 01 , if 0 . 05 ≤ x ≤ 0 . 15 , H ( x , 0) = 1 − b ( x ) , else . Initial velocity: v = 0

Small Perturbation Test Water surface w = H + b (for N = 2 , K = 46360 , α i = 5 N ∆ t h i , p = 1) t = 0 . 12 t = 0 . 24 t = 0 . 36 t = 0 . 48 Problem: boundary conditions (30 contour levels from 0 . 99 to 1 . 01)

Positivity Treatment Xing/Zhang/Shu-approach assumes Non-negative cell means ¯ H n i at time t n Positivity preserving flux (e.g. HLL, Lax-Friedrichs) Non-negative values of H at certain quadrature nodes Zhang, Xia, Shu ’11: Maximum-principle-satisfying and PP DG schemes on triangular grids (scalar & Euler equations) ⇒ Non-negative cell means at next time level t n +1 � i − ∆ t H n +1 ¯ = ¯ H n F WB ( u n , − , u n , + , n ) d σ 1 ∗ ∗ i | τ i | ∂τ i Under suitable CFL-like time step restriction

Positivity Treatment Xing/Zhang/Shu-approach assumes Non-negative cell means ¯ H n i at time t n Positivity preserving flux (e.g. HLL, Lax-Friedrichs) Non-negative values of H at certain quadrature nodes ⇒ Non-negative cell means at next time level t n +1 � i − ∆ t H n +1 ¯ = ¯ ( u n , − , u n , + H n F WB , n ) d σ 1 ∗ ∗ i | τ i | ∂τ i Under suitable CFL-like time step restriction � if H < 10 − 6 0 Comp. velocities via v = 2 H · ( H v ) else (Bryson et al.) H 2 +max { H 2 ,ǫ } � − s i α i η 2 p � u i , mod u i Filter modification: ˆ = exp · ˆ lm lm H i ) − 1 · min � 1000(5 N 4 + 1) ω res � s i = (¯ , 0 i

Oscillating Lake Paraboloidal vessel on Ω = [ − 2 , 2] 2 : b ( x ) = 0 . 1( x 2 + y 2 ) . 2D cut at y = 0 Periodic analytical solution known H ( x , t ) = max { 0 , 0 . 05(2 x cos( ω t ) + 2 y sin( ω t )) + 0 . 075 − b ( x )) } , v 1 ( x , t ) = − 0 . 1 ω sin( ω t ) , � v 2 ( x , t ) = 0 . 1 ω cos( ω t ) , ω = 0 . 2 g

Oscillating Lake Water surface w = H + b (for N = 2 , K = 23138 , α i = 10 N ∆ t h i , p = 1) t = T per / 6 t = T per / 3 t = T per / 2 t = T per

Positivity Treatment – Implicit Euler Case Extension of Xing/Zhang/Shu-approach Non-negative cell means ¯ H n i at time t n Positivity preserving flux (e.g. HLL, Lax-Friedrichs) Non-negative values of H at certain quadrature nodes ⇒ Non-negative cell means at next time level t n +1 � i − ∆ t H n +1 ¯ = ¯ ( u n +1 , − , u n +1 , + H n F WB , n ) d σ 1 ∗ ∗ i | τ i | ∂τ i No time step restriction! Implicit Euler is unconditionally SSP Higueras ’05