Discontinuous Galerkin and spectral element method for rotating shallow water equation on the sphere Praveen Chandrashekar praveen@math.tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore-560065, India http://cpraveen.github.io Computational Science Symposium Department of Computational and Data Sciences, IISc Bangalore 16-18 March 2017 Supported by Airbus Foundation Chair at TIFR-CAM, Bangalore http://math.tifrbng.res.in/airbus-chair 1 / 56
Shallow water model 1 , 2 • Shallow atmosphere/ocean ◮ 3 / 4 of total mass within 11 km ◮ Circumference = 40,000 km • Assumptions ◮ Incompressible flow ◮ Hydrostatic equilibrium in radial direction 1 Pedlosky, Geophysical Fluid Dynamics 2 Vallis, Atmospheric and Oceanic Fluid Dynamics 2 / 56
Rotation, spherical coordinates Fig. 2.3 The spherical coordinate sys- tem. The orthogonal unit vectors i , j and k point in the direction of in- creasing longitude λ , latitude ϑ , and altitude z . Locally, one may apply a Cartesian system with variables x , y and z measuring distances along i , j and k . from Vallis, Atmospheric and Oceanic Fluid Dynamics Radius: R = 6 . 37122 × 10 6 m , Rotation rate: Ω = 7 . 292 × 10 − 5 s − 1 Acc. due to gravity: g = 9 . 80616 m · s − 2 3 / 56
Rotating shallow water model v = velocity , H = height of atm. rel. to mean surface H s = height of ground rel. to mean surface , D = H − H s = depth Vector invariant form ∂ v ∂t + ∇ Φ + ( ω + f ) v ⊥ = 0 ∂D ∂t + ∇ · ( v D ) = 0 where K = 1 2 | v | 2 , Φ = gH + K, ω = k · ∇ × v v ⊥ = k × v , f = 2Ω sin θ (Note that k · v = 0 ) 4 / 56
Energy is conserved � ∂ v � � ∂D � ∂t + ∇ Φ + ( ω + f ) v ⊥ D v · + Φ ∂t + ∇ · ( v D ) = 0 leads to a conservation law ∂ � 1 2 D | v | 2 + 1 � �� gH + 1 � � 2 gH 2 2 | v | 2 + ∇ · v D = 0 ∂t Total energy is conserved � 1 � � 2 D | v | 2 + 1 2 gH 2 d s = const. S 5 / 56
Vorticity dynamics Absolute vorticity η := ω + f Vorticity equation: Take curl of the velocity equation ∂η � ∂t + ∇ · ( v η ) = 0 = ⇒ η d s = const. S Potential vorticity q := η D PV is advected by the flow ∂q ∂t + v · ∇ q = 0 q is constant along streamlines = ⇒ minimum/maximum principle S q ( x ′ , y ′ , z ′ , 0) ≤ q ( x, y, z, t ) ≤ max q ( x ′ , y ′ , z ′ , 0) min S 6 / 56
Potential enstrophy is conserved Another conservation law � η 2 � η 2 ∂ � � + ∇ · = 0 D v ∂t D Potential enstrophy is conserved η 2 � � Dq 2 d s = const. D d s = S S Infinite number of conserved quantities (Casimirs) � DF ( q )d s S for any function F , e.g., � Dq i d s = const. , i = 0 , 1 , 2 , . . . , S 7 / 56
Conservation form Momentum and mass conservation ∂ � D v ⊗ v + 1 � 2 gD 2 I ∂t ( D v ) + ∇ · = − f ( k × D v ) − gD ∇ H s ∂D ∂t + ∇ · ( D v ) = 0 Finite volume, SEM, DG 8 / 56
Vorticity-divergence form (Z grid) Tangential divergence: δ = ∇ · v ∂D ∂t + ∇ · ( v D ) = 0 ∂η ∂t + ∇ · ( v η ) = 0 ∂δ ∂t − ∇ × ( v η ) + ∇ 2 ( gH + | v | 2 / 2) = 0 −∇ 2 ψ = − ( η − f ) −∇ 2 χ = − δ ∇ χ + ∇ ⊥ ψ = v Usually solved with global pseudo-spectral methods 3 or finite volume method 3 Hack & Jakob (1992) 9 / 56
Hyperbolic property In Cartesian coordinates, with the z -axis along the axis of rotation ∂ U ∂ U ∂ U ∂y + ˜ ∂t + A 1 ∂x + A 2 S = 0 (1) where v 1 v 1 0 g v 2 0 0 , , U = v 2 A 1 = 0 v 1 0 A 2 = 0 v 2 g D D 0 v 1 0 D v 2 For any unit vector n = ( n 1 , n 2 ) −√ gn 1 √ gn 1 − n 2 v n 0 gn 1 −√ gn 2 √ gn 2 , A 1 n 1 + A 2 n 2 = 0 v n gn 2 R = n 1 √ √ Dn 1 Dn 2 v n 0 D D � � Eigenvalues: v n , v n − gD, v n + gD 10 / 56
Riemann solver based numerical flux ∂ U ∂t + ∂ F ∂x = 0 with initial condition � U l x < 0 U ( x, 0) = x > 0 U r Linearized problem (Roe, 1981) ∂ U ∂t + A ( U l , U r ) ∂ U A = R Λ R − 1 ∂x = 0 , Solve Riemann problem exactly and compute flux at x = 0 F lr = 1 2[ F ( U l ) + F ( U r )] − 1 2 R | Λ | R − 1 ( U r − U l ) 11 / 56
Previous works 4 • Finite volume Arakawa & Lamb (1981), Salmon (2007), Eldred & Randall (2016), Toy & Nair (2017) • DG, conservation form Giraldo et al. (2002) • DG, vector invariant form Nair et al. (2004) • SEM Taylor et al. (1997), Giraldo (2001), Taylor & Fournier (2010) • Compatible finite elements Natale et al. (2016) 4 See review article by Marras et al., Arch Comp Meth Eng, 2015 12 / 56
Cubed-sphere grid 5 ( X 1 , X 2 ) ∈ [ − a, + a ] 2 ( x 1 , x 2 , x 3 ) ∈ S 5 Sadourny (1972), Ronchi et al. (1996) 13 / 56
Cubed-sphere grid ( X 1 , X 2 ) ∈ [ − a, + a ] 2 ( x 1 , x 2 , x 3 ) ∈ S 14 / 56
Cubed-sphere grid ( χ 1 , χ 2 ) ( X 1 , X 2 ) ( x 1 , x 2 , x 3 ) χ i = arctan( X i /a ) ∈ [ − π/ 4 , π/ 4] , i = 1 , 2 X i = a tan χ i ∈ [ − a, + a ] 15 / 56
Cubed-sphere map From Nair et al., MWR, 133, 2005 Cube: [ − a, + a ] 3 √ a = R/ 3 � � On P 1 : ( x 1 , x 2 , x 3 ) = R x 2 x 1 , x 3 r ( a, X 1 , X 2 ) , ( X 1 , X 2 ) = a x 1 16 / 56
Finite element space • T h = cubed sphere grid • ˆ T = [0 , 1] × [0 , 1] is the reference cell • Each cell T ∈ T h obtained by mapping reference cell ˆ T F T : ˆ T → T, ( ξ 1 , ξ 2 ) → ( χ 1 , χ 2 ) → ( x 1 , x 2 , x 3 ) • Q N = 2-d tensor product polynomials of degree at most N • Space of broken polynomials for scalar and vector functions V N h = { ψ ∈ L 2 ( S ) : ψ | T ◦ F − 1 W N h = V N h × V N h × V N ∈ Q N } , T h • Lagrange polynomials using Gauss-Legendre nodes • Same nodes used for quadrature 17 / 56
Tangential gradient 6 Coordinates ξ = ( ξ 1 , ξ 2 ) ∈ ˆ T, x = ( x 1 , x 2 , x 3 ) ∈ S First fundamental form g ij = ∂ x · ∂ x , G = [ g ij ] ∂ξ i ∂ξ j and its inverse g ij = [ G − 1 ] ij Tangential gradient 2 2 g ij ∂x d ∂φ � � ( ∇ φ ) d = , d = 1 , 2 , 3 ∂ξ i ∂ξ j i =1 j =1 6 Dziuk & Elliott, Acta Numerica, 2013 18 / 56
DG for v − D model For each cell T ∈ T h d � � v h · ϕ h d s − Φ( v h , D h , H s ) ∇ · ϕ h d s d t T T � � ˆ F v · ϕ − ( ω h + f ) v ⊥ ∀ ϕ h ∈ W N + h d σ + h · ϕ h d s = 0 , h ∂T T d � � � ˆ F D ψ − ∀ ψ h ∈ V N D h ψ h d s − v h D h · ∇ ψ h d s + h d σ = 0 , h d t T T ∂T F v , ˆ ˆ F D are numerical fluxes obtained from Rusanov scheme 1 2(Φ − + Φ + ) n − − 1 2 λ ( v + − v − ) ˆ F v = 1 2( D − v − + D + v + ) n − − 1 2 λ ( D + − D − ) ˆ F D = where λ = max {| v − · n − | + gD − , | v + · n + | + � � gD + } 19 / 56
DG for v − D model Vorticity computed locally on each cell � � � h n − · ˆ ∇ ψ h · v ⊥ ψ − v ⊥ ∀ ψ h ∈ V N ω h ψ h d s = h d s − h d σ, h T T ∂T where v = 1 2( v − + v + ) ˆ 20 / 56
DG for v − D − η model For each cell T ∈ T h d � � v h · ϕ h d s − Φ( v h , D h , H s ) ∇ · ϕ h d s d t T T � � ˆ F v · ϕ − η h v ⊥ ∀ ϕ h ∈ W N + h d σ + h · ϕ h d s = 0 , h ∂T T d � � � ˆ F D ψ − ∀ ψ h ∈ V N D h ψ h d s − v h D h · ∇ ψ h d s + h d σ = 0 , h d t T T ∂T d � � � F η ψ − ˆ ∀ ψ h ∈ V N η h ψ h d s − v h η h · ∇ ψ h d s + h d σ = 0 , h d t T T ∂T F v , ˆ ˆ F D , ˆ F η are numerical fluxes obtained from a Riemann solver. 21 / 56
Numerical flux for v − D − η model In Cartesian coordinates, with the z -axis along the axis of rotation ∂ U ∂ U ∂ U ∂y + ˜ ∂t + A 1 ∂x + A 2 S = 0 (2) where v 1 v 1 0 g 0 v 2 0 0 0 v 2 0 v 1 0 0 0 v 2 g 0 U = , A 1 = , A 2 = D D 0 v 1 0 0 D v 2 0 η η 0 0 v 1 0 η 0 v 2 For any unit vector n = ( n 1 , n 2 ) v n 0 gn 1 0 0 v n gn 2 0 A 1 n 2 + A 2 n 2 = Dn 1 Dn 2 v n 0 ηn 1 ηn 2 0 v n 22 / 56
Numerical flux for v − D − η model � � Eigenvalues: v n , v n , v n − gD, v n + gD Corresponding linearly independent eigenvectors −√ gDn 1 √ gDn 1 0 − n 2 −√ gDn 2 √ gDn 2 0 n 1 R = 0 0 D D 1 0 η η Numerical flux (Φ − + Φ + ) n − v + 1 − v − 1 1 (Φ − + Φ + ) n − v + 2 − v − F = 1 − 1 ˆ 2 R | Λ | R − 1 2 2 ( D − v − + D + v + ) · n − D + − D − 2 ( η − v − + η + v + ) · n − η + − η − Rotate velocity flux to global ( x 1 , x 2 , x 3 ) coordinates. 23 / 56
Time integration: ˙ U = R ( U ) 3-stage, strong stability preserving, Runge-Kutta scheme U n − ∆ t · R ( U n ) U (1) = 3 4 U n + 1 4[ U n + ∆ t · R ( U (1) )] U (2) = 1 3 U n + 2 3[ U n + ∆ t · R ( U (2) )] U n +1 = After each stage, project the velocity to tangent plane at all GL points. CFL condition CFL h T ∆ t = 2 N 2 + 1 min �| v | + √ gD � L ∞ ( T ) T ∈T h 24 / 56
Recommend
More recommend
