Efficient High-Order Methods for Geophysical Fluid Dynamics Models * Frank Giraldo

1. Efficient High-Order Methods for Geophysical Fluid Dynamics Models* Frank Giraldo Department of Applied Mathematics Naval Postgraduate School,Monterey CA USA http://faculty.nps.edu/fxgirald Collaborators: Matthias Läuter (AWI, Potsdam), Marco Restelli (Max-Planck, Hamburg), Emil Constantinescu (Argonne NL, Chicago), Jim Kelly (NPS) *Funded by the Office of Naval Research: (1) Computational Mathematics and (2) Meteorology

2. Motivation for this Work The goal of this research is to construct the best possible climate and numerical weather prediction models for the atmosphere and ocean that: • Produce high-order accuracy solutions (or good quality); • Can use unstructured (adaptive) grids to better handle the physical geometry of the problem; • Are efficient to run on all types of computers; and • Scale well on parallel computers. • Producing good quality solutions efficiently requires different methods for different types of flow problems (e.g., continuous versus discontinuous methods). • Element-based Galerkin (EBG) methods (e.g., finite elements/volumes, SE, DG) allow for easy construction of the discrete derivatives on unstructured grids. • To construct efficient models on any type of computer requires the use of implicit, semi-implicit, or Lagrangian time-integrators (TI). • To scale well on parallel computers requires very good domain decomposition techniques. EBG methods are such methods but the matrix problem required in the TI must use better preconditioners to achieve linear scalability.

3. Talk Summary • Equation Sets being explored in this work • Hydrostatic Primitive Equations • Euler/Navier-Stokes Equations • Shallow Water Equations • Spatial Discretization • Continuous and Discontinuous Element-based Galerkin Methods • Time-Integrators • Explicit SSP Methods • Semi-implicit Methods • Fully-Implicit Methods IV. Results of the three models under development • Global Hydrostatic Atmospheric Model (NSEAM) • Mesoscale Atmospheric Models • Coastal Ocean Model

4. I. Equation Sets

5. Hydrostatic Primitive Equations(Euler Equations with no Vertical Acceleration) (Mass) (Momentum) (Energy/Entropy) (State)

6. Mesoscale Equations(Navier-Stokes Equations) (Mass) (Momentum) (Energy) (Equation of State) (Viscous Fluxes)

7. Coastal Ocean Equations(Shallow Water Equations) (Mass) (Momentum) (Geopotential)

8. II. Spatial Discretization

9. Galerkin Methods • Primitive Equations: • Approximate the solution as: • Interpolation O(N) • Write Primitive Equations as: • Weak Problem Statement: Find • such that • Integration O(2N) • Choice of  determines the method (e.g., global function defines spectral method, local function defines EBG) (CG) (DG)

10. Advantages of Element-based Galerkin Methods • Order of the method can be changed automatically (via Basis Functions) • For many problems, high-order is the most efficient way to reach a certain level of accuracy • Unstructured adaptive (conforming and non-conforming) grids can be used quite naturally (since the method is based on constructing discrete operators on elements/volumes • Excellent performance on MPP (especially in high-order mode) since the amount of on-processor work is huge while the communication footprint (stencil) is very small. Numerous Gorden Bell prizes at SC have been awarded to Spectral Element codes (at least 5 to my knowledge). • Easy to maintain and improve codes since basis functions can be whatever you wish. All one needs to do is to swap out the basis functions routines.

11. Comparison of EBG Methods Continuous Galerkin Methods • High order accurate yet local construction (via DSS) • globally conservative = good for hydrostatic primitive equations • Theoretically optimal for self-adjoint operators (elliptic equations) • Excellent for incompressible Navier-Stokes • also extremely good for non-self-adjoint operators (hyperbolic equations) • Simple to construct efficient semi-implicit time-integrators • Semi-implicit = nonlinear terms are explicit and linear terms are implicit Discontinuous Galerkin Methods • High order like Continuous Galerkin • Completely local in nature • no global assembly/DSS required as in CG (truly local) • High order generalization of the FV • locally and globally conservative = excellent choice for non-hydrostatic equations (mesoscale models) • upwinding and BCs implemented naturally (via Riemann solvers) • Designed for hyperbolic equations (i.e., shock waves) • Simple provision for elliptic equations (e.g., LDG) • Not so straightforward to construct efficient semi-implicit time-integrators, due to the inherent nonlinear nature of the numerical flux, until recently (Restelli-Giraldo, SIAM J. Sci. Comp. and Giraldo-Restelli, Int. J. Num. Meths Fluids)

12. Geometric Flexibility of EBG Methods Icosahedral Telescoping Banded Hexahedral Quadrilaterals Triangles Adaptive Adaptive Adaptive Icosahedral

13. Convergence Rate for EBG Methods(DG Spatial Operators) Note that the method achieves the expected convergence rate; that is, error =O(hN+1)

14. NSEAM Scalability of EBG Methods(Surface Values for T185 L26 during 0-30 days) Temperature Pressure 30 day simulation of a Baroclinic Instability for a fully 3D atmospheric model using 8th order polynomials

15. NSEAM Scalability of EBG Methods IBM SP4 (1.7 GHz) T498 (20 km) L60 DT=150 secs T249 (54 km) L30 DT=300 secs Note that the problem size from T249 L30 to T498 L60 has increased by a factor of 16 (22 hor x 2 vert x 2 time-steps). However, the Wallclock Time has only decreased by a factor of 8. Furthermore, NSEAM T498 L60 scales linearly with processor count! At T498 L60NSEAM can theoretically accommodate 70,000 processors!

16. III. Time-Integrators

17. Explicit Time-Integrators • Let’s rewrite the governing equations as • SSP-RK(2,3,4) For k=1,…,K • SSP-BDF2 Which we write as or in general, more compactly, as

18. Semi-Implicit Time-Integrator(Building an implicit method on top of an explicit one) • Let’s rewrite the governing equations as • Now, if we knew the linear operator L, then we could write • Discretizing by Kth order time-integrator yields Where and

19. Semi-Implicit (IMEX) Time-Integrators (BDFK)(SWE: Linear Kelvin Wave) Performance Accuracy DG with N=10

20. BDFK Time-Integrators(Stability Regions) Explicit Implicit

21. Additive Runge-Kutta Time-Integrators(Stability Regions) ARK3 ARK3 RK3 BDF3 RK2 Implicit Explicit

22. Semi-Implicit (IMEX) Time-Integrators(Constructing the Schur Complement) • The implicit problem that we have is: • Where the dims are: • For 2D Euler d=4, 3D d=5, etc. • For 2D Euler we solve a 16N2 system • For 2D SWE d=3, and solve a 9N2 • A better approach is to write out the system as follows (for 2D SWE): • Applying block LU decomposition: • Where the Schur Complement is: • And the dimensions are:

23. Semi-Implicit (IMEX) Time-Integrators(Important Properties of Schur Complement) • The original system is: • The Schur Complement system is: • This system is clearly smaller than the original system (NxN instead of 3Nx3N for 2D SWE). • Equally important is that the Schur System is better conditioned than the original system. This means that fewer iterations are required by an iterative solver to reach convergence. • Key Point: A semi-implicit method should be more efficient than the most efficient explicit methods, but with the addition of the Schur Complement system, the gains are much bigger. Thus, whenever possible we must strive to find such a system. This is not always very easy.

24. Semi-Implicit (IMEX) Time-Integrators(Which Methods offer a Schur Complement?) • It turns out that it can be done as long as one can write the method in the IMEX form: • Additive Runge-Kutta Methods can be written in the Schur Complement form. This now extends the order of accuracy in time as well as the maximum time-step allowed for efficiency (big gains in the explicit part for increasing K).

25. 4. Fully-Implicit Time-Integrators • Writing the governing equations as • We then discretize this implicitly using implicit methods (eg., BDFK, IRK) • Where the matrix problem is now • Which clearly requires a nonlinear implicit solution(eg., Newton-Krylov methods) • We write the problem as the functional: • And then solve the nonlinear problem: • With the resulting Linear System:

26. Comparison of Various Time-Integrators(Rising Thermal Bubble) z x

27. IV. Results of the Three Models

28. NSEAM IV. Results of the Three ModelsSEGlobal Hydrostatic ModelCollaborators: Y.J. Kim, Maria Flatau, Chi-Sann Liou, and Melinda Peng (NRL-Monterey)

29. Aqua-Planet Experiments Surface Temperature & Convective Precipitation Day 0~180 Control Control

30. Control 3.75° x 2.5° L19 w/ reduced Hyperviscosity (=5e5) 180 0 180 0E 0E 180 180 0W 0W Neale & Hoskins (2001) Convective Precipitation (-5°~5°) Control Control UKMO NSEAM E6P8L20 (~2.2°~T54) NICAM Cloud Resolving Model Tomita

31. IV. Results of the Three ModelsSE/DGNonHydrostatic Atmospheric ModelCollaborators: Matthias Läuter (AWI, Potsdam), Emil Constantinescu (Argonne NL, Chicago), Marco Restelli (MPI,Hamburg), Saša Gaberšek (NRL, Monterey), Jim Doyle (NRL, Monterey), Jim Kelly (NPS)

32. Rising Thermal Bubble(nelx=nely=20 N=10, 5 m res, Time=700 seconds) z x Profile along x=500 meters Potential Temperature

33. Inertia-Gravity Wave(nelx=120, nely=4 N=10, 250 m res, Time=3000 seconds) z x Potential Temperature Profile along z=5000 meters

34. Density Current(nelx=128, nely=32 N=8, 250 m res, Time=900 seconds) z x Profile along z=1200 meters Potential Temperature

35. Linear Hydrostatic Mountain(nelx=20, nely=10 N=10, 250 m res, Time=10 hours) z Vertical Velocity Momentum Flux Vertical Velocity x

36. Linear Nonhydrostatic Mountain Waves(360 x 310 meter resolution with 10th order polynomials) Vertical Velocity

37. Convective Storm Simulation(SE2NC-SIBDF2: 500 x 200 m res, N=10, Time=10 hours) Vapor Rain Clouds

38. IV. Results of the Three ModelsTriangular DGCoastal Ocean ModelCollaborators: Tim Warburton (Rice), Marco Restelli (MPI, Hamburg), Dimitrios Alevras (Hellenic Navy Hydrographic Service, Athens), Jörn Behrens (Univ. of Hamburg)

39. Propagation of the 2004 Indian Ocean Tsunami(Grid Dimensions: Np=66715, Ne=130444, N=1, K=2, J=1) y x Time evolution of the water surface height (grid data provided by J. Behrens, AWI-Bremerhaven, and data formatted by D. Alevras, NPS)

40. Propagation of the 2004 Indian Ocean Tsunami(Grid Dimensions: Np=66715, Ne=130444, N=1, K=2, J=1) y

41. Riemann Problem on Sphere(Discontinuous Galerkin Shallow Water) Height V-Velocity U-Velocity

42. Current Work and Future Directions • Continuous and Discontinuous EBG methods are a good choice for the development of next-generation GFD models • To make this new class of models ready for operational use requires the introduction of good Time-Integrators such as semi-implicit, fully-implicit and Lagrangian methods • Current Areas of Research are: • Time-Integration • Semi-Implicit Runge-Kutta Methods (and Rosenbrock Methods) • Fully-Implicit Methods • Multi-rate methods based on extrapolation • Preconditioners • Adaptivity • Conforming and non-conforming for triangles • Non-conforming for quadrilateral • Physical Parameterizations • Continue testing simple Kessler microphysics • Positivity-preserving schemes for moisture • Extension to more complicated physics

43. Semi-Implicit (IMEX) Time-Integrators(Additive RK Methods) • Starting with: • We can then write this in the IMEX form: • Discretizing by Kth order time-integrator yields For k=1,…,K Where and