Recent developments in Numerical Methods, incl. report from the PDEs on the sphere meeting

1. Recent developments in Numerical Methods, incl. report from the PDEs on the sphere meeting (Part II) WGNE-Meeting 18.-22. Oct. 2010, Tokyo Michael Baldauf (DWD)

2. Spatial discretization methods • Finite Volume schemes • advantages • 'automatically' have local conservation properties of the prognostic variables • designed to handle shocks and other discontinuities • can easily made to satisfy monotonicity and positivity constraints • standard method in many CFD fields • presentations: • Riemann solvers for 2D shallow water equations [Castro] • convergence proofs for FV methods on Riemannian manifolds [T. Müller] • higher order Riemann solvers, proper choice of the flux formulation (AUSM+-up is superior to the Rusanov flux) [Ullrich]

3. Spatial discretization methods • Riemann solver for shallow water equations (Tsunami-simulation) • [Castro] • ADER scheme • very accurate treatment of shocks and other discontinuities (solution of a generalized Riemann problem) • well-balancing ('lake at rest'): exactly by numerical reconstruction of the bathymetry • local time stepping for higher efficiency homogeneous Riemann-problem generalized Riemann-problem from Castro, Toro (2008) JCP

4. Mediterranean tsunami from: C. Castro, talk at 'PDE's on the sphere', Aug. 2010, Potsdam

5. Spatial discretization methods Discontinuous Galerkin Method Arbitrary system of PDE's: Weak formulation of the PDE-system with test functions v(r): Representation by base functions pl(r): This leads to an ODE-System for the coefficients Qj,l(k)(t) • To choose: • flux formulation for boundary integral (Lax-Friedrich, Rusanov, AUSM+, ...) • quadrature formula (Gauss-Lobatto, ...) • time integration scheme (Runge-Kutta, implicit, ...)

6. Spatial discretization methods • Discontinuous Galerkin Method • Seek weak solutions of a balance equation(correspondance to finite volume methods  conservation) • Expand solution into a sum of base functions on each grid cell (correspondance to finite element methods  higher order)(classical Finite-volume-scheme is recovered by base function = const.) • useable on arbitrary grids  suitable for complex geometries • discontinuous elements  mass matrix is block-diagonal • in combination with an explicit time integration scheme (e.g. Runge-Kutta  RKDG-methods)  highly parallelizable code • but: how to solve vertically expanding sound waves efficiently? • in general: max. allowable Courant numbers are quite small • Applicationto LES [A. Müller] • hp-adaptive methodforshallowwaterequations[Blaise]:DG well suitedfor adaptive meshrefinement (AMR) • shallow-waterequations on thespherewithcurvilinearcoordinates[Nair] 6 M. Baldauf (DWD) 20.04.2010

7. Spatial discretization methods Spectral elements (SE) method = continuous Galerkin method for non-hydro. compressible models (global, regional) 3D semi-implicit (GMRES) benefits: 1.) higher order accuracy 2.) flexibility on unstructured and adaptive grids 3.) efficient on massively parallel computers [Kelly] 'High Order Methods Modeling Environment' (HOMME) [Levy] HOMME-spectral element method needs quadrilateral grids  cubed sphere refinement Finite element methods more often used in oceanographic applications [Danilov], [Sidorenko]: horizontally highly unstructured and stretched grids

8. Equation systems for 3D atmospheric models • non-hydrostatic, compressible equation sets • anelastic equation sets • for global CSU-model development: [Heikes] • anelastic Z-grid model (i.e. vertical vorticity , horizontal divergence , vertical velocity w, Exner function ') • anelastic vector-vorticity-eqns.  (Jung, Arakawa (2008) MWR) • Lipps, Hemler (1982) eqns. [Smolarkiewicz], [Szmelter]

9. Mesh Refinement (AMR) Adaptive (AMR) - 'amatos' (software package (Behrend) use of space filling curves with DG discretization [A. Müller] - wavelets [Dubos] Static: - by non-conformal refinement in triangular grid in ICON model [Zängl] - HOMME - by stretching function of hexagon grid in MPAS model [Klemp]

10. Grid structure (schematic view) ICON (G. Zängl, DWD) Triangles are used as primal cells Mass points are in the circumcenter Velocity is defined at the edge midpoints Red cells refer to refined domain Boundary interpolation is needed from parent to child mass points and velocity points

11. Mesh refinement in ICON (G. Zängl, DWD) Temperature at lowest model level on day 14 70 km 35 km 70 km, nested nest, 35 km

12. staggered vs. unstaggered grids • staggered grids have: • (+) advantages in the dispersion relation (high wavenumbers) • (--) disadvantages towards unstaggered grids: • stronger timestep constraint • horizontal-vertical splitting a bit easier in unstagg. grids • in the vicinity of changes in the grid resolution, much stronger wave reflection than in unstaggered grids can occur [Ullrich], [Thuburn]important for AMR! • change of grid position in one time-step: staggered for gravity wave, unstaggered elsewhere: good dispersion relations and no Gibbs-phenomena in CCAM [McGregor] • probably the distinction between these two grids types becomes less important for higher order methods (?)

13. ICON-project [Baldauf] Linearised shallow water equations analytic dispersion relation ( ): Rossby-deformation radius: LR= c / f inertial-gravity wave (Poincaré wave) geostrophic mode • 'Standard' C-grid spatial discretisation: • divergence: sum of fluxes over the edges (Gauss theorem) • gradient: centered differences ( -> 2nd order ) • Coriolis term: Reconstruction of v by next 4 neighbours with RBF-vector reconstruction (with an 'idealised' RBF-fct. =1) • temporal discretisation: assumed exact

14. ICON-project Triangle grid [Baldauf] 5 branches: 1,2 =  fI-G(k), 3,4 =  fHF(k), 5=0 (exact in equilat. grid!) ‚high-frequent mode‘ interpretation 1: internal excitation of an elementary cell interpretation 2: (poor) resolution of shorter waves (until 'x') ( not artificial)  consequences in practice: reduction of time step 1st Brillouin-zone

15. ICON-project Hexagon grid [Baldauf]  analytic solution inertial- gravity wave f geostrophic mode ka 4 branches: 1,2 = ..., inertial-gravity-waves 3,4  0 false geostrophic mode with 4-point reconstruction (Nickovic et al., 2002) Thuburn (2008), Thuburn, Ringler, Skamarock, Klemp (2009): 8-point reconstruction of Coriolis term  correct geostrophic mode 3,4 =0

16. Advection schemes • flux-form Semi-Lagrangian scheme (Miura, 2007) [Reinert], [Mittal] 1. Rekonstruktion of the mass contained in this area 2. flux through the lower boundary from: D. Reinert (DWD) runs successfully in ICON

17. Advection schemes • Runge-Kutta based advection scheme [Skamarock]: • requirements: scheme for icosahedral hexagon grid, should match with the explicit Runge-Kutta time integration scheme • idea: higher order flux-formulation of Hundsdorfer et al. (1995) in 1D can be extended to unstructured grid by the reformulation least-square fit to reconstruct the 2nd derivative this produces much less phase/ampitude errors than 2nd order scheme dramatically reduce phase errors in baroclinic wave test case in MPAS and ICON

18. Thank you for attention

19. Homepage of PDE's 2010 'Workshop on the Solution of Partial Differential Equations on the Sphere' 24.-27. April 2010 in Potsdam, Germany organized by Alfred-Wegener-Institute (Bremerhaven, Potsdam ), Konrad-Zuse-Institute, Berlin Purposes: A fundamental issue in climate and weather modeling is the dependence of errors in general circulation models on their basic dynamical cores and their subgrid scale parameterizations. To advance the state-of-the-art for numerical weather prediction and climate simulation, this organization encourages research on numerical algorithms and solution methods for partial differential equations (PDEs) posed in a spherical geometry. We also foster the development of tests and diagnostics for atmospheric and ocean model dynamical cores. The development of new methods in spherical geometry and how these eliminate shortcomings of current methods include applications to the shallow water equations and to complete hydrostatic and non-hydrostatic baroclinic models. Topics of interests: 1. Advection schemes 2. Discretization methods and adaptive grids 3. Comparison study of methods 4. Computational performance 5. Discussion and proposal of test cases 6. Global hydrostatic and non-hydrostatic models

20. 'Workshop on the Solution of Partial Differential Equations on the Sphere' 24.-27. April 2010 in Potsdam, Germany Topics of interests: 1. Advection schemes: monotone advection schemes, TVD methods, flux correction methods, and Eulerian and semi-Lagrangian methods applied to three dimensional baroclinic equations, two dimensional shallow water equations, and pure transport in spherical geometry and/or plane geometry. 2. Discretization methods and adaptive grids: spectral transform methods, gridpoint methods using polyhedral grid systems such as icosahedral and cubed-sphere systems as well as the latitude longitude grids, spectral element methods, discontinuous Galerkin methods, and alternative discretizations on regular and unstructured adaptive grids. 3. Comparison study of methods: the presentation of results from current or new methods, and the comparison of these methods with standard schemes. 4. Computational performance: the performance of current methods of choice and new methods on parallel computers, including limitations associated with algorithm and computer design are appreciated; implementation details for distributed and shared memory architectures as well as discussion of attendant difficulties are of interest; characteristics for parallel computers that may affect computational algorithms, in particular, future design aspects could be considered as they might affect algorithm implementation. 5. Discussion and proposal of test cases: the development and evaluation of new test methods, in particular, the potential ability of test cases to evaluate dynamical cores of atmospheric and ocean general circulation models; summaries of experiences with existing test cases, including what has been learned and what can be learned from such tests, and the introduction of additional test cases and experiences with them. 6. Global hydrostatic and non-hydrostatic models: the development of global models based on the hydrostatic and non-hydrostatic equations; now that recent computer facilities enable us to perform very high-resolution simulation over the globe, we are in a period of transition of the choice of the equations system.

21. Time integration methods Semi-Implicit Semi-Lagrangian models UK MetOffice: lat-lon-grid Conformal-cubic atmosphere model (CCAM) [McGregor] cubed sphere; stretching Schmidt transformation Russia: SL-AV model [Tolstykh] France: ALADIN model multigrid solver facilitated by the recursive data structure of the icosahedral grid (used e.g. for anelastic eqns. in the CSU-model [Heikes])

22. Time integration methods • Split-Explicit models • distinguish slow (advection) and fast (sound, gravity waves) dynamical processes  different time steps • VCAM [McGregor] • COSMO • NICAM [Tomita] • Explicit (HE/VI) Models • ICON-model [Zängl], [Gassmann] • MPAS[Klemp] Exact integration by Laplace-transformation of fast processes, no orographic resonance even with coupling with a Semi-Lagrangian transport scheme [Clancy] Lynch (1985)