Implicit–explicit Runge–Kutta methods applied to atmospheric chemistry-transport modelling

Chemistry-transport calculations are highly stiff in terms of time-stepping. Because explicit ODE solvers require numerous short time steps in order to maintain stability, it seems that especially sparse implicit–explicit solvers are suited to improve the numerical efficiency for atmospheric chemistry applications. In the new version of our mesoscale chemistry-transport model MUSCAT [Knoth, O., Wolke, R., 1998a. An explicit–implicit numerical approach for atmospheric chemistry–transport modelling. Atmospheric Environment 32, 1785–1797.], implicit–explicit (IMEX) time integration schemes are implemented. Explicit second order Runge–Kutta methods for the integration of the horizontal advection are used. The stiff chemistry and all vertical transport processes (turbulent diffusion, advection, deposition) are integrated in an implicit and coupled manner utilizing the second order BDF method. The horizontal fluxes are treated as ‘artificial’ sources within the implicit integration. A change of the solution values as in conventional operator splitting is thus avoided.The aim of this paper is to investigate the interaction between the explicit Runge–Kutta scheme and the implicit integrator. The numerical behavior is discussed for a 1D test problem and 3D chemistry-transport simulations. The efficiency and accuracy of the algorithm are compared to results obtained using the Strang splitting approach. The numerical experiments indicate that our second order implicit–explicit Runge–Kutta methods are a valuable alternative to the conventional operator splitting approach for integrating atmospheric chemistry-transport-models. In mesoscale applications and in cases with stronger accuracy requirements the ‘source splitting’ approach shows a better performance than Strang splitting.     

A semi-implicit time-splitting scheme for a regional nonhydrostatic atmospheric model

In this study a semi-implicit finite difference scheme is proposed for the nonhydrostatic atmospheric model based on the Euler equations for compressible ideal gas. Fast acoustic and gravity waves are approximated implicitly, while slow inertial processes are treated explicitly. Such time approximation requires the solution of 3D elliptic equations at each time step. An efficient elliptic solver is based on decoupling in the vertical direction and splitting in the horizontal directions. Stability analysis of the scheme shows that the time step is restricted only by the maximum wind speed and does not depend on the propagation velocity of the fast waves. A specific approximation of the advection terms keeping the second order of accuracy and possessing extended stability is employed to achieve larger time steps. The performed numerical experiments show the computational efficiency of the designed scheme and accuracy of the predicted atmospheric fields.     

Direct Discretization Method for the Cahn–Hilliard Equation on an Evolving Surface

We propose a simple and efficient direct discretization scheme for solving the Cahn–Hilliard (CH) equation on an evolving surface. By using a conservation law and transport formulae, we derive the CH equation on evolving surfaces. An evolving surface is discretized using an unstructured triangular mesh. The discrete CH equation is defined on the surface mesh and its dual surface polygonal tessellation. The evolving triangular surfaces are then realized by moving the surface nodes according to a given velocity field. The proposed scheme is based on the Crank–Nicolson scheme and a linearly stabilized splitting scheme. The scheme is second-order accurate, with respect to both space and time. The resulting system of discrete equations is easy to implement, and is solved by using an efficient biconjugate gradient stabilized method. Several numerical experiments are presented to demonstrate the performance and effectiveness of the proposed numerical scheme.     

Variably saturated flow described with the anisotropic Lattice Boltzmann methods

This paper addresses the numerical solution of highly nonlinear parabolic equations with Lattice Boltzmann techniques. They are first developed for generic advection and anisotropic dispersion equations (AADE). Collision configurations handle the anisotropic diffusion forms by using either anisotropic eigenvalue sets or anisotropic equilibrium functions. The coordinate transformation from the orthorhombic (rectangular) discretization grid to the cuboid computational grid is equivalent for the AADE to the anisotropic rescaling of the convection/diffusion terms. The collision components (eigenvalues and/or equilibrium functions) become discontinuous on the boundaries of the computational sub-domains which have different space scaling factors. We focus on the analysis of the boundary continuity conditions by using anisotropic LB techniques. The developed schemes are applied to Richards’ equation for variably saturated flow. The anisotropy of the Richard’s equation originates from distinct soil conductivity values in both the vertical and horizontal directions. The method should on the interface between the heterogeneous layers maintain the continuity of the normal component of the Darcy’s velocity (total flux). Also the method should accommodate steep jumps of the moisture content variable (conserved quantity) resulting from the continuity of the pressure variable, a given non-linear function of the moisture content. The coupling between heterogeneity and the anisotropy is examined by using the distinct space steps in neighboring layers and tested against uniform grid solutions. Different formulations of the Richard’s equation illustrate the construction of distinct diffusion forms and their integral transforms via specification of the equilibrium components. Integral transforms are used to overcome the difficulties coming from the rapid change of the main variables on sharp fronts. The numerical assessment of the stability criteria and the interface boundary conditions extend the analysis of the Lattice Boltzmann schemes to nonlinear problems with discontinuous coefficients.     

Robust and accurate viscous discretization via upwind scheme – I: Basic principle

In this paper, we introduce a general principle for constructing robust and accurate viscous discretization, which is applicable to various discretization methods, including finite-volume, residual-distribution, discontinuous-Galerkin, and spectral-volume methods. The principle is based on a hyperbolic model for the viscous term. It is to discretize the hyperbolic system by an advection scheme, and then derive a viscous discretization from the result. A distinguished feature of the proposed principle is that it automatically introduces a damping term into the resulting viscous scheme, which is essential for effective high-frequency error damping and, in some cases, for consistency also. In this paper, we demonstrate the general principle for the diffusion equation on uniform grids in one dimension and unstructured grids in two dimensions, for node/cell-centered finite-volume, residual-distribution, discontinuous-Galerkin, and spectral-volume methods. Numerical results are presented to verify the accuracy of the derived diffusion schemes and to illustrate the importance of the damping term for highly-skewed typical viscous grids.     

Lattice Boltzmann method for variable density shallow water equations

A lattice Boltzmann method is developed for solution of a form of the shallow water equations that is suitable for flows which are fully mixed in the vertical direction but have variable density in the horizontal plane. In the present approach, double distribution functions are applied: one for the shallow water flows and the other for the mass transport. Direct coupling between the water flow and mass transport is achieved by updating the flow density from the concentration during simulation. Accuracy and applicability of the model are demonstrated by two numerical tests: the stationary hydrostatic equilibrium of liquid of variable density in a tank with non-uniform bed terrain, and the horizontal diffusion of species with an initial Gaussian distribution of concentration in a uniform flow field.     

Three time integration methods for incompressible flows with discontinuous Galerkin Boltzmann method

This paper presents three time integration methods for incompressible flows with finite element method in solving the lattice-BGK Boltzmann equation. The space discretization is performed using nodal discontinuous Galerkin method, which employs unstructured meshes with triangular elements and high order approximation degrees. The time discretization is performed using three different kinds of time integration methods, namely, direct, decoupling and splitting. From the storage cost, temporal accuracy, numerical stability and time consumption, we systematically compare three time integration methods. Then benchmark fluid flow simulations are performed to highlight efficient time integration methods. Numerical results are in good agreement with others or exact solutions.     

DIMEX Runge–Kutta finite volume methods for multidimensional hyperbolic systems

We propose a class of finite volume methods for the discretization of time-dependent multidimensional hyperbolic systems in divergence form on unstructured grids. We discretize the divergence of the flux function by a cell-centered finite volume method whose spatial accuracy is provided by including into the scheme non-oscillatory piecewise polynomial reconstructions. We assume that the numerical flux function can be decomposed in a convective term and a non-convective term. The convective term, which may be source of numerical stiffness in high-speed flow regions, is treated implicitly, while the non-convective term is always discretized explicitly. To this purpose, we use the diagonally implicit–explicit Runge–Kutta (DIMEX-RK) time-marching formulation. We analyze the structural properties of the matrix operators that result from coupling finite volumes and DIMEX-RK time-stepping schemes by using M-matrix theory. Finally, we show the behavior of these methods by some numerical examples.     

Numerical Strategies for Solving the Nonlinear Rational Expectations Commodity Market Model

In this paper, I compare the accuracy, efficiency and stability of different numerical strategies for computing approximate solutions to the nonlinear rational expectations commodity market model. I find that polynomial and spline function collocation methods are superior to the space discretization, linearization and least squares curve-fitting methods that have been preferred by economists in the past.     

High-order compact ADI (HOC-ADI) method for solving unsteady 2D Schrödinger equation

In this paper, a high-order compact (HOC) alternating direction implicit (ADI) method is proposed for the solution of the unsteady two-dimensional Schrödinger equation. The present method uses the fourth-order Padé compact difference approximation for the spatial discretization and the Crank-Nicolson scheme for the temporal discretization. The proposed HOC-ADI method has fourth-order accuracy in space and second-order accuracy in time. The resulting scheme in each ADI computation step corresponds to a tridiagonal system which can be solved by using the one-dimensional tridiagonal algorithm with a considerable saving in computing time. Numerical experiments are conducted to demonstrate its efficiency and accuracy and to compare it with analytic solutions and numerical results established by some other methods in the literature. The results show that the present HOC-ADI scheme gives highly accurate results with much better computational efficiency.     

A High Order Conservative Semi-Lagrangian Discontinuous Galerkin Method for Two-Dimensional Transport Simulations

In this paper, we develop a class of high order conservative semi-Lagrangian (SL) discontinuous Galerkin methods for solving multi-dimensional linear transport equations. The methods rely on a characteristic Galerkin weak formulation, leading to \(L^2\) stable discretizations for linear problems. Unlike many existing SL methods, the high order accuracy and mass conservation of the proposed methods are realized in a non-splitting manner. Thus, the detrimental splitting error, which is known to significantly contaminate long term transport simulations, will be not incurred. One key ingredient in the scheme formulation, borrowed from CSLAM (Lauritzen et al. in J Comput Phys 229(5):1401–1424, 2010), is the use of Green’s theorem which allows us to convert volume integrals into a set of line integrals. The resulting line integrals are much easier to approximate with high order accuracy, hence facilitating the implementation. Another novel ingredient is the construction of quadratic curves in approximating sides of upstream cell, leading to quadratic-curved quadrilateral upstream cells. Formal third order accuracy is obtained by such a construction. The desired positivity-preserving property is further attained by incorporating a high order bound-preserving filter. To assess the performance of the proposed methods, we test and compare the numerical schemes with a variety of configurations for solving several benchmark transport problems with both smooth and nonsmooth solutions. The efficiency and efficacy are numerically verified.     

Über die Integration von Differentialgleichungssystemen 1. Ordnung mit exponentiell angepaßten numerischen Methoden

A general theoretical basis for the design of exponentially fitted numerical integration methods will be derived for obtaining methods which, from accuracy and stability, allow large discretization intervals even for systems with large eigenvalues. The technique of transforming existing numerical methods to an exponentially fitted version will be exemplified and the influence of exponential fitting on discretization error and stability for one-step methods will be discussed.     

Solution of contaminant transport with equilibrium and non-equilibrium adsorption

In this paper we present a numerical approximation scheme for the solution of contaminant transport problems with diffusion and adsorption in equilibrium and non-equilibrium mode. The method is based on time stepping and operator splitting. The non-linear transport is solved semi-analytically via the multiple Riemann problem, the non-linear diffusion by a finite volume method and by Newton’s type of linearisation, and finally the reaction part, incorporating the non-equilibrium adsorption, is transformed to an integral equation which is solved numerically using time discretization. Various results of numerical experiments are shown, and the method is applied to the practical dual-well problem.     

Russian scatterometer: discussion of the concept and the numerical simulation of wind field retrieval

Orbital scatterometry is briefly overviewed and its trends are indicated. Two scatterometer concepts are currently considered for trade-offs: with fixed and rotating antenna systems. The concept with a rotating antenna system was selected and SeaWinds was chosen as the prototype for the first Russian scatterometer. The scatterometer concept was then further developed and instead of two pencil beams, a fan-beam antenna was proposed (about 1° × 6°). The fan-beam antenna allows successive measurements for horizontal and vertical polarization in each wind vector cell (WVC). This increases the number of observations of the WVC at different incidence and azimuth angles during flight. The scatterometer parameters required to implement the proposed measurement geometry for an orbit altitude of 650 km and a swath width of 1525 km are discussed. A numerical scatterometer model that accounts for both the specifications and the observation geometry is developed. The scatterometer performance, with subsequent formation of a swath and splitting into WVCs, is simulated. The procedure of wind vector retrieval includes two stages: 1) determining wind speed and wind direction in a single WVC; and 2) using the information from adjacent WVCs to correct wind direction. It is shown that the accuracy of wind direction retrieval by a WVC can be increased by simultaneous radar cross-section (RCS) measurements at vertical and horizontal polarization. The basic error in determining wind direction is due to a 180° wind direction ambiguity caused by the form of RCS azimuth dependence. Two-dimensional median filtering is commonly employed in scatterometry to increase the accuracy of wind direction retrieval. In this study, two-dimensional angular median filtering was employed and it is shown that the error in wind direction retrieval significantly decreased. The results of the research indicate that wind field can be retrieved by the new scatterometer with the level of precision required.     

A family of stable numerical solvers for the shallow water equations with source terms

In this work we introduce a multiparametric family of stable and accurate numerical schemes for 1D shallow water equations. These schemes are based upon the splitting of the discretization of the source term into centered and decentered parts. These schemes are specifically designed to fulfill the enhanced consistency condition of Bermúdez and Vázquez, necessary to obtain accurate solutions when source terms arise. Our general family of schemes contains as particular cases the extensions already known of Roe and Van Leer schemes, and as new contributions, extensions of Steger–Warming, Vijayasundaram, Lax–Friedrichs and Lax–Wendroff schemes with and without flux-limiters. We include some meaningful numerical tests, which show the good stability and consistency properties of several of the new methods proposed. We also include a linear stability analysis that sets natural sufficient conditions of stability for our general methods.     

A numerical method for the vorticity-velocity Navier-Stokes equations in two and three dimensions

This paper provides a numerical method for solving the steady-state vorticity-velocity Navier-Stokes equations in two and three dimensions. The vorticity transport equation is considered together with a Poisson equation for the velocity vector, the latter equation being parabolized in time according to the false transient approach. The two vector equations are discretized in time using the implicit Euler time stepping and the delta form of Beam and Warming. A staggered-grid spatial discretization is employed in conjunction with a deferred correction procedure. Second-order-accurate central differences are used to approximate the steady-state residuals, written in conservative form for accuracy reasons, whereas upwind differences are used for the advection terms in the implicit operator, to obtain diagonally-dominant tridiagonal systems. The discrete equations are solved sequentially by means of a robust alternating direction line-Gauss-Seidel iteration procedure combined with a simple multigrid strategy. For the model driven-cavity-flow problem in two and three dimensions, the method is found to be efficient and very accurate. For the first time, the three-dimensional discrete vorticity and velocity fields, computed using a Poisson equation for the velocity vector, are both solenoidal and satisfy their mutual relationship, exactly.     

Discretization of parabolic inequalities

We consider the time discretization with linear multistep methods of an abstract parabolic variational inequality. For such a discretization, we prove stability and convergence under suitable regularity assumptions. Moreover, we prove that, under suitable assumptions, the methods considered converge with order at least one. We performed numerical essays which lead to the conjecture that the order is higher.     

Fractional Jacobi Galerkin spectral schemes for multi-dimensional time fractional advection–diffusion–reaction equations

One of the ongoing issues with time fractional diffusion models is the design of efficient high-order numerical schemes for the solutions of limited regularity. We construct in this paper two efficient Galerkin spectral algorithms for solving multi-dimensional time fractional advection–diffusion–reaction equations with constant and variable coefficients. The model solution is discretized in time with a spectral expansion of fractional-order Jacobi orthogonal functions. For the space discretization, the proposed schemes accommodate high-order Jacobi Galerkin spectral discretization. The numerical schemes do not require imposition of artificial smoothness assumptions in time direction as is required for most methods based on polynomial interpolation. We illustrate the flexibility of the algorithms by comparing the standard Jacobi and the fractional Jacobi spectral methods for three numerical examples. The numerical results indicate that the global character of the fractional Jacobi functions makes them well-suited to time fractional diffusion equations because they naturally take the irregular behavior of the solution into account and thus preserve the singularity of the solution.

     

The role of the kinetic parameter in the stability of two-relaxation-time advection–diffusion lattice Boltzmann schemes

In general, explicit numerical schemes are only conditionally stable. A particularity of lattice Boltzmann multiple-relaxation-time (MRT) schemes is the presence of free (“kinetic”) relaxation parameters. They do not appear in the transport coefficients of the modelled second-order (macroscopic) equations but they have an impact on the effective accuracy and stability of the algorithm. The simplest uniform choice (the well known BGK/SRT model) is often inadequate, and therefore a compromise in the complexity of the model is sought. For this purpose, the von Neumann stability analysis is performed for the d1Q3 two-relaxation-time (TRT) advection–diffusion model. The extended optimal (EOTRT) model, which relates the two collision times such that the most stable scheme is set by a suitable choice of the equilibrium parameters, equal for any Peclet number, is then developed. This extends the very recently derived optimal subclass (OTRT) to larger combinations of “physical” and “kinetic” collision rates. Next, we provide the necessary and/or sufficient stability limits on the EOTRT subclass for a wide range of velocity sets, with and without numerical diffusion, and delineate the interesting choices of free equilibrium weights for the d2Q9 and d3Q15 models. The BGK/SRT model is without advanced advection properties; we prove (for minimal stencil schemes d1Q3, d2Q5 and d3Q7) that the non-negativity of the equilibrium distribution is necessary for its stability in the advection-dominated limit. Beyond the EOTRT and BGK/SRT subclasses of the TRT model, blind choices of the “ghost” collision number may result in quite unstable schemes, even for positive equilibrium. However, we find that the d1Q3 stability curves govern the advection properties of the multi-dimensional models and a fuller picture of the TRT stability properties begins to emerge.