Linear System Solutions of the Navier-Stokes Equations with Application to Flow over a Backward-Facing Step ()
1. Introduction
Numerical solution of the Navier-Stokes equations is a crucial problem in engineering and physical sciences. We consider the solution of the incompressible Navier-Stokes equations governing the flow of viscous Newtonian fluids whose form reads
(1)
where u is the velocity vector, p the pressure and f the field of external forces. The constant
is the kinematic viscosity. The first equation models the conservation of momentum fluid. The second equation models the conservation of mass. We consider the problem posed on a domain
of dimension 2 or 3 with boundary conditions defined by
(2)
where
denote normal vector and
. We consider the following notations
and Sobolev space is defined as follow:
We define solution and test spaces:
The variational formulation of such problems is the following (1) find
and
such that:
(3)
Let
et
. The linearized form of (3) is given as follow find
and
, such that:
(4)
where
et
. We use finite dimensional spaces. Let
et
, the discret form of (4) is defined as follow: find
and
(5)
for all
and
.
Then we use (5) as in function of the basis
and
, to find the following saddle point system
(6)
where
,
, A denote the laplacian matrix
B denote the divergence matrix
D is the nonlinear matrix
where
,
et
. The right-hand side is defined as follow
For finding the root, the Newton method is often implemented. It is an iterative algorithm consisting in assigning the root of the linearized equation at the current iterate to the next iterate to refine the solution until a convergence criterion is satisfied. Implementing the Newton method requires evaluating the Jacobian matrix. In the Newton method, algebraic linear systems (6) involving the Newton matrices must thus be solved at each iterate. We are thus interested in the following linear system:
(7)
where matrix
is generally nonsymmetric and vector
is the right-hand side at the current iterate. This linear system can be solved directly by decomposing the Newton matrix into the product of a lower triangular and upper triangular factor (LU factorization) and by solving the two resulting triangular systems through forward and backward substitutions. Whenever it is prohibitively costly to assemble Newton matrix
, the linear system (1) is to be solved using a matrix free method. It is indeed possible to resort to iterative linear methods, for instance a Krylov subspace projection technique. Implementing any iterative method in this context would require to repeatedly compute the application of the Newton matrix on the current iterate. Numerical differentiation via the difference quotient being prone to truncation and round-off errors, this inexact way of proceeding is usually referred to as the inexact Newton method. More general approximations of the Newton matrix can be used to perform the Newton iterations. Proceeding this way precisely amounts to implementing a modified Newton method. This means that an approximate inverse of the Newton matrix is employed to indicate the tangent direction. The exact Newton method resorting to the exact Newton matrix
, as opposed to the inexact or modified Newton methods, is seldom used because assembling the exact Newton matrix at each iteration is computationally too expensive in most practical problems of interest. In this work, we investigate three alternative strategies to perform the Newton iterations based on the factorization of an approximate Newton matrix.
The first strategy consists in implementing a modified Newton method and solving the linear systems (7) appearing at the Newton iterations using lower-upper factorization (LU). Is an efficient solver to perform the LU factorization of a sparse matrix and the triangular resolutions. Factoring the Newton matrix being computationally intensive, the task is done occasionally, not at every Newton iteration nor even time step. The previously updated Newton matrix is recycled several times. The Newton matrix is computed again when the number of Newton iterations to satisfy the converge criterion exceeds a certain parameter value. This is the default way of proceeding of the IFISS software [1] used in numerous incompressible flow simulation [2] . The second strategy is to solve the linear systems occurring during the Newton iterations using an iterative method. Any iterative solver consists in generating a sequence of improving approximate solutions by repeatedly applying the current approximate onto the matrix until a convergence criterion is satisfied. It can possibly be implemented within both the modified and inexact Newton methods, depending on whether the matrix-vector product involves an approximate Newton matrix or is computed inexactly from a different quotient approximating the exact Newton matrix. We apply the generalized minimum residual (GMRES) method [3] [4] [5] . This iterative method is based on the standard Arnoldi algorithm [6] , is a particular Krylov subspace projection method and is well adapted to nonsymmetric indefinite linear systems of equations, unlike the conjugate gradient method that rather deals with symmetric definite positive matrices. Krylov subspace projection methods are suited to solve linear systems in which the involved matrix is large and sparse. An important number of iterations are however necessary to obtain a reliable approximate of the solution when the condition number of the matrix is large. For ill-conditioned matrices, it is preferable to apply a preconditioner to accelerate the convergence of the associated GMRES method. In practice, the preconditionning matrix
should be chosen close to the inverse of the Newton matrix. Here, we consider the Newton matrix obtained at a previous time stem and compute its LU factorization [7] or an incomplete LU factorization [8] . Preconditioning then amounts to applying the inverted matrix on the current approximate by solving the two associated triangular linear systems by forward and backward substitutions.
The third strategy consists in solving the linear systems using a Schur approach. The goal is to take advantage of the particular structure of the saddle point matrix. N corresponds to the number of differential equations. The saddle point system to solve can be rewritten as follows:
(8)
where:
is a sparse and nonsymmetric matrix,
maybe a rank deficient matrix with
,
is a dense matrix,
and
are given vectors. Besides, we have
with
. The Schur approach consists in computing the Schur complement of the partitioned Newton matrix, in which the linear system associated with the large sparse block is handled by LU factorization. The paper is organized as follows. We next introduce the direct approach (Section 2.1), iterative approach (Section 2.2) and Schur approach (Section 2.3) used for solving the saddle point systems arising within the modified or inexact Newton algorithms. We present simulation results for the lid driven cavity and chanel flow systems obtained using the various solvers in Section 3. We finally conclude in Section 4 and give some recommendation on which approach to implement depending on the nature of the problem.
2. Implementation of Methods
With the exact Newton method, the saddle point matrix must be assembled at every iteration, which represents a costly task. A modified Newton method that is often computationally more efficient in practice consists in reusing the Newton matrix several times. This way of proceeding usually necessitates more Newton iterations, but these ones are performed faster. In all implementations of Newton methods, we allow at most 10 Newton iterations, but this limit can be changed by the user. The input constants used by the main solver (maximum number of Newton) as well as their meaning are described in detail in Ref [1] .
2.1. Direct Approach
Linear Equation (8) is solved directly by resorting to lower-upper factorization implemented [9] . Once the saddle point matrix has been factorized into
, the linear system (8) is solved via forward and backward substitution. This amounts to successively solving the two following triangular systems of equations
(9)
(10)
The solution Y of system (9) as obtained via forward substitution is to be injected into system (10), yielding X via backward substitution. Solving the two triangular linear systems amounts to performing a modified Newton iteration whose cost is much smaller than that of factoring the Newton matrix.
2.2. Preconditioned Iterative Approach
Iterative methods are ideally suited to solve high-dimensional sparse linear systems of equations of the form (7). Their advantage over direct methods is that they don’t require to factorize saddle point matrix
nor even evaluate it. This is the case for Krylov subspace projection methods and, among them the GMRES method implemented in the following. The only request is the ability to compute the application of the matrix on any vector v. In the present context, this task can be done without evaluating the saddle point matrix. Hence, vector
can be obtained through a finite difference along direction v.
The drawback of iterative methods is that they are inexact and may require a large number of iterations so that the iterate satisfies the tolerance condition. The remedy to this technical drawback is called preconditioning. Here, the preconditioner
applies a linear transformation to system (7) so as to reduce the condition number of the transformed matrix
. The preconditioned linear system to solve writes:
(11)
If the condition number of the transformed matrix
is smaller than that of matrix
then the number of iterations is generally reduced. Several ways of implementing the preconditioned iterative approach can be designed depending on the choice of the preconditioning linear transformation
and Newton matrix
. We furthermore consider the preconditioner built upon the incomplete lower-upper (ILU) factorization of the last-updated Newton matrix. Its potential advantage resides in the smaller amount of memory necessary to perform an incomplete factorization. This feature is useful when the model system is very large and the Newton matrix cannot be factored due to memory limitations. The three preconditioners are denoted by PD, PT and PR respectively, while notations PDGMRES, PTGMRES and PRGMRES stand for the corresponding preconditioned GMRES methods. When solving Equation (11) using any of the preconditioned GMRES methods described above, the iterations are stopped as soon as the Euclidean norm of the current residue is lower than a tolerance threshold (
):
(12)
where
denotes the current iterate and
is the threshold value. The maximum number of iteration of the PGMRES is 200.
2.3. Schur Approach
The Schur approach is an alternative direct method aiming at solving linear system (8). The goal is to take advantage of the structure of the saddle point matrix
. The preliminary steps for implementing the Schur approach consists in computing the Schur complement matrix
associated with matrix
as illustrated in Figure 1 and listed below:
Figure 1. Schematic diagram illustrating the block LDU decomposition of S and the computation of the Schur complement matrix.
1) Perform the block decomposition of
into four block matrices such that the diagonal block
and matrix D.
2) Formally perform the lower-diagonal-upper block decomposition (LDU) of
.
3) Formally define the Schur complement matrix
.
4) Solve the systems
, with multiple right-hand sides
using LU factorization.
5) Evaluate the Schur complement
by computing the product of
and
.
The next steps of the Schur approach consist in solving the two triangular systems of linear block equations using the block LDU factorization so as to obtain the solution of overall linear system (8). These additional steps are depicted in Figure 2 and detailed below:
1) The matrix, unknown and right-hand side of linear system (8) are first partitioned.
2) The solution part y is computed through forward block substitution; the dense system involving the previously computed Schur matrix
is solved through LU factorization.
3) The solution part x is computed through backward substitution and the sublinear system involving A is solved directly using LU.
4) The solution X of block linear system (8) is eventually obtained by concatenating partial solutions x and y.
2.4. Modified Newton Method
In the following applications, the root-finding algorithm will be a modified Newton method when used in combination with preconditioned iterative solver, also Schur approach and direct approach. In this situation, the Newton matrix
is computed at a previous time step and is fixed throughout the nonlinear iterations. When used in combination with an iterative approach, i.e. the preconditioned GMRES method, the root-finding algorithm is either a modified Newton method, depending on whether the matrix-vector products
is computed using a partially updated Newton matrix.
Figure 2. Schematic diagram of forward and backward block substitution in Schur approach. The main costs at each modified Newton iteration consist in solving the two sub-linear systems involving the sparse high-dimensional matrix Ae. The small and dense linear system involving matrix S is solved using LU factorization since the dimension of the matrix is too small.
3. Numerical Results
In this section, we carry out two numerical example simulations, to assess the relative efficiencies of the various algebraic numerical methods described in Sections 2.1, 2.2 and 2.3. For the direct, iterative and Schur approaches, In the meanwhile the feasibility and effectiveness of the approaches, from the point of view of the elapsed Central Process Units time for the solution of saddle point problem (8) (denoted by CPU) as well as the number of iterations of PGMRES. The parameters introduced in the lid driven cavity and flow over backward facing step are listed in Table 1, for more details see [1] . We use the IFISS software package developed by Elman et al. [1] to generate the linear systems corresponding to
,
and
. The IFISS software provides the matrices A, B, and the right-hand side f and g. Generic information of the test problems, including n and m, are provided in Table 1.
In practice, the preconditioners used for solving (6) are
,
and
, where
,
and
are given as follows:
(13)
Here S is a sparse approximation of the pressure Schur complement
, and Q is one of the matrices I or (S). The parameter of the
preconditioner is chosen as to implement the regularized preconditioner effeciently, we need to choose the parameters α appropriately since the analytic determination of the parameters which results in the fastest convergence of the preconditioned GMRES iteration appears to be quite a difficult problem. In all Tables, to implement the regularized preconditioner efficiently, we need to choose the parameters α appropriately since the analytic determination of the parameters which results in the fastest convergence of the preconditioned GMRES iteration appears to be quite a difficult problem. In the regularized preconditioner, the parameter α is taken as
, which balances the matrices A and
in the Euclidean norm.
Example: Flow over backward facing step
We consider the Flow over backward facing step. The simulation parameters used for Lid driven cavity model are summarized in Table 1, for more details see [1] . Numerical results for all approaches are presented in Tables 2-4.
In all Tables indicates that the PR preconditioner with
leads to much better numerical results than the PD and PT preconditioned GMRES methods less CPU times in all trials and iteration. the preconditioned PRGMRES method with the proper parameter α has a better performance than the preconditioned PDGMRES and the preconditioned PTGMRES methods in terms of the iterations and CPU times. In the following figures, we display a streamline plot for the velocity solution, and a plot of the pressure solution of the Flow in a symmetric step channel, over a plate and over a backward facing step.
Figure 3 illustrates that the flow and pressure rendering in a rectangular duct containing a sudden expansion. We set the conditions of Dirichlet at the entry of the rectangular conduit and the conditions of Neumann at the exit of the duct. The vertical speed is zero, which means that the lines of the fluid flow propagate in a parallel way to the rectangular duct walls. Rendering pressure becomes closer to zero at the exit of the duct.
Table 1. The size of the matrices A and B on
.
aFlow over backward facing step.
Table 2. The numerical results of direct, Iterative and schur approaches.
aDirect solver. bIterative solver. cSchur solver.
Table 3. The numerical results of direct, iterative and schur approaches.
aDirect solver. bIterative solver. cSchur solver.
Figure 3. The velocity and pressure solutions of flow over backward facing step.
Table 4. The numerical results of direct, iterative and schur approaches.
aDirect solver. bIterative solver. cSchur solver.
4. Conclusion
The numerical results reported in Section 3 show that the preconditioned GMRES methods are more efficient than LU and Schur approaches. The two latter approaches are more expensive computationally in term of CPU times outperforms the other methods. For the test presently illustrated, the speed-up factors increases with the size of the Jacobian matrix and reaches a value of three. The Flow over backward facing systems for which the PGMRES method preconditioning exhibits speed-ups that are even more important compared with the implementation of Schur and direct approaches. As a result, the system of ordinary differential equations is extremely stiff. The stiffness implies that the saddle point matrix has a large condition number, i.e. is ill-conditioned. In this situation, Modified-Newton methods based on the inverse of the approximate Jacobian matrix require an important number of iterations. With PRGMRES, the iterations are based on the exact Jacobian inverse and for this reason the method reverts to the standard Newton method that converges in fewer iterations. Future developments will focus on the deterministic/Navier-Stokes that is enabled in [1] code when the simulated system is too large. When this situation occurs, the linear solver being very sensitive to solve the linear system, other direct linear solver will be tested. Future research should be devoted to the construction of parallel iterative solvers robust with respect to higher order virtual element discretizations and non-symmetric problems, such as saddle point systems deriving from discretizations of the Navier-Stokes equations.