Problem and background: The idea of predicting the weather by solving fluid equations was developed many years ago. In 1922 Richardson published a seminal contribution on predicting the weather by a numerical process, which created the first formulation of algorithms for the approximate solution of meteorological equations, leading ultimately to the predictions based on primitive equations. The shallow water equations (SWE) have become a standard tool for the investigation of prototypes of the meteorological models. In particular, they provide a testing ground for accurate and stable time integration methods in designing Numerical Weather Prediction (NWP) models. The equations describing flow of a thin layer of fluid on sphere $S$ can be cast in the form
$\dfrac{\partial {\bf u} } {\partial t } = – ({\rm Curl}_n {\bf u} + f {\bf n}) \times {\bf u} – {\rm grad} \bigg(\dfrac{ {|{\bf u}|^2} }{2} +g(h+h_s)\bigg), \quad \dfrac{\partial h } {\partial t } + {\rm div} (h {\bf u}) = 0$
where ${\bf u}$ is the smooth velocity field on $S$, $h$ is the thickness of the fluid layer, $h_s$ is the height of the surface level << the radius of the sphere), $g$ is the gravitational acceleration, $f=2\Omega \sin(\theta)$ is the Coriolis parameter, $\Omega$ is the angular velocity of the rotation, and $\times$ denotes the vector product in $\mathbb{R}^3$.
The difficult task of solving the SWE equations is mainly due to the existence of vastly differing time-scales in atmospheric phenomena, ranging from a relatively slow advection to very fast gravity waves. This poses a major challenge for time integration techniques and has been a central issue in the numerical solution of meteorological models over the past 70 years. Explicit (classical) methods are often forced to use very short time-steps due to stability for gravity waves. Semi-implicit schemes (such as the SISL schemes, see ) are thus commonly used in many global NWP centers as they can take longer time steps without sacrificing accuracy compared to classical methods. However, when increasing the complexity of the models, e.g., by including smaller scale processes, these standard methods require more computational effort and introduce unrealistic representation of high frequencies. Recently, Pudykiewicz et al. (JCP 2016) suggested the use of a classical third-order exponential multistep scheme epi3, which shows an advantage compared to semi-implicit schemes in terms of accuracy and is comparable efficiency to them. As flows progress into the nonlinear regime, however, it turns out that this scheme does not capture well their resolution due to the lower-accuracy approximation of the nonlinearity.
Contribution:
In this context, I ask the following question: Is it possible to develop a new method that can retain all the positive characteristics of the established methods (outlined above) and is able to take even longer time-steps as well as offers higher-accuracy in approximating the nonlinearity? Clearly, if one exists, such a method can reduce computational expense while maintaining sufficient solution accuracy. Working on this, I came up with the idea of using exponential Rosenbrock methods that is based on dynamic linearization of the flow, leading to a much smaller nonlinearity to solve in each time step. I performed evaluations and identified a candidate set of my three recent exponential Rosenbrock schemes of orders 4 and 5 (exprb42, pexrb43, exprb53) for use on the SWE model. I, along with Dr. Pudykiewicz (Environment Canada) and Dr. Reynolds (SMU), then investigated the performance of these methods on a suite of four challenging test problems, from the relatively simple flows exhibited in the Lauter test to complex turbulent flows exhibited in the unstable jet test, comparing them against epi3 investigated previously. In all cases, our methods enable accurate solutions at much longer time-steps than epi3, proving considerably more efficient as either the desired solution accuracy decreases or as the nonlinearity increases. For more details, see Luan, Pudykiewicz, and Reynolds: Journal of Computational Physics, Vol. 376, 817-837 (2019).
Time history of error ($l_{\infty}$) in the height field simulated over 15 days for mountain test
Conservation errors for total mass and energy for the Rossby-Haurwitz test
The vorticity field (scaled by $10^5$) for the unstable jet test.