Parallel computation of unsteady three-dimensional incompressible viscous flow using an unstructured multigrid method

The development and validation of a parallel unstructured tetrahedral non-nested multigrid (MG) method for simulation of unsteady 3D incompressible viscous flow is presented. The Navier-Stokes solver is based on the artificial compressibility method (ACM) and a higher-order characteristics-based finite-volume scheme on unstructured MG. Unsteady flow is calculated with an implicit dual time stepping scheme. The parallelization of the solver is achieved by a MG domain decomposition approach (MG-DD), using the Single Program Multiple Data (SPMD) programming paradigm. The Message-Passing Interface (MPI) Library is used for communication of data and loop arrays are decomposed using the OpenMP standard. The parallel codes using single grid and MG are used to simulate steady and unsteady incompressible viscous flows for a 3D lid-driven cavity flow for validation and performance evaluation purposes. The speedups and efficiencies obtained by both the parallel single grid and MG solvers are reasonably good for all test cases, using up to 32 processors on the SGI Origin 3400. The parallel results obtained agree well with those of serial solvers and with numerical solutions obtained by other researchers, as well as experimental measurements.     

Domain decomposition based coupling between the lattice Boltzmann method and traditional CFD methods – Part II: Numerical solution to the backward facing step flow

The lattice Boltzmann method (LBM) and traditional finite difference methods have separate strengths when solving the incompressible Navier–Stokes equations. The LBM is an explicit method with a highly local computational nature that uses floating-point operations that involve only local data and thereby enables easy cache optimization and parallelization. However, because the LBM is an explicit method, smaller grid spacing requires smaller numerical time steps during both transient and steady state computations. Traditional implicit finite difference methods can take larger time steps as they are not limited by the CFL condition, but only by the need for time accuracy during transient computations. To take advantage of the strengths of both methods, a multiple solver, multiple grid block approach was implemented and validated for the 2-D Burgers’ equation in Part I of this work. Part II implements the multiple solver, multiple grid block approach for the 2-D backward step flow problem. The coupled LBM–VSM solver is found to be faster by a factor of 2.90 (2.87 and 2.93 for Re = 150 and Re = 500, respectively) on a single processor than the VSM for the 2-D backward step flow problem while maintaining similar accuracy.     

From multi- to single-grid CFD on massively parallel computers: Numerical experiments on lid-driven flow in a cube using pressure–velocity coupled formulation

A parallel implementation of a fully pressure–velocity coupled multigrid solver based on analytical solution accelerated coupled line Gauss Seidel (ASA-CLGS) smoother with grid partitioning is carried out. The parallelized algorithm is characterized by an enhanced scalability that results from a formulation enabling an intermediate analytical solution for the entire row (column) of control volumes. General strategies of applying single- or multigrid approach depending on flow characteristics are discussed. Performance of the parallelized algorithm is studied for up to 2048 processors. The developed approach is applied to analysis of a time-dependent three-dimensional incompressible lid-driven cavity flow. The steady state results of benchmark quality are reported for Re = 10³, 1.5 × 10³ and 1.9 × 10³. A new benchmark case of a fully 3D flow in a cubic cavity driven by the lid moving at 45° relatively to its lateral boundaries is proposed and the corresponding data is reported.     

A parallel multigrid solver for incompressible flows on computing architectures with accelerators

An efficient parallel multigrid pressure correction algorithm is proposed for the solution of the incompressible Navier–Stokes equations on computing architectures with acceleration devices. The pressure correction procedure is based on the numerical solution of a Poisson-type problem, which is discretized using a fourth-order finite difference compact scheme. Since this is the most time-consuming part of the solver, we propose a parallel pressure correction algorithm using an iterative method based on a block cyclic reduction solution method combined with a multigrid technique. The grid points are numbered with respect to the red–black ordering scheme for the parallel Gauss–Seidel smoother. These parallelization techniques allow the execution of the entire simulation computations on the acceleration device, minimizing memory communication costs. The realization is developed using the OpenACC API, and the numerical method is demonstrated for the solution of two classical incompressible flow test problems. The first is the two-dimensional lid-driven cavity problem over equal mesh sizes while the other is the Stokes boundary layer, which is a decent benchmark problem for unequal mesh spacing. The effect of several multigrid components on modern and legacy acceleration architectures is examined. Eventually the performance investigation demonstrates that the proposed parallel multigrid solver achieves an acceleration of more than 10\(\times \) over the sequential solver and more than 4\(\times \) over multi-core CPU only realizations for all tested accelerators.     

Numerical simulation of 3D fluid–structure interaction flow using an immersed object method with overlapping grids

The newly developed immersed object method (IOM) [Tai CH, Zhao Y, Liew KM. Parallel computation of unsteady incompressible viscous flows around moving rigid bodies using an immersed object method with overlapping grids. J Comput Phys 2005; 207(1): 151–72] is extended for 3D unsteady flow simulation with fluid–structure interaction (FSI), which is made possible by combining it with a parallel unstructured multigrid Navier–Stokes solver using a matrix-free implicit dual time stepping and finite volume method [Tai CH, Zhao Y, Liew KM. Parallel computation of unsteady three-dimensional incompressible viscous flow using an unstructured multigrid method. In: The second M.I.T. conference on computational fluid and solid mechanics, June 17–20, MIT, Cambridge, MA 02139, USA, 2003; Tai CH, Zhao Y, Liew KM. Parallel computation of unsteady three-dimensional incompressible viscous flow using an unstructured multigrid method, Special issue on “Preconditioning methods: algorithms, applications and software environments. Comput Struct 2004; 82(28): 2425–36]. This uniquely combined method is then employed to perform detailed study of 3D unsteady flows with complex FSI. In the IOM, a body force term F is introduced into the momentum equations during the artificial compressibility (AC) sub-iterations so that a desired velocity distribution V₀ can be obtained on and within the object boundary, which needs not coincide with the grid, by adopting the direct forcing method. An object mesh is immersed into the flow domain to define the boundary of the object. The advantage of this is that bodies of almost arbitrary shapes can be added without grid restructuring, a procedure which is often time-consuming and computationally expensive. It has enabled us to perform complex and detailed 3D unsteady blood flow and blood–leaflets interaction in a mechanical heart valve (MHV) under physiological conditions.     

A non-intrusive parallel-in-time adjoint solver with the XBraid library

In this paper, an adjoint solver for the multigrid-in-time software library XBraid is presented. XBraid provides a non-intrusive approach for simulating unsteady dynamics on multiple processors while parallelizing not only in space but also in the time domain (XBraid: Parallel multigrid in time, http://llnl.gov/casc/xbraid). It applies an iterative multigrid reduction in time algorithm to existing spatially parallel classical time propagators and computes the unsteady solution parallel in time. Techniques from Automatic Differentiation are used to develop a consistent discrete adjoint solver which provides sensitivity information of output quantities with respect to design parameter changes. The adjoint code runs backwards through the primal XBraid actions and accumulates gradient information parallel in time. It is highly non-intrusive as existing adjoint time propagators can easily be integrated through the adjoint interface. The adjoint code is validated on advection-dominated flow with periodic upstream boundary condition. It provides similar strong scaling results as the primal XBraid solver and offers great potential for speeding up the overall computational costs for sensitivity analysis using multiple processors.

     

Application of single-level,pointwise algebraic,and smoothed aggregation multigrid methods to direct numerical simulations of incompressible turbulent flows

Single- and multi-level iterative methods for sparse linear systems are applied to unsteady flow simulations via implementation into a direct numerical simulation solver for incompressible turbulent flows on unstructured meshes. The performance of these solution methods, implemented in the well-established SAMG and ML packages, are quantified in terms of computational speed and memory consumption, with a direct sparse LU solver (SuperLU) used as a reference. The classical test case of unsteady flow over a circular cylinder at low Reynolds numbers is considered, employing a series of increasingly fine anisotropic meshes. As expected, the memory consumption increases dramatically with the considered problem size for the direct solver. Surprisingly, however, the computation times remain reasonable. The speed and memory usage of pointwise algebraic and smoothed aggregation multigrid solvers are found to exhibit near-linear scaling. As an alternative to multi-level solvers, a single-level ILUT-preconditioned GMRES solver with low drop tolerance is also considered. This solver is found to perform sufficiently well only on small meshes. Even then, it is outperformed by pointwise algebraic multigrid on all counts. Finally, the effectiveness of pointwise algebraic multigrid is illustrated by considering a large three-dimensional direct numerical simulation case using a novel parallelization approach on a large distributed memory computing cluster.     

A new class of massively parallel direction splitting for the incompressible Navier–Stokes equations

We introduce in this paper a new direction splitting algorithm for solving the incompressible Navier–Stokes equations. The main originality of the method consists of using the operator (I ? ?_xx)(I ? ?_yy)(I ? ?_zz) for approximating the pressure correction instead of the Poisson operator as done in all the contemporary projection methods. The complexity of the proposed algorithm is significantly lower than that of projection methods, and it is shown the have the same stability properties as the Poisson-based pressure-correction techniques, either in standard or rotational form. The first-order (in time) version of the method is proved to have the same convergence properties as the classical first-order projection techniques. Numerical tests reveal that the second-order version of the method has the same convergence rate as its second-order projection counterpart as well. The method is suitable for parallel implementation and preliminary tests show excellent parallel performance on a distributed memory cluster of up to 1024 processors. The method has been validated on the three-dimensional lid-driven cavity flow using grids composed of up to 2 × 10⁹ points.     

Parallel preconditioned WENO scheme for three-dimensional flow simulation of NREL Phase VI Rotor

The preconditioned weighted essentially non-oscillatory (P-WENO) solver for viscous flows (Huang et al. (2009) [9]) is extended to non-inertial reference frames. In the present scheme, patched multi-block grid system is employed and parallel computing is adopted as well. With the present parallel P-WENO solver, three-dimensional flows of the Phase VI Rotor from National Renewable Energy Laboratory (NREL) can be simulated and analyzed for different wind speeds. Our simulation results show good agreement with the numerical predictions based on incompressible Navier–Stokes (N–S) equations as well as the available wind tunnel data from NREL. The flow phenomena, including separation and attachment line, can be captured by the present scheme. The parallel strategy adopted is a block-domain decomposition method for the patched multi-block grid system. To balance the load among different computing nodes, a Tabu search algorithm is adopted for the parallelization. The parallel efficiency of the parallel P-WENO scheme is examined for node numbers ranging from 1 to 64. It is found that the parallel efficiency is monotonically decreased as the node number adopted is increased; the parallel efficiency is retained over 90% for all cases of different node numbers. Due to the high parallel efficiency, our parallel P-WENO solver is validated for applying to practical fluid problems from compressible to incompressible limits.     

多重网格技术对SIMPLER算法的加速性能研究

将多重网格技术引入SIMPLER算法以加快其收敛速度,从而节约计算时间。通过计算不同雷诺数下的二维方腔顶盖驱动流,研究了多重网格方法中的V循环、W循环对SIMPLER算法的加速效果,并讨论了网格层数对加速性能的影响。研究结果表明,在不同雷诺数下,多重网格方法均可以起到良好的加速效果;在相同雷诺数和精度要求下,W循环方式的外迭代次数少于V循环方式的外迭代次数,而且网格层数对多重网格加速性能的影响并不显著。     

Towards a three-dimensional moving body incompressible flow solver with a linear deformable model

In this study, a three-dimensional fluid–structured parallelized solver is extended from the previous work (Niu et al., 2009 [1]) for moving body simulations. Based on the unified Eulerian and Lagrangian coordinate transformations, the unsteady three-dimensional incompressible Navier–Stokes equations with artificial compressibility (Chorin, 1967 [2]) in a dual-time stepping approach are first derived. To implement unsteady flow calculations, the dual-time stepping strategy including the LU decomposition method is used in the pseudo-time iteration and the second-order accurate backward difference is adopted to discretize the unsteady flow terms. Also, a third-order Roe type flux limited splitting is derived to evaluate the spatial difference of the convective fluxes. The original FORTRAN code is converted to the MPI code and tested on a 64-CPU IBM SP2. The parallel strategy here is based on the partitions of all do-loops in the original FORTRAN code and transferring the calculations inside the do-loop into different CPUs. The partition of the do-loop can be applied on the innermost loop, only or the last two inner loops depending on two-dimensional or three-dimensional problems. This kind of the parallel data partition of the loops is independent of what kind of the explicit or implicit type numerical algorithm used. Therefore, the current parallel approach can take advantage of the MPI language fully to transfer data efficiently among CPUs even for solving the governing equation implicitly. The test results show that a significant reduction of computing time in running the model and a near-linear speed up rate is achieved up to 32 CPUs at IBM SP2. The speed up rate is as high as 31 for using 64 IBM SP2 processors The test shows efficient parallel processing to provide prompt simulation of 3D cavity, unsteady dropping airfoil and blood flows in an aortic tube with a linear elastic modeling of wall motion is included here.     

Implementation of density-based solver for all speeds in the framework of OpenFOAM

In the framework of open source CFD code OpenFOAM, a density-based solver for all speeds flow field is developed. In this solver the preconditioned all speeds AUSM+(P) scheme is adopted and the dual time scheme is implemented to complete the unsteady process. Parallel computation could be implemented to accelerate the solving process. Different interface reconstruction algorithms are implemented, and their accuracy with respect to convection is compared. Three benchmark tests of lid-driven cavity flow, flow crossing over a bump, and flow over a forward-facing step are presented to show the accuracy of the AUSM+(P) solver for low-speed incompressible flow, transonic flow, and supersonic/hypersonic flow. Firstly, for the lid driven cavity flow, the computational results obtained by different interface reconstruction algorithms are compared. It is indicated that the one dimensional reconstruction scheme adopted in this solver possesses high accuracy and the solver developed in this paper can effectively catch the features of low incompressible flow. Then via the test cases regarding the flow crossing over bump and over forward step, the ability to capture characteristics of the transonic and supersonic/hypersonic flows are confirmed. The forward-facing step proves to be the most challenging for the preconditioned solvers with and without the dual time scheme. Nonetheless, the solvers described in this paper reproduce the main features of this flow, including the evolution of the initial transient.     

PetRBF — A parallel O(N) algorithm for radial basis function interpolation with Gaussians

We have developed a parallel algorithm for radial basis function (rbf) interpolation that exhibits O(N) complexity, requires O(N) storage, and scales excellently up to a thousand processes. The algorithm uses a gmres iterative solver with a restricted additive Schwarz method (rasm) as a preconditioner and a fast matrix-vector algorithm. Previous fast rbf methods — achieving at most O(NlogN) complexity — were developed using multiquadric and polyharmonic basis functions. In contrast, the present method uses Gaussians with a small variance with respect to the domain, but with sufficient overlap. This is a common choice in particle methods for fluid simulation, our main target application. The fast decay of the Gaussian basis function allows rapid convergence of the iterative solver even when the subdomains in the rasm are very small. At the same time we show that the accuracy of the interpolation can achieve machine precision. The present method was implemented in parallel using the petsc library (developer version). Numerical experiments demonstrate its capability in problems of rbf interpolation with more than 50 million data points, timing at 106 s (19 iterations for an error tolerance of 10^? 15) on 1024 processors of a Blue Gene/L (700 MHz PowerPC processors). The parallel code is freely available in the open-source model.     

Evaluation of multigrid acceleration for preconditioned time-accurate Navier-Stokes algorithms

The development of a two-dimensional time-accurate dual time step Navier-Stokes flow solver with time-derivative preconditioning and multigrid acceleration is described. The governing equations are integrated in time with both an explicit Runge-Kutta scheme and an implicit lower-upper symmetric-Gauss-Seidel scheme in a finite volume framework, yielding second-order accuracy in space and time. Issues concerning the implementation of multigrid for preconditioned, dual time step algorithms are discussed. Steady and unsteady computations were made of lid driven cavity flow, thermally driven cavity flow and pulsatile channel flow for a variety of conditions to validate the schemes and evaluate the effectiveness of multigrid for time-accurate simulations. Significant speedups were observed for steady and unsteady simulations. The speedups for unsteady simulations were problem dependent, a function of how rapidly the flow varied in time and the size of the allowable time step.     

Enhancement of an industrial finite-volume code for large-eddy-type simulation of incompressible high Reynolds number flow using near-wall modelling

We present a validation strategy for enhancement of an unstructured industrial finite-volume solver designed for steady RANS problems for large-eddy-type simulation with near-wall modelling of incompressible high Reynolds number flow. Different parts of the projection-based discretisation are investigated to ensure LES capability of the numerical method. Turbulence model parameters are calibrated by using a minimisation of least-squares functionals for first and second order statistics of the basic benchmark problems decaying homogeneous turbulence and turbulent channel flow. Then the method is applied to the flow over a backward facing step at Re_h = 37,500. Of special interest is the role of the spatial and temporal discretisation error for low order schemes. For wall-bounded flows, present results confirm existing best practice guidelines for mesh design. For free-shear layers, a sensor to quantify the resolution quality of the LES based on the resolved turbulent kinetic energy is presented and applied to the flow over a backward facing step at Re_h = 37,500.     

High-order accurate simulation of low-Mach laminar flow past two side-by-side cylinders using spectral difference method

This paper reports development of a two-dimensional solver for compressible viscous flow using spectral difference (SD) method and its applications on simulating laminar flow past two side-by-side cylinders at various spacings. The high-order spectral difference solver is based on unstructured quadrilateral grids. High-order curved wall boundary representation is developed for cylinders. Nine different spacings (center-to-center distance/diameter s = 1.1, 1.4, 1.5, 1.7, 2, 2.5, 3, 3.4 and 4) are investigated. The simulation results are compared to experimental results and other numerical results. As s increases, single bluff-body, flip-flopping, anti-symmetric and symmetric wake patterns are predicted.     

A block LU-SGS implicit unsteady incompressible flow solver on hybrid dynamic grids for 2D external bio-fluid simulations

A hybrid dynamic grid generation technique for two-dimensional (2D) morphing bodies and a block lower-upper symmetric Gauss-Seidel (BLU-SGS) implicit dual-time-stepping method for unsteady incompressible flows are presented for external bio-fluid simulations. To discretize the complicated computational domain around 2D morphing configurations such as fishes and insect/bird wings, the initial grids are generated by a hybrid grid strategy firstly. Body-fitted quadrilateral (quad) grids are generated first near solid bodies. An adaptive Cartesian mesh is then generated to cover the entire computational domain. Cartesian cells which overlap the quad grids are removed from the computational domain, and a gap is produced between the quad grids and the adaptive Cartesian grid. Finally triangular grids are used to fill this gap. During the unsteady movement of morphing bodies, the dynamic grids are generated by a coupling strategy of the interpolation method based on ‘Delaunay graph’ and local remeshing technique. With the motion of moving/morphing bodies, the grids are deformed according to the motion of morphing body boundaries firstly with the interpolation strategy based on ‘Delaunay graph’ proposed by Liu and Qin. Then the quality of deformed grids is checked. If the grids become too skewed, or even intersect each other, the grids are regenerated locally. After the local remeshing, the flow solution is interpolated from the old to the new grid. Based on the hybrid dynamic grid technique, an efficient implicit finite volume solver is set up also to solve the unsteady incompressible flows for external bio-fluid dynamics. The fully implicit equation is solved using a dual-time-stepping approach, coupling with the artificial compressibility method (ACM) for incompressible flows. In order to accelerate the convergence history in each sub-iteration, a block lower-upper symmetric Gauss-Seidel implicit method is introduced also into the solver. The hybrid dynamic grid generator is tested by a group of cases of morphing bodies, while the implicit unsteady solver is validated by typical unsteady incompressible flow case, and the results demonstrate the accuracy and efficiency of present solver. Finally, some applications for fish swimming and insect wing flapping are carried out to demonstrate the ability for 2D external bio-fluid simulations.     

An algebraic variational multiscale–multigrid method for large eddy simulation of turbulent flow

An algebraic variational multiscale–multigrid method is proposed for large eddy simulation of turbulent flow. Level-transfer operators from plain aggregation algebraic multigrid methods are employed for scale separation. In contrast to earlier approaches based on geometric multigrid methods, this purely algebraic strategy for scale separation obviates any coarse discretization besides the basic one. Operators based on plain aggregation algebraic multigrid provide a projective scale separation, enabling an efficient implementation of the proposed method. The application of the algebraic variational multiscale–multigrid method to turbulent flow in a channel produces results notably closer to reference (direct numerical simulation) results than other state-of-the-art methods both for mean streamwise and root-mean-square velocities. For predicting highly sensitive components of the Reynolds-stress tensor in the context of turbulent recirculating flow in a lid-driven cavity, the algebraic variational multiscale–multigrid method also shows a remarkably good performance in predicting reference results from experiment and direct numerical simulation compared to other methods.     

Boundary domain integral method for high Reynolds viscous fluid flows in complex planar geometries

The article presents new developments in boundary domain integral method (BDIM) for computation of viscous fluid flows, governed by the Navier–Stokes equations. The BDIM algorithm uses velocity–vorticity formulation and is based on Poisson velocity equation for flow kinematics. This results in accurate determination of boundary vorticity values, a crucial step in constructing an accurate numerical algorithm for computation of flows in complex geometries, i.e. geometries with sharp corners. The domain velocity computations are done by the segmentation technique using large segments. After solving the kinematics equation the vorticity transport equation is solved using macro-element approach. This enables the use of macro-element based diffusion–convection fundamental solution, a key factor in assuring accuracy of computations for high Reynolds value laminar flows. The versatility and accuracy of the proposed numerical algorithm is shown for several test problems, including the standard driven cavity together with the driven cavity flow in an L shaped cavity and flow in a Z shaped channel. The values of Reynolds number reach 10,000 for driven cavity and 7500 for L shaped driven cavity, whereas the Z shaped channel flow is computed up to Re = 400. The comparison of computational results shows that the developed algorithm is capable of accurate resolution of flow fields in complex geometries.     

Application of NASA general-purpose solver to large-scale computations in aeroacoustics

Of several iterative and direct equation solvers evaluated previously for computations in aeroacoustics, the most promising was the NASA-developed general-purpose solver (winner of NASA's 1999 software of the year award). This paper presents detailed, single-processor statistics of the performance of this solver, which has been tailored and optimized for large-scale aeroacoustic computations. The statistics compiled using an SGI ORIGIN 2000 computer with 12 Gb available memory (RAM) and eight available processors, are the central processing unit time, RAM requirements, and solution error. The equation solver is capable of solving 10 thousand complex unknowns in as little as 0.01 s using 0.02 Gb RAM, and 8.4 million complex unknowns in slightly less than 3 h using all 12 Gb. This latter solution is the largest aeroacoustics problem solved to date with this technique. The study was unable to detect any noticeable error in the solution, since noise levels predicted from these solution vectors are in excellent agreement with the noise levels computed from the exact solution. The equation solver provides a means for obtaining numerical solutions to aeroacoustics problems in three dimensions.