A scalable geometric multigrid method for non-symmetric linear systems arising from elliptic equations

Center for Turbulence Research Annual Research Briefs 2015 211 A scalable geometric multigrid method for non-symmetric linear systems arising from elliptic equations By M. Esmaily-Moghadam, L. Jofre, G. Iaccarino A ND A. Mani 1. Motivation and objectives Elliptic partial differential equations (PDE) are commonplace in physics. A few examples range from heat transfer in solids to the effect of pressure in incompressible flows with variable density. Solutions to these PDE systems, under various boundary conditions as well as geometries, are not often available analytically. As a result, various numerical methods have been adopted to discretize the problem in space and time and to reduce it to a system of linear equations. Solving these systems remains a challenge, because scientific and engineering applications are continuously growing in size, demanding more efficient solvers (Saad 2003; Shakib et al. 1989; Saad et al. 1998; Esmaily-Moghadam et al. 2013, 2015a,b). In a common PDE solver, a large portion of the computational cost is dedicated to solving the underlying linear system. These linear systems have a typical size of over hundreds of thousands of unknowns. Since direct methods are of O(Np) with Nbeing the number of unknowns and pgenerally greater than 1.5, they quickly become an unattractive option in terms of computational cost (Saad 2003). To offset the cost, iterative methods are generally favored in these applications. The high computational cost, however, remains a critical issue in solving linear systems obtained from the Poisson equation (a subset of elliptic PDEs). This high cost is caused by the dependency of the local solution on the solution of the entire domain. In this class of problems, the entire domain is affected by a perturbation at any arbitrary point in space. This strong coupling introduces several challenges in terms of ill-conditioning of the linear system as well as designing efficient parallel algorithms. The ill-conditioning becomes a pressing issue as the size of the system increases (Saad 2003). To explain this in simple terms, consider two fixed points in space that are at the location of two unknowns on the discrete setting. The number of unknowns between these two points varies proportional to the total number of unknowns. Considering the intrinsic characteristic of the iterative methods, in which at each iteration information is propagated from a point to its neighbors, the number of iterations required for propagation of information between two points increases proportionally with the total number of grid points, i.e., the size of the system. This effect, that is also referred to as ill-conditioning of the linear system, leads to a scaling versus size worse than O(N) in iterative solvers. Although techniques based on Fourier transformation, which are of O(N) log(N), do not encounter this challenge, their applicability is limited to only symmetric systems with very simple geometries and boundary conditions. Another effective technique to tackle this issue is the multigrid approach, which can be divided into two classes: algebraic and geometric (Wesseling 1995; Vanˇek et al. 1996; Ghia et al. 1982). The algebraic multigrid techniques, although applicable to a wider range of linear systems, are often less robust, 212 Esmaily-Moghadam et al. more complex to implement, and less efficient than the geometric counterpart. The geometric multigrid methods, on the other hand, are well suited for applications with simple geometry that are large enough in size to justify the added complexity. In this study, we target a class of systems that arise from a specific physical problem. The underlying physical problem is flow in a duct, which is a cuboid domain with walls on four sides and an inflow-outflow boundary conditions on the third direction. Point particles are dispersed in the flow and the entire apparatus is subjected to radiative heating, causing temperature and density variations. The Navier-Stokes equations are formulated based on the low Mach assumption, leading to a variable density, but incompressible governing equations. Using a fractional time-stepping method, the conservation of mass is satisfied by correcting the velocity field through solving a linear system (Jofre et al. 2014). This linear system, as a result of density variation and grid non-uniformity, is not symmetric. For accuracy and cost considerations, the grid is refined near the walls, leading to a wide spectrum of cell sizes that translates to a linear system with diagonal entries varying significantly. Therefore, the objective of this study is to design an efficient and scalable linear solver for a non-symmetric linear system obtained from an elliptic operator on a cuboid domain. We target systems with up to a billion unknowns, a typical size for direct numerical simulation. Based on the earlier discussion, the geometric multigrid method is well suited for this class of problems, considering the simplicity of cuboid geometry and the size of the linear system. The massive size of the linear systems under consideration mandates the use of parallel computations. The main challenge here is to design an algorithm that scales up to a large number of processors without significantly losing parallel efficiency. The design of such an algorithm is fully intertwined with the underlying iterative method, because a scalable algorithm minimizes processor-to-processor communications and maintains a good load balance on all processors. A multigrid technique that relies on a static partitioning of a stretched grid can encounter load imbalance on coarser levels. To address these challenges, we will discuss a partitioning approach for the introduced geometric multigrid technique. This paper is organized as follows. First, we present a geometric multigrid technique for stretched grids that relies on a mapper between non-uniform and uniform grids, and discuss a recursive implementation of this algorithm for optimal convergence. Then, we introduce our partitioning approach, designed for improved parallel scalability. To benchmark our method, we compare it against a multigrid technique from the Trilinos package (Hestenes & Stiefel 1952; Heroux et al. 2005) and conventional Krylov-based iterative methods. In Section 4 we present: (1) triply periodic isotropic turbulence on uniform grids, and (2) variable density duct simulations on non-uniform grids. We present weak and strong scaling results to establish the parallel performance of our method. Finally, we draw conclusions in Section 5. 2. A geometric multigrid method In this section, we discuss our geometric multigrid method, which includes restriction (fine-to-coarse) and interpolation (coarse-to-fine) operations, a recursive implementation of the multigrid method, and the main algorithm. In what follows, roman superscripts are used to denote variable names and italic subscripts are used as indices. Superscripts c and f are used to denote variables on coarse and fine grids, respectively. Indices i,j, and kare used for fine, and I,J,Kare used for coarse grid. The multigrid technique tackles the ill-conditioning issue by reducing the number of A scalable multigrid method for elliptic linear systems 213 grid points in all directions, thus reducing the number of iterations required for the propagation of information. Hence, given the linear system defined on a fine grid Afxf=bf,(2.1) the objective is to solve an equivalent system on a coarse grid as Acxc=bc,(2.2) such that the difference between xcand xfis minimal in a physical sense. To obtain Eq. (2.2) from Eq. (2.1), one needs to map bfto bcand define Acsuch that solutions to both systems are similar. To achieve this goal, we take the PDE that produces Eq. (2.1) and re-discretize it on a coarse grid. This process reduces to two core operations, which are mapping a field from the fine grid to the coarse grid and vice versa. We denote the former (restriction) by Cand the latter (interpolation) by F. Exploiting these two operators, the right-hand side (RHS) of Eq. (2.1) is mapped to the coarse grid using bc=C(bf),(2.3) and the solution to Eq. (2.2) is mapped back to the fine grid using xf=F(xc).(2.4) We denote the cardinality of a set or vector by |•|,F:R|xc|7→ R|xf|and C:R|bf|7→ R|bc|. The left-hand side matrices in Eqs. (2.1)-(2.2) are discrete non-symmetric elliptic operators. In a continuous form, they represent a class of PDEs as ∇ · (κ∇T) = q, (2.5) which fits to the heat equation in a solid with variable material property, κ, under volumetric heating, −q, and temperature, T. The heat equation is chosen to provide a physical intuition. It directly translates to the equation solved for satisfying the continuity equation in flow with variable density by replacing Twith pressure, κwith the inverse of density, and qwith the divergence of uncorrected velocity field plus the contribution from the energy equation. Denoting the discrete form of an elliptic operator on grid G for a given κby Dκ{G}, by definition Af=Dκf{Gf}.(2.6) Through numerical experiments, in which we considered various combination of calculating Acdirectly from Afand the underlying PDE itself, we found that Ac=Dκc{Gc}(2.7) provides a good estimate of an optimal Ac. Specifically, an optimal Ac, given arbitrary non-singular Afand bf, can be potentially defined as argminAc  (Af)−1bf− F (Ac)−1C(bf)  .(2.8) Before proceeding any further, it is necessary to exactly specify operators Cand F. To satisfy conservation of energy (or mass), one needs to ensure that the integral of q in Eq. (2.5) over the entire computational domain is preserved. In a discrete setting, this constraint translates to a similar condition on Cpertaining to Eq. (2.3). To ensure this constraint is satisfied on any interval, we enforce it on all grid points, i.e., bc Iwith I∈ {1,··· ,|bc|}. This is accomplished by summing bfover all the cells that overlap with bc I. This is schematically shown for a one-dimensional grid in Figure 1. 214 Esmaily-Moghadam et al. bi-1 fbi fbi+1 f bI c w I,i-1 w I,i+1 w I,i w I-1,i-1 li-1 fli fli+1 f lI c Figure 1. One-dimensional mapping of a field bbetween a fine (top/green) and a coarse (bottom/red) grid. Based on Figure 1, defining wI,i as the fraction of length of cell ithat overlaps with cell I,Coperator turns into a multiplication by matrix win one dimension, i.e., bc=wbf. Note that wI,i ≥0 and wef=lc, where lcis an array that contains the length of coarse grid cells and efis a vector with one at all its entries defined on Gf. Also, wTec=lf, in which lfis an array that contains length of coarse grid cells. These two properties are a direct consequence of the fact that the integral of a quantity is preserved by w. In three dimensions, wis extended to three matrices, wx,wy, and wzfor x, y, and z directions, respectively. Cresults into three summations as C(bf)IJK =bc IJK =X iX jX k wx I,iwy J,jwz K,kbf ijk.(2.9) In Eq. (2.9), each summation must be carried out on a limited number of overlapping cells; hence, the computational cost remains as O(|bf|). Operator Fis also defined and extended to three dimensions in a similar manner. The only difference is that wis normalized by lfin this case. Defining ui,I =wI,i lf i ,(2.10) in one dimension we have bf=ubc, and in three dimensions F(bc)ijk =bf ijk =X IX JX K ux i,I uy j,J uz k,K bc IJK,(2.11) in which ux,uy, and uzare the corresponding weights for x, y, and z directions, respectively. Note uec=ef, viz. a field on Gcwith constant value is mapped to a constant value on Gf. Conventionally, multigrid methods have a preset pattern of restriction and interpolation. The most common patterns are V cycle, i.e., mrestriction followed by minterpolation, and W cycle, i.e., mlevels of restriction followed by ˜m≤minterpolation, ˜m restriction, and minterpolation. The performance of each of these methods depends on the quality of the solution on the coarser levels. For one class of problems, one cycle of restriction might be sufficient, while for another, several cycles might be necessary. Therefore, to achieve an optimal performance, we do not preset a pattern, but rather allow it to be determined dynamically (see Figure 2). This is achieved by adopting a recursive implementation. The goal of this algorithm is to optimally solve a linear system at any given level, regardless of the number of previous restrictions. Achieving this A scalable multigrid method for elliptic linear systems 215 0 10 20 30 40 4 3 2 1 0 Level Accumulated iterations Uniform grid Non-uniform grid Figure 2. Dynamic pattern of multigrid restriction and interpolations for two representative cases. Blue is a uniform 1283grid and red is a non-uniform 120 ×80 ×80 grid. objective at any level (including the coarsest level), by recursion, the entire algorithm would be optimal. A second important factor for achieving an optimal performance is to ensure smoothing operations, viz. reducing residuals at any given level including the coarsest level, are also optimal. The Krylov-based iterative methods, namely the conjugate gradient (CG) for symmetric matrices and the bi-conjugate gradient (BCG) for non-symmetric matrices, are among the most efficient smoothing techniques. Therefore, we use the CG and BCG for symmetric and non-symmetric matrices, respectively. Experimenting with other smoothing techniques such as Gauss-Seidel technique has confirmed near optimality of our choice of smoothers. In general, the CG outperforms the BCG. Exploiting this to maximize the performance, we ensure all restriction operations lead to a symmetric matrix. This way, the BCG is only used to solve a given non-symmetric matrix at the finest level, i.e., smoothing of the original system. To ensure a restriction operation leads to a symmetric matrix, we build the coarse grid such that it is uniform in any given direction, i.e., lc i=lc, in which lc≡ |lc|−1X i lc i(2.12) is the average grid size. Note that |lc|and Pilc iare equal to the number of grid points and the total length of the domain, respectively. Through numerical experiments, we found that the best restriction strategy is to increase the average grid size at maximum by a factor of two in each direction such that the average grid size in all directions merges to the same value. On a highly stretched grid this translates to significant coarsening to no coarsening (or even refining) on different segments of the grid. This strategy, however, remains superior in comparison with a coarsening strategy that ensures similar degree of coarsening in the entire domain. Mathematically, given the average grid size of the fine grid in x, y, and z directions as, lf,x,lf,y, and lf,z, respectively, we define a target coarse grid size as ∆≡2 min lf,x,lf,y,lf,z.(2.13) Setting the size of the grid to ∆ in all directions can lead to refinement in the direction with the largest lf. Hence, lcin each direction is determined based on max ∆,lf. More 216 Esmaily-Moghadam et al. specifically, lc,x= max ˜ ∆x,lf,x,(2.14) in which ˜ ∆xis the nearest value to ∆ with Pilf,x i/˜ ∆x∈N. Equations similar to Eq. (2.14) can also be written for lf,yand lf,z. We next put together all these components to build a multigrid algorithm. The main problem is defined as: given A,b, and ǫ, find xsuch that kAx −bk ≤ ǫ. First, let A={A0,A1,···Ammax }be a set of matrices in which A0=A,Aiis the matrix obtained from iconsecutive restriction operations (Eq. (2.7)), and mmax is the maximum number of multigrid levels. Note Airelative to Ai+1 is analogous to Afrelative to Acwith their relationship fully described above. Additionally, we define function K(˜ A,˜ b, nmax) to be an iterative solver (the CG for a symmetric and the BCG for a non-symmetric ˜ A) that returns ˜ xas the solution to ˜ A˜ x=˜ bwith a tolerance ǫKand a maximum number of iterations nmax. Now, we define a recursive multigrid function MA,ˆ b, mthat returns ˆ xas the solution to Amˆ x=ˆ bwith a tolerance ǫ. With this definition, it is clear that x=M(A,b,0) is the solution to the original problem. In this setting, operator Mis defined in details in Algorithm 1 in Appendix. Note that the maximum number of iterations of the iterative solver at the coarsest level, when m=mmax, is ˜nmax as opposed to nmax with ˜nmax ≫nmax. This is to ensure the solution to the coarsest problem converges properly before return. In other words, two later calls to Kare smoothers while the first call is for solving a linear system. Also, note the recursive call to the next level of the multigrid solver inside the while loop. To prevent an infinite loop, one may add a counter that calls an exit command in case of excessive iterations. 3. Partitioning approach The scalability of an iterative algorithm is primarily a function of its communication overhead and load balance. The communication overhead is directly proportional to granularity of the partitioning, which is the ratio between the number of cells shared between processors to the total number of cells. The load balance relies on an equal distribution of the number of cells between processors. Partitioning the grid to equal-size cubes with an equal number of cells in each direction is an optimal choice with regard to both the communication overhead and load balance. However, this is not possible for a given grid with an arbitrary number of cells in combination with utilization of an arbitrary number of processors. Hence, a compromise must be made between these two factors. To improve flexibility and increase the number of partitioning combinations, we create a partition-grid †with each of its block having different sizes. This is to obtain an algorithm that remains optimal for a wider combination of grids and number of processors. As a one-dimensional illustration, consider a grid with 10 cells. A partitioning strategy that only allows for equal-size partitions, only accepts the number of processors to be 1, 2, 5, or 10. In three dimensions, such a restriction significantly limits the range of possible partitioning options. With this limitation lifted, the average number of grid points in a direction, denoted by ni, can be a positive real number, i.e., ni∈R+∀i∈ {1,2,3}, and not necessarily an integer. †A partition-grid is a rectilinear grid of partitions that is constructed by dividing the computational grid into a certain number of slices in each direction. A scalable multigrid method for elliptic linear systems 217 To calculate ni, we first decompose the number of partitions, denoted by np, to prime numbers np=Y i pki i,(3.1) in which p={2,3,5,7,11,13,17,19}is a truncated sequence of prime numbers, and ki∈Nis the calculated repetition associated with each prime number. Given pand k, and the total number of grid points in each direction, |l|,nis calculated via Algorithm 2 in Appendix. For example, a 300 ×200 ×100 grid partitioned to np= 30 will produce n={60,66.667,50}, viz. partitioning x-direction to 5, y-direction to 3, and z-direction to 2. In y-direction, this leads to slices with a thickness of 67, 67, and 66 grid points, hence committing a negligible load imbalance of 1.5%. Each of the mmax grids that are generated is separately partitioned according to the algorithm described above. These independent partitions, specifically for a non-uniform grid, can be significantly misaligned at each level of the multigrid method. This has a direct implication on the implementation of Cand Fin terms of processor-to-processor communication. To simplify such implementation, we perform each of those operations in three steps. In the first step, considering Cfor example, values from the fine grid are aggregated into a buffer according to Eq. (2.9). Second, the buffer is communicated between processors, and third, it is assigned to the coarse grid. Using a buffer simplifies treatment of boundaries, specifically when multiple cells, each from a different partition, contribute to a single coarse grid cell. Finally, we note that in the case of large mmax and np, there can be excessive number of processors for solving a small system at the coarsest grid. To prevent this issue, we impose a lower bound on n(10 in our calculations), hence ensuring that there is always a minimum number of grid points per processor as the grid becomes coarser. Therefore, some processors can be idle when calculations are performed on the coarser grids, yet this has minimal effect on the overall performance due to the negligible cost of coarse grid calculations. Among others, this lower bound depends on the communication latency of the hardware, hence must be calibrated for each cluster. However, in fixing this parameter in all our calculations on various machines, we have not observed a significant change in the overall performance. 4. Results The parameters of the multigrid method, described in Section 2, are kept fixed for all the results reported in this section. The maximum number of restriction levels is mmax = 4. Increasing this number to 3, 5, or 6 has a minimal effect on the results. The tolerance of the solution is ǫ= 10−7. Decreasing or increasing this value will only shift the absolute computational costs, hence our final conclusion is not affected by this value. The maximum number of iterations at the coarsest level is ˜nmax = 500. This number is much larger than the actual number of iterations, hence, as long as a large value is used, no effect on the result is expected. The maximum number of iterations of the smoother is nmax = 8. The reported results can be moderately affected by this parameter; however, the change in the performance is minimal for 4 ≤nmax ≤10. The tolerance of the smoother is ǫK= 0.15. This parameter can also moderately affect the results; however, 0.15 is nearly optimal for a wide range of tested cases. Two sets of cases are studied in this section. The first is the direct numerical simulations of isotropic homogeneous turbulence in a cubic box. In these cases, the grid is 218 Esmaily-Moghadam et al. uniform in all directions and density (κin Eq. (2.5)) is constant. These cases test most of the components of the algorithm under simplified conditions. The second is the direct numerical simulation of a radiated particle-laden flow in a duct with square cross section. In these cases, the grid is not homogeneous in yand z-directions, in which cells near the wall have a much smaller liin comparison with those at the duct center. max(ly)/min(ly)≈15, ly/min(ly)≈7, and lx/min(ly)≈25. The domain is cuboid with an aspect ratio of 6 ×1×1. The number of grid points of the finest grid in x-direction is 1.5 times of that of yor z-direction. For each of these two classes of problem, we ensure that the flow is sufficiently evolved in time such that all wave numbers are present in the RHS of the linear system. All of the results reported here are obtained from calculations on Titan. Briefly, Titan contains 18,688 compute nodes, each containing a 16-core 2.2GHz AMD Opteron 6274 processor and 32 GB of RAM. Two nodes share a Gemini high-speed interconnect router. To benchmark our results, we repeated some of the above calculations using the Trilinos software project (Heroux et al. 2005). For these cases, the same linear systems employed for testing our algorithm were loaded to the Trilinos. To set up the problem, pre-existing template files, found in the library’s documentation, are used. The pre-installed Trilinos library on Titan (version cray-trilinos/11.12.1.3) is used in our calculations. The multigrid and CG solvers from Prokopenko et al. (2014) and Bavier et al. (2012) packages, respectively, are used for benchmarking our results. The adjustable parameters are tuned based on the recommended criteria. In particular, we set the problem type to Poisson3D for the multigrid solver, and we followed performance indications for parallel runs. Namely, the Chebyshev polynomial is used as the smoother with a V cycle pattern, as recommended by Trilinos for Poisson linear systems. With these options, as will be shown in the following, the Trilinos multigrid solver scaled poorly on more than 4 nodes (64 processors). We note that these results are obtained based on the best combination of parameters that we were able to find by exploring the parameter space. This does not, however, suggest that there is no combination of parameters that could produce better results. The direct numerical simulation of isotropic turbulence on a uniform 2563grid is performed using the present multigrid method and the results are compared against the CG (Figure 3) and the Trilinos multigrid solver (Figure 4). Based on Figure 3, the presented multigrid method is over an order of magnitude faster than the two implementations of the CG. In a practical sense, performing the same simulation without use of the multigrid technique entails approximately 10 times increase in cost. On a larger grid, this ratio will increase. At a higher number of processors, communications latency becomes a significant portion of the cost for the multigrid method. This leads to a lower cost ratio at higher number of processors in Figure 3. The higher and lower cost of CG-I in comparison with CG-T at a lower and higher number of processors, respectively, shows that our in-house implementation scales better in comparison with Trilinos’s. The scalability of the present multigrid method relies heavily on the scalability of its underlying linear solver, i.e., the CG-I. Therefore, the superior scalability of the CG-I also explains the better scalability of the present multigrid in comparison with that of Trilinos, which has led to a rising trend in Figure 4. Employing the present multigrid instead of Trilinos’s for computations with an intermediate number of processors (50 to 100 k cells in each partition) entails approximately 3 times reduction in cost. The strong scalability of the present multigrid on the same 2563grid is shown in A scalable multigrid method for elliptic linear systems 219 5 10 15 20 0 2 1 2 2 2 3 2 4 2 5 2 6 2 7 2 Number of nodes Number of cells/core (k) 524 262 131 66 33 16 81048 Cost ratio CG-I CG-T Figure 3. The ratio between the cost, defined as the total CPU hours, of CG to multigrid. Our in-house CG is CG-I (blue circles) and Trilinos’s CG is CG-T (red squares). Results are related to the isotropic turbulence simulation on a 2563grid. 1 2 3 4 0 2 1 2 2 2 3 2 4 2 5 2 6 2 Number of nodes Number of cells/core (k) 524 262 131 66 33 161048 Cost ratio Figure 4. The ratio between the cost of the Trilinos multigrid to the present multigrid. Results are related to the isotropic turbulence simulation on a 2563grid. Figure 5. An almost linear speedup is obtained up to 1012 processors, where there are 16,000 cells per partition. Compared to the case with 16 processors, the overall speedup is approximately 38 times. To show scalability of the present multigrid method on a wider range of grid sizes, we performed a weak scaling study in which the number of cells in each partition is held fixed as the number of partitions is increased. These results, which are also based on the simulations of isotropic turbulence on uniform grids, are presented in Figure 6. Two sets of results are presented in which the number of cells per partition is 25,000 and 50,000. The parallel efficiency of the case with the higher number of cells per partition remains above 60%, while for the other case, it drops to 40% at 4,048 processors. This drop in efficiency is mainly attributable to the latency of collective communications, which increases proportional to the number of processors. However, partitioning the grid with more than 50,000 cells per processor leads to a satisfactory scaling of the problems with size less than a billion grid points. The grid associated with the simulation run on eight nodes is refined by a factor of two in all directions in comparison with that of a single node in Figure 6. However, the parallel efficiency remains above 90% for the case of 50 k, indicating that the number of iterations is almost independent from the grid size. This independence is the main advantage of the