{'fraction_non_alphanumeric': 0.02517686225551394, 'fraction_numerical': 0.006361096248736698, 'mean_word_length': 4.930371937246607, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A GPU cluster optimized multigrid scheme for computing unsteady incompressible fluid flow

A multigrid scheme is proposed for the pressure equation of the incompressible unsteady fluid flow equations, allowing efficient implementation on clusters of modern CPUs, many integrated core devices (MICs), and graphics processing units (GPUs). It is shown that the total number of the synchronization events can be significantly reduced when a deep, 2h grid hierarchy is replaced with a two-level scheme using 16h-32h restriction, fitting to the the width of the SIMD engine of modern CPUs and GPUs. In addition, optimal memory transfer is also ensured, since no strided memory access is required. We report increasing arithmetic intensity of the smoothing steps when compared to the conventional additive correction multigrid (ACM), however it is counterbalanced in runtime by the decreasing number of the expensive restriction steps. A systematic construction methodology for the coarse grid stencil is also presented that helps in moderating the excess arithmetic intensity associated with the aggressive coarsening. Our higher order interpolated stencil improves convergence rate via minimizing spurious interference between the coarse and the fine scale solutions. The method is demonstrated on solving the pressure equation for 2D incompressible fluid flow: The benchmark setups cover shear driven laminar flow in cavity, and direct numerical simulation (DNS) of a turbulent jet. We have compared our scheme to the ACM in terms of the arithmetic intensity of the iterations and the number of the synchronization calls required. Also the strong scaling is plotted for our scheme when using a hybrid OpenCl/MPI based parallelization.


Introduction
Numerical studies of incompressible fluid flow are important in both academic research and engineering applications. The incompressibility constraint and constant density represent a good approximation when fluid flow velocity is significantly smaller than the speed of sound in the media. For subsonic unsteady flows this approach usually allows larger time-steps, because, contrary to the compressible fluid models, density waves travelling at the speed of sound do not need to be resolved. However, solving for the pressure instead of density changes the continuity equation from parabolic to elliptic type. The major computational challenge in simulating incompressible unsteady fluid flow is to develop an effective, parallel pressure equation solver [1].
The demand for more detailed descriptions of convection phenomena, such as microscale description of multiphase systems, large eddy simulation (LES) [2,3] or direct numerical simulation (DNS) [4], drives evolution of mathematical and numerical tools along with the evolving computer architectures. Recently, the share of accelerator cards (i.e. GPUs), and multi-core coprocessors is rapidly increasing in the supercomputing scene. Alongside, programming paradigms are changing to fully exploit the fine grained parallelism available on these hardwares. Writing GPU optimized codes for solving parabolic PDEs using explicit time marching is straightforward [5], while computationally efficient treatment of elliptic problems such as the Poisson equation together with incompressibility requires more complex methods. The difficulty lies in the multi-level parallelism of the current supercomputing hardwares. While efficient methods exist for large scale elliptic problems on distributed memory hardware, it is not straightforward to optimize these algorithms for both fine-grain and coarse-grain parallelism coincidentally [6,7]. Although multiple compute speed is available on recent accelerator cards such as GPUs, they are accompanied with significant latency when copying data between such devices. A simple synchronization between cards involves a device to host, a host to host, and a host to device memory transfer, which at least triples synchronization latency. It is not uncommon, that the ratio of the computational time and the memory transfer latency increases by even 2 or-ders of magnitude compared to CPU clusters. Therefore, in the present work we pay special attention to optimize our scheme for the case when latency limits the parallel efficiency of the computations.
Multigrid (MG) methods [8] are developed for to accelerate the iterative solution of large algebraic equations, such as the discretized pressure equation. The number of iterations required to reach the convergence criteria is reduced by using a multi-resolution discretization that accelerates the relaxation of long-wavelength components of the residual. Compared to simple iterative methods, MG has a low memory footprint [6], which is a significant benefit on accelerator devices having limited onboard memory. The MG methods traditionally use a deep hierarchy of discretization, where each grid is 2h coarser than the previous level. The scheme development for this multiresolution hierarchy needs special attention, since parallel execution of the iterations on fine meshes are usually bandwidth limited, while the solution of the coarser level discrete equations can be latency limited. Besides, it is not straightforward to achieve optimal fine grain parallelism when using these schemes. It is difficult to optimize memory load performance, cache usage and SIMD utilization at the same time, even if the programming model allows control over these procedures. The ease of programming is also an issue when evaluating a numerical procedure. Here, not only field variables, but also boundary conditions must be consecutively transfered through the grid hierarchy ensuring proper padding to keep data alignment. Despite challenges, various multigrid schemes have been successfully implemented on GPUs [9,10], and GPU clusters [11,12].
Here we present an aggressive coarsening strategy, that results in a simple but effective two-level multigrid scheme. Our scheme bypasses some of the programming difficulties arising when implementing to many-core hardware architectures, and decrease the memory footprint of the solver. Our method is compared to the Additive Correction Multigrid (ACM) method for unsteady incompressible fluid flows. Simulation of a laminar flow in a shear driven cavity, and a direct numerical simulation of a turbulent jet have been used to benchmark against the ACM scheme. Finally we present strong scaling figures for a channel-flow setup, that mimics 2D turbulence in soap film experiments [13].
2. The fluid model and its discretization

The equations of the incompressible fluid
Three assumptions on fluid behavior yields the Navier-Stokes equation: (i) The viscous dissipation is a linear function of the strain rates, (ii) the fluid is isotropic (rotational invariance), and (iii) for a fluid at rest hydrostatic pressure applies. Assuming incompressibility and constant density, the Navier-Stokes equation and the continuity of mass read as follows: where equation (1) is solved for the fluid velocity v, and equation (2) is solved for the pressure p when substituting v [see Eq. (4)]. η is the kinematic viscosity.

The time stepping scheme
A simple variant of Chorin's projection method [14], the Goda's incremental pressure correction method [15] has been chosen to solve the governing equations. We predict velocity v * for the next time-step using a simple first order time integration of the explicitly known terms. The pressure can be decomposed as: p t+1 = p t + δp. Using p t the predicted velocity reads as: Substituting v t+1 = v * + ∇δp into Eq. (2) yields the following Poison equation: which has to be solved for the pressure change δp. Finally v * is corrected with the pressure-change term ∇δp to obtain the divergence free v t+1 .

The spatial discretization
A staggered grid arrangement [16] has been chosen for the velocity to avoid odd-even decoupling and the resulting spurious checkerboard patterns in the solution. Finite differences of second order accuracy were used to discretize the Navier-Stokes equation [17], and a five point Laplacian stencil was used to discretize the pressure equation. Alternatively, the discretization can be derived from finite volumes (FV) that gives a clear explanation of constructing flux preserving conservative stencils.

Simple methods
While the explicit time marching scheme for the Navier-Stokes equation can be solved line by line, discretizing the continuity equation leads to a large system of linear equations. Direct solution of such systems is computationally too expensive for most practical problems, and also difficult to perform in a parallel manner. The discretization error is usually larger than the accuracy of the computer arithmetic, thus an approximate solution of the discrete problem is satisfactory, and an iterative solution can be applied. In parallel environment, the Gauss-Sedel (GS) iteration is commonly applied. Concurrent computing of the well known 5 point stencil is feasible when using redblack reordering [18], while multicolor ordering [19] can be used for higher order stencils.

The multigrid method
The GS iterations are very effective in removing the high frequency components of the error, but the low frequency components decay with a very low rate. Convergence requires O(N 2 ) iterations, where N is the linear size. Multigrid methods are commonly used to accelerate convergence: coarser discretizations are used to eliminate the lower frequency components of the residual. The MG technique consists of a sequence of smoothing, restriction, and prolongation operators. Smoothing steps are the actual iterations using a solver (i.e. Gauss-Seidel) on a grid level. Restriction is a downsampling the residual to a coarser discretization, and prolongation is an interpolation from a coarser to finer discretization. A typical multigrid cycle starts with smoothing of the residual error on the finest level, then the approximate solution is transfered to a coarser level (restriction), and smoothened again. The procedure is repeated until the coarsest level is reached, then the way back. The optimal depth of the grid hierarchy depends on the actual problem, the discretized equations must be solved in a reasonable time on the coarsest level. A direct solver is often used at the coarsest grid.
When using a geometric multigrid (GMG) scheme, the original PDE is re-discretized on each grid, data structures must be constructed on each level. Note that the solutions on coarse grids does not necessarily approximate the solution on the fine grid, but it approximates the original problem. Possible inconsistency between fine and coarse grid solutions may results in extra iteration cycles before reaching convergence.
When using an algebraic multigrid (AMG) scheme, agglomeration of the coefficient matrix is performed. In the case of the widely used AMG procedure, the Additive Correction Multigrid scheme [20], the coarse grid equations are obtained by the summation of the finer grid equations. The results obtained on the coarser grid are simply added to the finer grid solution, no interpolation or extrapolation operators are required in this procedure. The real advantage of this method is that on each level it approximates the solution of the discrete problem on the finest grid. Also note that conserving property of the spatial discretization scheme is satisfied at all discretization level.

The construction of an interpolated stencil for the coarse mesh.
We extend the basic idea of the ACM scheme in our MG procedure. It is often useful to interpret the discretization of the conservation equations (either momentum or mass) in the finite volume manner. When constructing a FV scheme, both sides of the PDE are integrated over a control volume (CV). Applying the Gauss theorem on the homogeneous part of the Poisson equation results in a surface integral of the fluxes over the CV faces. Thus simplest discretization is obtained, when the inhomogeneous part is assumed to be the CV average, and the surface normals of the pressure gradients are assumed to be constant at all CV faces. This recovers the well-known five point Laplacian stencil. The overall accuracy of the respective FV scheme depends both on the accuracy of calculating the fluxes on the CV face (using finite differences), and the accuracy of the surface integral [21].
After restricting the fine grid values to a coarser grid cell (hatched area in Fig. 1.), the same finite volume discretization procedure applies when the traditional geometric multigrid (GMG) method is used. However, one can use many values that are already pre-calculated on the fine mesh to make the discretization more accurate and easier to compute. For example, instead of downsampling v * to the coarse mesh and calculate its divergence using the corresponding finite volume stencil, one can sum ∇v * already calculated for the fine grid discretization. In the latter case the integration for the coarse CV will be the same accuracy as for the fine grid. The procedure is identical to the agglomeration of the coefficient matrix in the ACM method. Naturally, the integration (summation) can be applied for larger volumes too (a 8h restriction is shown in Fig. 1). Moreover one can perform similar sum for the fluxes derived from applying the Gauss theorem for homogeneous part of the PDE. It is clear that fluxes not crossing the coarse grid cell faces cancel each other in the volumetric integral (summation), therefore, we have to account for the discretization only at the cell faces. When applying Coarse grid is drawn by solid black lines, some cell centers (P C ,P E ,P N and P N E ) are denoted with shaded circles. The actual center node for the coarse grid is hatched. The interpolation area for the north-east corner of the center node is indicated using dashed lines, and fine grid cells whose interpolated values used to construct the coarse grid discretized equations are shaded. White arrows show pressure gradients that are summed for the coarse grid discretization.
to domains of arbitrary size, an interpolation rule must be constructed to define fine grid pressure values p i,j near the boundaries. Thus the finite volume scheme for the coarse grid is constructed by summing fluxes calculated from the p i,j interpolated values. At this point the method differs from the ACM scheme, where a simple summation is performed that corresponds to assuming constant p i,j in the entire coarse CV. The neighboring coarse cells must use the same interpolated values to maintain the conservative character. This can be easily ensured, if these cells use a common domain to interpolate over. A natural choice is to interpolate between four neighboring cell centers (dashed rectangle in Fig. 1.). We have used a simple bilinear interpolation to calculate p i,j as the function of pressure values at the coarse cell centers P i,j . Naturally, the same interpolation was used to prolongate the pressure values to the fine grid. If the origin is (Q 11 ) the bilinear interpolation reads as: where ∆x W |E and ∆y S|N are the horizontal, and the vertical size of the interpolation rectangle, while x and y are the coordinates of the fine grid control volume centers.
The superscripts indicate that the size of the interpolation rectangle can differ on the west-east and the south north faces, because of padding, whereas the coarse mesh can be nonuniform. Q I,J are pressure values at the corners. Now, the fluxes across fine CV faces S i and S j (white arrows in Fig. 1.) can be computed as follows: ∂p(x, y) ∂x and ∂p(x, y) ∂y Summing up the above fluxes gives the discretization scheme for the coarse mesh. We note that the simple five point stencil used on the fine mesh extends to a 9 point stencil on the coarse mesh. In practice the proper handling of the sum on the padded domain requires some integer arithmetics. Note that the above summation leads to a nine point stencil.

Multigrid strategy
Multigrid strategy including cycling (i.e. V, W or Fcycling) and the coarsening strategy is a subtle issue, and often examined in term of the arithmetic complexity [22,23]. The effect of cycling strategy is also known to have effect on parallel performance of the scheme [24], however the coarsening strategy is rarely examined in this respect. It is known, that not all discretization levels are equally important in reducing residual errors, the optimal coarsening may not even uniform [25]. Herein we wish to investigate the effect of a more aggressive coarsening strategy, in terms of the computational complexity and the parallel scalability in the latency limit. Fig. 2 shows a comparison of a common 2h V cycle and its two level "sparse" counterpart. Note that in the latter case, the number of restriction and prolongation steps are greatly reduced. In the next section we are going to investigate the effect of substituting common V cycles with sparse ones.

Implementation and benchmarking
Our ACM solver was implemented using NVIDIA Cuda and supports only a single device. Compute tasks for all grid levels run largely on the graphics card, except when global reduction is needed (i.e. convergence check). We use the five point Laplacian stencil on the fine grid, and the coarser grid discrete equations are generated by summing fine grid equations over the coarser grid, that also gives a 5 point stencil. Iterations are executed parallel by using a red-black Gauss-Seidel (GS) method. Reductions are executed in two steps starting on graphics board, and partially reduced data are transferred to the main memory to finish reduction.
In our improved ISMG solver all computations on the fine grid are implemented using the OpenCL streaming language and executed on accelerator cards accordingly. Our ACM implementation may also use multiple graphics cards over MPI protocol. While the fine grid equations are solved on the graphics boards, coarse grid equations are implemented in standard C language, and run on the host machine(s). This adds less latency when iterating on the small sized coarse grid, and avoids poor utilization of the graphics card. We note, that the free time-slot on the GPU can be utilized for arithmetically intensive pre or post-processing, or for solving coupled scalar equations (i.e. convenction-diffusion equation) that is common in many fluid dynamics applications. As for the ACM solver we use parallel red-black Gauss-Seidel iterations on the fine grid, but we use a block iterative GS method [26] on coarse grid level. GS iteration blocks are aligned with the MPI parallel layout, and simplifies to the serial GS when running on a single node. The coarse grid solver was implemented as host code, using C language. We note however, that using parallel multicolor iteration schemes [19], are also an option, which may increase computing efficiency of the CPU device on the expense of extra synchronization, and may increase latency.
The OpenCL language allows an efficient use of SIMD engines via defining blocks of operations (work-groups) that are executed asynchronously. Work-groups can be either one, two or three-dimensional, thus explicit finite difference stencils are straightforward to implement; workgroups are paired with blocks of gridded data (data "tiles" in 2D) [27]. However, we have to note that the implementation in 2D and 3D may differ because of the limited number of registers, and out of order execution is not commonly available on graphics hardware. For the sake of simplicity and easier programming we have chosen a 2D implementation. Despite our results are in 2D, the extrapolation to estimate 3D performance can be done since convergence rate for GS iterations are depend on the linear size of the simulation domain. In our software code each workgroup corresponds to a coarse-grid cell, thus code complexity is significantly reduced. The best computing efficiency can be achieved, when the tiles/work-groups are fitted to the SIMD width of the compute device. Using an optimal tile size on GPUs, the amount of data to be processed on the coarse grid is greatly reduced (i.e., by a factor of 256 for a 16 × 16 tile size).
Code performance and our multigrid scheme was analyzed on various hardware platforms. We wish to emphasize that we have not put extra efforts in optimizing for specific platforms. Our aim was to measure the effi-ciency of the MG scheme on a more generic way, and give estimates on possible bottlenecks.

OpenCl kerel performance
Our OpenCl compute kernel implementations use global memory, except reduction kernels (i.e. restrictions steps in MG) which also use local address space. This means that the kernels are quasi-optimal for CPUs, and further speedup is possible on graphics processors. The reason for using global memory access is its simplicity. We wanted to use exactly the same compute kernels on each platform, and we have not found a generic way to utilize the full potential of the specific hardwares. Thus, our kernel timings can easily be achieved by a scientist or fluid dynamics expert, who is familiar with GPU programming at a basic level. Implementations of the finite difference stencils may also differ in flexibility of handling boundary conditions and grid size. Flexibility of the code often leads to conditional statements, which may dramatically affect code performance. In practice, conditional statements can mostly be replaced with integer, or bit-wise integer arithmetics, but leads to extra operations. We also note that architectures and platforms may significantly differ in performing integer operations that may result in differences in computing time for codes with flexible boundary conditions. Our kernel implementation uses automated padding if the grid size is not a multiple of the tile size, and can hold symmetric, Dirichlet and periodic boundary conditions. We use wall times of kernel execution as the measure of the performance. The results benchmarked on some recent hardware are summarized in Table 1. According to the benchmarks, execution times are decreasing as work-group size increases on all hardware and saturates at the SIMD width. Its variation can be significant on GPUs. We note that with increasing restriction length we have double benefit; not just the computational cost of solving an iteration on the coarse mesh is decreasing, but also iterating on the fine grid becomes more efficient.

Accuracy test
In order to test accuracy and the righteousness of our implementation we have chosen the well-known lid-driven cavity test first published by Ghia [8]. Our computations were executed until steady solution was reached. The extrema of the velocities along the centerlines are in reasonable agreement with the values reported in previous works (see table 2.). We also note that deviations among the reported values are within the 32 bit floating point precision (after approx. 40000 non-dimensional time), indicating that our implementation using single precision is reasonable.

Multigrid performance
We have measured the computational cost (arithmetic intensity) of the MG scheme using the average number of computed Laplacians corresponding to a fine grid control volume in a time increment (N Lap ). This can be used to compare different MG schemes in an ideal case, where the network communication and the memory bandwidth are not limiting the calculation.
Although the computational cost of performing an iteration is an essential feature in evaluating a numerical method, the number of iterations required to reach convergence criteria is more fundamental, when evaluating scalability on large distributed memory computers. While the raw computational power can be increased by stacking computers, the interconnect latency and bandwidth have strict limitations. When performing finite difference iterations, the communication can be hidden up to the point, when compute time decreases to the network communication time. Down to this limit, the compute power can be fully utilized: theoretically the speedup is proportional to the number of processors. When pushing to this limit (i.e. by increasing the number of the compute devices), the total compute time can be reduced only via decreasing the iteration count, if we maintain computing efficiency. We use the total GS iteration count (I t ) and the related total network communication count (synchronizations) per time step (N CC t ) as the measure of performance on distributed memory clusters. We distinguish iteration numbers on fine (I f ) and coarse grids (I c ) since they behave differently with varying numerical parameters, and network communication characteristics also varies with grid size. The actual choice for the solver also influences network communication requirements: the simple block iterative GS involves 1 synchronization per iteration that doubles when red-black reordering is applied. Besides, when compute devices are external cards, an extra, compute device to host transfer is also involved. In this section we use laminar flow in a rectangular cavity as a benchmark. To make tests more challenging, we put a twist into the commonly used lid-driven cavity by prescribing opposite velocities on the fixed walls. The setup and the numerical parameters are shown in Fig.3. When using advanced multigrid methods, often complicated patterns are used to restrict and prolongate among various grid-levels, such as v-cycles and w-cycles. Here, the simplest accommodative cycles strategy is used: Iteration is started on the coarsest grid, and when the absolute maximum of the residual is decreased below a limit, the solution is prolongated to the finest grid. We can vary this criterion to find a quick solution. Then, iteration continues on the fine grid until reaching the fine grid convergence criterion. In our benchmark we require the residual to decrease below 10 −6 on the fine grid, and vary the coarse grid criterion. Reaching this criterion requires significant efforts, but we are still comfortably far from the limits of the 32 bit float arithmetics.  We have compared our results to the ACM scheme [32] that is using the conventional 2h restriction length. On the coarsest mesh, we have applied the same convergence criterion as for the ISMG scheme. We have performed heuristic tests to find a quasi-optimal cycle: V cycles without pre-smoothing, and one post-smoothing per grid level was found to be a good candidate. The grid depth was varied, thus we have the ACM substitutes for our ISMG scheme using varying restriction length (e.g. tile size). In the case of our ACM implementation we denote only the coarsest mesh iterations as "coarse", and all the others are accounted as "fine". Iteration counts at each grid level can be calculated, since the number of the smoothing steps are the same, except for the coarse mesh.
We have measured scalability (N CCs) and arithmetic intensity (N Lap ) that are summarized in table A.3 for the ISMG scheme and in Table A.4 for the ACM scheme, respectively. The optimum values we have found are highlighted with bold characters. We have measured similar computational cost while our ISMG method requires six times less synchronizations. We note that the greatly reduced synchronization count also results in less kernel launches that is preferable in streaming computer languages. Besides, restriction and prolongation kernels are executed fewer times, however it is difficult to make a direct comparison on this basis, since it depends largely on the multigrid cycling pattern.
The low iteration counts we measured are partly due to the fact that in the ISMG scheme the coarse grid solution approximates the fine grid discrete solution. Fig.  4.b shows that in the moderate transient regime after prolongating to the fine grid the solution occasionally passing convergence criteria with no iterations. The same can be observed the way back, sometimes no coarse iteration is required after a converged solution of the previous time-step is applied in the predictor step. The previous test setup can be regarded as a synthetic benchmark, since the pressure we calculate is deterministic, therefore, one can tune the setup to provide unbiased benchmarks for the actual scheme (e.g. by varying the grid-size). In this section, we present benchmarks for Direct Numerical Simulation (DNS) of a turbulent jet. The actual setup, and the numerical parameters applied are shown in Fig. 5. Due to the chaotic nature of turbulence, it would be very difficult to tweak the test setup to show unbiased benchmarks, it can be accounted as a real-life application of our scheme. We have increased the grid size to 1000 × 2000. Here we have maximized the number of GS iterations at 20000 for both ACM and ISMG computations to prevent convergence stagnation, however we believe this can be bypassed by using a more advanced cycling strategy. The linear sizes we have used in this setup are close to the supercomputing scale 3D direct numerical simulations.
Due to the increase in the linear size the small grid depth computations for ACM scheme and the corresponding 4 × 4 and 8 × 8 tile sizes in the ISMG scheme the solution is far from optimal and requires excessive compute time, therefore results are not presented for these cases. The results for 5-6 grid levels of the ACM and the corresponding 16 × 16 and 32 × 32 tile size ISMG schemes are presented in Table A.5. For comparison we also have made a benchmark for this setup using an aggressive coarsening strategy, but with the same five point coarse grid stencil as is usual in geometric multigrid (GMG) schemes. Results (last two columns of Table A.5) provide insight how arithmetic intensity decreases when using our interpolated stencil methodology.
In the case of the ACM scheme the optimum has been found at 6 grid levels contrary to the ISMG scheme, where the optimum has been found at 16 × 16 tile size (corresponds to 5 grid levels of the ACM method). Comparing the optima of the two methods, the ISMG scheme maintained the sixfold advantage in scalability (N CC t ), but at the cost of approximately twofold increase in the computational need (N Lap ). We note that in case of the applied cycling strategy this computational overhead is compensated by reduced restriction and prolongation counts, and their more efficient execution. When comparing the interpolated stencil to the conventional five point Laplacian, we have measured 26% reduction in the arithmetic intensity of the iterations.

Benchmark setup III.: Strong scaling
Finally, we asses multi-GPU performance and efficiency of our algorithm for different problem sizes. The benchmark setup is a channel flow having equally spaced jets from a side. Each jets having width of 25h, (just as for benchmark II), and their spacing is 200h. Above construction allows a rough comparison of different problem sizes. However the solution is inherently stochastic, averaging the first 20000 time-step gives a reasonable basis for the comparison.
Our tests has been performed on two nodes totaling eight graphics processors, and connected with Mellanox SDR infinband. The detailed configuration for a nodes is: CPU: Dual Intel Xeon E5530 (2.4 GHz), 2 x Nvidia GTX 590 boards (4 x GF110 GPU per node). Our benchmarks were run under CentOS 6.5, and using Cuda 5.5 along with Nvidia driver version 319.37.
The benchmarks are summarized in Fig. 6. The largest benchmark (15360x8640) saturates the 1.5 GB memory per GPU, therefore strong scaling curves showe for 4-8 graphics cores. For large enough setups, we have measured up to 80% efficiency at 8 graphics cores. Losing efficiency comes partly from gathering communication when checking convergence, and partly from transferring the coarse grid problem forth and back between the GPU memory and the main memory. The latter can be hidden, if some coupled phenomenon (i.e. thermal convection and diffusion) are needed to solve along with the fluid flow. We also note, that in three dimensions, the efficiency curve is expected to be even shallower due to the more rapid decrease in the problem size to coarse meshes.

Concluding remarks
We have presented an interpolated stencil multigrid method that is efficient when 16-32h wide restriction step is applied. The stencil allows the replacement of a deep grid hierarchy (optimally 5-6 grid levels) with a two grid scheme, while keeping the number of iterations low. The scheme was benchmarked against an ACM scheme with deep grid hierarchy. It is shown that ISMG needs significantly less iterations compared to a conventional ACM scheme, since it involves less synchronizations (i.e. kernel calls and network communications). Also, at larger grid sizes the raw computational cost (the number of the computed Laplacians) is increased by a factor of two, however it is overcompensated by the reduced number of prolongation and restriction kernel calls. We have measured the parallel efficiency near 80% for a cluster includes 8 graphics cores. For further fine tuning our scheme, we plan to replace the simple one-step accommodative cycling strategy with more complex cycling patterns to further decrease computational cost. Using our method with higher order interpolation is also an issue that we wish to investigate. Since Gauss-Seidel has convergence rate of O(N 2 ), where N is the linear size of the grid, we expect even faster convergence in three dimensional cases. Besides, in 3D the coarse graining results in faster reduction of the coarse problem size, thus helping the more effective use of the parallel hardware and better strong scaling.