An accurate and efficient numerical framework for adaptive numerical weather prediction

An accurate and efficient numerical framework
for adaptive numerical weather prediction

Giovanni Tumolo,  Luca Bonaventura

We present an accurate and efficient discretization approach for the adaptive discretization of typical model equations employed in numerical weather prediction. A semi-Lagrangian approach is combined with the TR-BDF2 semi-implicit time discretization method and with a spatial discretization based on adaptive discontinuous finite elements. The resulting method has full second order accuracy in time and can employ polynomial bases of arbitrarily high degree in space, is unconditionally stable and can effectively adapt the number of degrees of freedom employed in each element, in order to balance accuracy and computational cost. The adaptivity approach employed does not require remeshing, therefore it is especially suitable for applications, such as numerical weather prediction, in which a large number of physical quantities are associated with a given mesh. Furthermore, although the proposed method can be implemented on arbitrary unstructured and nonconforming meshes, even its application on simple Cartesian meshes in spherical coordinates can cure effectively the pole problem by reducing the polynomial degree used in the polar elements. Numerical simulations of classical benchmarks for the shallow water and for the fully compressible Euler equations validate the method and demonstrate its capability to achieve accurate results also at large Courant numbers, with time steps up to 100 times larger than those of typical explicit discretizations of the same problems, while reducing the computational cost thanks to the adaptivity algorithm.

Earth System Physics Section

The Abdus Salam International Center for Theoretical Physics

Strada Costiera 11, 34151 Trieste, Italy

MOX – Modelling and Scientific Computing,

Dipartimento di Matematica “F. Brioschi”, Politecnico di Milano

Via Bonardi 9, 20133 Milano, Italy

Keywords: Discontinuous Galerkin methods, adaptive finite elements, semi-implicit discretizations, semi-Lagrangian discretizations, shallow water equations, Euler equations.

AMS Subject Classification: 35L02, 65M60, 65M25, 76U05, 86A10

1 Introduction

The Discontinuous Galerkin (DG) spatial discretization approach is currently being employed by an increasing number of environmental fluid dynamics models see e.g. [12], [17], [28],[39],[19],[25] and a more complete overview in [5]. This is motivated by the many attractive features of DG discretizations, such as high order accuracy, local mass conservation and ease of massively parallel implementation.

On the other hand, DG methods imply severe stability restrictions when coupled with explicit time discretizations. One traditional approach to overcome stability restrictions in low Mach number problems is the combination of semi - implicit (SI) and semi - Lagrangian (SL) techniques. In a series of papers [41], [42], [19], [14], [52] it has been shown that most of the computational gains traditionally achieved in finite difference models by the application of SI, SL and SISL discretization methods are also attainable in the framework of DG approaches. In particular, in [52] we have introduced a dynamically adaptive SISL-DG discretization approach for low Mach number problems, that is quite effective in achieving high order spatial accuracy, while reducing substantially the computational cost.

In this paper, we apply the technique of [52] to the shallow water equations in spherical geometry and to the the fully compressible Euler equations, in order to show its effectiveness for model problems typical of global and regional weather forecasting. The advective form of the equations of motion is employed and the semi-implicit time discretization is based on the TR-BDF2 method, see e.g. [21], [32]. This combination of two robust ODE solvers yields a second order accurate, A-stable and L-stable method (see e.g. [27]), that is effective in damping selectively high frequency modes. At the same time, it achieves full second order accuracy, while the off-centering in the trapezoidal rule, typically necessary for realistic applications to nonlinear problems (see e.g. [7], [11], [52]), limits the accuracy in time to first order. Numerical results presented in this paper show that the total computational cost of one TR-BDF2 step is analogous to that of one step of the off-centered trapezoidal rule, as well as the structure of the linear problems to be solved at each time step, thus allowing to extend naturally to this more accurate method any implementation based on the off-centered trapezoidal rule. Numerical simulations of the shallow water benchmarks proposed in [55], [29], [22] and of the non-hydrostatic benchmarks proposed in [47], [23] have been employed to validate the method and to demonstrate its capabilities. In particular, it will be shown that the present approach enables the use of time steps even 100 times larger than those allowed for DG models by standard explicit schemes, see e.g. the results in [38].

The method presented in this paper, just as its previous version in [52], can be applied in principle on arbitrarily unstructured and even nonconforming meshes. For example, a model based on this method could run on a non conforming mesh of rectangular elements built around the nodes of a reduced Gaussian grid [20]. For simplicity, however, no such implementation has been developed so far. Here, only a simple Cartesian mesh has been used. If no degree adaptivity is employed, this results in very high Courant numbers in the polar regions. These do not result in any special stability problems for the present SISL discretization approach, as it will be shown by the numerical results reported below. On the other hand, even with an implementation based on a simple Cartesian mesh in spherical coordinates, the flexibility of the DG space discretization allows to reduce the degree of the basis and test functions employed close to the poles, thus making the effective model resolution more uniform and solving the efficiency issues related to the pole problem by static adaptivity. This is especially advantageous because the conditioning of the linear system to be solved at each time step is greatly improved and, as a consequence, the number of iterations necessary for the linear solver is reduced by approximately , while at the same time no spurious reflections nor artificial error increases are observed.

Beyond these computational advantages, we believe that the present approach based on adaptivity is especially suitable for applications to numerical weather prediction, in contrast to adaptivity approaches (that is, local mesh coarsening or refinement in which the size of some elements changes in time). Indeed, in numerical weather prediction, information that is necessary to carry out realistic simulations (such as orography profiles, data on land use and soil type, land-sea masks) needs to be reconstructed on the computational mesh and has to be re-interpolated each time that the mesh is changed. Furthermore, many physical parameterizations are highly sensitive to the mesh size. Although devising better parameterizations that require less mesh-dependent tuning is an important research goal, more conventional parameterizations will still be in use for quite some time. As a consequence, it is useful to improve the accuracy locally by adding supplementary degrees of freedom where necessary, as done in a adaptive framework, without having to change the underlying computational mesh. In conclusion, the resulting modeling framework seems to be able to combine the efficiency and high order accuracy of traditional SISL pseudo-spectral methods with the locality and flexibility of more standard DG approaches.

In section 2, two examples of governing equations are introduced. In section 3, the TR-BDF2 method is reviewed. In section 4 the approach employed for the advection of vector fields in spherical geometry is described in detail. In section 5, we introduce the SISL-DG discretization approach for the shallow water equations in spherical geometry. In section 6, we outline its extension to the fully compressible Euler equations in a vertical plane. Numerical results are presented in section 7, while in section 8 we try to draw some conclusions and outline the path towards application of the concepts introduced here in the context of a non hydrostatic dynamical core.

2 Governing equations

We consider as a basic model problem the two-dimensional shallow water equations on a rotating sphere (see e.g. [15]). These equations are a standard test bed for numerical methods to be applied to the full equations of motion of atmospheric or oceanic circulation models, see e.g. [55]. Among their possible solutions, they admit Rossby and inertial gravity waves, as well as the response of such waves to orographic forcing. We will use the advective, vector form of the shallow water equations:


Here represents the fluid depth, the bathymetry elevation, the Coriolis parameter, the unit vector locally normal to the Earth’s surface and the gravity force per unit mass on the Earth’s surface. Assuming that are orthogonal curvilinear coordinates on the sphere (or on a portion of it), we denote by and the components of the (diagonal) metric tensor. Furthermore, we set where and are the contravariant components of the velocity vector in the coordinate direction and respectively, multiplied by the corresponding metric tensor components. We also denote by the Lagrangian derivative

so that . In particular, in this paper standard spherical coordinates will be employed.

As an example of a more complete model, we will also consider the fully compressible, non hydrostatic equations of motion. Following e.g. [10], [4], [11], they can be written as

where, being a reference pressure value, is the potential temperature, is the Exner pressure, while are the constant pressure and constant volume specific heats and the gas constant of dry air respectively. Here the Coriolis force is omitted for simplicity. Notice also that, by a slight abuse of notation, in the three-dimensional case denotes the three dimensional velocity field and the operators are also three-dimensional, while we will assume in the description of two dimensional, vertical slice models. It is customary to rewrite such equations in terms of perturbations with respect to a steady hydrostatic reference profile, so that assuming with one obtains for a vertical plane


It can be observed that equations (3)-(6) are isomorphic to equations (1)-(2), which will allow to extend almost automatically the discretization approach proposed for the former to the more general model.

3 Review of the TR-BDF2 method

We review here some properties of the so called TR-BDF2 method, which was first introduced in [2]. Given a Cauchy problem


and considering a time discretization employing a constant time step the TR-BDF2 method is defined by the two following implicit stages:


Here is an implicitness parameter and

It is immediate that the first stage of (3) is simply the application of the trapezoidal rule (or Crank-Nicolson method) over the interval It could also be substituted by an off centered Crank-Nicolson step without reducing the overall accuracy of the method. The outcome of this stage is then used to turn the two step BDF2 method into a single step, two stages method. This combination of two robust stiff solvers yields a method with several interesting accuracy and stability properties, that were analyzed in detail in [21]. As shown in this paper, this analysis is most easily carried out by rewriting the method as


In this formulation, the TR-BDF2 method is clearly a Singly Diagonal Implicit Runge Kutta (SDIRK) method, so that one can rely on the theory for this class of methods to derive stability and accuracy results (see e.g. [27]). Notice that the same method has been rediscovered in [6] and has been analyzed and applied also in [18], to treat the implicit terms in the framework of an Additive Runge Kutta approach (see e.g. [26]). As shown in [21], the TR-BDF2 method is second order accurate and A-stable for any value of Written as in (9), the method can also be proven to constitute a (2,3) embedded Runge-Kutta pair, with companion coefficients given by

provided that no off centering is employed in the first stage of (3). This equips the method with an extremely efficient estimator of the time discretization error. Furthermore, for it is also L-stable. Therefore, with this coefficient value it can be safely applied to problems with eigenvalues whose imaginary part is large, such as typically arise from the discretization of hyperbolic problems. This is not the case for the standard trapezoidal rule (or Crank-Nicolson) implicit method, whose linear stability region is exactly bounded by the imaginary axis. As a consequence, it is common to apply the trapezoidal rule with off centering, see e.g. [7], [11] as well as [52], which results in a first order time discretization. TR-BDF2 appears therefore to be an interesting one step alternative to maintain full second order accuracy, especially considering that, if formulated as (3), it is equivalent to performing two Crank-Nicolson steps with slightly modified coefficients. In order to highlight the advantages of the proposed method in terms of accuracy with respect to other common robust stiff solvers, we plot in figure 1 the contour levels of the absolute value of the linear stability function of the TR-BDF2 method without off centering in the first stage, compared to the analogous contours of the off centered Crank-Nicolson method with averaging parameter in figures 2, 3, respectively, and to those of the BDF2 method in figure 4. It is immediate to see that TR-BDF2 introduces less damping around the imaginary axis for moderate values of the time step. On the other hand, TR-BDF2 is more selective in damping very large eigenvalues, as clearly displayed in figure 5, where the absolute values of the linear stability functions of the same methods (with the exception of BDF2, for which an explicit representation of the stability function is not available) are plotted along the imaginary axis.

Figure 1: Contour levels of the absolute value of the stability function of the TR-BDF2 method without off centering in the first stage. Contour spacing is from to
Figure 2: Contour levels of the absolute value of the stability function of the off centered Crank-Nicolson method with averaging parameter (equivalent to an off centering parameter valued ). Contour spacing is from to
Figure 3: Contour levels of the absolute value of the stability function of the off centered Crank-Nicolson method with averaging parameter (equivalent to an off centering parameter valued ). Contour spacing is from to
Figure 4: Contour levels of the absolute value of the stability function of the BDF2 method. Contour spacing is from to
Figure 5: Graph of the absolute value of the stability functions of several L-stable methods along the imaginary axis.

4 Review of Semi-Lagrangian evolution operators for vector fields on the sphere

The semi-Lagrangian method can be described introducing the concept of evolution operator, along the lines of [37, 35]. Indeed, let denote a generic function of space and time that is the solution of

To approximate this solution on the time interval a numerical evolution operator is introduced, that approximates the exact evolution operator associated to the frozen velocity field that may coincide with the velocity field at time level or with an extrapolation derived from more previous time levels. More precisely, if denotes the solution of


with initial datum at time , then the expression denotes a numerical approximation of where and the notation is used. Since is nothing but the position at time of the fluid parcel reaching location at time , according to standard terminology, it is called the departure point associated with the arrival point . Different methods can be employed to approximate ; in this paper, for simplicity, the method proposed in [34] has been employed in spherical geometry. Furthermore, to guarantee an accuracy compatible with that of the semi-implicit time-discretization, an extrapolation of the velocity field at the intermediate time level was used as in (10). On the other hand, in the application to Cartesian geometry (for the vertical slice discretization), a simple first order Euler method with sub-stepping was employed, see e.g. [16], [45].

In case of the advection of a vector field

as in the momentum equation (2), the extension of this approach has to take into account the curvature of the spherical manifold. More specifically, unit basis vectors at the departure point are not in general aligned with those at the arrival point, i.e., if represent a unit vector triad, in general

To deal with this issue two approaches are available. The first, intrinsically Eulerian, consists in the introduction of the Christoffel symbols in the covariant derivatives definition, giving rise to the well known metric terms, before the SISL discretization, and then in the approximation along the trajectories of those metric terms. This approach has been shown to be source of instabilities in a semi-Lagrangian frame, see e.g. [44, 8, 9, 13] and therefore is not adopted in this work. The second approach, more suitable for semi-Lagrangian discretizations, takes into account the curvature of the manifold only at discrete level, i.e. after the SISL discretization has been performed. Many variations of this idea have been proposed, see e.g. [44, 8, 9, 3, 49]. In [48], they have all been derived in a unified way by the introduction of a proper rotation matrix that transforms vector components in the departure-point unit vector triad into vector components in the arrival-point unit vector triad . To see how this rotation matrix comes into play, it is sufficient to consider the action of the evolution operator on a given vector valued function of space and time , defined as an approximation of


and to write this equation componentwise. is known through its components in the departure point unit vector triad:


Therefore, via (11), the components of in the unit vector triad at the same point are given by projection of (12) along :

i.e., in matrix notation:

Under the shallow atmosphere approximation [51], can be reduced to the rotation matrix


where, as shown in [48], Therefore, in the following the evolution operator for vector fields will be defined componentwise as


5 A novel SISL time integration approach for the shallow water equations on the sphere

The SISL discretization of equations.(1)-(2) based on (3) is then obtained by performing the two stages in (3) after reinterpretation of the intermediate values in a semi-Lagrangian fashion. Furthermore, in order to avoid the solution of a nonlinear system, the dependency on in is linearized in time, as common in semi-implicit discretizations based on the trapezoidal rule, see e.g. [7],[52]. Numerical experiments reported in the following show that this does not prevent to achieve second order accuracy in the regimes of interest for numerical weather prediction. The TR stage of the SISL time semi-discretization of the equations in vector form (1)-(2) is given by


The TR stage is then followed by the BDF2 stage:


For each of the two stages, the spatial discretization can be performed along the lines described in [52], allowing for variable polynomial order to locally represent the solution in each element. The spatial discretization approach considered is independent of the nature of the mesh and could also be implemented for fully unstructured and even non conforming meshes. For simplicity, however, in this paper only an implementation on a structured mesh in longitude-latitude coordinates has been developed. In principle, either Lagrangian or hierarchical Legendre bases could be employed. We will work almost exclusively with hierarchical bases, because they provide a natural environment for the implementation of a adaptation algorithm, see for example [56]. A central issue in finite element formulations for fluid problems is the choice of appropriate approximation spaces for the velocity and pressure variables (in the context of SWE, the role of the pressure is played by the free surface elevation). An inconsistent choice of the two approximation spaces indeed may result in a solution that is polluted by spurious modes, for the specific case of SWE see for example [31, 53, 54] as well as the more recent and comprehensive analysis in [30]. Here, we have not investigated this issue in depth, but the model implementation allows for approximations of higher polynomial degree for the velocity fields than for the height field. Even though no systematic study was performed, no significant differences were noticed between results obtained with equal or unequal degrees. In the following, only results with unequal degrees are reported, with the exception of an empirical convergence test for a steady geostrophic flow.

All the integrals appearing in the elemental equations are evaluated by means of Gaussian numerical quadrature formulae, with a number of quadrature nodes consistent with the local polynomial degree being used. In particular, notice that integrals of terms in the image of the evolution operator i.e. of functions evaluated at the departure points of the trajectories arriving at the quadrature nodes, cannot be computed exactly (see e.g. [36, 40]), since such functions are not polynomials. Therefore a sufficiently accurate approximation of these integrals is needed, which may entail the need to employ numerical quadrature formulae with more nodes than the minimal requirement implied by the local polynomial degree. This overhead is actually compensated by the fact that, for each Gauss node, the computation of the departure point is only to be executed once for all the quantities to be interpolated.

After spatial discretization has been performed, the discrete degrees of freedom representing velocity unknowns can be replaced in the respective discrete height equations, yielding in each case a linear system whose structure is entirely analogous to that obtained in [52]. The non-symmetric linear systems obtained from the TR-BDF2 stages are solved in our implementation by the GMRES method [46]. A classical stopping criterion based on a relative error tolerance of was employed (see e.g. [24]). For the GMRES solver, so far, only a block diagonal preconditioning was employed. As it will be shown in section 7, the condition number of the systems to be solved can be greatly reduced if lower degree elements are employed close to the poles. In any case, the total computational cost of one TR-BDF2 step is entirely analogous to that of one step of the standard off centered trapezoidal rule employed in [52], since the structure of the systems is the same but for each stage only a fraction of the time step is being computed. Once has been computed by solving this linear system, then can be recovered by back substituting into the momentum equation.

6 Extension of the time integration approach to the Euler equations

In this section, we show that the previously proposed method can be extended seamlessly to the fully compressible Euler equations as formulated in equations (3) - (6). For simplicity, only the application to the two dimensional vertical slice case is presented, but the extension to three dimensions is straightforward. Again, in order to avoid the solution of a nonlinear system, the dependency on in and the dependency on in are linearized in time, as common in semi-implicit discretizations based on the trapezoidal rule, see e.g. [10], [4].

The semi-Lagrangian counterpart of the TR substep of (3) is first applied to to (3) - (6), so as to obtain:


Following [10] the time semi-discrete energy equation (22) can be inserted into the time semi-discrete vertical momentum equation (6), in order to decouple the momentum and the energy equations as follows


Equations (6), (6) and (6) are a set of three equations in three unknowns only, namely and that can be compared with equations (15), (5) with and (Cartesian geometry). From the comparison it is clear that the two formulations are isomorphic under correspondence

We can then consider the semi-Lagrangian counterpart of the BDF2 substep of (3) applied to (3) - (6) to obtain:


Again, following [10], the time semi-discrete energy equation (27) can be inserted into the time semi-discrete vertical momentum equation (26), in order to decouple the momentum and the energy equations:


Now equations (24), (25) and (6) are a set of three equations in three unknowns only, namely and that can be compared with equations (17), (5) with and (Cartesian geometry). Again, it is easy to see that also in this case exactly the same structure results as in equations (17)-(5) with the correspondence , so that the approach (and code) proposed for the shallow water equations can be extended to the fully compressible Euler equation in a straightforward way.

7 Numerical experiments

The numerical method introduced in section 5 has been implemented and tested on a number of relevant test cases using different initial conditions and bathymetry profiles, in order to assess its accuracy and stability properties and to analyze the impact of the adaptivity strategy. Whenever a reference solution was available, the relative errors were computed in the and norms at the final time of the simulation according to [55] as:


where denotes the reference solution for a model variable and is a discrete approximation of the global integral


computed by an appropriate numerical quadrature rule, consistent with the numerical approximation being tested, and the maximum is computed over all nodal values.

The test cases considered for the shallow water equations in spherical geometry are

  • a steady-state geostrophic flow: in particular, we have analyzed results in test case 2 of [55] in the configuration least favorable for methods employing longitude-latitude meshes;

  • the unsteady flow with exact analytical solution described in [29];

  • the polar rotating low-high, introduced in [33], aimed at showing that no problems arise even in the case of strong cross polar flows;

  • zonal flow over an isolated mountain and Rossby-Haurwitz wave of wavenumber 4, corresponding respectively to test cases 5 and 6 in [55].

For the first two tests, analytic solutions are available and empirical convergence tests can be performed. The test cases considered for the discretization of equations (3)-(6) are

  • inertia gravity waves involving the evolution of a potential temperature perturbation in a channel with periodic boundary conditions and uniformly stratified environment with constant Brunt-Wäisälä frequency, as described in [47];

  • a rising thermal bubble given by the evolution of a warm bubble in a constant potential temperature environment, as described in [23].

In all the numerical experiments performed for this paper, neither spectral filtering nor explicit diffusion of any kind were employed, the only numerical diffusion being implicit in the time discretization approach. We have not yet investigated to which extent the quality of the solutions is affected by this choice, but this should be taken into account when comparing quantitatively the results of the present method to those of reference models, such as the one described in [22], in which explicit numerical diffusion is added. Sensitivity of the comparison results to the amount of numerical diffusion has been highlighted in several model validation exercises, see e.g. [43].

Since semi-implicit, semi-Lagrangian methods are most efficient for low Froude number flows, where the typical velocity is much smaller than that of the fastest propagating waves, all the tests considered fall in this hydrodynamical regime. Therefore, in order to assess the method efficiency, a distinction has been made between the maximum Courant number based on the velocity, on one hand, and, on the other hand, the maximum Courant number based on the celerity, or the maximum Courant number based on the sound speed, defined respectively as

where is to be interpreted as generic value of the meshsize in either coordinate direction. For the tests in which adaptivity was employed, if denotes the local polynomial degree used at timestep to represent a model variable inside the element of the mesh, while is the maximum local polynomial degree considered, the efficiency of the method in reducing the computational effort has been measured by monitoring the evolution of the quantities

where is the total number of elements, denotes the total number of GMRES iterations at time step for the adapted local degrees configuration and the total number of GMRES iterations at time step for the configuration with maximum degree in all elements, respectively. Average values of these indicators over the simulations performed are reported in the following, denoted by and respectively. The error between the adaptive solution and the corresponding one obtained with uniform maximum polynomial degree everywhere has been measured in terms of (29). Finally, in some cases conservation of global invariants has been monitored by evaluating at each time step the following global integral quantities:


where has been defined in (32) and is the density associated to each global invariant. According to the choice of , following invariants are considered: mass, i.e. , total energy, i.e. ), and potential enstrophy, i.e.

7.1 Steady-state geostrophic flow

We first consider the test case 2 of [55], where the solution is a steady state flow with velocity field corresponding to a zonal solid body rotation and field obtained from the velocity ones through geostrophic balance. All the parameter values are taken as in [55]. The flow orientation parameter has been chosen here as making the test more challenging on a longitute-latitude mesh. Error norms associated to the solution obtained on a mesh of elements for different polynomial degrees are shown in tables 1, 2 and 3 for and respectively. All the results have been computed at days at fixed maximum Courant numbers so that different values of have been employed for different polynomial order. We remark that the resulting time steps are significantly larger than those allowed by typical explicit time discretizations for analogous DG space discretizations, see e.g. the results in [38]. The spectral decay in the error norms can be clearly observed, until the time error becomes dominant. For better comparison with the results in [38], we consider again the configuration with on elements, which corresponds to the same resolution in space as for the grid used in [38]. While s is used in [38] giving a the proposed SISLDG formulation can be run with s, in which case and the average number of iterations required by the linear solver is 1 for the TR substep and 4 for the BDF2 substep.

Table 1: Relative errors on for different polynomial degrees, SWE test case 2 with at time days.
Table 2: Relative errors on for different polynomial degrees, SWE test case 2 with at time days.
Table 3: Relative errors on for different polynomial degrees, SWE test case 2 with at time days.

Another convergence test was performed for increasing the number of elements and correspondingly decreasing the value of the time step. In this case, the maximum Courant numbers vary because of the mesh inhomogeneity, so that The results are reported in tables 4, 5 and 6 for and respectively. The empirical convergence order based on the norm errors has also been estimated, showing that in this stationary test convergence rates above the second order of the time discretization can be achieved.

Table 4: Relative errors on for different number of elements, SWE test case 2 with at time days.
Table 5: Relative errors on for different number of elements, SWE test case 2 with at time days.
Table 6: Relative errors on for different number of elements, SWE test case 2 with at time days.

7.2 Unsteady flow with analytic solution

In a second, time dependent test, the analytic solution of (1)-(2) derived in [29] has been employed to assess the performance of the proposed discretization. More specifically, the analytic solution defined in formula (23) of [29] was used. Since the exact solution is periodic, the initial profiles also correspond to the exact solution an integer number of days later. The proposed SISLDG scheme has been integrated up to days with and on meshes with increasing number of elements, while the time step has been decreased accordingly. In this case, the maximum Courant numbers vary because of the mesh dishomogeneity, so that