Elsevier

Journal of Computational Physics

Volume 227, Issue 24, 20 December 2008, Pages 10226-10242
Journal of Computational Physics

A discontinuous Galerkin method for the shallow water equations in spherical triangular coordinates

https://doi.org/10.1016/j.jcp.2008.08.019Get rights and content

Abstract

A global model of the atmosphere is presented governed by the shallow water equations and discretized by a Runge–Kutta discontinuous Galerkin method on an unstructured triangular grid. The shallow water equations on the sphere, a two-dimensional surface in R3, are locally represented in terms of spherical triangular coordinates, the appropriate local coordinate mappings on triangles. On every triangular grid element, this leads to a two-dimensional representation of tangential momentum and therefore only two discrete momentum equations.

The discontinuous Galerkin method consists of an integral formulation which requires both area (elements) and line (element faces) integrals. Here, we use a Rusanov numerical flux to resolve the discontinuous fluxes at the element faces. A strong stability-preserving third-order Runge–Kutta method is applied for the time discretization. The polynomial space of order k on each curved triangle of the grid is characterized by a Lagrange basis and requires high-order quadature rules for the integration over elements and element faces. For the presented method no mass matrix inversion is necessary, except in a preprocessing step.

The validation of the atmospheric model has been done considering standard tests from Williamson et al. [D.L. Williamson, J.B. Drake, J.J. Hack, R. Jakob, P.N. Swarztrauber, A standard test set for numerical approximations to the shallow water equations in spherical geometry, J. Comput. Phys. 102 (1992) 211–224], unsteady analytical solutions of the nonlinear shallow water equations and a barotropic instability caused by an initial perturbation of a jet stream. A convergence rate of O(Δxk+1) was observed in the model experiments. Furthermore, a numerical experiment is presented, for which the third-order time-integration method limits the model error. Thus, the time step Δt is restricted by both the CFL-condition and accuracy demands. Conservation of mass was shown up to machine precision and energy conservation converges for both increasing grid resolution and increasing polynomial order k.

Introduction

Modeling atmospheric flows for climate simulations as well as for weather prediction is a complex problem, due to the nonlinear structure of the dynamical and physical phenomena on widely varying spatial and temporal scales and their multi-scale interaction processes. Depending on the complexity of an atmospheric model the governing equations are the fundamental atmospheric conservation laws for mass, momentum, and energy or appropriate simplifications of them. If the regarded equation set is a hyperbolic system, energetic shocks can develop theoretically. Although this is usually not the case in atmospheric models, the discretization should represent regions of scale collapse and breaking waves generating discontinuities in the velocity field; the discrete conservation properties of the discontinuous Galerkin (DG) method are appropriate for this task.

The shallow water equations (SWE), valid for a homogeneous atmosphere with small vertical velocities and horizontal velocities independent in the vertical direction, constitute a hyperbolic system of conservation laws. It is one of the simplest nonlinear hyperbolic systems, covering important planetary atmospheric features, like the Rossby wave formation.

For the spherical SWE the spatial domain is the sphere S, a two-dimensional surface in R3. In a regional or mesoscale SWE model the momentum is a two-dimensional vector. In contrast, the Cartesian formulation of the spherical case in [43], [8], [44] represents the tangential momentum of the flow as a three-dimensional vector and includes a Lagrangian multiplier to constrain the momentum to be tangential. Applying this form to a numerical model usually leads to three momentum equations and requires a correction step to satisfy the constraint discretely, see e.g. [17]. Models in standard spherical coordinates satisfy a two-dimensional momentum representation but have to pay additional attention to phenomena near the poles due to singularities of the coordinate mapping, see e.g. [27].

The idea to avoid the drawbacks of the Cartesian and the spherical coordinates formulation is to represent the spherical SWE in terms of local coordinate transformations. On a cubed-sphere grid, a spherical quadrilateral grid, the following works [37], [33], [31], [45], [35] achieved a two-dimensional momentum representation avoiding any pole problem. Our new model achieves the same flexibility but on unstructured triangular grids using spherical triangular coordinates.

Numerous models on spherical triangles have been proposed in the last three decades. For example, the early work by Sadourny at al. [38] and Williamson [49] introduced to the atmospheric community the use of triangular grids based on the icosahedron to develop the underlying grid for the construction of finite difference operators. Work on triangular grids based on the icosahedron for discretizing the sphere lay dormant for another 15 years until the work was resumed by Baumgardner and Frederickson [1]. Ten years later work on these grids was resurrected by Heikes and Randall [20] and followed by Giraldo [13], Thuburn [47], Stuhne and Peltier [42], Tomita et al. [48], and Heinze and Hense [21]. All of these models rely directly on either the triangular grid being derived from the icosahedron or on a linear representation of the discrete operators; this way there is an easily computable dual grid which is based on hexagons (icosahedron) or the operators can be constructed using the vertices of an element which are co-planar (linear representation). However, for high-order operators on general triangular grids, one needs to construct the discrete spatial operators directly on the curved manifold which then requires the derivation of the Christoffel symbols from differential geometry (see [35] for a summary of the use of differential geometry for atmospheric flow).

One of the contributions of this manuscript is to show how to use these ideas for constructing high-order spatial operators on general triangulations on the sphere. While these ideas have been used extensively for quadrilateral-based grids, they have not been used at all (or not often); to our knowledge, there are currently no existing models on the sphere which use high-order discretizations on the triangle with the exception being those developed by the present authors. Läuter et al. [25] developed a shallow water model on the sphere using second order finite elements on dynamically adaptive triangular grids, while Giraldo and Warburton [17] and Giraldo [15] developed shallow water models on the sphere using up to 15 order spectral element and discontinuous Galerkin operators on unstructured triangular grids. Extracting the highlights of these three separate works results in a robust, accurate, and efficient model. However, in order to achieve this aim requires writing the equations directly on the parametric space defined by the curved spherical triangles. We consider the spherical SWE in flux form on the surface S using coordinate independent differential operators on S known from differential geometry. By using an appropriate local coordinate mapping γE on the curved spherical triangle E we can then obtain a two-dimensional representation of the tangential momentum vectors. This then allows us full grid independence such that our discrete operators can be constructed to arbitrarily high-order while doing so on generalized unstructured triangulations on the sphere (or any other curved manifold, for that matter). This approach allows us the same flexibility enjoyed with the Cartesian methods discussed previously while now only requiring tangential momentum equations.

Numerous numerical methods have been proposed for next generation global atmospheric models including finite volumes [27], [34], spectral elements [45], [11], [9], [17], and DG [16], [30], [15] methods. We have selected the DG method for our model because it allows us to achieve high-order accuracy as in spectral elements while conserving all quantities both locally and globally as in finite volumes, see the review in [5]. Furthermore, our use of unstructured triangular grids allows for much flexibility in future work on adaptivity.

The organization of this article is as follows. In Section 2 the governing spherical SWE are given using surface differential operators. Section 3 describes the numerical discretization by a Runge–Kutta discontinuous Galerkin method applying spherical triangular coordinates. In Section 4 the atmospheric model based on the discretization is validated in terms of standard tests from Williamson et al. [50], steady-state and unsteady analytical solutions and a barotropic instability generated by a small initial perturbation.

Section snippets

Spherical shallow water equations

The spherical SWE are a system of conservation laws for the geopotential layer depth (mass) and the flow momentum. Because the integration domain of the SWE is the sphere, a two-dimensional surface in R3, the system can be formulated in the surrounding Cartesian space R3, see Côté [8]. Côté’s formulation is equivalent to a conservative form of the SWE on the surface S, which is the formulation used further below.

Let us consider the sphere S={xR3||x|=a} with the Earth’s radius a=6.371×106m, the

Discontinuous Galerkin method

The DG method is applied to the conservative form (1) of the SWE on the surface S. On each curved triangle (element) E of the grid tesselation spherical triangular coordinates are introduced, which are local coordinate mappings γE on E. The polynomial representation on each grid element uses high-order Lagrange polynomials based on specially chosen Lagrange points (see Section 3.2). This approach leads to the local representation of the tangential momentum fields by two components only. An

Numerical results

The RK-DG method described in Section 3 has been used to implement an atmospheric model which has been validated performing numerical experiments. The validation process is carried out in three steps. At first, standard model tests (flow over an isolated mountain and a Rossby–Haurwitz wave) from Williamson [50] are performed. After that, a convergence study considering steady-state and unsteady analytical solutions of the nonlinear SWE is done. Finally, a barotropic instability in a localized

Summary and outlook

A global model of the atmosphere has been developed based on the spherical shallow water equations and discretized by a Runge–Kutta discontinuous Galerkin method. The equations are formulated as a hyperbolic conservation law on the sphere, a two-dimensional surface in R3, using coordinate independent differential operators.

Spherical triangular coordinates, which are local coordinate mappings γE, and high-order polynomial spaces are defined on each curved element ES of a given spherical

Acknowledgements

The authors gratefully acknowledge the support of the Office of Naval Research Global through the Grant N00014-07-1-4038. The second author also gratefully acknowledges the support of the Office of Naval Research through PE-0602435N.

References (50)

  • M. Läuter et al.

    Unsteady analytical solutions of the spherical shallow water equations

    J. Comput. Phys.

    (2005)
  • M. Läuter et al.

    A parallel adaptive barotropic model of the atmosphere

    J. Comput. Phys.

    (2007)
  • C. Ronchi et al.

    The cubed sphere: a new method for the solution of partial differential equations in spherical geometry

    J. Comput. Phys.

    (1996)
  • J.A. Rossmanith

    A wave propagation method for hyperbolic systems on the sphere

    J. Comput. Phys.

    (2006)
  • J.A. Rossmanith et al.

    A wave propagation algorithm for hyperbolic systems on curved manifolds

    J. Comput. Phys.

    (2004)
  • C.-W. Shu et al.

    Efficient implementation of essentially non-oscillatory shock-capturing schemes

    J. Comput. Phys.

    (1988)
  • G.R. Stuhne et al.

    New icosahedral grid-point discretizations of the shallow water equations on the sphere

    J. Comput. Phys.

    (1999)
  • M. Taylor et al.

    The spectral element method for the shallow water equations on the sphere

    J. Comput. Phys.

    (1997)
  • H. Tomita et al.

    Shallow water model on a modified icosahedral geodesic grid by using spring dynamics

    J. Comput. Phys.

    (2001)
  • D.L. Williamson et al.

    A standard test set for numerical approximations to the shallow water equations in spherical geometry

    J. Comput. Phys.

    (1992)
  • J.R. Baumgardner et al.

    Icosahedral discretization of the two-sphere

    SIAM J. Numer. Anal.

    (1985)
  • B. Cockburn et al.

    The Runge–Kutta local projection p1-discontinuous-Galerkin finite element method for scalar conservation laws

    Math. Mod. Num. Anal.

    (1991)
  • B. Cockburn et al.

    Runge–Kutta discontinuous Galerkin methods for convection-dominated problems

    J. Sci. Comput.

    (2001)
  • J. Côté

    A Lagrange multiplier approach for the metric terms of semi-Lagrangian models on the sphere

    Q.J.R. Meteorol. Soc.

    (1988)
  • J. Dennis et al.

    High-resolution mesh convergence porperties and parallel efficiency of a spectral element atmospheric dynamical core

    Int. J. High Perfor. Comput. Appl.

    (2005)
  • Cited by (47)

    • Dispersion analysis of the P<inf>n</inf>−P<inf>n−1</inf><sup>DG</sup> mixed finite element pair for atmospheric modelling

      2018, Journal of Computational Physics
      Citation Excerpt :

      However as noted by Adcroft et al. [1], Le Roux et al. [28] and others, the C-grid supports a stationary Coriolis mode due to the spatial averaging necessary to reconstruct the tangential velocity components which can become problematic for poorly resolved Rossby radius. In recent years high-order methods such as spectral elements [21,38], and discontinuous Galerkin (DG) based methods [18,20,32] have increased in popularity, due in part to the high level of computational efficiency obtainable on massively parallel supercomputers [32] as well as providing a high order of accuracy on non-uniform, non-orthogonal grids. This gives a greater flexibility in the choice of grids, thus allowing a more uniform tiling of the sphere or adaptive/variable resolution grids [24,7].

    View all citing articles on Scopus

    This work relates to Department of the Navy Grant N00014-07-1-4038 issued by the Office of Naval Research Global. The United States Government has a royalty-free licence throughout the world in all copyrightable material contained herein.

    View full text