'''Abstract:''' Among the main problems in the field of Computational Fluid Dynamics (CFD), there is the two-dimensional laminar flow in a transient regime of an incompressible fluid modeled by Navier-Stokes equations. Among the decoupled solutions for this equation system, that is, solutions for velocity regardless to pressure, there are the Projection methods, which separate the solution in three parts, solved in each time steps applying an iterative process regardless to the iterative process for the other variables. It may result, according to the Projection method applied, in two Reaction-Diffusion equations and one Poisson equation to be solved in each time step. This paper sought to develop an algorithm to solve the Navier-Stokes equation, applying the Finite Volume Method (FVM) with second order approximation scheme (CDS), beside a Projection method with incremental pressure-correction scheme, so that each Reaction-Diffusion and the Poisson equation are solved efficiently. Therefore, several solvers were tested for each equation, resulting in an algorithm with the combination that achieved the best result for each equation, with the preconditioned Conjugate Gradient method (PCG) with the Multigrid method (MG) and ILU solver (Incomplete LU factorization) being the methodology used in the whole problem solving process. The geometric Multigrid with V cycle, the correction scheme (CS), the full weighting restriction, the prolongation through bilinear interpolation and the maximum number of levels for the studied cases were utilized. The results achieved were satisfactory, since the proposed methodology accelerated the iterative process considerably in relation to the classical methods available in the literature.

'''Keywords:''' Navier-Stokes, Multigrid, Conjugated Gradient, Projection methods, ILU factorization.

=1 Introduction=

In computational mathematics, determining the solution for linear or nonlinear equation systems is an important problem, since they often occur in the discretization process of mathematical models in Computational Fluid Dynamics (CFD).  In order to achieve such solutions, iterative methods are widely applied. They can be categorized as linear, nonlinear, stationary and nonstationary. These approaches may result in difficulties related to the slow convergence of the iterative process applied. [1]. In order to improve the convergence of the iterative methods, the Multigrid method has been proving very efficient [2-6]. Its philosophy is based on the application of several grids with several degrees of refinement, which are gone through during the iterative process. Nevertheless, the conditioning of the matrix for the linear system's coefficients may slow the iterative method's convergence. It may be accelerated by modifying this coefficient matrix in order to improve its conditioning. This process is called matrix preconditioning, [7]. Among the preconditioning technics, one can mention the incomplete LU factorization (ILU), which consists in changing a matrix ''A ''into a product of a lower triangular matrix ''L ''and an upper triangular matrix ''U'', attached to a matrix called residual ''R''.

=2 Mathematical model and solution method=

The mathematical model in two dimensional Cartesian coordinates for the laminar flow of an incompressible fluid, isothermal and in a transient regime, given by the Navier-Stokes equations [16], in this paper, is represented by the continuity equation that depicts the physical principle of the conservation of mass, and by the equations of dimensionless movement in the directions ''x'' and ''y '':

 

where ''x ''and ''y'' are the coordinates; ''t'' the temporal coordinate, ''u ''and ''v'' the velocity components of the vector '''u''' in the directions ''x ''and ''y'', respectively; ''p ''the static fluid pressure and Re'' ''the Reynolds number associated to the viscous and the inertia forces. According to [17], when Re << 1, the viscous forces prevail and when Re >> 1 the inertia forces prevail.

\end{array}</math> ,

where ''T_i'' is the initial time and ''t ''> 0

According to [19], Eqs. (5) and (6) are Helmholtz equations

| <math>{\nabla }^2u=f</math> .

=3 Numerical models=

==3.1 Utilized grid==

In this paper the Finite Volume method (FVM) with staggered grids was applied, as described in [21], also applied by [12], with the velocities placed on the faces and the pressure on the center of the volumes. Applying this strategy avoids the numerical instabilities on the pressure solution [6, 22, 23]. Such grid and its ordering are represented in Fig. 1 (4 to 4 volumes).

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image9.png|246px]] </div>

==3.2 Spatial discretization by finite volumes==

Equations (1) a (3) can be integrated in each control volume described in Fig. 2:

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

| [[Image:Draft_Anunciacao_160657346-image10-c.jpeg|174px]]

| [[Image:Draft_Anunciacao_160657346-image10-c1.jpeg|center|204px]]

| <math>\underset{S}{\oint }u.n\mbox{ }dS</math> ,

| <math>\frac{\partial }{\partial t}\int_Vu=-\underset{S}{\oint }uu.n\mbox{ }dS-</math><math>\underset{S}{\oint }pn_x\mbox{ }dS+\frac{1}{\mbox{Re}}\underset{S}{\oint }\nabla u.n\mbox{ }dS</math>

| <math>\frac{\partial }{\partial t}\int_Vv=-\underset{S}{\oint }vu.n\mbox{ }dS-</math><math>\underset{S}{\oint }pn_y\mbox{ }dS+\frac{1}{\mbox{Re}}\underset{S}{\oint }\nabla v.n\mbox{ }dS\mbox{,}</math>

The terms of Eqs. (8), (9) and (10) were discretized applying  the FVM with staggered grids, where the integrals are calculated over the control volumes defined in Fig. 2, where for a matter of simplicity, the same spatial refinement was applied in all directions, that is, ''h_x = h_y = h ''(uniform grids).

In order to discretized the mass conservation, a volume centered in ''p'' was applied, Fig. 2 (b), where there is:

| <math>\underset{S}{\oint }u.n\mbox{ }dS\approx u_{i+\frac{1}{2},j}^{n+1}-</math><math>u_{i-\frac{1}{2},j}^{n+1}+v_{i,j+\frac{1}{2}}^{n+1}-</math><math>v_{i,j-\frac{1}{2}}^{n+1}=0</math> .

The pressure terms were discretized applying the volumes centered in the velocity components defined in Fig. 2 (a and c), generating:

| <math>\underset{S}{\oint }pn_x\mbox{ }dS\approx \left(p_{i+1,j}-\right. </math><math>\left. p_{i,j}\right)h</math>

| <math>\underset{S}{\oint }pn_y\mbox{ }dS\approx \left(p_{i,j+1}-\right. </math><math>\left. p_{i,j}\right)h</math> .

The same control volumes defined in Fig. 2 (a and c) were applied in the discretization of the diffusive terms, where one obtains:

| <math>\underset{S}{\oint }\nabla u.n\mbox{ }dS\approx u_{i+\frac{3}{2},j}^n+</math><math>u_{i-\frac{1}{2},j}^n+u_{i+\frac{1}{2},j+1}^n+u_{i+\frac{1}{2},j-1}^n-</math><math>4u_{i+\frac{1}{2},j}^n</math>

| <math>\underset{S}{\oint }\nabla v.n\mbox{ }dS\approx v_{i,j+\frac{3}{2}}^n+</math><math>v_{i,j-\frac{1}{2}}^n+v_{i+1,j+\frac{1}{2}}^n+v_{i-1,j+\frac{1}{2}}^n-</math><math>4v_{i,j+\frac{1}{2}}^n</math> .

| <math>\underset{S}{\oint }uu.n\mbox{ }dS\approx \left[{\left(u^2\right)}_{i+1,j}^n+\right. </math><math>\left. {\left(u^2\right)}_{i,j}^n+{\left(uv\right)}_{i+\frac{1}{2},j+\frac{1}{2}}^n-\right. </math><math>\left. {\left(uv\right)}_{i+\frac{1}{2},j-\frac{1}{2}}^n\right]h</math>

| <math>\underset{S}{\oint }vu.n\mbox{ }dS\approx \left[{\left(v^2\right)}_{i,j+1}^n+\right. </math><math>\left. {\left(v^2\right)}_{i,j}^n+{\left(uv\right)}_{i+\frac{1}{2},j+\frac{1}{2}}^n-\right. </math><math>\left. {\left(uv\right)}_{i-\frac{1}{2},j+\frac{1}{2}}^n\right]h</math> .

| <math>{\left(u^2\right)}_{i,j}^n={\left[\frac{1}{2}\left(u_{i+\frac{1}{2},j}^n+u_{i-\frac{1}{2},j}^n\right)\right]}^2</math> ,

| <math>{\left(v^2\right)}_{i,j}^n={\left[\frac{1}{2}\left(v_{i,j+\frac{1}{2}}^n+v_{i,j-\frac{1}{2}}^n\right)\right]}^2</math> .

==3.3 Projection methods==

The Projection methods were introduced by [24,25], and they are characterized by the Navier-Stokes solution in two steps. Initially, an auxiliary velocity '''u''' ''^t''is calculated ignoring the incompressibility condition given by Eq. (1) and the pressure term (considering ''p ''= 0 in Eqs. (2) and (3)). In this step, a Reaction-Diffusion equation is solved by temporal discretization.  In the second step, pressure is calculated by solving a Poisson equation. The most attractive property in the Projection methods is the fact that in each step a sequence of decoupled elliptical equations are solved for pressure and velocities, enabling large scale numerical simulations to be executed efficiently; however, it is not trivial to develop and analyze high order Projection methods [14].

There are several formulations for the Projection methods, differing in the way to calculate the auxiliary velocity '''u''' ''^t''and the projection step. Each of these comes from the method of [24,25], which is a Projection method with a simpler correction scheme for the pressure, which is not yet a first order method (regarding time) for pressure on norm ''l''₂ [14] due to the absence of the pressure gradient on the first step of the method. [26] observed that by adding the pressure gradient, Chorin's algorithm presents an improvement in accuracy, an idea applied by [27] to develop a incremental pressure-correction scheme of second order, which was improved by [28] through a modification on the boundary conditions for pressure from the method mentioned previously. This algorithm is known [29] as “incremental correction for pressure on rotational form” and it will be adopted in this paper for the solution of the Navier-Stokes equations. Such algorithm will be given below and it will be designated as Algorithm 1.

Furthermore, [29] show that the incremental version on rotational form is a second order scheme for velocity on ''l''₂-, ''l₁''- and ''l_∞''-norms. For pressure the scheme is a second order scheme for norms  ''l''₂- and ''l₁''-norms  and 3/2 order for ''l_∞''-norms. The method consists on the following steps:

\frac{3u^t-4u^n+u^{n-1}}{2\nu h_t}=\frac{\beta g}{\nu }+{\nabla }^2u^t+\frac{\nabla P^n}{\nu }\\

\end{array}</math> ,

|}

 

where '''u''' ''^t'' is the auxiliary velocities field, '''u''' ''ⁿ'' is the velocities field on the time step ''n'', ''h_t'' is the temporal refinement, ''Pⁿ'' is pressure on the time step ''n'' and ''bⁿ'' the boundary conditions on time step ''n'',  <math>\beta g=u\cdot \nabla u</math> '' ''with constant  <math>\beta </math> and  <math>\nu =1\mbox{/Re}</math> .

Rearranging the terms on Eq. (21) , there is:

-{\nabla }^2u^t+\frac{3\mbox{Re}}{2h_t}u^t=\mbox{Re}\left(\beta g+\nabla P^n+\frac{4u^n-u^{n-1}}{2h_t}\right)\\

\end{array}</math> ,

which has the form of a Reaction-Diffusion equation (Eq. (6)) for the variable '''u''' ''^t'', with:

| <math>\eta =-\frac{3\mbox{Re}}{2h_t}</math>

| <math>f=\mbox{Re}\left(\beta g+\nabla P^n+\frac{4u^n-u^{n-1}}{2h_t}\right)</math> .

{\nabla }^2{\phi }^{n+1}=\frac{3}{2h_t}\nabla \cdot u^t\\

u^{n+1}=u^t-\frac{2h_t}{3}\nabla {\phi }^{n+1}\\

\end{array}</math> ,

where  <math>{\phi }^{n+1}</math> is the pressure correction.

 

==Algorithm 1: Projection method==

 

 

Estimate initial conditions for ''u'' and ''v'':  <math>u^{-1}=u^0</math>

 

Estimate initial conditions for ''P'': <math>P^0</math>

 

 

For ''n'' = 0, 1, 2, ...,  <math>N_t</math> – 1, do:

 

Step 1: given the velocities  <math>u^n\mbox{, }u^{n-1}</math> and the pressure  <math>P^n</math> , calculate '''u''' ''^t'', solving the Reaction-Diffusion equation:

 

<math>-{\nabla }^2u^t+\frac{3\mbox{Re}}{2h_t}u^t=\mbox{Re}\left(\beta g+\right. </math><math>\left. \nabla P^n+\frac{4u^n-u^{n-1}}{2h_t}\right)</math>

 

 

 

 

\end{array}</math>

End

If the time step ''h_t''is too large in relation to spatial refinement (''h_x'' and ''h_y''), the simulations can become unstable [30]. The criterion applied in this paper to avoid instabilities is given by [12,30]:

| <math>h_t<\frac{1}{2\mbox{Re}}\frac{h_x^2h_y^2}{h_x^2+h_y^2}</math> .

Since this paper considers ''h_x'' = ''h_y'' = ''h'', the criterion can be rewritten as

| <math>h_t<\frac{h^2}{4\mbox{Re}}</math> .

| <math>{\int }_{\Omega }{\nabla }^2{\phi }^{n+1}={\int }_{\Omega }\frac{1}{h_t}\nabla \cdot u^{n+1}=</math><math>0</math> .

=4 Method for linear systems solution=

| <math>Ab=f</math> ,

 

==4.1 Gauss-Seidel method (GS)==

''A'' = ''P'' – ''N'', (30)

| <math>Pb^{\left(k+1\right)}=Nb^{\left(k\right)}+f</math> ,

| <math>b^{\left(k+1\right)}=Bb^{\left(k\right)}+P^{-1}f</math> ,

 

If the diagonal entries of ''A'' are non-zeros, the unknown corresponding value can be isolated in each equation, resulting in the Gauss-Seidel method, given an initial estimate ''b''⁽⁰⁾:

| <math>b_i^{\left(k+1\right)}=\frac{1}{a_{ii}}\left[f_i-\sum_{j=1}^na_{ij}b_j^{\left(k+1\right)}-\right. </math><math>\left. \sum_{j=1+1}^na_{ij}b_j^{\left(k\right)}\right],\mbox{     }i=</math><math>0,...,n</math> ,

==4.2 Preconditioned Conjugated Gradient (PCG)==

Descent methods are based on obtaining a solution ''b'' for a system of the type given by Eq. (29) by the minimization property of an appropriate norm for the error. From an initial estimate ''b''⁽⁰⁾, these methods determine a solution for the problem by diminishing the error.

The vectorial norm being given by  <math>{\Vert b\Vert }_2=\langle b,Sb\rangle </math> , where  <math>S\in R^{nxn}</math> is a positive defined symmetrical matrix. Considering the residual equation ''Ae''⁽''^k''⁾ = ''r''⁽''^k''⁾, it can be deduced that:

| <math>{\Vert e^{\left(k\right)}\Vert }_2=\langle e^{\left(k\right)},Se^{\left(k\right)}\rangle =</math><math>\langle r^{\left(k\right)},Rr^{\left(k\right)}\rangle </math> ,

where  <math>R=A^{-T}SA^{-1}</math> , which means that minimizing the error equals minimizing the residue in the same norm [7].

 

One of the ways to choose the best direction in each point ''b''⁽''^k''⁾ would be to choose the direction of'' p''⁽''^k''⁾ which is lined up with the direction of the function gradient, but in the opposite orientation (Gradients method).

| <math>b^{\left(k+1\right)}=G_kb^{\left(k\right)}+c_k</math> ,

where <math>G_k=I-{\alpha }_kS</math> and  <math>c_k={\alpha }_kSA^{-1}f</math> ,  <math>{\alpha }_k</math> being a constant.

Conjugate Gradient method (CG) is an iterative descent method applied in cases where the iteration matrix ''A ''is a symmetric positive definite matrix. This method demands a more careful choice of the descent directions ''p''⁽''^k''⁾.

In order to exemplify it, consider a two-dimensional problem. The lines <math>{\Vert e^{\left(k\right)}\Vert }_2={\Vert r^{\left(k\right)}\Vert }_2=</math><math>\mbox{constante}</math> , are concentric ellipses with the center in ''b'', as illustrated in Fig.3.

 

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image65-c.png|282px]] </div>

Supposing there is an initial estimate ''b⁽''⁰⁾, according to the Gradient method, one can start from this point and research the minimum of  <math>{\Vert e^{\left(k\right)}\Vert }_2</math> , across the direction ''r''⁽⁰⁾, obtaining ''b''⁽¹⁾. Observing Fig. 3, it can be verified that the best progression of ''b''⁽¹⁾would not be the more inclined direction, but the direction ''p''⁽¹⁾ which points to the center of the ellipse.

In fact, the new minimizing point ''b''⁽²⁾ would coincide with the exact solution of ''b'', and thus the iterative process would end with just two iterations. Therefore a direction ''p''⁽¹⁾ can be chosen as:

| <math>p^{\left(1\right)}=b-b^{\left(1\right)}</math> .

| <math>0=\langle r^{\left(1\right)},r^{\left(0\right)}\rangle =</math><math>\langle r^{\left(1\right)},p^{\left(0\right)}\rangle =</math><math>\langle f-Ab^{\left(1\right)},p^{\left(0\right)}\rangle =</math><math>\langle A(b-b^{\left(1\right)}),p^{\left(0\right)}\rangle </math> ,

| <math>\langle p^{\left(1\right)},Ap^{\left(0\right)}\rangle =</math><math>0</math> .

| <math>p^{\left(k+1\right)}=r^{\left(k+1\right)}+{\omega }_kp_k</math> ,

| <math>{\omega }_k=\frac{\langle r^{\left(k+1\right)},r^{\left(k+1\right)}\rangle }{\langle r^{\left(k\right)},r^{\left(k\right)}\rangle }</math> .

| <math>\Vert e^{\left(k\right)}\Vert \leq \left(\frac{{\left(cond_2A\right)}^{\frac{1}{2}}-1}{{\left(cond_2A\right)}^{\frac{1}{2}}+1}\right)\Vert e^{\left(0\right)}\Vert </math> ,

where  <math>cond_2A={\Vert A\Vert }_2\cdot {\Vert A^{-1}\Vert }_2</math> .

 

One way to make this iterative method faster consists in transforming the system  <math>Ab=f</math> into an equivalent system, where the matrix has a more favorable conditioning number. According to Eq. (41), if it is possible to reduce  <math>cond_2A</math> , consequently a fast convergence is achieved.

 

Considering that in the CG method the transformed matrix has to keep being a symmetric positive definite matrix, the preconditioning will also have to preserve this property. In fact, it can be obtained making ''P'' = ''S'' and ''Q'' = ''P^T'' = ''S^T'' , and in this manner, it can be concluded that  <math>\overline{A}\overline{b}=\overline{f},</math>  <math>\overline{A}=SAS^T</math> ,  <math>\overline{b}=S^{-T}b</math> and  <math>\overline{f}=Sf</math> .

| <math>W={\left(S^TS\right)}^{-1}</math> .

| <math>{\alpha }_0=\langle {\overline{r}}^{\left(0\right)},r^{\left(0\right)}\rangle /\langle p^{\left(0\right)},Ap^{\left(0\right)}\rangle </math> ,

| <math>b^{\left(1\right)}=b^{\left(0\right)}+{\alpha }_0p^{\left(0\right)}</math> ,

| <math>r^{\left(1\right)}=r^{\left(0\right)}-{\alpha }_0Ap^{\left(0\right)}</math> ,

| <math>p^{\left(1\right)}=r^{\left(1\right)}+{\beta }_0p^{\left(0\right)}</math> ,

| <math>W{\overline{r}}^{\left(j+1\right)}=r^{\left(j\right)}</math> ,

==Algorithm 2. Preconditioned Conjugated Gradient method (PCG)==

 

 

 

Estimate ''b''⁽⁰⁾

 

 

 

 

 

 

 

 

 

 

 

 

 

 

End

==4.3 Incomplete LU Decomposition (ILU)==

==4.4 Multigrid method (MG)==

| <math>v^H={\left[I\right]}_{\mbox{ }h}^{\mbox{ }H}v^h</math> ,

| <math>v^h={\left[I\right]}_{\mbox{ }H}^{\mbox{ }h}v^H</math> .

Figure 4 illustrates the cycle-V with CS for five grids. The superscripts ''h'', 2''h'' , 4''h'', ... indicate the grid where the vectors and matrixes are defined, the smoothing numbers (S) being achieved in the restriction process (R) (pre-smoothing) and prolongation (P) (post-smoothing) being represented, respectively, by  <math>{\nu }_1</math> and  <math>{\nu }_2</math> .

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image122-c.jpeg|372px]] </div>

=5. Results=

==5.1 Methodology and computational details==

:1) Multigrid with ILU solver (MG-ILU).

:2) Preconditioned Conjugated Gradient with Multigrid and Gauss-Seidel solver (PCG-MG-GS).

:3) Preconditioned Conjugated Gradient with Multigrid and ILU solver (PCG-MG-ILU).

In all simulations and for all the variables of interest, the stop criterion applied was the nondimensionalized residue ''l₂''-norm based on the initial estimate, given by:

 

In this paper the null vector was utilized as the initial estimate for all the variables of interest and with the use of double precision. The standard coarsening ratio was utilized, that is , ''r'' = 2 [5], where the Multigrid method started from the finest grid and got to the coarsest grid possible, that is, the grid with ''N'' = 2x2 volumes. The number of unknowns in the finest grids was:'' N'' = 16x16, 32x32, 64x64, 128x128, 256x256, 512x512 and 1024x1024. The number of pre- and post-smoothing was equals to  <math>\nu ={\nu }_1={\nu }_2=3</math> .

:* ''l_∞''-norm for numerical error (||''E''||''_∞'');

:* l₂-norm for nondimensionalized residue based on the initial estimate (Eq. (54));

| <math>q^{\mbox{(}k\mbox{)}}=\frac{{\Vert r^{\left(k\right)}\Vert }_2}{{\Vert r^{\left(k-1\right)}\Vert }_2}</math>

==5.2 Analysis of the pressure==

The behavior of the ''l_∞''-norm was initially evaluated for the four studied methods, where it could be noticed that all of them presented a reduction in the numerical error norm with the increase in the size of the problem; a desirable characteristic in a convergent numerical method (see Fig. 5).

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image126.jpeg|432px]] </div>

Figure 6 presents the number of iterations necessary in order to reach the stop criterion for the problem with ''N'' = 1024X1024. Based on this figure, it can be noticed that PCG-MG-ILU reached the stop criterion in a lower number of iterations.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image127.jpeg|456px]] </div>

Figure 7 illustrates the behavior of the empirical convergence factor with the increase of the size of the problem. Based on this figure, it can be claimed again that the PCG-MG-ILU methodology presented the best results, since it is the one that generated the best values for the convergence factors, that is, values near to zero [5,34]:

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image128.jpeg|456px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

for ''p''.</div>

 

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image129.jpeg|432px]] </div>

Table 1 shows the Speed-up of the three methods MG-GS, MG-ILU and PCG-MG-GS in relation to the PCG-MG-ILU method.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

'''Table 1: '''Speed-up values for MG-GS, MG-ILU and PCG-MG-GS in relation to PCG-MG-ILU</div>

{| style="width: 62%;margin: 1em auto 0.1em auto;border-collapse: collapse;"  

|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|''N''

|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|MG-GS

|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|MG-ILU

|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|1

|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|1

|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|1

|  style="text-align: center;vertical-align: top;"|32x32

|  style="text-align: center;vertical-align: top;"|64x64

|  style="text-align: center;vertical-align: top;"|128x128

|  style="text-align: center;vertical-align: top;"|256x256

|  style="text-align: center;vertical-align: top;"|512x512

|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|1024x1024

|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|1.074

|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|1.194

|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|1.167

==5.3 Analysis of the velocities==

Figure 9 shows the behavior of the numerical error ''versus'' the size of the problem, for the several solution methods. It can be noticed that the four studied methods presented an expected characteristic for the convergent iterative methods, that is, the decrease in the numerical error with the increase in the size of the problem. It can also be noticed that the presented error have values which are very close among themselves.

Figure 10 displays the behavior of ''l''₂-norm of the nondimensionalized residue by ''l''₂-normof the initial estimate in function of the number of iterations for the problem with ''N'' = 1024X1024 volumes. Based on this figure, it can be noticed that the PCG-MG-ILU and the MG-ILU methods were the ones which achieved the stop criterion in the lowest number of iterations.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

'''Figure 9:''' ''l_∞''-norm for the numerical error (||''E''||''_∞'') ''versus'' size of the problem (''N'') for ''u''.</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image132.jpeg|456px]] </div>

Figure 11 illustrates the behavior of the empirical convergence factor in function of the size of the problem. Based on this figure, it can be concluded that the PCG-MG-ILU method presented the best results, since it generated the best values for the convergence factors (closest values to zero), besides decreasing its value with the increase of the size of the problem.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image133.jpeg|456px]] </div>

Figure 12 shows the CPU times (in seconds) of the simulations with the increase of the size of the problem. Based on this figure, it can be noticed that the computational times were very close, with a slight advantage for the PCG-MG-ILU methodology, which can be confirmed in the results presented in the sequence.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image134.jpeg|432px]] </div>

Table 2 presents the Speed-up (Eq.(56)) of the MG-GS, MG-ILU and PCG-MG-GS methods in relation to the PCG-MG-ILU method. In this table, it can be confirmed that the CPU times are close among themselves, but there is an advantage with the use of PCG-MG-ILU, that is, the method that most accelerated the convergence of the iterative process for the presented problem.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

'''Table 2: '''Values of Speed-up for MG-GS, MG-ILU and PCG-MG-GS in relation to PCG-MG-ILU</div>

{| style="width: 62%;margin: 1em auto 0.1em auto;border-collapse: collapse;"  

|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|''N''

|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|MG-GS

|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|MG-ILU

|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|1

|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|1

|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|1

|  style="text-align: center;vertical-align: top;"|32x32

|  style="text-align: center;vertical-align: top;"|64x64

|  style="text-align: center;vertical-align: top;"|128x128

|  style="text-align: center;vertical-align: top;"|256x256

|  style="text-align: center;vertical-align: top;"|512x512

|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|1024x1024

|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|1.337

|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|1.228

|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|1.535

==5.4 Algorithm==

After the tests of the two previous phases (sections 5.2 and 5.3), whose objective was to identify the best methods for the solution of the equations for the variables ''p'' (Poisson equation) and ''u'' and ''v'' (Reaction-Diffusion equations), we get to the point where the best methods were applied simultaneously to solve the Navier-Stokes equations.

The results presented below are compared the classic Gauss-Seidel (GS) method in its singlegrid version (only grid method) so that the gains can be measured when using the proposed algorithm (adaptation of the Projection method presented in Algorithm 1), that is, using the PCG-MG-ILU solver in the solution of the linear systems given in steps 1 and 2 of Algorithm 1:

Figure 13 presents the behavior of ''l''₂-norm of the nondimensionalized residue by the initial estimate in function of the number of iterations for the problem with ''N'' = 128X128 volumes for the variables ''u'', ''v'' and ''p''. Based on this figure, it can be noticed that the PCG-MG-ILU method reached the stop criterion ε = 10^–10 in an extremely lower number of iterations, both for ''p ''and for ''u'' and ''v''.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image135.jpeg|456px]] </div>

Figure 14 illustrates the behavior of the empirical convergence factor with the increase of the size of the problem. Based on this figure, it can be concluded that the PCG-MG-ILU method presented much better results for the convergence factors for the three variables of interest, comparing with the classic GS method.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image136.jpeg|456px]] </div>

Figure 15 contains the values of the CPU times (in seconds) of the simulations with the increase of the size of the problem. Based on this figure, it can be noticed that the PCG-MG-ILU methodology was the one that converged faster, being quite evident the difference in CPU time between the proposed methodology and the classic GS method.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

[[Image:Draft_Anunciacao_160657346-image137.jpeg|432px]] </div>

Table 3 presents the Speed-up (Eq.(56)) on the GS method in relation to the PCG-MG-ILU method. According to the data on this table, the methodology that uses the PCG-MG-ILU method proved much more efficient than the methodology that utilizes the GS method, reaching a 105 times higher velocity for the biggest problem compared in this table. It can also be noticed that the Speed-up increases as the size of the problem increases, which is a highly desirable property.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

'''Table 3: '''Values of GS Speed-up in relation to PCG-MG-ILU</div>

{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;"  

|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|''N''

|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|''S_P''

|  style="text-align: center;vertical-align: top;"|16x16

|  style="text-align: center;vertical-align: top;"|32x32

|  style="text-align: center;vertical-align: top;"|64x64

|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|128x128

|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|105.242