1 Introduction
Mathematical models of continuum mechanics Batchelor ; LandauLifshic1986
describe the transport of scalar and vector quantities due to advection. In particular, the basic equation of hydrodynamics is the continuity equation. Advective transfer causes the fulfillment of conservation laws
Godunov ; LeVeque . In addition, there are properties of the positivity and monotonicity of the solution. Such important properties of the differential problem must be inherited when passing to the discrete problem Anderson ; Wesseling .For spatially approximation, conservative approximations are constructed on basis of using the conservative (divergent) form of the advection equation. Most naturally such technology is implemented when using the integrointerpolation method (balance method) on regular and irregular grids
Samarskii1989 , in the control volume method LeVeque ; Versteeg . The construction of monotonic approximations is discussed in many papers (see, for example, kulikovskii2000mathematical ; HundsdorferVerwer2003 ; Kuzmin ). In MortonKellogg1996 ; SamarskiiVabischevich1999a standard linear approximations are considered for convectiondiffusion problems.Currently, the main computing technology for solving applied problems is the finite element method Guermond ; Larson . It is widely used in computational fluid dynamics Donea ; Zienkiewicz . Monotonization of the solution is achieved by using various linear and nonlinear variants of stabilization techniques. It should be noted that the standard formulation of the equations of continuous medium mechanics in the conservative or nonconservative form is poorly suited for applying finite element approximations, for which the Hilbert spaces are natural.
Separate attention deserves the problems of constructing and investigating approximations in time. When solving boundary value problems for partial differential equations, twolevel schemes (
method, schemes with weights) are traditionally widely used HundsdorferVerwer2003 ; Ascher2008 ; LeVeque2007 . Research of this schemes can be based on the general theory of stability (wellposedness) of operatordifference schemes Samarskii1989 ; SamarskiiMatusVabischevich2002. In particular, unimprovable (coinciding necessary and sufficient) stability conditions can be used, which are formulated as operator inequalities in a finitedimensional Hilbert space. The achieved level of theoretical stability research allows us to abandon the heuristic methods widely used in computational fluid dynamics to study the stability of difference schemes: the von Neumann method for stability analysis, Fourier analysis, the principle of frozen coefficients, the consideration of the problem without taking into account the boundary conditions.
In this paper, the standard finiteelement approximation in space is used for the nonstationary equation. The advection equation is written in the socalled symmetric form SamarskiiVabischevich1999a , when the advection operator is the halfsum of the advection operators in the conservative (divergent) and nonconservative (characteristic) forms. Thus the continuum mechanics equations are written using the SD (Square root from Density) variables VabForm ; form1 . In this case the corresponding conservation laws are a direct consequence of the skewsymmetry of the advection operator. The conservativeness property is related to the preservation of the norm of the solution of the nonstationary advection equation, with stability with respect to the initial data. This property takes place not only for the solution, but also for some of its transformations. In this case, we are talking about the property of multiconservativeness. The principal point is related with the fact that the most important skewsymmetry property of the advection operator is inherited for finite element approximations.
For the advection equation an explicit twolevel scheme as well as all schemes with a weight lower than 0.5, are aboslutly unstable. At the same time, extensive computational practice aims us to use the explicit schemes in these problems. In such schemes, the conditional stability (the CourantFriedrichsLewy (CFL) condition) with a timestep limited by the Courant number is provided in fact by refusing the skewsymmetry of the advection operator — the use of dissipative approximations. In addition, the standard explicit schemes are not conservative.
We construct conditionally stable schemes for the advection equation based on the principle of regularization of operatordifference schemes Samarskii1967 ; Vabishchevich2014 , when stability is provided by a small perturbation. The explicit secondorder RungeKutta scheme Butcher2008 ; HairerWanner2010 is considered. In this context, the classical explicit LaxWendroff scheme lax1964difference is also considered as regularization of the secondorder RungeKutta scheme — the perturbation of the advection operator squared. Using the stability criteria of twolevel operatordifference schemes, the stability conditions of regularized schemes and the explicit LaxWendroff scheme are obtained. Effective computational implementation of explicit schemes using finite element approximation in space is provided by using the diagonalization procedures of the mass matrix (masslumping procedure) Thomee2006 ; Cohen .
Of greatest interest are implicit twolevel schemes for the advection equation, which belong to the unconditionally stable class. The classical CrankNicolson scheme has a secondorder of accuracy, is unconditionally stable and multiconservative. For the advection problems under the consideration, we can use a scheme of the fourthorder of accuracy, which is also unconditionally stable and multiconservative. A certain drawback of this scheme is associated with the need to use the lumping procedure. An implicit version of the LaxWendroff scheme is proposed, which is conditionally stable, but has a higher accuracy than the explicit LaxWendroff scheme and does not require the diagonalization of the mass matrix.
The paper is organized as follows. A model twodimensional problem for the advection equation is formulated in Section 2. Approximation in space is constructed using Lagrangian finite elements, the main properties of the problem solution are noted. In Section 3, we consider known and new explicit difference schemes for the advection equations, and investigate the stability conditions. Central for this work is Section 4. Implicit schemes, their stability and conservatism are studied here. In Section 5, numerical experiments on the accuracy of explicit and implicit schemes are discussed for the model IBV problem. The results of the work are summarized in Section 6.
2 Problem statement
In a bounded twodimensional domain , we consider the advection equation written in the symmetric form. A standard finite element approximation in space is used. The problem of constructing approximations in time is formulated in such a way that the approximate solution inherits the basic properties of the solution of the differential problem.
2.1 Differential problem
The Cauchy problem is considered in the domain ()
(1) 
(2) 
using notation . The operator of advection (convective transport) is assumed to be constant and is written in the symmetric form:
(3) 
Thus, we take the halfsum of the transfer operator in the divergent (conservative part with ) and nondivergent (characteristic part with ) forms SamarskiiVabischevich1999a ; zang1991rotation . The convective transport is determined by the velocity of the medium , and the nopermeability condition on the boundary of the domain
(4) 
where is the outer normal to the boundary .
The standard continuity equation has a divergent form:
(5) 
where is the density. From (4), (5), by direct integration, we obtain the law of mass conservation:
When orienting to finite elemental approximations in space, it is convenient to use of SD (Square root from Density) variables for the equations of hydrodynamics VabForm . Let , then equation (5) is written in the form (1), (3).
In the Hilbert space , we define the scalar product and norm in the standard way:
Under constraint (4), the convective transfer operator is skewsymmetric:
(6) 
and therefore, in particular,
Using the scalar product of equation (1) and , taking into account (2), (6), we obtain
(7) 
This relation with respect to the continuity equation (5) expresses the law of mass conservation. The property (7) for the solution of the the Cauchy problem (1), (2) is associated with the nondissipativity of the system. The equation (1) for (6) has many other conservation laws. If the constant operator is permutable with , then we have the equality
(8) 
In particular, equality (8) holds for .
2.2 Finite element approximation in space
To solve numerically the problem (1)–(3), we employ finite element approximations in space (see, e.g., Thomee2006 ; brenner2008mathematical ). For (3), we define the bilinear form
By (4), we have
Define the subspace of finite elements and the discrete operator as
The operator acts on the finite dimensional space and, similarly to (8), is skewsymmetric:
(9) 
In finite element approximation, the main property of the advection operator, namely, its skewsymmetry, is inherited and, as a consequence, there is energy neutrality (the equality ).
The Cauchy problem (1), (2) is associated with the problem
(10) 
(11) 
for , where with denoting projection onto . Similarly (7), for the solution of the problem (10), (11), the conservative property is established:
(12) 
The multiconservative property (see (8)) is associated with the equality
(13) 
provided that
When the problem (10), (11) is approximated in time, it is necessary to focus on the fulfillment of the properties (12) and (13) at separate time levels.
2.3 Approximation in time
Let, for simplicity, be a step of a uniform grid in time such that , . To solve numerically the problem (10), (11), we use explicit and implicit twolevel schemes, when the solution at the new level is determined by the previously found solution , .
When approximating in time, we focus, first of all, on the scheme stability. With respect to our problem, stability will be ensured by the following estimate of the solution at each time level:
(14) 
The most favorable situation is associated with the fulfillment of the equality
(15) 
In this case, not only stability is ensured, but also the conservatism of the solution takes place (see (12)). Multiconservatism (see (13)) is due to the fulfillment of the equality
(16) 
for some operators .
Conservative approximations are related to the property of time reversibility for the problem (9)–(11). The computational algorithm is reversible in time if we calculate the solution at the level and then change the transfer velocity () to the opposite (by ), then at the next step in time we get the solution that coincides with the solution at the level .
3 Explicit schemes
The explicit scheme for the problem (9)–(11) has the firstorder error in time and is absolutely unstable. A conditionally stable scheme can be constructed on the basis of the regularization principle of difference schemes. We also consider a scheme of secondorder approximation, namely, the LaxWendroff scheme.
3.1 Necessary and sufficient conditions for the stability of schemes with weights
Strict results on the stability of difference schemes in finitedimensional Hilbert spaces are obtained in the theory of stability (wellposedness) of operatordifference schemes Samarskii1989 ; SamarskiiMatusVabischevich2002 . With respect to the subject of our investigation, we give the best possible conditions for the stability of a standard twolevel scheme with weights (method) to be unimprovable (matching necessary and sufficient) for the solution of the problem (10), (11).
When passing from the level to the level , the approximate solution of the problem (10), (11) is determined from the equation
(17) 
with some constant operator . In the simplest case, . The initial condition is
(18) 
If , we have then explicit scheme, for , we obtain the fully implicit scheme and for , the CrankNicolson scheme appears. The stability criterion is formulated as follows.
Lemma 1
The proof of this statement is given, for example, in the book SamarskiiGulin1973 . The condition (19) can be written in the form of the operator inequality
3.2 Explicit schemes of the first and second order of approximation
It is natural to start with the simplest explicit scheme for the problem (10), (11). In this case, in (17), i.e.
(20) 
For the explicit scheme, the stability criterion (19) takes the form
(21) 
In the case of a skewsymmetric operator , inequality (21) is not satisfied for any . Therefore, the explicit scheme (18), (20) is absolutely unstable.
Among twolevel schemes, explicit schemes of the secondorder approximation are deserved separate consideration. Instead of (20) we will use the explicit secondorder RungeKutta scheme:
(22) 
In this case, we have
Taking into account (9), we get
In view of this, the stability condition (21) takes the form
Again, we cannot specify such that the stability of the scheme (22) holds.
Thus, we can formulate the following result.
3.3 Conditionally stable nonstandard scheme
We note the possibility of constructing a conditionally stable scheme based on the correction of the explicit scheme (18), (20). We consider more weak stability requirements, allowing the growth of the solution norm in accordance with (14) with the choice of some positive constant .
The construction of the considering approximation in time is based (see, for example, afanas2013unconditionally ; VabExact ) on the solution representation of the problem (10), (11) in the form
(23) 
For , from (10), (11), (23), we have the Cauchy problem
(24) 
(25) 
where is the identity operator.
For the approximate solution of the problem (24), (25), the explicit difference scheme is used
(26) 
(27) 
Taking into account relation
scheme (26), (27) corresponds to the use of
(28) 
with the initial condition (18). Such schemes belong to the class of nonstandard mickens1994nonstandard ; mickens2002nonstandard . The levelwise estimate
for the solution of the difference scheme (26), (27) corresponds to the estimate (14) for the scheme (18).
Theorem 3
3.4 Regularized schemes
The regularization principle for difference schemes Samarskii1989 ; Vabishchevich2014 provides great opportunities in constructing difference schemes of a prescribed quality. The standard approach to the construction of stable schemes on the basis of the regularization principle is associated with the introduction of additional terms (regularizers) in operators of an origional (generating) difference scheme.
Instead of the explicit scheme (20), we will use the regularized scheme
(30) 
with the regularizer . Let us formulate the stability conditions for this scheme.
For the scheme (30), we have
Taking into account the skewsymmetry of the operator , we obtain
The stability condition (21) gives
(31) 
Therefore, we can rely on the conditional stability of the scheme (30) for .
Let
be the maximum eigenvalue of the spectral problem
(32) 
and . Then the inequality (31) will be satisfied for
(33) 
Thus, we can formulate the stability conditions.
Theorem 4
We highlight some interesting possibilities for the choice of the regularizer . The simplest version is related to the regularizer
(34) 
The conditions (32), (33) give the following restrictions on the time step:
(35) 
Thus, the explicit scheme (18), (30), (34) is conditionally stable under the time step constraints of Courant type (see (35)). In the special case, for , we have
The regularized scheme (18), (30), (34) has the firstorder approximation in time. It is natural to expect a higher accuracy when choosing , i.e. when the scheme (30) is close to the explicit secondorder RungeKutta scheme (22). Here we can indicate two variants of such schemes.
The first variant is related with choosing the regularizer in the form
(36) 
The stability condition (21) for this case takes the form
It will be satisfied for
(37) 
The restrictions (37) for the scheme (18), (30), (36) are substantially more strong than the restrictions (35) for the scheme (18), (30), (34). This disadvantage is compensated by the fact that we have the secondorder approximation in time instead of the first order.
Let us now consider the second variant of the regularized scheme of the secondorder approximation. The most wellknown modification of the absolutely unstable explicit secondorder RungeKutta scheme (18), (22) is the LaxWendroff scheme lax1964difference ; leveque2002finite . This scheme is actually based on replacing the operator in the explicit secondorder RungeKutta scheme by a selfadjoint nonnegative operator close to it.
3.5 Computational implementation of explicit schemes
The explicit schemes under the consideration are explicit only in form, according to the time approximation. The computational implementation of explicit finiteelement approximations, unfortunately, is connected with the solution of systems of linear equations on a new level in time. The standard approach is related to the correction of approximations based on mass lumping (see, for example, Thomee2006 ).
The use of certain schemes for the Cauchy problem (10), (11
) is based on the solution of matrix problems for finding an approximate solution at a new time level. We associate with the differentialoperator equation the corresponding system of ordinary differential equations.
We consider a standard quasiuniform triangulation of the domain into triangles (or tetrahedra in 3D). Let be vertices of this triangulation. We introduce the finite dimensional space of continuous functions that are liner over each finite element, see, e.g. Thomee2006 . As a nodal basis we take the standard hat function . Then for , we have the representation
where .
From (10), we obtain the equation
(39) 
where is the vector of unknowns . Here is the mass matrix and is the stiffness matrix, and
If we set , then equation (39) is written in the form (9), (10) for
(40) 
Thus, we can construct approximations in time for the equation (39) on the basis of approximations for equation (10), taking into account that and (40).
For example, the explicit scheme (20) corresponds to the scheme
(41) 
It is necessary to solve the problem
at a new level in time. Because of this, the scheme (41) is explicit in terms of time approximation, but implicit in terms of the computational implementation, since the symmetric matrix is offdiagonal.
In order to provide an explicit computational implementation, various diagonalization procedures for the mass matrix are used:
In the simplest mass lumping procedure Thomee2006 we have
Instead of (39), we seek an approximate solution of the Cauchy problem for the equation
This equation corresponds to setting
(42) 
4 Implicit schemes
Unconditionally stable schemes for the advection equation are built on the basis of implicit approximations. It seems reasonable to employ the CrankNicolson scheme, which has the second order of accuracy in time and has conservative properties. In addition, schemes of higher accuracy are outlined, which have the fourthorder approximation. The implicit version of the LaxWendroff scheme is highlighted.
4.1 The CrankNicolson scheme
Among the implicit schemes for the problem (10), (11), the most important is the symmetric scheme whith in (17). In particular, it has the best SM (Spectral Mimetic) properties in the class of schemes with the skewsymmetric operator Vabischevich2011 . In this case, the solution is determined from the equation
(43) 
By virtue of (19), the scheme (18), (43) is absolutely stable and approximates the problems (10), (11) with the second order in .
The CrankNicolson scheme is nondissipative. To show this fact, it is sufficient to multiply the scalar equation (43) by . Taking into account the skewsymmetry of the operator , we have
Taking into account (18), we arrive at (15). Also there is the multiconservative property, see (13). For any constant operator , which commutes with , (16) is satisfied. We immediately see that the scheme (18), (43) generates a computational algorithm that is time reversible. The result of our consideration is the following theorem.
4.2 Scheme of higher accuracy order
The CrankNicolson scheme (43) corresponds to the use of the following Pade approximant for the exponential:
When considering more accurate approximations, we can consider
The corresponding difference scheme () has the form
(44) 
We set
and apply scalar multiplication of equation (44) by . This leads to the equality
For , this equality, taking into account , implies the stability estimate, the conservativeness and multiconservative properties. The positivity of the operator is satisfied for at not very large steps in time. Really,
Thus for
(45) 
The restrictions (45) can be removed by modifying the proof. We write the scheme (44) in the form
(46) 
with notation
(47) 
and . The first question is related with the proof of the existence . Necessary and sufficient condition for the existence of the inverse operator is (see, for example, (Anderson, ; Lusternik, )) the fulfillment of the inequality
wherein .
Lemma 6
For the operator , which is defined according to (47), the inequality is satisfied
(48) 
[Proof.] The inequality (48) is equivalent to the inequality
We have
taking into account the skewsymmetry of the operator .
Taking into account lemma 6, we can write (46) as follows
(49) 
For the transition operator from one level in time to another, the following statement holds.
Lemma 7
The operator under the conditions (47) is unitary, i.e.
(50) 
[Proof.] We have
Taking into account the permutability of the operator factors on the righthand side, we get the equality (50).
The noted unitarity property (50) of the operator allows (49) to come to the conservativity property (15). Similarly, the multiconservative property is established. The result of our consideration is the following statement.
Theorem 8
In the computational implementation, we solve the Cauchy problem for equation (39). Because of this, the representation (40) is used in equation (10). The scheme (44) corresponds to the scheme
(51) 
Thus, for the transition to a new time level, it is necessary to solve equation with the operator
As in the case of explicit schemes, the diagonalizing procedure of the mass matrix saves the situation. In this case, instead of (40), we use the representation (42), and the scheme (45) is replaced by the scheme
(52) 
The scheme (52) is stable under constraint (45) taking into account the fact that the representation (42) holds.
4.3 Implicit LaxWendroff scheme
On the base of the implicit scheme (44), we can construct an implicit version of the LaxWendroff scheme. Taking into account the notation introduced above, instead of (44), we will use the scheme
(53) 
As in the case of the explicit scheme, the construction is based on replacing the operator by the operator , which is defined by (38).
We rewrite the scheme (53) in the form
where
The operator explicitly highlights the approximation error in space due to the replacement of by .
For a positive constant selfadjoint operator , we define the Hilbert space , where the scalar product and norm are defined as follows
Analogously to theorem 8, the following result is formulated.
Theorem 9
5 Numerical experiments
The possibilities of the explicit and implicit schemes under the consideration for solving the advection equation are illustrated by the numerical results for the model problem.
5.1 Model problem
We will consider the problem (1)(3) in the unit square:
The initial condition is taken in the form
The components of the velocity are defined by the stream function
so that
The calculations are performed for .
The computational code was implemented using the FEniCS fenics numerical framework. The finite element approximation in space is based on the use of continuous Lagrange element, namely, piecewiselinear elements. A uniform grid is used for spatial domain. The grid with the step is shown in Fig. 1.
The accuracy of different approximations in time will be estimated by a reference solution. It was obtained using the scheme under the consideration with an essentially small time step: . The timeevolution of the solution on the grid with (the main grid in space) is illustrated in Fig. 2. For the relative error of the approximate solution, we have
where is the reference solution.
5.2 Explicit schemes
Prevously the explicit regularized schemes (18), (30) have been outlined. For the choice of a regularizer in the form (34), stability takes place with constraints (35) on the time step. It should be noted that we must orient on schemes with the operator (42).
We present the results of calculations for the regularized scheme with . The constraints for the step are related with the norm of the operator . Taking into account its skewsymmetry property, we have
where are the eigenvalues of the operator :
Taking into account the representation (40), this spectral problem corresponds to the spectral problem
Similarly, using the procedure of mass lumping with allowance for (42), we get the spectral problem
To solve the spectral problems with symmetrical matrices, we use the SLEPc library (Scalable Library for Eigenvalue Problem Computations) hernandez2005slepc ). We apply the KrylovSchur algorithm, a variation of the Arnoldi method, proposed by stewart2002krylov . The calculation results of the norm of the operator according to (40) and (42) on different grids are presented in the table 1. It should be noted that and the norm of the operator decreases approximately by two times when using mass lumping, because the maximum of permissible time step (see (35)) increases approximately twice.
space grid ()  for  for 

0.02  1.05288993e+02  5.59579462e+01 
0.01  2.16001186e+02  1.14622718e+02 
0.005  4.37491174e+02  2.31964151e+02 
The timehistory of the error for various time steps for the regularized scheme (18), (30), (34) with (42) and are shown in the Fig. 3. Convergence with first order in is observed.
Using the mass lumping procedure and taking into account (38), the explicit LaxWendroff scheme is written as
(56) 
It corresponds to the regularized scheme (30), whereh the operator is defined according to (42), and for the operator , we have
(57) 
The stability conditions of the scheme (30) are formulated in theorem 4. The inequality will be satisfied if the constant in the inequality . We consider the spectral problem
then . Taking into account (42), (57), this spectral problem corresponds the spectral problem
(58) 
The constraints on time step (33) is related to the spectral problem (32). In the case (42), (56), it is equivalent to the spectral problem
(59) 
The calculation results for the determination of and from the numerical solution of the spectral problems (58), (59) are given in table 2. These data demonstrate the conditional stability of the LaxWendroff scheme, moreover .
space grid ()  

0.02  1.00098795e+00  1.73477111e02 
0.01  1.00025320e+00  8.47323207e03 
0.005  1.00006414e+00  4.17705891e03 
We write the explicit LaxWendroff scheme (56) in the form
Here the term with
distinguishes this scheme from the explicit RungeKutta scheme of the secondorder accuracy in time. In our problem, the convergence of the approximate solution with the second order in time is manifested (see Fig. 4), first of all, in the initial time interval, at which the smoothness of the solution is large enough and the influence of the term is insignificant. The first order of accuracy due to the term begins to appear at .
5.3 Implicit schemes
In the class of implicit schemes, the CrankNicolson scheme (18), (43) seems to be the basic one. Its accuracy in solving the model problem is illustrated in Fig. 5. In comparison with the explicit regularized scheme (see Fig. 3), a higher accuracy of the approximate solution is observed, convergence is approximately of the second order in . In addition, this scheme stands out among all those considered schemes due to the fact that it is unconditionally stable.
Similar results for the scheme of higher order of accuracy are shown in Fig. 6. The approximate solution using the lumping procedure is found from the equation (52). This scheme is absolutly stable and demonstrates the high accuracy of the approximate solution for substantially large time steps.
The accuracy of the implicit version of the LaxWendroff scheme is illustrated in Fig. 7. As for the explicit scheme (see Fig. 4), there is a higher accuracy for the initial time, when the solution is smooth in space. Further, the effect of occures that results in the difference between the implicit LaxWendroff scheme (53) and the fourthorder accuracy scheme (44). The scheme is stable for , where corresponds to the right side (54). Numerical results for and , which are obtained from the solution of the spectral problem for the operator , are given in table 3 and show that and .
space grid ()  

0.02  3.22933843e+04  1.92767512e02 
0.01  1.33745164e+05  9.47221570e03 
0.005  5.44513748e+05  4.69446600e03 
6 Conclusions

In numerical simulation of advection processes, it is necessary to focus on writing the equation in a symmetric form (the halfsum of the advection operator in the divergent (conservative) form and the advection operator in the nondivergent (characteristic) form). In this case the advection operator is skewsymmetric. For the solution of the nonstationary advection equation, the multiconservative property holds, which is associated with the fulfillment of the set of conservation laws.

The stability of known and new twolevel difference schemes is investigated for the approximate solution of the Cauchy problem for the advection equation. Investigation of stability is carried out on the basis of the general theory of stability (wellposedness) of operatordifference schemes. The standard finiteelement approximation is used in space, which ensures the skewsymmetry of the discrete advection operator.

The standard explicit scheme for the advection equation belongs to the class of absolutely unstable ones. The explicit regularized scheme of the first order of accuracy is proposed, and the stability conditions are formulated. The stability conditions are established for the finiteelement version of the classical explicit LaxWendroff scheme. For finiteelement approximation in space, it is necessary to use various diagonalization algorithms of the mass matrix.

For solving advection problems, the CrankNicolson scheme seems to be very attractive. It demonstrates the second order of accuracy and belongs to the class of unconditionally stable schemes. In addition, it demonstrates many conservation laws. The possibilities of using schemes with higher order of accuracy, which is also unconditionally stable and multiconservative, are also considered. The implicit version of the LaxWendroff scheme is proposed, and conditions for its stability are formulated.

The properties of the explicit schemes are illustrated by the solution of the model twodimensional problem for the advection equation. The focus is on studing the accuracy of various approximations in time.
Acknowledgements
This work is supported by the megagrant of the Russian Federation Government (# 14.Y26.31.0013).
References
 (1) G. K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, 2000.
 (2) L. D. Landau, E. Lifshitz, Fluid Mechanics, ButterworthHeinemann, 1987.
 (3) S. K. Godunov, E. I. Romenskii, Elements of Continuum Mechanics and Conservation Laws, Springer, 2003.
 (4) R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002.
 (5) J. Anderson, Computational Fluid Dynamics: The Basics with Applications, McGrawHill, 1995.
 (6) P. Wesseling, Principles of Computational Fluid Dynamics, Springer, 2010.
 (7) A. A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, New York, 2001.
 (8) H. Versteeg, W. Malalasekra, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Prentice Hall, 2007.
 (9) A. G. Kulikovskii, N. V. Pogorelov, A. Y. Semenov, Mathematical Aspects of Numerical Solution of Hyperbolic Systems, Taylor & Francis, 2000.
 (10) W. H. Hundsdorfer, J. G. Verwer, Numerical Solution of TimeDependent AdvectionDiffusionReaction Equations, Springer Verlag, 2003.
 (11) D. Kuzmin, A Guide to Numerical Methods for Transport Equations, University ErlangenNuremberg, 2010.
 (12) K. W. Morton, R. B. Kellogg, Numerical Solution of ConvectionDiffusion Problems, Chapman & Hall London, 1996.
 (13) A. A. Samarskii, P. N. Vabishchevich, Numerical Methods for Solving ConvectionDiffusion, URSS, Moscow, 1999.
 (14) A. Ern, J.L. Guermond, Theory and Practice of Finite Elements, Springer, 2004.
 (15) M. G. Larson, F. Bengzon, The Finite Element Method: Theory, Implementation and Applications, Springer, 2013.
 (16) J. Donea, A. Huerta, Finite Element Methods for Flow Problems, Wiley, 2003.
 (17) O. C. Zienkiewicz, R. L. Taylor, N. P., The Finite Element Method for Fluid Dynamics, ButterworthHeinemann, 2013.
 (18) U. M. Ascher, Numerical Methods for Evolutionary Differential Equations, Society for Industrial Mathematics, 2008.
 (19) R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: SteadyState and TimeDependent Problems, Society for Industrial Mathematics, 2007.
 (20) A. A. Samarskii, P. P. Matus, P. N. Vabishchevich, Difference Schemes with Operator Factors, Kluwer Academic, Dordrecht, 2002.
 (21) P. N. Vabishchevich, On the form of the hydrodynamics equations, in: WestEast High Speed Flow Field conference, 1922, November 2007, Moscow, Russia, Moscow, 2007, pp. 1–9.
 (22) A. Churbanov, P. Vabishchevich, Numetical methods for solving convectiondiffusion problems, in: C. Zhao (Ed.), Focus on Porous Media Research, Nova Science Publishers, 2013, pp. 1–83.
 (23) A. A. Samarskii, Regularization of difference schemes, Zh. Vychisl. Mat. Mat. Fiz. 7 (1967) 62–93, in Russian.
 (24) P. N. Vabishchevich, Additive OperatorDifference Schemes: Splitting Schemes, de Gruyter, Berlin, 2014.
 (25) J. C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, 2008.
 (26) E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems, Springer Verlag, 2010.
 (27) P. D. Lax, B. Wendroff, Difference schemes for hyperbolic equations with high order of accuracy, Communications on Pure and Applied Mathematics 17 (3) (1964) 381–398.
 (28) V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Springer Verlag, Berlin, 2006.
 (29) G. C. Cohen, HigherOrder Numerical Methods for Transient Wave Equations, Springer, 2002.
 (30) T. A. Zang, On the rotation and skewsymmetric forms for incompressible flow simulations, Applied Numerical Mathematics 7 (1) (1991) 27–40.
 (31) S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer, New York, 2008.
 (32) A. A. Samarskii, A. V. Gulin, Stability of Difference Schemes, Nauka, Moscow, 1973, in Russian.
 (33) N. M. Afanas’eva, P. N. Vabishchevich, M. V. Vasil’eva, Unconditionally stable schemes for convectiondiffusion problems, Russian Mathematics 57 (3) (2013) 1–11.
 (34) P. N. Vabishchevich, Fundamental mode exact schemes for nonstationary problems, arXiv (arxiv.org/abs/1705.07010) (2017) 1–17.
 (35) R. E. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, 1994.
 (36) R. E. Mickens, Nonstandard finite difference schemes for differential equations, The Journal of Difference Equations and Applications 8 (9) (2002) 823–847.
 (37) R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002.
 (38) P. N. Vabishchevich, Twolevel schemes of higher approximation order for timedependent problems with skewsymmetric operators, Computational Mathematics and Mathematical Physics 51 (6) (2011) 1050–1060.
 (39) L. Lusternik, V. Sobolev, Elements of Functional Analysis, Wiley, 1975.
 (40) A. Logg, K. A. Mardal, G. Wells, Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Springer, 2012.
 (41) V. Hernandez, J. E. Roman, V. Vidal, SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems, ACM Transactions on Mathematical Software (TOMS) 31 (3) (2005) 351–362.
 (42) G. W. Stewart, A Krylov–Schur algorithm for large eigenproblems, SIAM Journal on Matrix Analysis and Applications 23 (3) (2001) 601–614.