Updating incomplete factorization preconditioners for shifted linear systems arising in a wind model

Updating incomplete factorization preconditioners for shifted linear systems arising in a wind model A. Suárez, H. Sarmiento, E. Flórez, M.D. García, G. Montero University Institute for Intelligent Systems and Numerical Applications in Engineering, University of Las Palmas de Gran Canaria. Edificio Central del Parque Científico y Tecnológico, Campus Universitario de Tafira, 35017 Las Palmas de Gran Canaria, Spain Abstract The efficiency of a finite element mass-consistent model for wind field adjustment depends on the stability parameter αwhich allows from a strictly horizontal wind adjustment to a pure vertical one. Each simulation with the wind model leads to the resolution of a linear system of equations, the matrix of which depends on a function ε(α), i.e., (M+εN)xε=bε, where Mand Nare constant, symmetric and positive definite matrices with the same sparsity pattern for a given level of discretization. The estimation of this parameter may be carried out by using genetic algorithms. This procedure requires the evaluation of a fitness function for each individual of the population defined in the searching space of α, that is, the resolution of one linear system of equations for each value of α. Preconditioned Conjugate Gradient algorithm (PCG) is usually applied for the resolution of this type of linear systems due to its good convergence results. In order to solve this set of linear systems, we could either construct a different preconditioner for each of them or use a single preconditioner constructed from the first value of εto solve all the systems. In this paper, an intermediate approach is proposed. An incomplete Cholesky factorization of matrix Aεis constructed for the first linear system and it is updated for each ε at a low computational cost. Numerical experiments related to realistic wind field are presented in order to show the performance of the proposed preconditioning strategy. Key words: Incomplete factorization, Shifted linear systems, Preconditioning, Conjugate gradient, Wind modelling, Genetic algorithms Preprint submitted to Elsevier 1 November 2009 1 INTRODUCTION The application of discretization techniques to problems defined by partial differential equations that model physical phenomena leads to linear systems of equations of which the matrix is sometimes given as a function of a parameter. This is specifically true in the numerical simulation of wind fields with mass consistent models [1,2]. In general, these problems are defined over regions with complex terrain, therefore a suitable discretization of the studied zone is necessary. Here, we have used the technique of Montenegro et al. [3] for constructing tetrahedral meshes such that they are adapted to the terrain orography and have a higher density of nodes near the terrain surface. Moreover some regions may need an additional refinement due, for example, to more accurate approximation in those zones. On the other hand, the combination of the model with a Gaussian plume approach makes the refinement along the trajectory of the plume necessary. So, in general, we are going to work with meshes including elements of very different size. This fact affects the conditioning of the linear systems of equations that arises from this type of discretization in this problem, i.e., Aεxε=bε(1) Thus, a suitable preconditioning technique should be applied for an efficient conjugate gradient iteration. We are interested on the preconditioning of (1) in the particular case Aε=M+εN (2) where Mand N, are symmetric and positive definite matrices and remain constant along the process for a given discretization level. These linear systems appear in each step of the model. Two different cases exist where a sequence of linear systems like (1) must be solved. On the one hand, we have applied genetic algorithms to estimate εfor each given set of station measurements [2,4]. Genetic algorithms (GAs) are optimisation tools based on the natural evolution mechanism. They produce successive trials that have an increasing probability of obtaining a global optimum. The main aspects of GAs are the construction of an initial population, the evaluation of each individual of a population with a fitness function, the selection of the parents of the next generation, the crossover of those parents to create the children and the mutation to increase diversity. In practice, the initial population is randomly generated and we use iteration limit exceeded as the stopping criterion. The selection phase allocates an intermediate population on the basis of the evaluation of the fitness function. This evaluation, which is in general related to the difference between the observed and computed wind at the stations, leads us to solve a complete wind simulation for each individual of the population. In other words, we have to solve a linear system of equations like (1) for each 2 value of ε. On the other hand, for a given sequence of εvalues corresponding to a time interval of a simulation episode, a set of linear systems must also be solved. Two extreme strategies for preconditioning such linear systems may be applied. Indeed, we could either construct a different preconditioner for each of them and improve the convergence of Preconditioned Conjugate Gradient (PCG) [5] with the consequent high computational cost related to the construction of each preconditioner, or use a single preconditioner constructed with a given value of εfor all the systems. In this latter case the convergence will be getting worse as the value of the parameter moves away from the initial value. In the particular case of shifted linear systems, Meurant [6] and Benzi et al. [7] propose an intermediate solution by using two type of preconditioners, respectively, which are constructed once at the beginning of the process and updated for each εat a low computational cost. These strategies also lead to an intermediate rate of convergence between the above extreme options. On the one hand, Meurant studied the case of a shifted matrix of the form Aε= M+εD, with Dbeing a diagonal matrix, and updated the preconditioner from an incomplete Cholesky factorization of M. Here, we generalize Meurant’s algorithm to the case of a generalised shifted matrix Aε=M+εN, with Mand NSPD matrices with the same sparsity pattern arising from a finite element discretization of mass consistent mathematical model. On the other hand, Benzi et al. develop the study of factorised approximate inverses using the SAINV algorithm [8] for the special case of Aε=M+εI, with Ibeing the unit matrix. This approach is oriented to parallel computating, since the construction of the preconditioner may be carried out in parallel. A similar updating method for a generalised shifted matrix using SAINV algorithm is proposed in [9]. The organisation of the paper is as follows. In Section 2, the mass consistent model is presented. It generates a wind velocity field for an incompressible fluid which adjusts to an initial field obtained from experimental measurements and physical considerations. The construction of the initial field may be found in [2] and is not directly involved in the updating of matrix Aε. Nevertheless, we remark the study of the stability parameter αof the wind model which is directly related to εand, thus, is the generator of the set of linear systems. The construction of the preconditioner is carried out in terms of an incomplete Cholesky factorization of matrix Aε=M+εN, with the corresponding simplifications that allow it to be update at a reasonable cost. This is described in Section 3. Section 4 is devoted to illustrating the performance of this preconditioner in same numerical experiments. Finally, our conclusions are presented in Section 5. 3 2 WIND MODEL This model [1] is based on the continuity equation for an incompressible flow where the air density is constant in the domain Ωand no-flow-through conditions on Γb(terrain and top) are considered ~ ∇ · ~u = 0 in Ω(3) ~n ·~u = 0 on Γb(4) The problem is formulated as a least-square approach in Ω, with ~u(e u, e v, e w)to be adjusted E(~u) = ZΩhα2 1(e u−u0)2+ (e v−v0)2+α2 2(e w−w0)2idΩ(5) where the interpolated wind ~v0= (u0, v0, w0)is obtained from experimental measurements and physical considerations, and α1,α2are the Gauss precision moduli. In practice, we use the so called stability parameter of the wind model, α=α1 α2 (6) since the minimum of the functional given by (5) is the same if we divide it by α2 2. So, if α >> 1, flow adjustment in the vertical direction predominates. However if α << 1, flow adjustment occurs primarily in the horizontal plane. Thus, the selection of αallows the air to go over a terrain barrier or around it, respectively. Moreover, the behaviour of mass consistent models in many numerical experiments has shown that they are very sensitive to the value chosen for α. In [4,2] a brief discussion about the selection of αby several authors is presented. Solving (5) constrained by (3) and (4) is equivalent to find a saddle point (~v, φ)of the Lagrangian E(~v) = min ~u∈KE(~u) + ZΩ φ~ ∇ · ~u dΩ(7) where ~v = (u, v, w),φthe Lagrange multiplier and Kthe set of admissible functions. The Lagrange multipliers technique is used to minimise the problem (7), whose minimum comes to form the Euler-Lagrange equations u=u0+1 2α2 1 ∂φ ∂x, v =v0+1 2α2 1 ∂φ ∂y , w =w0+1 2α2 2 ∂φ ∂z (8) 4 Since α1and α2are constant in Ω, the variational approach results in an elliptic problem substituting (8) in (3) ∂2φ ∂x2+∂2φ ∂y2+α2∂2φ ∂z2=−2α2 1 ∂u0 ∂x +∂v0 ∂y +∂w0 ∂z !in Ω(9) We consider Dirichlet condition for flow-through boundaries and Neumann condition for terrain and top φ= 0 on Γa(10) ~n ·T~ ∇µ=−~n ·~v0on Γb(11) with T=diag h1 2α2 1 ,1 2α2 1 ,1 2α2 2i. The problem given by (9)-(11), is solved using tetrahedral finite elements (see [4,2]), which leads to a set of 4×4elemental matrices and 4×1elemental vectors related to element Ωe, with ˆ ψibeing the form function of the i-th node, i= 1,2,3,4, defined in the reference element ˆ Ωeand |J|the Jacobian of the transformation from Ωeto ˆ Ωe, {Ae}ij =Zˆ Ωe((∂ˆ ψi ∂ξ ∂ξ ∂x +∂ˆ ψi ∂η ∂η ∂x +∂ˆ ψi ∂ϕ ∂ϕ ∂x )(∂ˆ ψj ∂ξ ∂ξ ∂x +∂ˆ ψj ∂η ∂η ∂x +∂ˆ ψj ∂ϕ ∂ϕ ∂x ) + (∂ˆ ψi ∂ξ ∂ξ ∂y +∂ˆ ψi ∂η ∂η ∂y +∂ˆ ψi ∂ϕ ∂ϕ ∂y )(∂ˆ ψj ∂ξ ∂ξ ∂y +∂ˆ ψj ∂η ∂η ∂y +∂ˆ ψj ∂ϕ ∂ϕ ∂y )(12) +α2(∂ˆ ψi ∂ξ ∂ξ ∂z +∂ˆ ψi ∂η ∂η ∂z +∂ˆ ψi ∂ϕ ∂ϕ ∂z )(∂ˆ ψj ∂ξ ∂ξ ∂z +∂ˆ ψj ∂η ∂η ∂z +∂ˆ ψj ∂ϕ ∂ϕ ∂z )) · |J|dξ dη dϕ {be}i=Zˆ Ωe −2α2 1(u0(∂ˆ ψi ∂ξ ∂ξ ∂x +∂ˆ ψi ∂η ∂η ∂x +∂ˆ ψi ∂ϕ ∂ϕ ∂x ) +v0(∂ˆ ψi ∂ξ ∂ξ ∂y +∂ˆ ψi ∂η ∂η ∂y +∂ˆ ψi ∂ϕ ∂ϕ ∂y )(13) +w0(∂ˆ ψi ∂ξ ∂ξ ∂z +∂ˆ ψi ∂η ∂η ∂z +∂ˆ ψi ∂ϕ ∂ϕ ∂z ))· |J|dξ dη dϕ Note that if we set ε=α2, the elemental matrix may be written as {Ae}ij ={Me}ij +ε{Ne}ij (14) The assembling of such matrices yields a symmetric matrix Agiven by equation (2). Here we propose to solve the corresponding linear system of equations 5 (1) by using PCG algorithm. Once we have obtained φ, the resulting wind field is computed using equation (8). 3 UPDATING OF THE INCOMPLETE CHOLESKY FACTORIZATION We will generalize the incomplete factorization proposed by Meurant [6] for the case of matrices Aε=M+εD, with Dbeing diagonal, to matrices Aε= M+εN, with Mand Nbeing two n×nsymmetric positive definite matrices that have the same sparsity pattern in this case. We can write Aεas follows, Aε= (mij) + ε(nij) =    m11 +εn11 (f1M+εf1N)T f1M+εf1NM2+εN2    where f1M, f1Nrepresent (n−1) ×1column matrices and M2, N2,(n−1) × (n−1) matrices. A factorization of the first row and column of Aεis carried out, Aε=    m11 +εn11 0 l1M+εl1NI    (m11 +εn11)−10 0C2      m11 +εn11 (l1M+εl1N)T 0 I   = L1Z1LT 1 with l1M=f1Mand l1N=f1N. Then, identifying term by term, we obtain for matrix C2, C2=M2+εN2−1 m11 +εn11 (l1M+εl1N) (l1M+εl1N)T(15) In order to build the preconditioner, if we only consider the diagonal entries of Nas first approximation, then equation (15) is simplified since l1N= 0, C2=εD2+M2−1 m11 +εn11 l1MlT 1M, An order 0algorithm is derived from C2=εD2+M2−1 m11 l1MlT 1M 6 and the entries of C2are computed by adding εD2to what we would have obtained for the incomplete decomposition of M. Another approximation consists of considering all the entries in N2and neglecting the products εl1Nin (15). So, the successive computations of matrices Cido not involve εand those matrices may be obtained easily from the M decompositions, C2=εN2+M2−1 m11 l1MlT 1M and thus, in matrix form, C2=εN2+   m(2) 22 fT 2M f2MM3   =   m(2) 22 +εn22 (f2M+εl2N) f2M+εl2NM3+εN3    Only the entries of f2Mcorresponding to non null entries of Mare computed in order to avoid the fill-in, obtaining l2M. So the decomposition of C2results in C2≈    m(2) 22 +εn22 0 l2M+εl2NI    m(2) 22 +εn22−10 0C3      m(2) 22 +εn22 (l2M+εl2N)T 0 I    where, identifying C3=M3+εN3−1 m(2) 22 +εn22 (l2M+εl2N) (l2M+εl2N)T Similarly, with the same simplifications, we have, C3=M3+εN3−1 m(2) 22 l2MlT 2M=   m33 +εn33 (f3M+εl3N)T f3M+εl3NM4+εN4    which is constructed following the same procedure of C2. In this way, once all the matrices Ciare constructed, the incomplete decomposition of Aεresults in, Aε≈L1Z1LT 1=L1L2Z2LT 2LT 1= (L1L2···Ln)Zn(L1L2···Ln)T(16) 7 Zbeing the diagonal matrix,               (m11 +εn11)−1. m(2) 22 +εn22−1. m(3) 33 +εn33−1. . . . . . .m(n) nn +εnnn−1               The diagonal entries of the lower triangular matrix L1L2···Lnare m(i) ii +εnii. The respective columns below the diagonal entries are defined by (n−i)×1 matrices ljM +εljN . 4 NUMERICAL EXPERIMENTS In this section we present the results obtained using PCG with the proposed preconditioners for solving the linear systems of equations arising from the elliptic equation related to a 3-D mass consistent model for wind field adjustment [1,2]. All the experiments were carried out on a XEON Precision 530 with Fortran Double Precision. In the resolution, we always started from the null vector and stopped if krkk2≤10−10 kr0k2or if the number of iterations was greater than 10000. Two cases are presented and two sets of values of parameter αwere generated for each of them. The initial population corresponding to the first set with values between 0and 20, is obtained with the random generator of the GAs package of Levine [10]. It should be as diverse as possible in order to obtain satisfactory results with GAs [11]. For the set with values between 0and 100, the initial population is generated with a normal distribution. The first case is related to a wind simulation in a region of Gran Canaria Island and we have used an adaptive mesh to produce linear systems of 98999 equations. The second case is also related to a wind simulation in the whole of Gran Canaria Island. We have used another adaptive mesh to produce linear systems of 100643 equations. Both matrices Mand Ncorresponding to the above problems, are SPD. Figures 1-4 show the behavior of the above mentioned preconditioners. For a wide range of values of α, the timings for reaching convergence are represented in each case. ICHOLDand ICHOLNare the ICHOL preconditioners obtained with the two approaches developed in section 3, respectively. These preconditioners are compared with full-ICHOL of matrix Aε, that is, computing a new 8 ICHOL decomposition for each ε, and with the use of a single preconditioner, ICHOL(Aε0), for all the linear systems. 0 5 10 15 20 0 50 100 150 200 250 300 α CPU time (s) icholA icholD icholN fullichol Fig. 1. Windfield.98999: Convergence of PCG with several ICHOL preconditioners for different values of parameter αrandomly calculated. 0 20 40 60 80 100 0 50 100 150 200 250 300 α CPU time (s) icholA icholD icholN fullichol Fig. 2. Windfield.98999: Convergence of PCG with several ICHOL preconditioners for different values of parameter αusing normal distribution. In all the cases, the proposed ICHOLNpreconditioner led to the fastest convergence. For low values of α, the full-ICHOL preconditioner displayed better behavior. However, for higher values of α, it got worse very rapidly and did not even allow convergence to be reached. This result may be explained since we can not always obtain an incomplete factorization ICHOL from the definite positive matrix (M+εN)which preserves symmetry and positivity. 9