Recent progresses in the development of a new generation adaptive DG - PowerPoint PPT Presentation

Recent progresses in the development of a new generation adaptive DG dynamical core for RegCM Giovanni Tumolo The Abdus Salam ICTP and OGS - Trieste, < gtumolo@ictp.it > Trieste, May 24, 2016

◮ In collaboration with Luca Bonaventura (MOX-Politecnico di Milano) ◮ with thanks to ◮ Filippo Giorgi (ICTP Abdus Salam) ◮ Graziano Giuliani (ICTP Abdus Salam) ◮ Marco Restelli (MPI for Plasma Physics) ◮ and acknowledgements for funding from ◮ The Abdus Salam International Centre for Theoretical Physics ◮ Istituto Nazionale di Oceanografia e Geofisica Sperimentale (whithin the HPC-TRES framework) ◮ The INdAM-GNCS ◮ The Royal Meteorological Society

Outline ◮ Motivation and introduction to the p-SISLDG formulation. ◮ Review of a novel SISL time integration approach. ◮ 1st extension: mass conservative mixed Eulerian/semi-Lagrangian variant. ◮ 2nd extension: meshes of deformed quadrilaterals on the sphere and on a vertical plane. ◮ p − adaptivity. ◮ Numerical validation: ◮ horizontal: ◮ efficiency gain by TR-BDF2: unsteady flow with analytic solution ◮ efficiency gain by p − adaptivity: Williamson’s test 6 ◮ effects of the mesh deformation on the solution (mass conservative version): Williamson’s test 2 and unsteady flow with analytic solution ◮ p − adaptive tracers transport: ◮ Solid body rotation ◮ Deformational flow ◮ Coupling with SWE solver: advection by Rossby-Haurwitz wave ◮ vertical: ◮ Interacting bubbles test ◮ Linear hydrostatic lee waves ◮ Nonlinear nonhydrostatic lee waves ◮ Where is the 3D dycore?? Current status and HPC requirements for using p-SISLDG as a three dimensional dynamical core. Further plans and conclusions.

Motivation ◮ Goal: use DG methods for the design of a new generation dynamical core for the regional climate modelling system RegCM of ICTP . ◮ This is challenging for: stability restrictions with explicit time stepping: ◮ “The RKDG algorithm is stable provided the following condition holds: u ∆ t 1 < 2 p + 1 h where p is the polynomial degree; (for the linear case this implies a CFL limit 1 3 )” Cockburn-Shu, Math. Comp. 1989 computational cost : DG requires more d.o.f. per element than CG . ◮ ◮ How to increase computational efficiency of DG ? coupling DG to semi-implicit semi-Lagrangian (SI-SL) technique (no CFL) ◮ introduction of p − adaptivity (flexible degrees of freedom) ◮ ◮ = ⇒ p-SISLDG method.

p-SISLDG: main features Main novel features of the proposed p-SISLDG formulation: ◮ is the first unconditionally stable DG formulation for the shallow water and for the Euler equations, ◮ is based on the first fully second order two-time-level SISL time integrator, ◮ is the first extensive application of p − adaptivity strategies in NWP , ◮ employs a unified discretization approach for the horizontal and vertical.

A novel approach for SISL time integration: TR-BDF2 Given a Cauchy problem for a system of ODEs: y ′ = f ( y , t ) , y ( 0 ) = y 0 , the TR-BDF2 method is defined by the two following implicit stages (Bank et al. IEEE trans. 1985): u n + 2 γ − γ ∆ t f ( u n + 2 γ , t n + 2 γ ∆ t ) u n + γ ∆ t f ( u n , t n ) , = u n + 1 − γ 2 ∆ t f ( u n + 1 , t n + 1 ) ( 1 − γ 3 ) u n + γ 3 u n + 2 γ , = with γ ∈ ( 0 , 1 / 2 ] fixed implicitness parameter and 1 − 2 γ γ 3 = 1 − γ 2 γ 2 = 2 ( 1 − γ ) , . 2 γ

SL reinterpretation of TR-BDF2 If suitable semi-Lagrangian approximate evolution operators for scalar and vector valued functions are introduced: [ E ( t n , ∆ t ) G ]( x ) = G ( t n , x D ) � t n + 1 u n � X ( t ; t n + 1 , x ) dt and X ( t ; t n + 1 , x ) is the solution of: � where x D = x − t n dt X ( t ; t n + 1 , x ) = u n � � � d X ( t ; t n + 1 , x ) , X ( t n + 1 ; t n + 1 , x ) = x x t n +1 X ( t n ; t n +1 , ( x )) X ( t n ; t n +1 , · ) t n i.e. two steps are required to compute [ E ( t n , ∆ t ) G ] ( x ) : 1. departure point x D computation ( e.g. McGregor, Mon. Wea. Rev.,1993); 2. interpolation of G n at departure point.

SL reinterpretation of TR-BDF2 ... and if governing equations in advective form are to be solved: (being D Dt the Lagrangian derivative operator) (SWE) Shallow Water Eqs. (no (VSE) Euler eqs. (no Coriolis force) on a Vertical Slice ( ∂ Coriolis force): ∂ y = 0): D Π � c p � Dt + c v − 1 Π ∇ · u = 0 , Dh Dt + h ∇ · u = 0 , Du Dt + c p Θ ∂π ∂ x = 0 , D u Dw Dt + c p Θ ∂π ∂ z − g θ Dt + g ∇ h = − g ∇ b , θ ∗ = 0 , D θ Dt + w d θ ∗ dz = 0 . � p � p � − R / c p , Π = � R / c p , with h , u = ( u , v ) T and b being fluid with Θ = T p 0 p 0 p , T , u = ( u , w ) T , pressure, temperature and depth, horizontal velocity and vertical velocity, c p , c v , R specific heats and gas bathymetry elevation respectively, constant of dry air, and Π( x , y , z , t ) = π ∗ ( z ) + π ( x , y , z , t ) , Θ( x , y , z , t ) = θ ∗ ( z ) + θ ( x , y , z , t ) ,

... then SISL-TR steps for SWE and VSE are isomorphic π n + 2 γ + γ ∆ t ( c p / c v − 1 ) Π n ∇ · u n + 2 γ = − π ∗ h n + 2 γ + γ ∆ t h n ∇ · u n + 2 γ = t n , 2 γ ∆ t � � t n , 2 γ ∆ t E [ h − γ ∆ t h ∇ · u ] , � � + E [Π − γ ∆ t ( c p / c v − 1 ) Π ∇ · u ] , u n + 2 γ + γ ∆ t c p Θ n ∂π n + 2 γ = ∂ x � u − γ ∆ t c p Θ ∂π � E ( t n , 2 γ ∆ t ) , ∂ x u n + 2 γ + γ ∆ t g ∇ h n + 2 γ = − γ ∆ t g ∇ b t n , 2 γ ∆ t � � + E { u − γ ∆ t [ g ( ∇ h + ∇ b )] } . d θ ∗ n + 2 γ � 1 + ( γ ∆ t ) 2 g � w n + 2 γ + γ ∆ tc p Θ n ∂π = θ ∗ dz ∂ z c p Θ ∂π ∂ z − g θ � � �� E ( t n , 2 γ ∆ t ) w − γ ∆ t θ ∗ θ − γ ∆ t d θ ∗ + γ ∆ t g � � θ ∗ E ( t n , 2 γ ∆ t ) dz w . h ← → π, u ← → u , v ← → w .

then SISL-BDF2 steps for SWE and VSE are isomorphic π n + 1 + γ 2 ∆ t ( c p / c v − 1 ) Π n + 2 γ ∇ · u n + 1 = h n + 1 + γ 2 ∆ t h n + 2 γ ∇ · u n + 1 = − π ∗ + ( 1 − γ 3 )[ E t n , ∆ t � � Π] t n , ∆ t � � � � 1 − γ 3 E h t n + 2 γ ∆ t , ( 1 − 2 γ )∆ t t n + 2 γ ∆ t , ( 1 − 2 γ )∆ t � � + γ 3 [ E Π] , + γ 3 E � � h , n + 1 u n + 1 + γ 2 ∆ t c p Θ n + 2 γ ∂π = ∂ x t n , ∆ t � � ( 1 − γ 3 )[ E u ] u n + 1 + γ 2 ∆ t g ∇ h n + 1 = t n + 2 γ ∆ t , ( 1 − 2 γ )∆ t + γ 3 [ E � � u ] , − γ 2 ∆ t g ∇ b t n , ∆ t + � 1 − γ 3 � E � � u d θ ∗ n + 1 t n + 2 γ ∆ t , ( 1 − 2 γ )∆ t 1 + ( γ 2 ∆ t ) 2 g w n + 1 + γ 2 ∆ t c p Θ n + 2 γ ∂π � � � � + γ 3 E u . = θ ∗ dz ∂ z t n + 2 γ ∆ t , ( 1 − 2 γ )∆ t t n , ∆ t ( 1 − γ 3 )[ E � � w ] + γ 3 [ E � � w ] + γ 2 ∆ t g t n + 2 γ ∆ t , ( 1 − 2 γ )∆ t t n , ∆ t � ( 1 − γ 3 )[ E � � θ ] + γ 3 [ E � � θ ] � θ ∗ h ← → π, u ← → u , v ← → w .

1st Extension: mass conservation, SISL-TR-BDF2 time discretization Considering the continuity equation in Eulerian flux form, while the momentum one in advective vector form: ∂η ∂ t = − ∇ · ( h u ) , D u Dt = − g ∇ η − f k × u , then, the TR stage of the SISL time discretization of previous equations is: η n + 2 γ + γ ∆ t ∇ · = η n − γ ∆ t ∇ · ( h n u n ) , � h n u n + 2 γ � u n + 2 γ + γ ∆ t g ∇ η n + 2 γ + f k × u n + 2 γ � � t n , 2 γ ∆ t � � = E [ u − γ ∆ t ( g ∇ η + f k × u )] . The TR stage is then followed by the BDF2 stage: η n + 1 + γ 2 ∆ t ∇ · ( h n + 2 γ u n + 1 ) = η n + γ 3 η n + 2 γ , � � 1 − γ 3 u n + 1 + γ 2 ∆ t g ∇ η n + 1 + f k × u n + 1 � � t n + 2 γ ∆ t , ( 1 − 2 γ )∆ t t n , ∆ t � � � � � � = 1 − γ 3 E u + γ 3 E u .

DG space discretization ◮ Defined a tassellation T h = { K I } N I = 1 of domain Ω and chosen ∀ K I ∈ T h three integers p π I ≥ 0, p θ I ≥ 0, p u I ≥ 0, at each time level t n , we are looking for approximate solution s.t. � � f ∈ L 2 (Ω) : f | K I ∈ Q p π h n , π n ∈ P h := I ( K I ) � � f ∈ L 2 (Ω) : f | K I ∈ Q p θ θ n ∈ T h := I ( K I ) � � g ∈ L 2 (Ω) : g | K I ∈ Q p u u n , v n , w n ∈ V h := I ( K I ) , ◮ modal bases are used to span P h , T h , V h , ◮ L 2 projection against test functions (chosen equal to the basis functions), ◮ introduction of (centered) numerical fluxes, ◮ substitution of velocity d.o.f. from momentum eqs. into the continuity eq., ◮ give raise, at each SI step, to a discrete (vector) Helmholtz equation in the fluid depth / pressure unknown only, i.e. a sparse block structured nonsymmetric linear system is solved by GMRES with block diagonal (for the moment) preconditioning.

Potential of p − adaptivity for atmospheric modelling applications ◮ No remeshing required of many physical quantities like orography profiles, data on land use and soil type, land-sea masks. ◮ Completely independent resolution for each single model variable. ◮ Easier coupling with SL technique, especially on unstructured meshes (no need to store two meshes). ◮ Possibility also of static p-adaptation: e.g. reduced p as counterpart of reduced grid, i.e. locally imposed p controlling the local Courant number ( = ⇒ significant # gmres-iterations reduction). ◮ Main potential problem: dynamic load balancing is mandatory for massively parallel implementations.