04, 2016, Reading Unified Framework for Discrete Integrations of Compressible and Soundproof PDEs of Atmospheric Dynamics 2016 Slide 1 ECMWF Part I: Analytic Formulations: Introduction; Eulerian & Lagrangian reference frames; Atmospheric PDEs; Nonhydrostatic (all-scale) PDEs; Perturbation forms; Unified framework;

Conservation law forms; Generalised coordinates; Examples Part II: Integration Schemes: Forward-in-time integrators, E-L congruence; Semi-implicit algorithms; Elliptic boundary value problems; Variational Krylov solvers; Preconditioning; Boundary conditions 2016 Slide 2 ECMWF

Part I: Analytic Formulations Meteorology has a large portfolio of diverse analytic formulations of the equations of motion, which employ variety of simplifying assumptions while focusing on different aspects of atmospheric dynamics. Examples include: shallow water equations, isosteric/isentropic models, hydrostatic primitive equations, incompressible Boussinesq equations, anelastic systems, pseudo-incompressible equations, unified equations, and fully compressible Euler equations. Many of these equations can be written optionally in Eulerian or Lagrangian reference frame and in terms of various dependent variables; vorticity, velocity or momentum for dynamics, and total energy, internal energy or entropy for thermodynamics. However, with increasing computational power the non-hydrostatic (i.e., allscale) systems come into focus, thus reducing the plethora of options. 2016 Slide 3 ECMWF

Two reference frames Eulerian Lagrangian (the archetype problem, AP) The laws for fluid flow --- conservation of mass, Newtons 2nd law, conservation of energy, and 2nd principle of thermodynamics --- are independent on reference frames the two descriptions must be equivalent, somehow. 2016 Slide 4

ECMWF Fundamentals: Newton, 1642-1727 physics (re measurement) math (re Taylor series) physics (relating observations in the two reference frames) Taylor , 1685-1731 2016 Slide 5 ECMWF

More math: parcels volume evolution; 0 < J < , for the flow to be topologically realizable flow divergence, definition flow Jacobian Leibniz, 1646-1716 Euler expansion formula, and the rest is easy Euler, 1707-1783

2016 Slide 6 ECMWF mass continuity key tools for deriving conservation laws Lebesque1875-1941 2016 Slide 7 ECMWF Digression: an alternative semiheuristic, popular Eulerian argument in terms of flow through a fixed volume:

mass flux through the bounding surface instantaneous total mass Gauss, 1777-1855 2016 Slide 8 ECMWF Elementary examples: Shallow-water equations

(Szmelter & Smolarkiewicz, JCP, 2010) anelastic system See: Wedi & Smolarkiewicz, QJR, 2009, for discussion; and a special issue of JCP, 2008, Predicting Weather, Climate and Extreme Events for an overview of computational meteorology 2016 Slide 9 ECMWF All leading weather and climate codes are based on the compressible Euler equations, yet much of knowledge about non-hydrostatic atmospheric dynamics derives from the soundproof equations descendants of the classical, reduced incompressible Boussinesq equations Euler, 1707-1783

2016 Slide 10 Boussinesq, 1842-1929 ECMWF Why bother? Handling unresolved acoustic modes, while insisting on large time steps relative to speed of sound, makes numerics of non-hydrostatic atmospheric models based on the compressible Euler equations demanding pressure and density solid lines, entropy long dashes, velocity short dashes 2016 Slide 11 ECMWF

From compressible Euler equations to incompressible Boussinesq equations 3D momentum equations under gravity mass continuity and adiabatic entropy equations perturbation about static reference (base) state: momentum equation, perturbation form: Helmholtz, 1821-1894 for problems with small vertical scales and density perturbations: incompressible Boussinesq equations 2016 Slide 12

ECMWF Perturbation forms in the context of initial & boundary conditions Take incompressible Boussinesq equations: which also require initial conditions for pressure and density perturbations. Then consider an unperturbed ambient state, a particular solution to the same equations subtracting the latter from the former gives the form that takes homogeneous initial conditions for the perturbations about the environment ! 2016 Slide 13 ECMWF Perturbation forms in terms of potential temperature and Exner function

compressible Euler equations F.M. Exner, 1876-1930 2016 Slide 14 ECMWF The incompressible Boussinesq system is the simplest nonhydrostatic soundproof system. It describes small scale atmospheric dynamics of planetary boundary layers, flows past complex terrain and shallow gravity waves, thermal convection and fair weather clouds. Its extensions include the anelastic equations of Lipps & Hemler (1982, 1990) and the pseudo-incompressible equations of Durran (1989, 2008).

In the anelastic system the base state density is a function of altitude; in the pseudo-incompressible system the base state density is a (different) function of altitude, and the pressure gradient term is unabbreviated. In order to design a common approach for consistent integrations of soundproof and compressible nonhydrostatic PDEs for all-scale atmospheric dynamics, we manipulate the three governing systems into a single form convenient for discrete integrations: 2016 Slide 15 ECMWF Unified Framework, combined symbolic equations: (Smolarkiewicz, Khnlein & Wedi, JCP, 2014) gas law

conservation-law forms 2016 Slide 16 ECMWF Combined equations, conservation form: specific vs. density variables Accounting for curvilinear coordinates: Riemann, 1826-1863

2016 Slide 17 recall the archetype problem, AP ECMWF example integration schemes Global baroclinic instability (Smolarkiewicz, Khnlein & Wedi , JCP, 2014) 8 days, surface , 128x64x48 lon-lat grid, 128 PE of Power7 IBM CMP, 2880 dt=300 s, wallclock time=2.0 mns

PSI, 2880 dt=300 s, wallclock time=2.3 mns, ANL, 2880 dt=300 s, wallclock time=2.1 mns, 2016 Slide 18 ECMWF The role of baroclinicity anelastic 2016 Slide 19 pseudoincompressible

compressible ECMWF 1.5h, surface ln, 320x160 Gal-Chen grid, domain 120 km x 60 km ``soundproof dt=5 s ``acoustic dt=0.5 s 320 PE of Power7 IBM CPS PSI

ANL 2016 Slide 20 ECMWF Part II: Integration Schemes A) Forward-in-time (FT) non-oscillatory (NFT) integrators for all-scale flows, Cauchy, 1789-1857 2016 Slide 21

Kowalevski, 1850-1891 Lax, 1926- Wendroff, 1930- Robert, 1929-1993 ECMWF Generalised forward-in-time (FT) nonoscillatory (NFT) integrators for the AP EUlerian/LAGrangian congruence

Eulerian Lagrangian (semi) 2016 Slide 22 ECMWF Motivation for Lagrangian integrals 2016 Slide 23 ECMWF Motivation for Eulerian integrals

forward-in-time temporal discretization: Second order Taylor expansion about t=nt & Cauchy-Kowalewski procedure Compensating 1st error term on the rhs is a responsibility of an FT advection scheme (e.g. MPDATA). The 2nd error term depends on the implementation of an FT scheme 2016 Slide 24 ECMWF Given availability of a 2nd order FT algorithm for the homogeneous problem (R0), a 2nd order-accurate solution for an inhomogeneous problem with arbitrary R is: Banach principle, an important tool for systems with nonlinear right-hand-sides:

1892-1945 Eulerian semi-implicit compressible algorithms 2016 Slide 25 ECMWF Semi-implicit formulations (solar MHD example) thermodynamic/elliptic problems for pressures 2016 Slide 26 ECMWF in some detail for compressible Euler PDEs

of all-scale atmospheric dynamics semi-implicit ``acoustic scheme: (RE: Banach principle) (RE: thermodynamic pressure) 2016 Slide 27 ECMWF simple but computationally unaffordable; example 8 days, surface , 128x64x48 lon-lat grid,

128 PE of Power7 IBM CPI2, 2880 dt=300 s, wallclock time=2.0 mns CPEX, 432000 dt=2 s, wallclock time=178.9 mns This huge computational-efficiency gain comes at the cost of increased mathematical/numerical complexity 2016 Slide 28 ECMWF (RE: elliptic pressure)

elliptic boundary value problems (BVPs): Poisson problem in soundproof models relies on the mass continuity equation and diagonally preconditioned Poisson problem for pressure perturbation 2016 Slide 29 ECMWF Helmholtz problems for large-time-step compressible models also rely on mass continuity equation:

combine the evolutionary form of the gas law & mass continuity in the AP for pressure perturbation, to then derive the Helmholtz problem And how does one solve this thing ? 2016 Slide 30 ECMWF Part II: Integration Schemes B) Elliptic solvers for boundary value problems (BVP) in atmospheric models Poisson,

1781-1840 2016 Slide 31 Helmholtz, 1821-1894 Krylov, 1863-1945 Schur, 1875-1941 Richardson, 1881-1953

ECMWF Taxonomy: Direct methods (e.g., spectral, Gaussian elimination, conjugate-gradients (CG)) vs. Iterative methods (e.g., Gauss-Seidel, Richardson, multigrid, CG) Matrix inversion vs. matrix-free methods Approximate vs. Exact projection user-friendly libraries vs. bespoke solvers Multiple terminologies & classifications; common grounds; the state-of-the-art Basic tools and concepts: Banach principle, Neumann series, Gaussian elimination, Thomas 3-diagonal algorithm, Fourier transformation, calculus of variations, multigrid Physical analogies: Heat equation, damped oscillation equation, energy minimisation

2016 Slide 32 ECMWF Notion of variational Krylov-subspace solvers: i) basic concepts and definitions symbolism: linear BVP pseudo-time augmentation

where as is the solution error, so where is the domain integral gives the exact solution to provided

:= negative definiteness (comments on dissipativity, semi-definiteness and null spaces) 2016 Slide 33 ECMWF i) basic concepts and definitions, cnt. Next, the energy functional given := self-adjointness or symmetry in matrix representation,

a common property of Laplacian, since (comment on suitable boundary conditions) 2016 Slide 34 ECMWF i) basic concepts and definitions, cnt. 2 only for exact solution, otherwise it defines the residual error can be rewritten as

Richardson iteration ( ) (comments on stability vs convergence, and spectral implications) 2016 Slide 35 ECMWF Notion of variational Krylov-subspace solvers: ii) canonical schemes Steepest descent and minimum residual By the same arguments like applied to continuous equations , Richardson iteration implies and For self adjoint operators

implies Because the exact solution minimises the energy functional, one way to assure the optimal convergence is to minimise And from self adjointness Steepest descent 2016 Slide 36 ECMWF Digression, orthogonality of subsequent iterates

Minimum residual: Self adjointness can be difficult to achieve in practical models, the minimum residual circumvents this by minimising instead Steepest descent and minimum residual are important for understanding, but otherwise uncompetitive. The tyru foundation is provided by conjugate gradients and residuals

3 term recurrence formula a.k.a 2nd order Richardson, due to Frankel 1950 2016 Slide 37 ECMWF Conjugte gradient (CG) and conjugate residual (CR): the coefficients of which could be determined via norms minimisation, analogous to steepest descent and minimum residual, or instead CG 2016 Slide 38

CR ECMWF variational Krylov-subspace solvers: iii) operator preconditioning The best asymptotic convergence rate one can get from plain CG methods is in the inverse proportionality to (condition number )1/2 of the problem at hand P (left preconditioner) approximates L but is easier to invert. Preconditioned conjugate residual: replacing 2016 Slide 39

ECMWF Preconditioners, e P -1(r), examples: 1) 2) 3) 2016 Slide 40

ECMWF Non-symmetric preconditioned generalized conjugate residual scheme GCR(k): 2016 Slide 41 ECMWF A few remarks on boundary conditions: LH 2016 Slide 42 CE

ECMWF Principal conclusions Soundproof and compressible nonhydrostatic models form complementary elements of a general theoretical-numerical framework that underlies non-oscillatory forward-in-time (NFT) flow solvers The respective PDEs are integrated using essentially the same numerics The resulting flow solvers can be available in compatible Eulerian and semi-Lagrangian variants The flux-form flow solvers readily extend to unstructured-meshes P.K. Smolarkiewicz, C. Khnlein, N.P. Wedi, A consistent framework for discrete integrations of soundproof and compressible PDEs of atmospheric dynamics, J. Comput. Phys. 263 (2014) 185-205 P.K. Smolarkiewicz, W. Deconinck, M. Hamrud, G. Mozdzynski, C. Khnlein, J. Szmelter, N. Wedi, A finite-volume module for simulating global all-scale atmospheric flows, J. Comput. Phys. 314 (2016) 287-304. C. Khnlein, P.K. Smolarkiewicz, An unstructured-mesh finite-volume MPDATA for compressible atmospheric dynamics, J. Comput. Phys. 334 (2017) 16-3. P.K. Smolarkiewicz, L. Margolin, Variational methods for elliptic problems in fluid models, in: ECMWF Proceedings,

Workshop on Developments in Numerical Methods for Very High Resolution Global Models, 2000, pp.137159, Reading, UK. The research leading to these results has received funding from the European Research Council under the European Unions Seventh Framework Programme (FP7/2012/ERC Grant agreement no. 320375) 2016 Slide 43 ECMWF