{'fraction_non_alphanumeric': 0.037055769744627244, 'fraction_numerical': 0.020345841126846585, 'mean_word_length': 4.098973285971486, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Application of Stabilized Explicit Runge-Kutta Methods to the Incompressible Navier-Stokes Equations by means of a Projection Method and a Differential Algebraic Approach

In this master thesis we have compared different second order stabilized explicit Runge-Kutta methods when applied to the incompressible Navier-Stokes equations by means of a projection method and a differential algebraic approach. We explored the stability and accuracy properties of the RKC, ROCK2 and PIROCK schemes when coupled with the projection and the differential algebraic approach. PIROCK has shown unexpected instabilities, ROCK2 resulted to be the most efficient and versatile Runge-Kutta method taken into account. The differential algebraic approach sounds computationally costly but it exhibits better accuracy and a larger stability region. These properties make it more efficient than the projection method. The theory presented in the first chapters is supported by numerical experiments.


Introduction
The aim of this master thesis is to compare the three explicit Runge-Kutta methods RKC, ROCK2 and PIROCK when solving the incompressible Navier-Stokes equations by means of a projection method or a differential algebraic approach. The dimensionless form the Navier-Stokes equations is with Dirichlet boundary condition and initial condition u(t 0 , x) = u 0 (x) x ∈ Ω. (1.1d) Here Ω is a bounded domain in R 2 , t 0 ∈ R, u = (u, v) is the velocity, p is the pressure and ν = 1/Re where Re is the Reynolds number. No boundary and initial conditions exists for the pressure. In fact it is given by a hidden constraint and imposing additional conditions would over-determine the system. The Navier-Stokes equations are the mathematical formulation of many physical phenomena. They may be used to model the weather, ocean currents, water flow in a pipe and air flow around a wing. They help with the design of aircraft and cars, the study of blood flow, the design of power stations, the analysis of pollution, and many other things. Coupled with Maxwell's equations they can be used to model and study magnetohydrodynamics. For this broad range of applications finding effective numerical methods for the approximation of the incompressible Navier-Stokes equations is of significant importance.
The spatial discretization of equation (1.1a) leads to a mildly-stiff system of ordinary differential equations. For this kind of systems it is not worth to use a computationally expensive implicit method. Since the system is not extremely stiff stabilized explicit methods are able to handle the stiffness of the system, in the following we will consider the three methods RKC, ROCK2 and PIROCK.
The above mentioned Runge-Kutta methods are of order two for ordinary differential equations. Moreover they are adaptive in the number of stages s and have a stability region which increases in size quadratically with s. RKC [28] and ROCK2 [3] are explicit, PIROCK [4] is also explicit when applied to the Navier-Stokes equations. The mentioned properties make these methods very attractive for the solution of mildly stiff problems of large dimension, like the one given by equation (1.1a) after spatial discretization.
However one cannot simply integrate equation (1.1a) for two reasons. The first one is that we need the value of the pressure p, the second is that the incompressibility constraint (1.1b) has to be taken into account. It turns out that the pressure is just a Lagrange multiplier [14] and its value is so that the velocity is divergence free. Hence the pressure value and the incompressibility constraint are handled together. In practice two different approaches are used to deal with the incompressibility constraint.
The first one is the use of projection methods. Projection methods have been proposed by Chorin [14] and Temam [30] independently years ago. They use a fractional step approach in which an intermediate velocity (often called virtual velocity) is obtained by solving the momentum equation (1.1a) keeping a constant pressure and disregarding the incompressibility constraint (1.1b). Then the virtual velocity is projected into the manifold of divergence free fields by solving a Poisson equation for the pressure. Due to they decoupled nature projection methods are much more efficient than fully coupled techniques. The price been paid, as we will see below, is that it introduces a numerical boundary layer on the velocity field and the velocity-pressure decoupling adversely affects the temporal accuracy of the numerical scheme. After the work of Chorin and Temam many new projection methods appeared [12,27,23,21,32], in the following we will concentrate on the method given in [33,16]. This is a second order accurate method for the velocity and first order for the pressure. An additional projection for the acceleration gives a second order accurate pressure when needed.
The second approach used to handle the incompressibility constraint is based on the differential algebraic nature of the incompressible Navier-Stokes equations. When discretized in space these equations are differential algebraic of differentiation index 2, the algebraic constraint comes from incompressibility. When Runge-Kutta methods are applied to this kind of equations [18,27] the algebraic constraint has to be satisfied in each stage and new order conditions arise. In the following we will see that for the particular case of the Navier-Stokes equations the new order conditions affect only the pressure accuracy and moreover they can be circumvented. In fact we will not consider the new conditions but instead other methods that allow to bring the pressure to the same order of accuracy as the velocity [27].
The goal of this master thesis is to compare RKC, ROCK2 and PIROCK when coupled with the projection method or the differential algebraic approach and see which one of these couplings is the most efficient. RKC and ROCK2 are well suited for diffusion dominated flows because their stability domain is in a neighborhood of the negative real axis. PIROCK instead is more versatile and behaves very well also for advection dominated flows. Unfortunately his partitioned structure has shown an unexpected instability when the projection method has been applied. We did not had time to investigate this issue, hence we focused on the RKC and ROCK2 methods.
For the description of the algorithm we will restrict our attention to homogeneous Dirichlet boundary conditions and assume that there are no external forces. We will make use of two dimensional test cases, which capture much of the computational difficulties for incompressible flow calculations. We will also assume that the mesh spacing is uniform in the x and y directions. These restrictions are not inherent limitations of the method but they have been adopted here for clarity of exposition.
The present report is organized as follows. In the first section we introduce the stabilized explicit Runge-Kutta methods and the projection method. In the second section we give and introduction to differential algebraic equations, we define the Runge-Kutta methods for semiexplicit differential algebraic equation of index 2 with emphasis on the Navier-Stokes equations and at the end we apply these methods to the stabilized Runge-Kutta methods. The fourth section contains the description of the spatial discretization and the method used for solving the Poisson problem. Finally in section five we show and discuss the results of the numerical experiments.

Stabilized explicit Projection method for the Navier-Stokes equations
The first difficulty in solving the incompressible Navier-Stokes equations arises from the coupling of the velocity and the pressure terms. In the 1960s, Chorin [14] proposed a projection method in which the approximation of the velocity and the pressure is done through a predictor-corrector procedure legitimated by the Helmholtz-Hodge decomposition. In the first step an intermediate velocity field (virtual velocity) is computed by solving the momentum equation (1.1a) ignoring the incompressibility constraint (1.1b) by taking a constant pressure. In the second step the intermediate velocity is projected into the space of divergence free vector fields and the corrector term is used to obtain the new pressure as well. Due to the decoupled approach this method is much more efficient than fully coupled techniques, this advantage makes the projection method attractive. Many improved methods based on this methodology have been seen in the past years, in this chapter we will present two such methods. Another difficulty arises from the spatial discretization of the incompressible Navier-Stokes equations. In the projection method the computation of the intermediate velocity is done integrating equation (1.1a) discretized in space. This semi-discrete system of equations is mildly stiff thus an integrator with a large stability region has to be used, moreover in order to avoid the solution of non-linear systems in must be explicit. Stabilized explicit Runge-Kutta methods are explicit methods with a large stability region along the negative real axis. Thanks to these schemes one handles the stiffness of the equations without solving non-linear systems.
In this chapter firstly we give a short introduction to Runge-Kutta methods, then we present the three stabilized explicit methods RKC, ROCK2 and PIROCK pointing out the major similarities and differences between them. In the second section we introduce two projection methods. The first one will be used in our numerical experiments in chapter 5, the second one is less general and could not be applied in our case. At the end we point out a minor issue of projection methods in general.

Stabilized Explicit Runge-Kutta methods
Historically, the goal when constructing a Runge-Kutta formula was to achieve the highest order possible with a given number of stages. In the case where the constructed method is explicit this approach leads to a relatively small stability domain which permits the integration of non stiff systems of equations only. A stiff problem would demand a severe restriction on the time step which prohibits the use of explicit methods. For this reason implicit methods with large stability domains were employed for the solution of stiff equations, requiring the solution of non Stabilized explicit Runge-Kutta methods are different in the sense that the main goal is to get a stability region which is as large as possible, in a sense that depends on the application. Thanks to this property one can explicitly integrate mildly stiff problems.
In this section we will consider three explicit methods which have a large stability domain along the negative real axis. This kind of schemes are very suitable for mildly stiff problems of large dimension where the eigenvalues are close to the negative real axis, like the ones arising from the spatial discretization of parabolic equations. In the next section we give an introduction to Runge-Kutta methods. In the following we present RKC and ROCK2 which are very similar, then we describe PIROCK, which strongly depends on ROCK2.

Introduction to Runge-Kutta methods
In this section an introduction to Runge-Kutta (RK) methods is given. We will restrict the presentation to the concepts needed to understand the rest of the discussion, for an exhaustive presentation of the subject we refer to [17,18].

Definition of the method
Let us consider the Cauchy problem y ′ (t) = f (t, y(t)) t > t 0 , y(t 0 ) = y 0 , (2.1) for which to we wish numerically approximate the solution. Here y : Let ∆t > 0 be the time step and t n = t 0 + n∆t for n ∈ N. Given y n an approximation of y(t n ) the next approximation y n+1 of y(t n+1 ) is given by the following definition. Setting k i := f (t n + c i ∆t, g i ) an equivalent definition is k i = f (t n + c i ∆t, y n + ∆t s j=1 a ij k j ) i = 1, . . . , s, For general a ij the computation of y n+1 requires the solution of a non-linear system, this should be avoided whenever possible. It can be done setting a ij = 0 for j ≥ i, in this case the RK method is called explicit, otherwise it is called implicit.
In the following two sections we will talk about two very important concepts in the framework of the RK methods, stability and accuracy. After we will explain how the time step ∆t can be chosen dynamically.

Stability
Let us consider the special case where f does not depend on t and is linear in y, i.e. f (t, y) = λy with λ ∈ C and y(t) = y 0 e λt . Applying the RK method it yields y 1 = R(z)y 0 , where z = λ∆t and R(z) is given by a ij and b i . R(z) is called the stability function of the RK method, it is rational if the method is implicit and a polynomial if it is explicit.
Suppose now that Re(λ) < 0, it follows that lim t→∞ y(t) = 0 and one would like to preserve this property when y(t) is numerically approximated by the RK method, i.e. lim n→∞ y n = 0. Since it is necessary that |R(z)| < 1. If |R(z)| ≤ 1 then y n remains bounded and the method is called stable, otherwise y n starts oscillating and explodes as n → ∞. The stability domain S of a RK method is defined by Let λ ∈ C − , for some RK methods λ∆t ∈ S for all ∆t > 0 while for others a condition ∆t < δt(λ) exists, these methods are called conditionally stable. All the explicit methods are conditionally stable. In the special case λ ∈ R − the condition becomes ∆t < l/|λ|, where The concept of stability is extended to non linear functions f in [18,IV.2]. A motivation is to preserve the Lyapunov stability of the system, which has significant importance when the computations are done in finite arithmetic. The parameter λ is replaced by ρ(t), the spectral radius of ∂f ∂y (t, y(t)).

Order conditions
Another important property is the accuracy of the numerical solution, this depends strictly on the order of the method, which is defined by the next definition.
where y(t n ) = y n .
A Runge-Kutta method must satisfy some conditions on a ij , c i , b i so that a certain order p is achieved. The order conditions for p ≤ 2 are given in (2.4).
For p = 2 the order conditions arise from the linearisation of f , which implies that R(z) = 1 + z + 1 2 z 2 + O z 3 is enough to have second order accuracy. The number of order conditions grows rapidly, for p = 6 there is 37 order conditions and for p = 10 they are 1205. For more details on how to find these conditions we refer to [17,II.2].
Also the internal stages g i have an order, it is given by the following definition.

Time step adaptivity with embedded formula
As we saw in the previous sections the size of ∆t has an effect on the stability and the accuracy of the solution. The parameter ρ(t) depends on time thus also the stability condition ρ(t)∆t ∈ S is time dependent. Likewise since the local error depends on the derivatives of f it depends on t, which means that for different ts different ∆ts are needed to achieve a required accuracy. These are the reasons for why it is important to dynamically choose ∆t, in the following we will explain how to do that. Consider two Runge-Kutta methods defined by a ij , b i and a ij ,b i of order p andp respectively, wherep < p. Let ∆t n+1 = t n+1 − t n . We integrate from t n to t n+1 computing the stages k i , which are the same for both methods, and we define the error at step t n+1 as which is the difference between the solutions given by the two RK methods. If this difference increases then probably y n+1 − y(t n+1 ) is increasing as well. The new time step ∆t new is given by and checking ∆t new ≤ l/|ρ(t n )|, where T ol is a user chosen parameter. If err n+1 > T ol the n+1th step is recomputed defining ∆t n+1 = ∆t new , otherwise y n+2 is computed setting ∆t n+2 = ∆t new . In practice there are some safety factors and checks in the implementation. For more details about time step adaptivity we refer to [18,IV.8].

The RKC and ROCK2 methods
RKC (Runge-Kutta-Chebyshev) [28] and ROCK2 (second order Orthogonal-Runge-Kutta-Chebyshev) [3] are two second order stabilized explicit Runge-Kutta methods intended for the time integration of parabolic partial differential equations. They are presented together in this section since they are very similar and share a lot of properties.

A brief description
The optimal stability polynomial of degree s of a second order stabilized explicit method is characterized by with l s as large as possible.
For every s these polynomials exist and are unique [26]. An analytic expression exists in terms of elliptic integrals (see [22]) and the stability region's size is l s = 0.821842s 2 . The practical computation of such polynomials is done numerically. However the realization of these optimal polynomials as Runge-Kutta methods suffer from internal instabilities [5] which could be avoided if one uses a recursion formula [5,31]. Both RKC and ROCK2 make use of Chebychev polynomials and recursion formulas in order to obtain a stability polynomial which is close the the optimal stability polynomials. Because of the recurrence relations the coefficients a ij , b i of the method are not explicitly given, nonetheless they can be computed recursively.
In order to include a strip of non zero width in the stability region a damping parameter η is used. Condition (2.7b) is replaced by |R s (z)| ≤ η for z ∈ [−l s , 0] with l s as large as possible, (2.7c) where η ∈]0, 1[. The smaller η is, the wider the strip will be. On the other hand this will decrease l s thus a compromise has to be found. Usually η = 0.95. One of the most remarkable properties of these methods is that they are adaptive in the number of stages s and the stability domain's size l s increases quadratically with s. Such schemes allow unrestricted integration steps as far as stability is concerned by simply taking s large enough. Thus one has to choose the step size ∆t taking into account only the accuracy requirements. Another good property is that thanks to recursive properties of the stability polynomial they need only a few storage vectors which do not depend on the number of stages s.
But RKC and ROCK2 are not the same method, so they have differences. The first one is that the stability polynomial of RKC is available analytically, while it is computed numerically for ROCK2. The second difference is that the stability bound l s of RKC increases as 0.653s 2 , while it increases as 0.811s 2 for ROCK2. The internal stages of RKC are of order 2, they are of order 1 for ROCK2. The last difference is about time step adaptivity. In ROCK2 the local error is estimated using an embedded formula. In RKC an expression of the local error's leading term is analytically computed, this expression is then numerically approximated in order to obtain the estimation of the local error.
We will now define the RKC and ROCK2 methods.

Definition of the RKC method
Here we will define the RKC method following [18,28]. The stability polynomial of RKC is realized by a three term recursion formula defining the internal stages. It uses scaled and shifted Chebychev polynomials. The Chebychev polynomial T s (x) of degree s is defined by the recursion T 0 (x) = 1, (2.8) and the stability polynomial of RKC is defined as where w 0 = 1 + ε/s 2 and ε ≪ 1 is a parameter defining the damping η = a s + b s ≈ 1 − ε/3. The other parameters w 1 , a s , b s are chosen so that the second order conditions are satisfied, this gives Thanks to the recursive properties (2.8) of Chebychev polynomials in the case of a non linear initial value problem (2.1) the scheme defined by realizes R s (z) defined in (2.9) as stability polynomial. In (2.12) the recursion coefficients are 14)  for j = 2, . . . , s and In Figure 2.1 we show the stability domain and polynomial of RKC with s = 15 stages and a damping η ≈ 0.95.

Definition of the ROCK2 method
The definition of ROCK2 and an exhaustive explanation about its derivation is given in [3], here we will only give its definition.
The construction of the stability polynomial R s (z) of ROCK2 strongly relies on the fact that the optimal stability polynomial can be written as R s (z) = w(z)P s−2 (z), where w(z) is a positive second degree polynomial and P s−2 (z) is a polynomial of degree s − 2 with s − 2 real roots (see [1]). In ROCK2 P s−2 (z) is approximated by a linear combination P s−2 (z) of Chebychev orthogonal polynomials. P s−2 (z) and the approximation w(z) of w(z) are computed numerically by a fixed point algorithm described in [3]. This algorithm computes a set of parameters which recursively define R s (z) = w(z)P s−2 (z), this parameters are stored in table which is used by ROCK2. The recursive scheme is defined by where µ j , ν j , κ j , σ, τ are computed by the fixed point algorithm and depend on s. When applied to the test problem y ′ = λy we have g j = P j (z)y 0 for j = 0, . . . , s − 2 and g s = w(z)P s−2 (z)y 0 . In Figure 2.2 we show the stability domain and polynomial of ROCK2 with s = 15 stages and a damping η = 0.95.

Time step and stages adaptivity
Both RKC and ROCK2 use formula (2.6) in order to choose the new time step ∆t new . Once that we have ∆t new the number of stages s is chosen so that ∆t new ≤ l s /ρ(t). The difference between the methods is in how err n+1 is computed.

RKC
Let us suppose that y n = y(t n ) and we compute y n+1 with the RKC method. The error is defined as The approximation in (2.17) is obtained by a Taylor expansion of y n+1 and y(t n+1 ) around t n (both depend on ∆t.).

ROCK2
ROCK2 uses an embedded formula in order to estimate err n+1 . We have that g * s and g s of (2.16) are first and second order approximations of y n+1 respectively. The error is defined as

The PIROCK method
In this section we will give a short description and define the PIROCK method, for more details we refer to [4].

A brief description
PIROCK is a partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev method for the time integration of diffusion-advection-reaction problems of the kind where F D , F A , F R are the space discretization of the diffusion, advection and reaction operators respectively. For each operator a well suited Runge-Kutta method is chosen and the coupling is done so that order conditions for additive Runge-Kutta methods are achieved. A symmetric diffusion operator has eigenvalues which lies in an interval of size [−O ∆x −2 , 0], where ∆x is the grid size. For this mildly-stiff problem the ROCK2 [3] explicit stabilized second order Chebyshev method is used. The advection operator has imaginary eigenvalues lying in [−iO ∆x −1 , iO ∆x −1 ] thus an explicit third order Runge-Kutta method is used. The reaction term can be very stiff thus an implicit second order integrator is used. This requires the solution of non-linear systems but since the reaction is local these systems are small and can be solved efficiently by LU decomposition and quasi-Newton methods. The case of a nonsymmetric diffusion operator is handled by decomposing F D in its symmetric and asymmetric parts F Ds and F Da respectively. Then F D is replaced by F Ds and F Da is absorbed in F A with some minor modifications in the method. The resulting method is second order accurate and is adaptive in time and in the number of stages s (see section 2.1.2). The number of evaluations of F A and F R is independent of the number of stages s needed to handle the stiffness of F D . Moreover PIROCK is stable even for non symmetric diffusion operators with eigenvalues located in a sector S π/4 = {−ρe iτ : ρ ≥ 0, −π/4 ≤ τ ≤ π/4}. This property also makes PIROCK well suited for advection dominated flows.

Definition of the PIROCK method
PIROCK is defined by the following algorithm for s ≥ 3: Diffusion integration -ROCK2 method

Advection-Reaction integration and coupling
Computation of y 1 ∂y (K s ) and α = 1 for l = 2 or α = 1/(2P ′ s−1 (0)) for l = 1. For l = 2 we obtain the original ROCK2 method. For l = 1 we have a larger damping thus the stability domain will be shorter but with a wider extension on the imaginary axis.

Time step and stages adaptivity
Three embedded methods y e,D , y e,A , y e,R are used for time step adaptivity. We obtain three error estimators The error is estimated as The factor 2/3 is due to the fact that err A = O ∆t 3 while err D and err R are O ∆t 2 . The new time step ∆t new is chosen with formula (2.6) and the new number of stages s and l = 1, 2 are chosen in order to keep the method stable.

Projection method based on stabilized Runge-Kutta methods
As we already said projection methods use a fractional-step technique. Firstly they compute an intermediate virtual velocity and then it is projected into the manifold of divergence-free fields by solving a Poisson equation, the pressure is updated during the projection step as well. The velocity-pressure decoupling of the method makes it cheap but affects the temporal accuracy of the numerical scheme. It has been observed numerically and analytically that second order accuracy for the velocity is readily obtained, while second order accuracy for the pressure is more involved, see [12,33] and chapter 3. Moreover the virtual velocity has no physical meaning and artificial boundary conditions should be used. If the same boundary conditions for the virtual and the physical velocity are used the tangential component will loose accuracy near the boundary. The methods presented in sections 2.2.1 and 2.2.2 do not fix this issue thus the solution is affected by a numerical pollution within a boundary layer. The size of the boundary layer has been studied in [16] for a particular projection method and results to be O ∆t/Re . In section 2.2.3 we propose a method that will be used in section 5 in order to motivate the employment of different boundary conditions for the virtual velocity.
In this chapter we will present two projection methods for the solution of the incompressible Navier-Stokes equations (1.1). The first one is very versatile and can be applied to any second or higher order time integrator. The second one has a limitation which rely on the advection term, in our case we could use it only for the Stokes equation. At the end we will present a fix which will be used in section 5 in order to point out a minor issue very common in projection methods [16]. We stress on the fact that this fix is not applicable to real life problems.

Method 1
In this section we will present the projection method (PM1) which will be used in our numerical experiments of section 5.
This numerical scheme is locally second order accurate for the velocity and first order accurate for the pressure. When necessary an additional projection allows to obtain for the pressure the same order of accuracy as for the velocity. The advantage of PM1 is that the projection step is independent from the integration step thus any second order time integrator can be used. We could apply this method to the three second order schemes RKC, ROCK2 and PIROCK in exactly the same way.
In the following we will define the method, then we do a local error analysis and give a remark about time step adaptivity. Finally we define a variant which is also used in the numerical experiments. In the definition of PM1 and in the local error analysis we follow [33] but PM1 is also present in [16]. Numerical experiments are shown in section 5. In section 5 we refer to this method as PM1.

Definition of the method
Our presentation differs from the one given in [33] in that here equations (1.1) are continuous in space. The following numerical scheme returns a solution u n+1 and p n+1 at time t n+1 given u n and p n at t n .
1. Compute u * n+1 a second order approximation of u * (t n+1 ) the exact solution of for φ 1 .
4. If at time t m a second order approximation for the pressure p m is needed an additional projection for the acceleration is performed.
for φ 2 and update p m = p m + φ 2 .
The motivation of the first projection comes from the constraint ∇ · u = 0 while the second projection arises from the hidden constraint ∇ · ∂u ∂t = 0. With this approach the normal component of the velocity at the boundary satisfies the boundary conditions since But ∂φ 1 ∂τ = 0 thus the tangential component of u n+1 will be affected by a numerical pollution, indeed The theoretical background of the projection u n+1 = u * n+1 − ∇φ 1 is the decomposition theorem of Helmholtz-Hodge [6, 1.16]. The theorem states that a vector field u on a simply connected domain can be uniquely decomposed into a divergence free part u div and an irrotational part u irrot . Thus u = u div + u irrot = u div + ∇φ (2.27) for some potential function φ. Taking the divergence of equation (2.27) yelds ∇ · u = ∇ 2 φ. If the vector field u is known the potential φ is found by solving the Poisson equation and the divergence free part of u can be extracted with the relation u div = u − ∇φ.

Local error analysis
In this section we will give the local error analysis of the above projection method as it is done in [33]. We fix x ∈ Ω and suppose that the exact solution is given at t n . In the following the explicit x-dependence will be left out of notation. Suppose u n = u(t n ) and p n = p(t n ), doing a formal Taylor expansion of equations (1.1a) and (2.21) around t n we get Since in (2.21) a second order accurate method is used we have Summing equations (2.28) and (2.29) we get From equation (2.30) one sees that setting we obtain u n+1 − u(t n+1 ) = O ∆t 3 , i.e. a second order accurate method for the velocity. In order to compute the correction term we take the divergence of equation (2.30) which leads to gives φ 1 = ∆t 2 2 ∂p ∂t (t n ) + O ∆t 3 up to a constant and setting u n+1 = u * n+1 − ∇φ 1 is equivalent to equation (2.31). Doing a Taylor expansion of the pressure p(t) around t n we get thus setting p n+1 = p n + 2 ∆t φ 1 leads to a first order approximation of the pressure. A second projection on the acceleration gives a second order accurate pressure. Supposing that the approximation of the velocity is pth order accurate, u m − u(t m ) = O ∆t p+1 , we can write and taking the divergence of (2.36) gives It follows that setting p m = p m + φ 2 brings the pressure to the same order of accuracy of the velocity. This is done only for output purposes since it is not necessary for maintaining second order velocity. In order to motivate this last assertion let us neglect the advection term, since u * n+1 has second order accuracy we can write and formally Thus the approximation of the pressure p n affects only the second and higher order terms of u n+1 . Because of equation (2.31) this projection method cannot have convergence order higher than 2, even if for computing u * n+1 one uses a third or higher order numerical integrator.

Time step adaptivity
The estimation of the error err n+1 is done before the projection, thus using u * n+1 and not u n+1 . This means that the real local error u n+1 − u(t n+1 ) is overestimated because we are taking into account also the non divergence free part of err n+1 .

Variant
As it is done in [33] we have implemented a variant of PM1 in RKC and ROCK2. It consists in projecting each stage g i during the computation of u * n+1 in (2.21). In section 5 we will refer to this method as PM1V. We will define this variant for the RKC method because the recursive formulas are simpler. In order to apply the variant to ROCK2 one just replaces the formulas.
As we said each stage g i is projected, we call g * i the stage before being projected, φ i the solution of steps 1, 2, 3 in 2.2.1 are replaced by the following algorithm: and finally set 6. u n+1 = g s and p n+1 = p n + 2 ∆t φ s . (2.40f) As you will see in chapter 3 this variant is closer to a Runge-Kutta method for differential algebraic equations than to a projection method. Nonetheless it exists a subtle difference which will be explained in section 3.3.1.

Method 2
The following method (PM2) has been presented in [23] showing good convergence results up to the boundary. In [23] the advection term is discretized explicitly at the beginning of the time step, therefore it can be put into a forcing term. It follows that the analysis can be done on the Stokes equations without loss of generality. However in the case of RKC and ROCK2 the diffusion and advection are advanced together in time, hence the latter cannot be embedded into the forcing term. And in PIROCK neither since we advance the solution in time using the diffusion and then we use these results in order to advance the advection. Thereby the following scheme when applied to RKC, ROCK2 or PIROCK holds only for the Stokes equations. In our numerical experiments the results of this method are very similar to the ones of PM1. In section 5 we will present only the results obtained with PM1.

Definition of the method
The method is very similar to the previous one, the main difference is the consistent pressure update. Consider the equation where q is an approximation of the pressure p which as to be defined. As in section 2.2.1 equation (2.42) is advanced in time from t n to t n+1 and a virtual velocity u * n+1 is computed. Afterward the divergence free velocity u n+1 is given as in section 2.2.1 by the Poisson equation for φ. Plugging the velocity u n+1 = u * n+1 − ∇φ into equation (2.42) gives and comparing equation (2.41a) with equation (2.43) it is clear that p n+1 has to be computed by In [23] it has been shown that setting second order accuracy is achieved for velocity and pressure. Observe that the boundary conditions are not satisfied in the tangential component of the velocity.

Errors at the boundary
As we already said in section 2.2.1 and 2.2.2 the projection methods PM1 and PM2 give a solution which do not satisfy the boundary conditions exactly. We will see in section 5.2 that this issue do not affects the convergence order of the methods but on the other hand we will also see that the errors are coming mostly from the boundary. Inspired by the numerical experiments done with PM1 we propose the following fix (PM3) which needs the derivatives of the exact solution at the boundary. Since it requires an additional information about the exact solution it can be used only to prove, numerically, that a better choice for the virtual velocity boundary conditions could decrease the error. In sections 2.2.1 and 2.2.2 we saw that the tangential component of the velocity do not satisfies the boundary conditions since we set

Equation (2.47) suggests the following boundary condition for the virtual velocity
One could approximate ∂ ∂n τ · u b n+1 by a finite difference formula using the inner points of the mesh but a problem arises when it is embedded in the spatial discretization of the diffusion. At this point the truncation errors become non negligible since they are divided by ∆x 2 , where ∆x is the grid size. We tried to make an approximation with a third order finite difference scheme but the boundary conditions are no more satisfied because of truncation errors. Moreover the errors propagate inside the domain as the solution advance in time.
Nevertheless we can use this fix to show that the boundary conditions for u * impacts the solution accuracy. Numerical experiments with PM1 and PM3 will be discussed in section 5.

Stabilized explicit Runge-Kutta methods for differential algebraic equations with application to the Navier-Stokes equations
Due to the incompressibility constraint (1.1b) the Navier-Stokes equations (1.1) are not ordinary differential equations (ODE) and one cannot simply apply a Runge-Kutta method for ODEs. Instead they are part of a class of equations called differential algebraic equations (DAE) of index 2. In this chapter we will consider Runge-Kutta methods for DAEs of index 2 with focus on the Navier-Stokes equations [27]. We will see that for general DAEs of index 2 new order conditions arise for the differential and algebraic variables (velocity and pressure) but for the specific case of Navier-Stokes most of them vanish. Moreover one can use workarounds in order to obtain high accuracy without satisfying the new order conditions.
In the first section of this chapter we follow [18, VII.6] giving a short introduction to differential algebraic equations and we show that the semi-discrete form of the Navier-Stokes equations has differentiation index 2. In the second section we define Runge-Kutta methods for DAEs of index 2 and apply them to the Navier-Stokes equations as it is done in [27]. At the end we will apply the methods defined in the second section to the stabilized Runge-Kutta methods RKC and ROCK2.

Differential algebraic equations and the Navier-Stokes differentiation index
In this section we will first give an introduction to differential algebraic equations and then prove that the semi-discrete form of the incompressible Navier-Stokes equations has differentiation index 2.

Introduction to differential algebraic equations
A differential algebraic equation is a system of differential equations The difference from an ordinary differential equation is that usually the system is not solvable for all the components of u ′ .

INTRODUCTION TO DIFFERENTIAL ALGEBRAIC EQUATIONS
If the system can be rewritten in the form it is called a semi-explicit DAE. The y variable is called differential, z is called algebraic. In the following we will restrict ourselves to semi-explicit DAEs.
Assuming that ∂g ∂z is invertible in a neighborhood of the solution we can transform the above DAE in an ODE by differentiating equation (3.2b) with respect to t. Indeed which gives The following definition (taken from [18, VII.1]) generalizes this idea.
such that equations (3.5) allow us to extract by algebraic manipulations an explicit ODE u ′ = ϕ(u) (which is called the underlying ODE).
Let us take a look at systems of index 1 and 2. Systems of Index 1. From (3.3) it follows that in the case where ∂g ∂z is invertible in the neighborhood of the solution the semi-explicit DAE (3.2) is of index 1. Systems of Index 2. Let us consider the system Differentiating one time the constraint (3.6b) gives is a hidden constraint of the system. Differentiating (3.6b) twice gives an expression for z ′ is found under the assumption that ∂g ∂y ∂f ∂z is invertible in a neighborhood of the solution. In this case the DAE system (3.6) has index 2.

Semi-discrete Incompressible Navier-Stokes equations differentiation index
The purpose of this section is to show that the spatially discretized incompressible Navier-Stokes equations have differentiation index 2. In the DAE context the velocity u is the differential variable y and the pressure p is the algebraic variable z. Let us discretize equations (1.1) in space, for this we follow [27]. Most of the spatial discretization techniques will give a semi-discrete problem which can be written where u(t) ∈ R Nu , p(t) ∈ R Np are the unknowns vectors for the velocities and pressure respectively, in the remainder their t-dependence will be left out of notation. M , C, D and G represent the discrete divergence, advection, diffusion and gradient operators respectively, Ω is a diagonal invertible matrix. r 1 is a vector with boundary conditions for the continuity equation, r 2 contains the boundary conditions and forcing terms of the momentum equation. Equation (3.8) is a DAE where the differential variable y is u and the algebraic variable z is p. Identifying (3.6) with (3.8) gives and Ω −1 has been absorbed in the definition of C, D, G and r 2 . In section 3.1.1 we saw that if is invertible in a neighborhood of the solution then the DAE system has index 2. The matrix L = M G is the discretized Laplacian operator. Let us derive the expression for p ′ . Deriving equation (3.9b) with respect to time we get and inserting u ′ = f (u, p, t) into (3.11) we obtain the hidden constraint for the pressure Differentiating (3.12) gives Because of the Dirichlet boundary conditions used in the matrices M and G any solution q ∈ R Np of Lq = f will satisfy Neumann boundary conditions, where f ∈ R Np . It follows that L is singular since any solution q plus a constant stills a solution. If one changes the last row of L with ones the matrix becomes regular since the mean value of q is imposed. Using this regularized matrix one can solve (3.13) and find an ODE for the algebraic variable p. Observe that in the case of a time varying mesh in equation (3.12) appears an additional term depending on M ′ , nonetheless this does not change the following result.
With two differentiations with respect to time we found an expression for p ′ , hence the differentiation index of the incompressible Navier-Stokes equations is 2.

Runge-Kutta methods for differential algebraic equations of index 2
We saw that given a semi-explicit differential algebraic equation of index 2 one can compute the underlying ODE u ′ = ϕ(u) (see definition 3.1.1). Hence we could solve the system (3.6) using a common ODE solver applied to u ′ = ϕ(u). This is the bad approach. The process used to obtain the underlying ODE is called index reduction by differentiation. Differentiating the constraint (3.6b) effectively lowers the index of the system but upon discretization lower index systems satisfy only the derived constraints as dg dt (y, t) = 0 and d 2 g dt 2 (y, t) = 0, not the original constraint g(y, t) = 0. This phenomenon produced by the local error is called drift-off, see [18,VII.2] for more details. In the following we will consider Runge-Kutta methods that satisfy the constraint g(y, t) = 0, the solutions will satisfy all the derived lower index systems as well.
We showed in section 3.1.2 that the incompressible Navier-Stokes equations are DAE of index 2. In this section we will firstly define the Runge-Kutta methods for DAEs of index 2 and apply them to the Navier-Stokes equations. Then we will talk about the new order conditions arising from the algebraic constraint and show how they can be circumvented. At the end a short remark about time step adaptivity is given.

Definition of the method
An explicit Runge-Kutta method for a semi-explicit index 2 system (3.6) is defined by where t i = t n + ∆tc i . In order to apply method (3.14) to the Navier-Stokes equations (1.1) we consider their space discretization, hence the functions f and g of (3.14) are given by (3.9). Applying method (3.14) to equations (3.9) gives where P i , U i and u n are approximations to p(t i ), u(t i ) and u(t n ) respectively and F j = F (U j , t j ). In order to simplify the notation we set a s+1,j = b j and c s+1 = 1, then method (3.15) can be rewritten as For solving this system we eliminate the pressure from equation (3.16a), this is done using the hidden constraint (3.12). Replacing the pressure p from (3.12) into (3.16a) gives If the internal stage order of the Runge-Kutta method is q we have that which suggests to insert the exact integral (3.18) into (3.17). Finally (3.17) becomes We have M (I − GL −1 M ) = 0 (even if a row of L has been modified, see [33]), thus supposing that and φ 1 = 0 we can write in such a way φ s+1 is a first order approximation of the pressure. Indeed, let U i be the stages of (3.17) but with u n replaced by where q is the stages order. It follows Let us set the vectors φ = (φ 2 , . . . , φ s+1 ) and ψ = (ψ 2 , . . . , ψ s+1 ) defined as and set u n+1 = U s+1 . Observe that at each stage the solution of a Poisson problem has to be computed in (3.24b).

Order conditions
The classical order conditions for ODEs are not enough to guarantee the expected order of accuracy of both the velocity and the pressure. This is because here f depends also on p and there is the algebraic constraint, clearly these two facts were not taken into account when developing order conditions for ODEs. The following theorem taken from [9] gives the order of convergence of the global error in function of the local error, allowing us to concentrate on local error only.
and that the local error δu(t) satisfies the method is convergent of order p, i.e.
The vector φ used in equation (3.23) depends only on u n and not on p n , thus p n+1 = ψ s+1 depends only on u n as well. Consequently the global error is p n − p(t n ) = O (∆t q ) provided that the local error is δp(t) = O (∆t q ) and u n − u(t n ) = O (∆t q ). This is because the pressure has an instantaneous character, its value is such that the velocity is divergence free and it is independent of the pressure at previous time. In the following we will concentrate on the order conditions for the local error.
As for the ODEs the order conditions are found by Taylor expansion of the exact and the numerical solution and comparing the coefficients of the differentials. For index 2 DAEs one additional condition for the velocity and two conditions for the pressure appear for order two. For order three there is four new conditions for the velocity and four for the pressure. As you can see the number of new conditions grows rapidly. Hopefully in the case of semi-explicit index 2 DAE some differentials vanish making some conditions trivially satisfied. For the even more specific case of the incompressible Navier-Stokes equations there is no additional order conditions for the velocity at least up to order five (see [27]). For the pressure all the conditions remain, for order two they are where ω ij are the coefficients ofÃ −1 .

Circumventing order conditions
The method presented in 3.2.1 needs to satisfy a considerable number of order conditions for the pressure so that the desired accuracy is achieved. If a Runge-Kutta method is built from scratch especially for the Navier-Stokes equations these conditions are taken into account and the method has the desired order. But if one wants to apply an existing Runge-Kutta method to the Navier-Stokes equations probably it has to modify the coefficients a ij so that the conditions for the pressure are satisfied. Moreover it is very likely that at least one stage must be added, therefore one more Poisson problem has to be solved. For these reasons in the following two sections we will show two methods taken from [27] which give higher order pressures without satisfying the additional order conditions. In the second method the stage order of the methods limits the accuracy of the pressure, in the last section we show how this can be avoided. These methods differ from the one in 3.2.1 only in the way the pressure is computed, the velocity is exactly the same and is given by algorithm (3.24).

Approach 1
This method (AP1) integrates the velocity as described in (3.24) without taking care of order conditions. If the Runge-Kutta method has order p then the velocity will converge with order p since no additional order conditions appear. Given the pth order accurate velocity u n+1 a pth order accurate pressure is computed solving Indeed, using (3.12) This approach is very simple but on the other hand one more Poisson problem must be solved. Moreover it requires that r 1 can be differentiated, something that is not required in the computation of u and φ. In some practical computations, for example involving a prescribed turbulent inflow, r ′ 1 might not be available. The approach presented in the next section does not need the derivative of r 1 .
In section 5 we will refer to this method with AP1.

Approach 2
In section 3.2.1 we saw that, under the assumptions that order conditions are satisfied, a higher order pressure is given by a linear combination of the pressure-like variables φ i (see (3.23)). Here we will give a different linear combination of the pressure-like variables which gives a rth order pressure assuming that the Runge-Kutta method has at least r stages and the order of these is at least r. This approach (AP2) is more involved than Approach 1 but it does not require any additional Poisson problem nor the availability of r ′ 1 . Let us denote the average of (3.29) For notation purposes we set φ 1 = 0 and φ i := φ(t i ) for i = 2, . . . , s + 1, yielding The goal is to find an accurate point value p n+1 from the average values φ i , this process is called reconstruction.
Let P (t) be the primitive function of p(t). We construct H(t) the polynomial interpolating P (t) at the points t k 1 , . . . , t km , where the t k j are distinct, k 1 = 1 and k m = s + 1. Then the derivative h(t) of H(t) has the same integral of p(t), indeed (3.31) Using the Lagrange basis and We set However the average values φ i are not known in practice and must be approximated. For a pth order method with qth order internal stages we have, using (3.12) and U i as in section 3.2.1, Finally, to obtain a rth order accurate pressure the Runge-Kutta method must have r stages and they must have rth order accuracy at least. In section 5 we will refer to this method with AP2.

Adapting Approach 2 to order one internal stages
In the previous section we saw that one can obtain a rth order approximation of the pressure by using rth order approximations of the φ i s. It follows that if only first order approximations of the φ i s are available one cannot obtain a second order accurate pressure. In this section we will show how to obtain a second order approximation to φ j using three first order approximations Consequently a second order pressure can be obtained even if only first order approximations to the φ i s are available.
Even if in the following we restrict ourselves to r = 2 with some effort this procedure can be generalized to any r. Using ). Let us consider equation (3.19) and write explicitly his dependence on ∆t, we have Using (3.38) we want to compute the error's leading term of φ i − φ i . Remember that we have (see (3.36) and (3.20))

39)
Let us estimate Lφ i − Lφ i . We call u(t) the solution of the Navier-Stokes equations but with initial condition u(t n ) = u n , it gives Observe that if the internal stages are of order two, i.e.
Remember that if we have two second order approximations of φ i we can compute a second order approximation of the pressure. For i = s + 1 we have φ s+1 − φ s+1 = O ∆t 2 because of the second order conditions of the method. It rests to find one more second order approximation to one of the φ i s, for this we will use (3.43).

Time step adaptivity
In section 2.2.1 we saw that for the projection method the local error is overestimated because also its non divergence-free component is taken into account. Oppositely if a Runge-Kutta method for DAEs is used each stage is projected into the manifold of divergence-free fields.
Then the local error computed with an embedded formula uses only divergence-free stages and a better estimation is given.

Stabilized Runge-Kutta methods for the Navier-Stokes equations: differential algebraic approach
Applying the Runge-Kutta method for differential algebraic equations (see section 3.2.1) to RKC or ROCK2 requiring that the order conditions are satisfied is not straightforward. Firstly because the coefficients a ij are not explicitly given and, more important, they depend on the number of stages s. One can recursively compute them but the new order conditions depends on w ij , which is s dependent as well. Second, all the a ij coefficients for i > j are non zero thus one stage must be added in order to gain degrees of freedom and satisfy the additional conditions. Finally, even if the conditions are satisfied one Poisson problem has been added (because of the new stage) and at this point Approach 1 is better since one does not have to modify the method. Approach 1 can be applied successfully to RKC and ROCK2 without any difficulty when r ′ 1 is available. Approach 2 can be applied to RKC obtaining a second order pressure. This is not true for ROCK2 since it has internal stages of order one preventing higher order pressures, thus the workaround for order one internal stages explained in section 3.2.3 has to be used. In the following we will explain how a Runge-Kutta method for differential algebraic equations is applied to a method which is realized by means of recursive formulas like RKC and ROCK2. Then we show how Approach 1 and 2 are applied to RKC and ROCK2. Finally we observe that in RKC time step adaptivity cannot be enabled when the solution is projected at each stage as in algorithm (3.24).

Recursive formulas and Runge-Kutta methods for differential algebraic equations of index 2
In this section we will see how a Runge-Kutta method for differential algebraic equations given in (3.24) can be applied to a stabilized Runge-Kutta method which uses recursive formulas instead of the coefficients a ij , b i . Again we will explain the method for the RKC scheme for simplicity, for ROCK2 one just replaces the recursive formulas. We write down the recursive formulas (2.12) of RKC when applied to Navier-Stokes without taking into account the pressure (as in (3.24)): for i = 3 for example we have i.e. a 3,1 = µ 2 κ 1 − κ 2 a 1 and a 3,2 = κ 2 . Consider now the method given in (3.24) where the following predictor-corrector procedure is used: where φ i is given in (3.24b). A naïve application to RKC would be and for i = 3, . . . , s + 1 do but (3.51d) is not equivalent to (3.50a), in fact for i = 3 we obtain which is wrong since the last term µ 2 c 2 ∆tφ 2 does not appear in (3.50a). This additional term comes out because in (3.51d) we used U i instead of U * i . The right realization of of RKC as a Runge-Kutta method for DAEs is given by the following algorithm: Observe that the evaluation of F is done on the projected stage U i , otherwise the non projected stage U * i is used. Formulation (3.52) is equivalent to (3.24). Algorithm (2.40) given in section 2.2.1 is equivalent to algorithm (3.51) but written with a different notation. The only difference is the different pressure update but the results do not change. We confirm this assertion by the numerical experiments of section 5.2.3.
In section 5 we will see that algorithm (2.40) (i.e. (3.51)) behave very similarly to (3.52). But when we use Approach 2 of section 3.2.3 (also its variant for order one internal stages) it is very important that the realization of the Runge-Kutta method with recursive formulas is consistent with its realization using the a ij coefficients.

Application to RKC and ROCK2
In order to advance the solution in time one does not have to compute the pressure at each time step so approaches 1 and 2 are used only when one computes the pressure.
It is straightforward to apply Approach 1 to RKC and ROCK2 by simply using algorithm (3.52) and solving the last Poisson problem (3.28).
In order to apply Approach 2 to RKC we need at least two stages of order two, in our implementation we use U s and U s+1 . Since U 1 and U 2 have stage order less than two, when the second order pressure is computed the minimal number of stages used in (3.52) is s = 3, so that U 1 , U 2 , U s and U s+1 are distinct. When the pressure is not computed s = 2 is allowed.
ROCK2 has order one internal stages thus to apply Approach 2 we must use the method explained at the end of section 3.2.3. In order to do that we need to compute the coefficients a ij , c i , b i of ROCK2. Since the method is explicit it follows that a ij = 0 if i ≤ j. From (2.16) it is not difficult to find the other coefficients recursively, we get and c i = s j=1 a ij . The workaround given in 3.2.3 for order one internal stages needs three stages U i , U j , U k of order one in order to compute a second order approximation of φ j with j = s + 1. In our implementation we use U 2 , U 3 , U 4 and we get a second order approximation of φ 3 . Since ROCK2 is a second order method the last stage U s+1 has order two, so we can use φ s+1 as second order approximation of φ s+1 . Hence we need the stages U 2 , U 3 , U 4 for φ 3 and U s+1 for φ s+1 . The stages U 4 and U s+1 can coincide meaning that the minimal number of stages is s = 3, as in the original ROCK2 method.
When using Approach 2 in RKC the computation of the pressure can affect also the velocity since the minimal number of stages changes, this happens only when very small ∆ts are used and it is not true for ROCK2.

Time step adaptivity issue in RKC
We saw in section 2.1.2 that in RKC the estimation of the local error err n+1 is done by an approximation of y ′′′ (t n ). This approximation is computed under the assumption that the equation being solved is an ODE, so there is no dependence on p(t). In the projection method explained in section 2.2.1 the pressure is kept constant during the integration, in such a way y ′′′ (t n ) does not depends on p(t) and the approximation holds. On the other hand when the differential algebraic approach of section 3.2.1 is used the velocity u is projected at each stage. This is equivalent to advance the pressure p in time since its value is so that the velocity is divergence free. So, the approximation of y ′′′ (t n ) given in (2.17) is not valid anymore. One can give another approximation taking into account the dependence on p(t) but this is not done here. This discussion is valid also for the variant of the projection method given in 2.2.1.
Observe that an embedded formula for time step adaptivity works also when the velocity is projected at each stage. This is because the theory of embedded formulas uses only the order conditions of the velocity and these are the same for ODEs and the Navier-Stokes equations, hence projecting the stages or not does not affects the correctness of the time step adaptivity procedure.

Chapter 4 Space discretization and Poisson solver
In this chapter we will describe the method used for the spatial discretization and explain how the Poisson problem is solved in our implementation. In chapter 3 we made use of matrices for the discretized differential operators in order to describe the theory, however these matrices are not explicitly built in the code. Instead the operators are implemented as functions which given the operand return the matrix multiplication. The boundary conditions and forcing terms are embedded in these functions as well. Consequently the Poisson problem is not solved by means of matrices neither.
For the spatial discretization we use a grid where the unknowns are staggered and in order to solve the Poisson problem we use a fast cosine transform well suited for this particular grid.

The MAC method
The Marker and Cell (MAC) method has been proposed in 1965 in the Los Alamos laboratory [19]. Recently, it has been shown in [24] that the MAC method is very competitive and particularly well suited for high Reynolds numbers and free surface problems.
In this chapter we will motivate the usage of the MAC method and apply it to the incompressible Navier-Stokes equations.

Motivating the use of the MAC method's staggered grid
The MAC method uses a staggered grid, in this arrangement the scalar variables (pressure, density,...) are stored in the cell centers and the momentum variables (velocity) are located at the cell faces. This is different from a collocated grid, where all the variables are stored at the same position. A staggered grid avoids the odd-even decoupling of pressure and velocity (see equation (4.2)). Odd-even decoupling introduces a discretization error which leads to checkerboard patterns in the solution.
Let us motivate the previous considerations with an example. Suppose that Ω = [0, 1]×[0, 1] is divided into a number of uniform cells as in Figure 4.1 and the equation to be solved is ∂u ∂t In a collocation grid the discretization of equation (4.1) with u = (u, v) ⊤ and f = (f 1 , f 2 ) ⊤ leads to

APPLICATION OF THE MAC METHOD TO THE NAVIER-STOKES EQUATIONS
where ∆x is the grid size. With this choice of positioning four independent pressure modes arises, as shown in Figure 4.1. Hence it is easy to build a functionp non constant but for which the discretized gradient vanish, makingp invisible to the momentum equation (4.1). In the staggered grid the pressure unknowns are all coupled together.

Application of the MAC method to the Navier-Stokes equations
In the presentation of the MAC method we will focus on the incompressible Navier-Stokes equations with Dirichlet boundary conditions. As in a collocation grid, in a staggered grid finite differences are employed but in this approach the velocity components u, v and the pressure p are not discretized at the same points of the domain.
and after discretization in space on a staggered grid by centered finite differences we obtain where v i,j−1/2 is computed by a simple mean without affecting the second order convergence of the overall scheme. Looking at Figure 4.2 we see that all the pressure unknowns are coupled together and only one pressure mode is allowed by this positioning of the unknowns. Near a boundary centered finite differences cannot be used since the grid spacing is not uniform, see Figure 4.2. A non centered finite difference scheme can be easily found using Taylor series, for x ∈ R and a function f ∈ C 4 ([x − ∆x/2, x + 2∆x]) we found the schemes given be the relations where · ∞ is the maximum norm on the interval [x − ∆x/2, x + 2∆x]. When we apply (4.5a) and (4.5b) to (4.3) at a point x i,1/2 we need the value u i,0 which is given by the Dirichlet boundary conditions. For the second component of equation (1.1a) we obtain a similar expression. We have thus obtained a system of 2N (N − 1) ODEs.
The equation (1.1b) is discretized at the pressure points x i−1/2,j−1/2 obtaining For both the projection method in section 2.2 and the differential algebraic approach in section 3.2 we need to solve a Poisson problem with Neumann boundary conditions (see for example step 2 of section 2.2.1 and equation (3.20)) discretized at the pressure points x i−1/2,j−1/2 . We decided to solve the Poisson problem with the discrete cosine transform. This is a direct method which automatically imposes Neumann boundary conditions.

Solving the Poisson problem by means of the discrete cosine transform
In this chapter we will present the method used in order to solve the Poisson equation on Ω = [0, 1] × [0, 1] with Neumann boundary conditions. We decided to solve the Poisson equation using a discrete cosine transform (DCT). The DCT is a special case of the discrete Fourier transform (DFT) which automatically imposes Neumann boundary conditions. Moreover, like the DFT, the differential equation can be easily solved by algebraic manipulations. Many variants of the DCT exists in the literature (see [11] for a survey) and the choice depends on the type of discretization of the domain. In the case of staggered grids the most appropriate DCT is called DCT-II, it will be presented later on.
In the following we will derive the DFT and its inverse from the Fourier transforms. Then the DCT-II and its inverse (called DCT-III) will be derived from the DFT. After we define

THE DISCRETE FOURIER TRANSFORM
the two dimensional DCT (2d-DCT) and explain how the Poisson problem can be solved with the 2d-DCT. Finally two fast algorithms for the computation of DCT-II and DCT-III will be derived and compared by numerical experiments.
For simplicity in the following we will refer to DCT-II with DCT and to DCT-III with IDCT.

The discrete Fourier transform
There is more than one way to obtain the discrete Fourier transform. In this section we will derive it from the Fourier transforms as in [10]. But it can be obtained also from Fourier series, trigonometric polynomials and using the delta function (see [10]).

Deriving the DFT from the Fourier transform
Let f ∈ L 1 (R), the Fourier transform of f is defined bŷ and if |f | ∈ L 1 (R) the inverse is defined by In the following we will suppose that f andf are sufficiently regular. For more details about the Fourier transform we refer to [29]. ω is called the frequency variable while x is called the spatial variable. Equation We wish approximate this integral numerically. Let N ∈ N be even, we define the N + 1 grid points x n = n∆x, where ∆x = A/N and n = −N/2, . . . , N/2. Using the composite trapezoidal rule we getf where The approximation (4.11) can be evaluated at any frequency ω but we would like to approximate also the integral in (4.8) so must we choose a discretization for the frequency domain. Imagine a function f which modes have an integer number of periods in the domain [−A/2, A/2], so that it fits perfectly in the interval. Consider the wave with one period, it has wave-length A and frequency 1/A units. The other waves will have frequency k/A for k ∈ N. It follows that the frequency unit has to be chosen as ∆ω = 1/A and the grid points ω k = k∆ω. In order to uniquely determinef (ω k ) from f (x n ) and vice-versa we choose also N points ω k in the frequency domain. The length of the latter is Ω = N ∆ω = N/A and ω k is defined for k = −N/2 + 1, . . . ., N/2. The following defines the DCT and its inverse IDFT.
Approximations tof (ω k ) are given by AF k . The IDFT is another sequence of N complex numbers given by It remains to show that the DFT and the IDFT are really inverse one of the other, i.e. D −1 (D(f n ) k ) n = f n and D(D −1 (F k ) n ) k = F k . This is easy to verify using the identity (4.14). (4.14)

The discrete cosine transform
The DCT is derived from the DFT using a sequence f n of sampled points which is even and real. Different sequences give rise of different DCTs, here the sequence f n is sampled on a staggered grid so that the arising DCT is appropriate for this kind of spatial discretization. Firstly we derive the one dimensional DCT, then we define the two dimensional DCT.
We define (g n ) 2N n=−2N +1 ∈ R as the sequence of 4N real numbers given by if n ≥ 0 and n is odd, f −n/2 if n ≤ 0 and n is odd.
Applying the DFT given by (4.12) to the above sequence we obtain kn .
Using g n = g −n and g n = 0 for n even it follows Inserting f n+1/2 in equation (4.15) the coefficients F k of the DCT are given by It is easy to see that F k = F −k and F k = −F 2N −k . Thanks to these symmetries the inverse DCT is given by As expected g n = 0 if n is even. For odd indexed g n we have the coefficients f n+1/2 of the IDCT: One sees that a sequence f n+1/2 given by (4.17) is symmetric in f 0 and f N , i.e. f −1/2 = f 1/2 and f N −1/2 = f N +1/2 . It is very common to merge together the scaling factors, the final form of the DCT is given by the following definition.

Definition 4.2.2.
Given a sequence f n+1/2 for n = 0, . . . , N − 1 the DCT coefficients are given by for k = 0, . . . , N − 1. The IDCT is given by for n = 0, . . . , N − 1. Here ′ is a sum in which the first term is weighted by one half.

The two dimensional DCT
With the same techniques used in 4.2.1 and 4.2.2 one can derive the two dimensional DFT (2d-DFT) from the two dimensional Fourier transform and the two dimensional DCT (2d-DCT) from the 2d-DFT. Here we will skip the derivation of the 2d-DCT giving only his definition.  for m, n = 0, . . . , N − 1.

Solving the Poisson problem
In this section we show how the Poisson problem with Neumann boundary conditions can be solved using the 2d-DCT. Consider the Poisson problem Let u m+1/2,n+1/2 and f m+1/2,n+1/2 be the u and f functions sampled at the staggered points x m+1/2,n+1/2 for m, n = 0, . . . , N − 1. We note the 2d-DCT coefficients of u and f by U j,k and F j,k respectively. The five points Laplacian in x m+1/2,n+1/2 gives Making a comparison term by term gives the equation for U j,k : The procedure for solving the Poisson equation with the 2d-DCT is: Notice that a solution of equation (4.22) for j = k = 0 fails to exists unless ∆x 2 F 0,0 = 0. Hopefully it is the case, we have Consequently it is valid to suppose ∆x 2 F 0,0 = 0. In this way any value of U 0,0 will satisfy equation (4.22), we set U 0,0 = 0 which means i.e. setting U 0,0 = 0 imposes the mean value of u.

Fast algorithms for the discrete cosine transform
In this section we will explain how the 2d-DCT is computed in our implementation. The 2d-DCT in (4.20) can be written Hence the 2d-DCT reduces to 2N DCTs. In the following we will present 2 algorithms for the computation of the DCT.

Iterative algorithm
This algorithm has been taken from [13]. It is an iterative algorithm for the computations of the DCT and the IDCT. Its computational complexity is O N 2 but the operations involved are mostly add and multiply. So the algorithm is much faster than a naïve computation of (4.18) and (4.19), which needs a lot of cosine computations.
DCT Let f n ∈ R for n = 0, . . . , N − 1 be a real sequence with N even. We set w k n = f n + (−1) k f N −1−n , equation (4.18) can be written Using recursive properties of sine and cosine we can compute G j (k) and H j (k) in a recursive way as well. Setting θ k = kπ/N we have IDCT When computing the IDCT we pre multiply F 0 by one half, so that we can write where Q j (n) = j k=0 F 2(j−k)+1 sin (2n+1) 2N (2k + 1)π . Using (4.24a) and (4.24b) we have As for G j , H j also for P j , Q j it exists a recursive definition. Setting θ n = (2n + 1)π/N we can show that

Recursive algorithm
The following algorithm has been presented in [20]. It is recursive and thanks to this property its computational complexity is O (N log 2 (N )). DCT Let f n ∈ R for n = 0, . . . , N − 1 be a real sequence with N = 2 m , m ∈ N. We set for n = 0, . . . , N/2 − 1. The DCT is given by

From (4.25)(c,d) we have
Multiplying the right hand side of (4.26) by cos 2m+1 N kπ and summing it over k gives doing the same on the left hand side we obtain (4.27) and in a very similar fashion we obtain Using (4.27) and (4.28) we have In this way the N point IDCT has been decomposed in two N/2 IDCTs. The DCT and IDCT recursion stops when N = 2, at this point the transform is computed using (4.18) or (4.19).

Comparison of the two algorithms
In

Numerical experiments
In this chapter we will compare the three schemes RKC, ROCK2 and PIROCK when coupled with the different approaches for solving the incompressible Navier-Stokes equations explained in chapters 2 and 3. Each numerical experiment is executed several times but using different methods or different test problems, for this reason we will describe each experiment once for all and tell which methods are used before showing the results. After we will consider three test problems, for each one of these we show and comment the results of the previously described experiments.

The codes used and description of the numerical experiments
In this section we will firstly explain how we name the different couplings between the Runge-Kutta methods and the different approaches for solving the Navier-Stokes equations, we also give a few remarks about these couplings. Then we will describe the numerical experiments and tell which methods are used to perform the experiment.

Names for the different methods and remarks
For each numerical method we give a name composed by three parts. The first part is the name of the numerical integrator, it can be RKC, ROCK2 or PIROCK. The second component tells which approach has been used to solve the incompressible Navier-Stokes equations. It is PM1 if the projection method of section 2.2.1 is used, PM1V if we use the variant of PM1, i.e. the method explained at the end of section 2.2.1, and PM3 if the method of section 2.2.3 is used. If we use approach 1 or 2 of section 3.2.3 the name will be AP1 or AP2 respectively, if the workaround of Approach 2 for order one stages is used it will be AP2W. The third component of the name tells if the second order pressure is computed at each time step or only in the last one, it will be CP1 or CP0 respectively (to be read Compute Pressure true or false).
For example the name 'ROCK2, AP1, CP0' means that we use ROCK2 coupled with AP1 and the pressure is computed only in the last time step. Now we give a few remarks about the possible couplings. After the first experiments PIROCK, PM1, CP0 has shown stability issues which we did not had time to investigate, for this reason we did not implemented PM1V, AP1 and AP2 in PIROCK. If RKC is used time step adaptivity can be enabled only when coupled with PM1, otherwise the local error estimation gives wrong results, the reason is given in 3.3.3. If PM1 or PM1V are used a first order pressure is computed after each time step, disregarding about CP0 or CP1. If AP1, AP2 or AP2W are used with CP0 no pressure is computed until the last time step. Observe that AP1, AP2 and AP2W in general give the same results for the velocity. The only exception is when in RKC a small time step is used and the pressure is computed, in this case the minimal allowed number of stages differs between AP1 and AP2, giving different velocities. Remember that with ROCK2 the approach AP2 cannot be used and we use AP2W, with RKC we use AP2. CP0 or CP1 do not change the way the pressure is computed but just the frequency, the method for the pressure computation is given by PM1, PM1V, AP1, AP2 or AP2W.

Convergence order
The convergence order test is done to check that the effective convergence order matches the theoretical one. In numerical experiments the time and space errors are mixed, in order to split these errors the following method is used.
Following [23] the numerical solution of a temporally pth order and spatially rth order accurate scheme can be expressed as where η n ij and η(x i , y j , t n ) are the numerical and analytical solutions respectively. β n ij (∆t) p and γ n ij (∆x) r are the errors corresponding to the time and spatial discretization and ǫ is the round-off error. Thanks to expression (5.1) we can carry out time and space convergence tests separately. For time convergence we compute a numerical reference solutionη n ij using a very small time step so that we can writẽ η n ij ≈ η(x i , y j , t n ) + γ n ij (∆x) r + ǫ, then various numerical solutions η n ij are computed with different time steps on the same grid giving the relation which gives the temporal convergence order. Similarly for the spatial convergence we compute a reference solutionη n ij with a very fine grid, yieldinĝ η n ij ≈ η(x i , y j , t n ) + β n ij (∆t) p + ǫ. Computing various numerical solutions η n ij with a fixed time step and different grid size we obtain which gives the spatial convergence order.
In These convergence tests are not done for CP1 but only for CP0. Also, the results of AP1, AP2 and AP2W are put in the same figure since the velocities are the same, only the pressure differs. We will now list the methods used in the convergence tests: • RKC with PM1, PM1V, AP1, AP2 and CP0, • ROCK2 with PM1, PM1V, AP1, AP2W and CP0, • PIROCK with PM1 and CP0.

Stability tests
The stability tests are done to investigate the size of the stability domain of the different methods. Two different tests are done, in order to confirm the generality of the results. As we already said PIROCK has stability issues which did not permit to perform this test, hence we did it for RKC and ROCK2.
The first test is done by taking a set of stage numbers and searching for the maximal stable time step when solving the Navier-Stokes equations with ∆x = 1/128 and Re = 5. The goal is to check if the stability domain size's is the same as the one for ODEs. We know that in the ODE context the stability domain grows as 0.653s 2 and 0.811s 2 along the negative real axis for RKC and ROCK2 respectively, where s it the number of stages. So the theoretical maximal stable time step is given by ∆t = 0.653s 2 /ρ and ∆t = 0.811s 2 /ρ for RKC and ROCK2 respectively, where ρ is the spectral radius of the diffusion operator. In the code we estimate ρ using Gershgorin's theorem.
The first test is done for Re = 5, in order to show that the results do not depend on Re we did a second test. For this test we neglected the advection term since we are interested on the Reynolds number effect. Fixing the time step size at ∆t = 10 −2 we take a set of Reynolds numbers and for each one of these we search for the minimal number of stages s which gives a stable solution. The theoretical minimal s is given by s = ∆tρ/0.653 and s = ∆tρ/0.811 for RKC and ROCK2 respectively, where ρ depends on Re.
In both the tests we integrate from t = 0 to t = 1. CP0 and CP1 do not affect the result thus these tests are done for CP0. Also, AP1, AP2 and AP2W give the same results since the velocity is the same for both. We will again list the methods used in these test: • RKC with PM1, PM1V, AP1, AP2 and CP0, • ROCK2 with PM1, PM1V, AP1, AP2W and CP0.

Errors at the boundaries, codes profiling, accuracy improvement by projections
This experiment has two goals. The first one is to motivate the employment of different boundary conditions for the virtual velocity u * of PM1. The second goal is to show that projecting the velocity after each stage gives a much higher accuracy being just a bit slower.
Using the ROCK2 scheme we will compare the virtual and physical velocities u * and u given by PM1, PM1V and PM3 (method of section 2.2.3) in the first two time steps when ∆t = 10 −1 .
The tests are done without the advection term thus the PIROCK scheme gives the same results of ROCK2 and the ones given by RKC are very similar. After we will integrate until the end of the interval using PM1, PM1V, AP1, AP2W showing the code profile and the velocity and pressure errors. In particular we will see if the time spent in projecting each stage (PM1V, AP1, AP2W) is worthwhile compared to the increased accuracy with respect to PM1. Finally the codes used for this test are • ROCK2 with PM1, PM1V, PM3, AP1, AP2W and CP0.

Numerical efficiency
These tests compare the different methods in order to find out which one gives the most accurate solution in less time. Following [18,II.10] we have compared the different methods with time step adaptivity enabled and different relative and absolute tolerances where m = 2, 3, . . . , 12. The tests are done on an Intel®Core ™ i7-2600S processor. The Fortran code has been compiled with the GNU compiler and the -Ofast option. We compare the solutions given by the different tolerances against a reference solution. We used a 128 × 128 grid, Re = 100, the starting time is t = 0 and the end time is t = 1, the starting time step for ROCK2 and PIROCK is ∆t = 10 −3 , for RKC it is automatically chosen. The reference solution has been computed with the fourth order RK4 method (see [17,II.1]) using ∆t = 10 −7 and compensated summation. In the figures the symbols represent the different tolerances.
For RKC we did the test only for PM1 since for the others methods we cannot enable time step adaptivity. For PIROCK also we did the test only for PM1 since the others are not implemented. For ROCK2 all the possible couplings are tested. Moreover this test is done for CP0 and CP1. To summarize the methods used are: • RKC with PM1 and CP0, CP1, • ROCK2 with PM1, PM1V, AP1, AP2W and CP0, CP1, • PIROCK with PM1 and CP0, CP1.

Different Reynolds numbers behavior
In this test we are interested in how the different methods behave with different Reynolds numbers. Using a fixed ∆x = 1/128 and time step adaptivity we consider different Reynolds numbers. For each Reynolds number we integrate from t = 0 to t = 1 and look at the following quantities: the velocity errors against the exact solution, the computational time, the average number of stages per time step, the total number of stages used, the number of time steps and the number of rejected time steps.
PIROCK has stability issues thus this test is done for RKC and ROCK2 only. Again for RKC we use only PM1 because of time step adaptivity. For large Reynolds numbers RKC and ROCK2 showed instabilities arising from the advection term, for this reason it has been neglected. The test is done under the same environment of the efficiency tests (same processor, compiler and options, see 5.1.5). The methods used are: • RKC with PM1 and CP0, • ROCK2 with PM1, PM1V, AP1, AP2W and CP0.

Comparing our results with an established reference
As a last test we compare our results with a well established reference in order to check the correctness of the implementation. We have compared the solution of the methods that allow time step adaptivity (except CP1) with tabular results given in [15] at the stationary point. The methods used are: • RKC with PM1 and CP0, • ROCK2 with PM1, PM1V, AP1, AP2W and CP0, • PIROCK with PM1 and CP0.
In the following sections we consider three test problems and show the results of the numerical experiments above described.

Convergence order
In Figures 5.1 and 5.2 we see the convergence results of the methods for the forced flow problem. The results of the spatial convergence are visually the same for all the methods and are summarized in Figures 5.1(a,b), we see that second order space convergence is achieved for all the methods for both the velocity and the pressure. Figures 5.2 show the results of the time convergence tests. In Figure 5.2(a) we see that PIROCK, PM1, CP0 achieve the expected time order of convergence for all the quantities but we do not have results for ∆t ≥ 10 −2 , this is because the method is not stable. In Figures 5.2(b-g) we show the results of RKC and ROCK2 for PM1, PM1V, AP1, AP2 and AP2W. We observe that in general the methods that project the velocity after each stage (PM1V, AP1, AP2, AP2W) have more regular results. Looking at Figures 5.2(f, g) we remark that the pressure computed with AP1 converges faster than the one computed with AP2 and AP2W. For RKC, PM1V, CP0 the pressure p 1 has second order of convergence, even if only first order was expected.

Space convergence order
(a) All the PM1 and PM1V methods.

Space convergence order
All the AP1, AP2 and AP2W methods.  Time convergence order Time convergence order Time convergence order Time convergence order Time convergence order Time convergence order Time convergence order

Stability tests
In Figure 5.3 we present the results of the stability test done on the forced flow problem. From Figures 5.3(a,b) we see that the methods that project the velocity after each stage (PM1V, AP1, AP2 and AP2W) maintain the same stability domain of the ODEs. On the other hand the stability domain grows only as 0.54s 2 and 0.67s 2 for RKC and ROCK2 respectively when we use PM1, so when we do just one projection per time step. Figures 5.3(c,d) show that these properties does not depend on the Reynolds number.

Errors at the boundaries, codes profiling, accuracy improvement by projections
In this test we show the error of the tangential component of the velocity for different methods, show the codes profile and give a deeper analysis on the errors. In Figure 5.4 u * n is the virtual velocity at time step t n of ROCK2, PM1 and ROCK2, PM3. g * s is the last stage of ROCK2, PM1V before being projected (see (2.40)). u n is the physical velocity and u e n is the exact velocity given by equations (5.2). The quantities in Figure 5 In Figure 5.4(a) we display the difference between the virtual velocity and the exact solution after one time step for PM1 and PM3. For PM1V we show the difference between the last stage g * s before projection and the exact solution. We plot the correction term ∂φ 1 ∂x as well for each method PM1, PM1V and PM3. For PM1 the imposition of the boundary conditions for u * nullify the error at the boundary y = 0 and y = 1. Moreover we see that the boundary conditions affect also the interior of the domain since several stages are computed without corrections. For PM1V the error also vanish in y = 0 and y = 1 but since the solution is projected after each stage there is no propagation inside the domain. For PM3 the error at the boundary satisfies ∂u * 1 −u e 1 ∂y = 0 which corresponds to equations (2.48). We see that PM1V naturally satisfies equations (2.48), without being imposed (it is the case for PM3). Observe that the correction term ∂φ 1 ∂x for PM1 and PM3 is much closer to the error u * 1 − u e of PM3 while for PM1 there is a larger gap between ∂φ 1 ∂x and u * 1 − u e near the boundary. For PM1V the error and the correction term have an excellent match. In Figure 5.4(b) we see the error |u 1 − u e 1 | between the physical velocity u 1 = u * 1 − ∂φ 1 ∂x and the exact solution. We see that the error of the physical velocity for PM1 is much larger at the boundary than the one of PM1V and PM3. In order to compare the errors of PM3 and PM1V we plot them in Figure 5.4(c) using a different scale. We see that the error of PM1V is much smaller than the one of PM3 and that PM3 still have a boundary layer, something that is not present at all for PM1V since the maximal error is inside the domain. Figure 5.4(d) shows the errors in the second time step, we see that the large error at the boundary of PM1 has propagated inside the domain in the second time step. The same is not true for PM1V and PM3. Comparing virtual velocities error  Analyzing the first two time steps we saw that changing the boundary conditions of PM1 (thus using PM3) can increase the accuracy of the solution. But more important we saw that projecting the velocity after each stage (so using PM1V) creates an excellent match between the last stage g * s and the correction term, as a consequence the solution is much more accurate. We did not plot the internal stages g * i for i < s but we suppose that they also match the correction term as well as g * s . We have not done this test for AP1, AP2 and AP2W but we are confident that the results are the same of PM1V, this is confirmed by the rest of this section where we will analyze the overhead of the extra projections in PM1V, AP1, AP2W and see if it is worthwhile compared to the increased accuracy.
In Table 5.1 we show part of the PM1, PM1V, AP1, AP2W codes profile when integrating from t = 0 to t = 1 again with ∆t = 10 −1 , so doing then time steps. We would like to stress on the fact that the timings given in Table 5.1 are adversely affected by the profiler tool (Score-P), in reality the code is much faster (about 4-5 times). What is interesting here is to see the relative difference between the methods, this quantity is not affected by the profiler. The first line of Table 5.1 shows that PM1 is about 1.3 times faster than the other methods. PM1 is the method which spends more time in computing the right hand side f , this is because it requires 15 stages per step while the other methods just 13 since they have a longer stability domain. On the other hand its time spent in projecting the velocity in negligible. For the other methods it is considerably high, looking at the third line we see that it takes about one third of the total time. The time spent in computing the pressure is negligible for all the methods since we do it in last time step only for PM1, PM1V and AP1. For AP2W it comes for free since no Poisson problems must be solved for the pressure.  In Table 5.2 we write the errors of the methods when compared against the reference solution computed with the RK4 method (see 5.1.5) at t = 1. We see that PM1V, AP1, AP2W give exactly the same velocity errors and PM1V and AP1 give the same pressure error, in fact they compute it in the same way.

PM1
PM1V AP1 AP2W Velocity error 1.88 · 10 −1 2.92 · 10 −4 2.92 · 10 −4 2.92 · 10 −4 Pressure error 2.22 3.5 · 10 −1 3.5 · 10 −1 3.19 · 10 −2 Table 5.2: Errors of the PM1, PM1V, AP1, AP2W methods with ∆t = 10 −1 at t = 1. Table 5.3 shows the ratio between the errors of PM1 and the errors of the other methods. In Figure 5.4 we saw that after the first time step PM1V gives an error that is much smaller than the error of PM1. With this table we want to quantify how much it is smaller, taking into account also AP1 and AP2W. We see in the first line that for the velocity PM1V, AP1, AP2W are 644 times more accurate than PM1. PM1V and AP1 are about 6 times more accurate than PM1 for the pressure, AP2W is about 70 times more accurate. PM1/PM1V PM1/AP1 PM1/AP2W Velocity error ratio 644× 644× 644× Pressure error ratio 6.34× 6.34× 69.6× Table 5.3: Ratio between the PM1 and other methods errors.
We saw that projecting the velocity after each stage makes the methods 1.3 times slower but on the other hand they are 644 times more accurate. Hence the additional projections are absolutely worthwhile.

Numerical efficiency
In Figure 5.5 we compare the numerical efficiency of the methods for the forced flow problem. In Figures 5.5(a,c) we see that among the PM1 and PM1V methods with CP0 the best one is ROCK2, PM1V while PIROCK is the less efficient. The efficiency lack in PIROCK can be associated to its stability issues. We also see that depending on the tolerance RKC, PM1, CP0 can beat ROCK2, PM1, CP0, especially for the pressure (see Figure 5.5(c)) but its accuracy increases slowly. When computing the second order pressure at each time step (see Figures  5.5(b,d)) ROCK2, PM1V is again the best method and ROCK2, PM1, CP1 is better than RKC, PM1, CP1 for all the tolerances. In Figure 5.5(e-h) we show the efficiency of ROCK2 with AP1 and AP2W. We see that for the largest tolerance there is an efficiency decrease. We think that it is not really an efficiency decreasing but instead it is the reference solution that is not accurate enough. This is even more visible for the pressure of AP2W since the pressure of the reference solution has been computed with AP1. However comparing Figures 5.5(e,f) we see that switching from CP0 to CP1 does not change the computation time of AP2W while for AP1 there is a small increase because of the additional Poisson problem for the pressure. Looking at the pressure ( Figures 5.5(g,h)) we see that AP1 is more efficient than AP2W, specially for CP0. For CP1 the difference is smaller, it would be interesting to see what happens with a more accurate reference solution. When using CP0 AP1 is better than AP2W, when using CP1 it depends if the user is interested in the velocity or pressure accuracy.
As we already saw in section 5.2.3 projecting the velocity at each stage (PM1V, AP1, AP2, AP2W) is more efficient than projecting only once per time step. The overhead of the extra projections is compensated by the increased accuracy (see 5.2.3), the better estimation of the local error (see 3.2.4) and the larger stability domain (see 5.2.2). Yet, in this example a lot of computational time is spent computing the right hand side f (see Table 5.1) because of the forcing terms, thus the time spent in projections is small in percentage. This means that in a problem without forcing term the overhead of the extra projections is more important and it could cause a decay in efficiency when projecting at each stage, we will investigate that in section 5.3.2. The projection time depends strongly on the grid size ∆x, it would be good to see what happens when a smaller grid size is used but we did not had time to do that. We are particularly interested in see how the results in Figures 5.5(f,h) would change. It is also interesting to see that except for RKC, PM1, CP0 the computational time increases only as the square root of the required accuracy, while it increases linearly for RKC, PM1, CP0. Switching from CP0 to CP1 completely changes the behavior of RKC, PM1 (see Figures 5.5(a-d)), it is much more sensible to the pressure accuracy than the other methods. In fact, comparing the first column of figures in 5.5 with the second column we remark that for ROCK2 and PIROCK computing a second order pressure after each time step slightly increases the computational time without improving accuracy at the final time step at t = 1 (where the errors are measured). Hence, for ROCK2 and PIROCK one should use CP1 only when accurate pressures are needed at intermediate time steps. While for RKC, PM1 there's a large difference between CP0 and CP1.
With a side by side comparison ROCK2, AP1 is the best method when CP0 is used. For CP1 the choice is between ROCK2, AP1 and ROCK2, AP2W, it depends if the user is interested on velocity or pressure accuracy. If it is interest in velocity it choses AP2W because it is faster, if it is interested in pressure it choses AP1 because of its accuracy. ROCK2, PM1V performs almost as well as ROCK2, AP1. Except the pressure update the difference between these methods is the recurrence formulation. PM1V uses algorithm (3.51) while AP1 uses (3.52) (see section 3.3.1). For this reason we prefer AP1 since its realization with recursive formulas is consistent

Different Reynolds numbers behavior
In this test we compare the different methods when changing the Reynolds number, the results are in Figure 5.6. Figure 5.6(a) shows the error of the different methods against the exact solution, we see that again projecting the velocity after each stage (PM1V, AP1, AP2W) is much more accurate than PM1. Figure 5.6(b) shows that RKC, PM1 is not only badly accurate but also the most computationally expensive method. ROCK2, PM1 is less accurate than ROCK2 with PM1V, AP1, AP2W but at least it is faster. However we saw in 5.2.3 and 5.2.4 that this velocity is not worthwhile. We see that the computational time of ROCK2, PM1V is higher than ROCK2 with AP1 and AP2W. Figure 5.6(c) shows the average number of stages used in each time step, we see that RKC, PM1 uses much less stages per time step but on the other hand in Figure 5.6(d) we see that in total it uses more stages than the other methods. In fact it has a smaller average ∆t, see Figure 5.6(e) which displays the total number of time steps. From this figure we see that PM1 and PM1V require more steps than AP1 and AP2W without gains in accuracy (see Figure 5.6(a)), proving that algorithm (3.52) is better than (3.51). Figure  5.6(f) shows that RKC has a considerable number of rejected time steps while ROCK2 do not have any rejected step.

(5.3)
This test has been taken from [27], here no forcing terms are needed. In the following we will present the convergence results and compare the numerical efficiency of the different methods.

Convergence order
In Figures 5.7 and 5.8 we see the space and time convergence results for the Green-Taylor vortex. It is shown in Figure 5.7 that velocity and pressure achieves second order convergence in space for all the methods. For this problem we had to compute a reference solution with ∆x = 1/1024 in order to see a good convergence. In Figure 5.8(a) we see that PIROCK, PM1, CP0 achieve the expected order of convergence for all the quantities but again it is not stable for ∆t ≥ 10 −2 . In Figures 5.8(b-g) we show the results of RKC and ROCK2 for PM1, PM1V, AP1, AP2 and AP2W. All of them have the expected order of convergence and as for the forced flow problem convergence results (see section 5.2.1) RKC, PM1V, CP0 exhibits second order of convergence for the pressure p 1 even if first order was expected. Looking at Figures 5.8(f, g) we notice that the pressure computed with AP1 is much more related to the velocity accuracy than the one computed with AP2 or AP2W. We believe that the reason is that using AP2 or AP2W the pressure is computed using the average values φ i s, hence large but very local errors are damped. Using AP1 a point value (the velocity) is used, so no damping occurs.

Numerical efficiency
In Figure 5.9 we compare the numerical efficiency of the methods for the Green-Taylor vortex. From Figures 5.9(a-d) we see that again among the PM1, PM1V methods ROCK2, PM1V is

Space convergence order
(a) PIROCK, PM1, CP0. Space convergence order Space convergence order Space convergence order Space convergence order       Time convergence order

Space convergence order
(a) PIROCK, PM1, CP0. Time convergence order Time convergence order Time convergence order Time convergence order (e) ROCK2, PM1V, CP0. Time convergence order (f) RKC, AP2 and AP2, CP0. Time convergence order (g) ROCK2, AP1 and AP2W, CP0. the best one. RKC, PM1, CP0 is better than ROCK2, PM1, CP0 only for a few tolerances, less than in the case of the forced flow. The gain in efficiency when using PM1V is less than in the forced flow problem (see Figures 5.5(a-d)) because in this example there is no forcing term, thus the overhead of the extra projections is more significant. Nonetheless PM1V remains more efficient than PM1. Comparing Figures 5.9(a,c) with Figures 5.9(b,d) we remark that again RKC has a much better performance when a second order pressure is computed after each step while ROCK2 and PIROCK are less sensible to this. Figures 5.9(g,h) show that again AP1 is more reliable than AP2W for the pressure but on the other hand Figure 5.9(f) shows that AP2W is faster than AP1 when the pressure is computed after each time step. Comparing the first column of figures in 5.9 with the second column we can conclude that for ROCK2 and PIROCK one should use CP1 only when accurate pressures are needed at intermediate time steps because at t = 1 there is no gain in accuracy using CP1 instead of CP0, as for the forced flow.
Again the best methods are ROCK2, AP1 and ROCK2, AP2W. The choice depends on the preference between the velocity and the pressure accuracy.

The Lid-driven cavity problem
The lid-driven cavity problem is considered as the classical test case for the validation of new Navier-Stokes simulation codes and algorithms. Its geometry and boundary conditions are very simple. It consists of a square domain with Dirichlet boundary conditions on all sides. Three sides are stationary while the last one is moving. The mathematical formulation is equations It is known [7] that this problem presents singularities at the corners (0, 1), (1,1) where the velocity is discontinuous. This property makes difficult to properly evaluate the accuracy of the numerical results, mainly in the neighborhood of these points where pressure and vorticity are not finite. Since the space accuracy is associated to the smoothness and boundedness of the of the solution derivatives it is then completely lost at these corners. This is the reason why one would use the regularized driven cavity problem [8] in order to evaluate the solution accuracy. In this regularized problem the velocity is smoothed, so that the above problems are avoided. On the other hand these difficulties makes the driven cavity problem a widely used test case for the evaluation of incompressible flow solvers.
In the following we will show the time convergence results and compare our results with tabular data given in [15].

Convergence order
As we said the solution is non smooth at the corners and it follows that the expected space convergence order is not achieved. However time convergence is not affected by these singularities. In Figure 5.11 we show the time convergence of the lid-driven cavity problem for the methods listed in section 5.1.2. We see that all the methods have the expected convergence order. In this problem not only PIROCK but also RKC and ROCK2 are not stable for ∆t ≥ 10 −2 . From Figures 5.11(f,g) we see that RKC is more stable than ROCK2, we suspect that it is because its stability region covers a larger area in the complex plane. Looking at Figures 5.11(d,f) we see that one more time algorithm (3.52) performs better than (3.51), in this case from the stability point of view.

Comparing our results with an established reference
Since this problem has been solved many times there is a great deal of data to compare with in the literature. An established reference is given by Ghia et al. in [15] since it contains tabular results for various Reynolds numbers. We have compared the results given by the methods that allow time step adaptivity at the stationary point t = 500 with the data given in [15] for Re = 1000. Like in [15] we did the simulation with a 128 × 128 grid. In our simulation we used an adaptive time step with rtol = atol = 10 −3 .
In the Figures the velocity u is taken at the vertical centerline of the cavity, the velocity v at the horizontal centerline. We did the test for the methods listed in section 5.1.7 and the results look the same, so they are summarized in Figure 5.10 under the label 'Our solution'. We see a very good matching between our solutions on the ones obtained by Ghia et al. in [15]

Conclusions
In the previous chapters we have presented the stabilized explicit Runge-Kutta methods RKC, ROCK2 and PIROCK and two approaches for handling the Navier-Stokes equations incompressibility constraint: the projection method and the differential algebraic approach. Then we described the method used for the spatial discretization and for the solution of the Poisson equation on this particular grid. At the end we showed the numerical results with discussions. Already in the first chapters we noticed that RKC has a limitation concerning time step adaptivity which is not present in ROCK2 and PIROCK. This issue can be solved easily but it was not our purpose. Due to this limit when using RKC we could enable time step adaptivity only with the projection method.
The main results obtained with numerical experiments are the following. We showed that all the methods work and attain the expected second order of convergence, in space and time.
The stability tests showed that the projection method has a stability domain shorter than the one obtained for ordinary differential equations. On the contrary the differential algebraic approach has exactly the same stability domain as the ordinary differential equations case. We did also a very close comparison of the projection method and the differential algebraic approach using ROCK2. It turned out that with the used time step the projection method has a non negligible boundary layer produced by the boundary conditions for the virtual velocity. And, more important, the boundary layer is not present at all when the differential algebraic approach is used. Using a fixed time step we saw that the differential algebraic approach was 1.3 times slower but on the other hand the solution was 644 times more accurate! In fact from the numerical efficiency tests it turned out that the differential algebraic approach is much more efficient than the projection method.
The comparison between RKC, ROCK2 and PIROCK showed that ROCK2 is the most efficient and PIROCK is braked by his stability issues. We showed that the accuracy of the pressure does not affects the velocity, so accurate pressures should be computed only if needed. When comparing the methods with different Reynolds numbers again ROCK2 is way more efficient than RKC and the differential algebraic approach is much better than the projection method.
For the differential algebraic approach two different methods have been used for the pressure computation: one uses a Poisson problem the other an interpolation with Lagrange polynomials. The latter is limited by the stages' order of the method but we showed how this limit can be circumvented. It turned out that the method which uses the Poisson problem is slightly more expensive but also more accurate and reliable. The one which uses the interpolation computes the pressure with a negligible cost. So the choice between them depends on how much the user is interested in the pressure accuracy. Anyway both of them give a second order accurate pressure.
As a last experiment we compared our results with an established reference and we found a very good match.
The differential algebraic approach resulted to be highly superior than the projection method: it his more accurate, it has a longer stability region and it permits to have a better estimation of the local error. Moreover the projection method cannot provide solutions of order higher than two. Differently, the differential algebraic approach accuracy depends only on the order of the Runge-Kutta method. It follows that it can be applied to higher order stabilized explicit methods, as ROCK4 [2], in a straightforward way.
In future we would like to compare the different methods with a finer grid, where the Poisson problems are more expensive. This might change the difference in efficiency between the projection method and the differential algebraic approach. Also the difference between the two methods for computing the pressure might change.
The main limit of RKC and ROCK2 is that their stability region is in the neighborhood of the negative real axis, hence they are unstable when the Navier-Stokes equations are advection dominated. In the future our main goal is to apply the differential algebraic approach to PIROCK (it has not been done here). We expect that it behaves better than the projection method. We saw that the differential algebraic approach has better stability properties, the same of ordinary differential equations, so it could be that the stability problems are also fixed. We expect that PIROCK behaves similarly to ROCK2 for diffusion dominated flows. For advection dominated flows we expect from it a better behavior due to its large stability domain on the complex plane.