A Finite-Element Coarse-Grid Projection Method: A Dual Acceleration/Mesh Refinement Tool for Incompressible Flows

A Finite-Element Coarse-Grid Projection Method: A Dual Acceleration/Mesh Refinement Tool for Incompressible Flows

A. Kashefi A. Kashefi Engineering Science and Mechanics Program, Virginia Tech, Blacksburg, VA 24061, USA
A. E. Staples Faculty of Engineering, Engineering Science and Mechanics Program, Virginia Tech, Blacksburg, VA 24061, USA
   A. E. Staples A. Kashefi Engineering Science and Mechanics Program, Virginia Tech, Blacksburg, VA 24061, USA
A. E. Staples Faculty of Engineering, Engineering Science and Mechanics Program, Virginia Tech, Blacksburg, VA 24061, USA
2email: kashefi@vt.edu
3email: aestaples@vt.edu
Received: date / Accepted: date

Coarse grid projection (CGP) methodology is a novel multigrid method for systems involving decoupled nonlinear evolution equations and linear elliptic Poisson equations. The nonlinear equations are solved on a fine grid and the linear equations are solved on a corresponding coarsened grid. Mapping operators execute data transfer between the grids. The CGP framework is constructed upon spatial and temporal discretization schemes. This framework has been established for finite volume/difference discretizations as well as explicit time integration methods. In this article we present for the first time a version of CGP for finite element discretizations, which uses a semi-implicit time integration scheme. The mapping functions correspond to the finite-element shape functions. With the novel data structure introduced, the mapping computational cost becomes insignificant. We apply CGP to pressure-correction schemes used for the incompressible Navier-Stokes flow computations. This version is validated on standard test cases with realistic boundary conditions using unstructured triangular meshes. We also pioneer investigations of the effects of CGP on the accuracy of the pressure field. It is found that although CGP reduces the pressure field accuracy, it preserves the accuracy of the pressure gradient and thus the velocity field, while achieving speedup factors ranging from approximately 2 to 30. Exploring the influence of boundary conditions on CGP, the minimum speedup occurs for velocity Dirichlet boundary conditions, while the maximum speedup occurs for open boundary conditions. We discuss the CGP method as a guide for partial mesh refinement of incompressible flow computations and show its application for simulations of flow over a backward-facing step and flow past a cylinder.

Pressure-correction schemes Geometric multigrid methods Unstructured grids Coarse-grid projection Finite elements Semi-implicit time integration

1 Introduction and motivation

Since the late 1960s, projection methods [1–3] have been widely used for the numerical simulations of transient incompressible Navier-Stokes equations in extensive scientific areas [4, 5] and industrial fields [6, 7]. The high popularity of these methods is due to the scheme’s abilities to execute the incompressibility constrains using a set of two decoupled elliptic equations: a nonlinear advection-diffusion equation, and subsequently, a linear pressure Poisson equation. Although each of the cascade equations imposes computational costs to the system, Poisson’s equation tends to be the most time-consuming component of the flow simulation in complex geometries [8] as well as in computations with open boundary conditions [6, 9].

To substantially lessen the computational expenses, a common approach is to design robust multigrid (MG) solvers, which are solely devoted to one of the existing elliptic systems (see e.g., Refs. [10–12] for MG solvers for the advection-diffusion equation, and Refs. [13, 14] for those associated with Poisson’s equation). In contrast with the methods cited above, a novel family of multigrid procedures, the so-called Coarse Grid Projection (CGP) methodology by the authors [15–18], involves both the nonlinear and linear equations to efficiently accelerate the computations. In the CGP methodology, the nonlinear momentum is balanced on a fine grid, and the linear Poisson’s equation is performed on a corresponding coarsened grid. Mapping functions carry out data transfer between the velocity and pressure fields. In this sense, the CGP methods not only effectively relieve the stiff behavior of the primary Poisson problem, but can also take advantage of any desired fast elliptic solver to achieve large speedups while maintaining excellent to reasonable accuracy. The CGP methodology is provided in detail in Sect. 2.2.

The CGP framework was originally proposed by Lentine et al. in 2010 [15] for visual simulations of inviscid flows for video game applications. The uniform grid finite-volume and explicit forward Euler methods were chosen, respectively, for spatial and temporal discretizations in their numerical simulations. This gridding format led to dealing with a volume-weighted Poisson equation [19] in the presence of a solid object (e.g., a sphere) inside the flow field. Aside from the added complication, the grid format would have caused considerable reductions in visual fidelities, if the authors had not employed the complex mapping operators that they did. This makes the CGP technique less practical for curved boundaries. Furthermore, the explicit time integration restricted their algorithm to low CFL (Courant-Friedrichs-Lewy) numbers, and therefore long run times. In 2014 Jin and Chen [18] implemented the same CGP scheme [15] in the fast fluid dynamics (FFD) models to calibrate the effect of this computational accelerator tool on simulating building airflows. A decrease in spurious fluctuations of ventilation rates was observed, but the maximum achieved speedup was only a factor of approximately 1.50.

In 2013 San and Staples [16] expanded the CGP technique (labeled “CGPRK3”) to the vorticity-stream function formulation of the incompressible Navier-Stokes equations and demonstrated speedup factors ranging from approximately 2 to 42 in their numerical studies. Additionally, they extended the method in order to dramatically lower computational costs associated with elliptic equations of potential vorticity in quasigeostrophic ocean models [17]. Notwithstanding these successes, the CGPRK3 strategy has four main shortcomings. First, the nine-term full weighting restriction [20] and bilinear prolongation [20] operators used can be exclusively utilized in equally spaced grids. Consequently, one must analytically reformulate all steps according to generalized curvilinear coordinates except in uniform rectangular domains. Second, CGPRK3 applies the third-order Runge-Kutta [21] temporal integration to the non-incremental pressure correction method [1], while the splitting error of this specific projection scheme is irreducibly first-order in time and higher-order temporal integration methods do not improve the overall accuracy [1]. Third, the third-order time stepping scheme unnecessarily forces the CGPRK3 scheme to run the Poisson solver three times at each time step, whereas the primary goal of CGP is to reduce the computational effort that arises from the Poisson equation. Fourth, their suggested mapping procedure becomes more costly than that of the Poisson equation solver for high coarsening levels.

To obviate the aforementioned problems, a semi-Implicit-time-integration Finite-Element version of the CGP method (IFE-CGP) is presented in the current study. The incorporation of a semi-implicit backward time integration results in a simple five-step CGP algorithm with nearly zero cost for the data restriction /prolongation operators. It typically enables flow simulations at large time steps and thus improves speedup [22]. The triangular finite element meshes improve the CGP method in the following ways. First, they enhance the scheme’s capacity for the solution of fluid problems defined on complicated domains, where irregular grids and realistic boundary conditions are unavoidable. Second, they facilitate the design of the required mapping modules so that the restriction/prolongation operators can be optionally equivalent to the shape functions approximating multilevel nested spaces of velocity and pressure fields. Third, the generation of the Laplacian and divergence operators on a coarsened grid is expedited by means of available geometric/algebraic multigrid (GMG/AMG) tools for the finite element method (see e.g., [23–26]). This feature is particularly important for obtaining a sufficiently accurate solution of Poisson’s equation in modeling flow over obstacles.

This article’s objective is to present a simple, elegant version of the CGP method for finite element discretizations of incompressible fluid flow problems. As accelerating incompressible flow computations is the major application of CGP, speedup rates of the computations and the corresponding reduction in the accuracy of velocity and pressure fields are calculated. Besides this main application, we explore mesh refinement usages of the CGP method for the first time. As a next concern, because the accuracy of the pressure field in projection methods suffers from the negative influence of artificial homogeneous Dirichlet and Neumann boundary conditions on formations of boundary layers [6, 9], the possibility of thickening the layers by the CGP prolongation operator is investigated. Lastly, because the CGP procedure reduces the degree of freedom for the pressure component, a greater reduction in the integrity of the pressure field appears in comparison with the velocity field [16]. On the other hand, the pressure gradient (instead of simply pressure) is applied to a velocity correction step of pressure projection schemes [27]. With these two hypotheses in mind, we examine the effect of the CGP process on variations of both the pressure and its gradient. All the above mentioned numerical challenges are investigated through several representative benchmark problems: the Taylor-Green vortex in a non-trivial geometry, flow over a backward-facing step, and finally flow past a circular cylinder.

This article structured as follows. Sect. 2.1 provides the governing equations for incompressible viscous flows and their finite element formulations. Details on the proposed CGP algorithm are presented in Sect. 2.2. The computational implementation and a computational cost analysis are described, respectively, in Sect. 2.3 and 2.4. Numerical results and their interpretations are collected in Sect. 3. Conclusions and notes for extensions of the work are given in Sect. 4.

2 Problem formulation

2.1 Governing equations

Mass and momentum conservation of an incompressible isothermal flow of a Newtonian fluid are governed by the Navier-Stokes and continuity equations, with boundary conditions


where u and respectively denote the velocity vector and the absolute pressure in the fluid domain . The external force and stress vectors are represented by f and , respectively. is the fluid density and is the dynamic viscosity. The boundary of the domain consists of two non-overlapping subsets of Dirichlet and Neumann boundaries, where n indicates the outward unit vector normal to them. The system of equations is temporally integrated by the semi-implicit first-order backward differentiation formula [28] with time increment , and takes the form:


In the next stage, the non-incremental pressure-correction method [1] decouples the solution of the velocity and pressure variables. Based on it, at each time step , the output of the momentum equation is a non-divergent vector field called the intermediate velocity . Next, the divergence of the intermediate velocity is fed to the source term of the pressure Poisson equation. Finally, the intermediate velocity is corrected using the obtained pressure such that the end-of-step velocity vector satisfies the incompressibility constraint. The procedure yields two elliptic problems along with one correction equation, expressed as


Notice that the Poisson equation is restricted to unrealistic homogenous Neumann boundary conditions, when the velocity is subject to Dirichlet types. Contrarily, natural Neumann conditions for the momentum equation result in boundaries with spurious homogenous Dirichlet pressures.

The velocity and pressure spaces are approximated by the piecewise linear basis functions (P/P) of standard Galerkin spectral elements [29]. In this way, the resulting finite element form of Eqs. (9)–(15) is


where M and N indicate, respectively, the velocity mass and the nonlinear convection matrices. L, L, D, and G denote the matrices associated, respectively, to the velocity laplacian, the pressure laplacian, the divergence, and the gradient operators. The vectors Ũ, U, F, and P contain, respectively, the nodal values of the intermediate velocity, the end-of-step velocity, the forcing term and the pressure at time . The desired boundary conditions are implicitly enforced in the discrete operators. Further details of the elemental matrices are available in references [29, 30].

Remark 2.1 (On the discrete Brezzi-Babuska condition[31, 32]). Because the projection method overcomes the well-known saddle-point issue of Eqs. (1)–(2), the discrete Brezzi-Babuska condition [31, 32] can be ignored [27] and therefore the degree of the polynomials over the triangular mesh elements is chosen to be equal for the velocity and pressure. In addition, the identical resolutions allow the comparison of computational effectiveness between the current approach and that described in previous works [15–18].

2.2 Coarse grid projection methodology

The CGP methodology provides a multiresolution framework for accelerating pressure-correction technique computations by performing part of the computations on a coarsened grid. Alternatively, CGP can be thought as a guide to mesh refinement for pressure-correction schemes. Figure 1 gives a schematic illustration of the IFE-CGP algorithm for triangular finite element meshes. As shown in Fig. 1, at each time step , the intermediate velocity field data Ũ obtained on a fine grid is restricted to a coarsened grid. The divergence of Ũ, the intermediate velocity field data on the coarsened grid, plays the role of the source term in solving the Poisson equation on the coarse grid to detrimine the pressure P on the coarse grid; then, the resulting pressure data P is prolonged to the fine grid and becomes of P.

Figure 1: Scheme of coarse grid projection methodology involving the restriction and prolongation of the intermediate velocity and pressure data.

Either GMG or AMG techniques are key numerical tools to conveniently formulate the grid transfer, laplacian, and divergence operators in the IFE-CGP algorithm. In the present work, GMG methods are used, because using AMG methods in the derivation of the divergence operator is challenging for anisotropic triangular grids [33]. Here we generate hierarchical meshes in which each triangular element of a coarse mesh is conformingly subdivided into four new triangles. As a consequence, if the available finest mesh (with -element resolution) is generated after refinement levels, a relatively coarse mesh (with a resolution of elements) and its basis functions are available. Generally, one can obtain the corresponding coarsened grid using the space decomposition algorithm discussed for nonuniform triangular grids in the literature (see e.g., Refs. [34, 35]).

In principle, there can be an arbitrary number, , of nested spaces of progressively coarsened grids such that: . In practice, however, reasonable levels of accuracy can only be obtained for a maximum of three levels of coarsening. We restrict our attention to such cases. Consider a sequence of four nested spaces, , wherein if () corresponds to a fine grid, characterizes the space of the next coarsest grid. Recall that for any arbitrary element on , there is a piecewise linear shape function that is capable of estimating the space of the four sub-elements over . With this in mind, linear interpolation is implemented to construct the prolongation operator and its matrix representation P. Mathematically, the transpose of the prolongation operator is a feasible option for the definition of the restriction operator [26]. In spite of this fact, the restriction matrix is designed so that fine data of is directly injected throughout the coarse grid only if the data belongs to a common node between the two spaces; otherwise, it cannot be accessed by the coarse grid. Unlike , is able to directly connect two non-nested spaces. It is worth noting that the CGP method is compatible with other advanced data interpolation/extrapolation architectures (see e.g., [36]); however, even the simple mapping techniques introduced here are adequate.

In the GMG method used, a relatively coarse mesh, , with as a basis is directly accessible. Hence, the relevant laplacian () and divergence () operators are derived by taking the inner products of a coarse-grid finite element shape function on . Note that the spaces of velocity and pressure variables are not physically the same, although the same notation is used here for the sake of simplicity.

Eqs. (19)–(23) summarize the explanations in the preceding paragraphs and depict the IFE-CGP scheme for executing one time interval of the simulation.
1. Calculate Ũ on by solving


2. Map Ũ onto and obtain Ũ via


3. Calculate P on by solving


4. Remap P onto and obtain P via


5. Calculate U via


2.3 Computational implementation

A ++ object oriented code is developed according to the concepts addressed in Ref. [37]. To accelerate linear algebra executions and minimize memory requirements, the standard compressed sparse row (SR) format [38] is employed for sparse matrix-vector multiplication (SpMV) except in data transfer applications. The preconditioned GMRES() algorithm [39, 40] is chosen as an iterative linear solver. Absolute and relative tolerances are set to . The public unstructured finite element mesh generator Gmsh [41] is utilized. Calculations are performed on a single Intel(R) Xeon(R) processor with a 2.66 GHz clock rate and 64 Gigabytes of RAM.

2.4 Computational cost analysis of IFE-CGP

A finite element analysis traditionally consists of three major portions: preprocessing, processing, and postprocessing. In the simulations undertaken here, the postprocessing stage mostly involves writing output files without a significant effect on the CPU time. As a result, the numerical cost, , of the pressure correction approach on a given coarse grid is estimated by


where is the preprocessing cost, whereas and comprise, respectively, the cost of the Poisson equation (see, Eq. (17)) and the remaining algorithms in the processing portion. For transient problems with a significant number of time steps, the processing cost always overcomes the preprocessing price. By -level uniform refinements of a two dimensional domain, a high-resolution simulation takes time , roughly expressed as


with two factors and . The factor depends on the computational resources and global matrix assembling algorithms. Because by quadrupling the element numbers, the global node numbers are doubled at minimum, is lower bound for the increment in the node numbers. Consequently, represents a lower bound associated with the matrix size enhancement. Additionally, note that is the lowest possible factor for cost scaling of Eqs. (16)–(17), because the cost of a matrix inversion is not linearly proportion to its size. Taking the advantages of the IFE-CGP method into account, the Poisson solution is performed on the coarse grid and its cost is not scaled. Therefore, is reduced to the IFE-CGP cost , given by


where is a new factor for the preprocessing and indicates the mapping expenses. The coefficient is not necessarily equal to the factor and the relation is variable. For instance, because the assembling process of rather than D is more cost-effective, it might be concluded that . But if one takes the transpose of G to establish D, that conclusion is questionable. Besides, the preparation of and involves an extra cost for IFE-CGP. A numerical comparison in Sect. 3.2 clarifies this point furthur. is shown to be negligible in comparison with the other three terms of Eq. (26) in Remark 2.2.

From a mesh refinement application point of view, the cost increment factor of the computational IFE-CGP tool () is approximated by


similarly, this factor for a regular triangulation refinement () is conjectured to be


Based on the above discussion, is greater than . This is mainly due to the factor of that multiples in Eq. (25). These results imply that mesh refinement using the CGP idea is more cost effective than the standard technique.

Remark 2.2 (On implementing the mapping operators). The implementation of in the CSR format is not possible because the matrix contains null vectors. Additionally, if is constructed in the CSR format, the data prolongation cost of each IFE-CGP loop is approximated by where is proportional to the number of non zero elements per row of (Depending on the mesh nonuniformity, it varies between 2 and 10 in the current grids.) and denotes the number of pressure unknowns on . Though CSR is inexpensive, an easier-to-implement method is introduced next. Consider two data structures and including three () and two () integer indices, respectively. An array of each data structure and , with equal to the node numbers on , is created as follows:
For find a pair of indices that satisfies


and for find an index that satisfies


where P is the pressure value at the th node of the spanned space , while Ũ is the restricted intermediate velocity of the th node of the discretized space . P, P, and Ũ are similarly defined. With respect to the prolongation function, this trick improves the performance by reducing the computational effort to , with only moderately increased memory usage. Likewise, the numerical expense of the suggested injection operator is of order .

3 Results and discussion

To evaluate the various aspects of the IFE-CGP method, three standard test cases are studied: the Taylor-Green decaying vortex problem, the flow over a backward-facing step, and the flow around a circular cylinder. Here, the grid resolution of a numerical simulation is denoted by , where indicates the element numbers of a fine grid used for the advection-diffusion solver, and demonstrates the element numbers of a corresponding coarsened grid for the Poisson solver. When is equal to , the standard, non-IFE-CGP, algorithm is recovered. indicates the level of mesh coarsening used in the IFE-CGP method, and is equal to zero for the standard algorithm. If function is assumed to be a finite element approximation of an exact solution, , on the domain , with elements, the difference between these two functions, , is defined as:


The discrete norms are defined as:


where is the discrete space of th element, , of . Note that when an exact solution is not available, the norms are measured with reference to the standard algorithm ().

3.1 Taylor-Green vortex in a non trivial geometry

The Taylor-Green vortex problem [42] is a widely-used benchmark problem which is an exact analytic solution of the unsteady incompressible Navier-Stokes equations in two-dimensions (see e.g., [6]). The velocity field is given by:


And the pressure field is given by:


A density value of kg/m and a viscosity of Pas are used. One goal of the Taylor-Green test case is an examination of the IFE-CGP method capability in complex geometries. For this purpose, a square domain with a circular hole is chosen such that

The geometry is depicted in Fig. 2 and details of the meshes are given. A similar computational domain (a rectangular hole instead of a circular one) has been implemented by J. M. Alam et al. [43] to perform this test case. A second goal of this section is an investigation of the effects of velocity Dirichlet boundary condition (and consequently artificial Neumann boundary conditions) on the rate of accelerating computations by the IFE-CGP algorithm. Therefore, the velocity domain boundaries are set to the exact solution of Eqs. (34)–(35). These types of boundary conditions have been previously applied to the Taylor-Green vortex problem in the literature [6, 27]. San and Staples [16] have also studied this problem to validate CGPRK3 performance, but using periodic boundary conditions. In this way, an opportunity for comparison is provided. As a last concern, the effects of the prolongation operator on the thickening of artificial boundary layers are investigated. The simulations are run with a constant CFL number. As a result, a time step of 0.01 s is selected for the 3384:3384 case, and this is halved for each quadrupling of the advection-diffusion solver grid.

Figure 2: Representation of the triangular finite element meshes used for solving Poisson’s equation in the simulation of Taylor-Green vortex. a After one level coarsening , 27520 nodes and 54144 elements; b After two levels coarsening , 6992 nodes and 13536 elements; c After three levels coarsening , 1804 nodes and 3384 elements.

The discrete norms of the velocity field for different mesh resolutions are tabulated in Table 1. Considering all the cases, for one and two levels (, and ) of the Poisson grid coarsening, the minimum and maximum of the error percentage relative to the finest mesh () are, respectively, 0.30% and 3.61%. However, a considerable reduction in the velocity accuracy is found for three levels of coarsening. For instance, the infinity norm computed for the velocity field obtained on the 216576:3384 grid resolution indicates a 99% reduction in the accuracy level, but it is still two orders of magnitude more accurate in comparison with the resulting data captured from the full coarse scale simulation performed on the 3384:3384 grid resolution. The speedup factors achieved range from 1.601 to 2.337. The velocity and the pressure field magnitudes for the result with three levels of coarsening are depicted in Fig. 3. The flow fields have reasonable levels of accuracy, however, dampened flows can be observed near the boundaries.

Resolution CPU time/t Speedup
0 216576:216576 1.59938E9 8.34306E10 91644000 1.000
1 216576:54144 1.59453E9 8.38567E10 57223280 1.601
2 216576:13536 1.63584E9 8.64473E10 40771280 2.247
3 216576:3384 3.18962E9 1.15333E9 40339600 2.272
0 54144:54144 1.29418E8 6.78120E9 1546544 1.000
1 54144:13536 1.29274E8 6.85697E9 761040 2.032
2 54144:3384 1.38146E8 7.48644E9 661788 2.337
0 13536:13536 1.05619E7 5.58118E8 30944 1.000
1 13536:3384 1.08807E7 5.73166E8 13598 2.275
0 3384:3384 8.73658E7 4.69697E7 661 1.000
Table 1: Error norms, CPU times per time step, and relative speedups for different grid resolutions of the Taylor-Green vortex simulation at s. Resolution in the form of represents the grid resolution of the advection-diffusion solver, elements, and Poisson’s equation, elements.
Figure 3: The solution of Taylor-Green vortex in a non trivial geometry at s obtained using the CGP method. The resolution of the nonlinear and linear equations is, respectively, 54144 and 3384 elements. a Velocity magnitude contour and streamlines; b Pressure field.

In the non-incremental pressure correction methods, for a non-smooth domain with artificial Neumann boundary conditions, the pressure error norms are not sensitive to the spatial resolution and they are mostly not improved by increasing the node numbers [27]. Contrarily, the pressure gradient error norms decrease by an increase in the number of degree of freedom [27]. The data collected in Table 2 demonstrates these findings. Regarding the standard computations (), the discrete norms of the pressure are only reduced by one order of magnitude with three levels of mesh refinement of both the pressure and the velocity spaces, whereas these quantities for the pressure gradient have a reduction of three orders of magnitude with the same grid resolution increase. By looking at the resulting pressure error norms using IFE-CGP, we find almost the same data obtained using the standard algorithm. But this does not mean that the IFE-CGP procedure conserves the accuracy of the pressure field, because the pressure field from the full fine scale simulations itself has a high level of error. In fact, it can be only concluded that the prolongation operator does not make the results worse. In contrast with the pressure norms, we observe similar trends between the velocity and pressure gradient norms for the IFE-CGP results. For instance, for one level () of coarsening, the minimum and maximum of the error percentage with reference to the finest mesh () are, respectively, 5.37% and 16.83%. As a note, the discrete norms calculated for the pressure gradient are relatively small in comparison with the velocity and pressure norms, and it is because we use the discrete gradient operator G in order to compute the gradient of both the exact solution and the numerical result. Recall that at the velocity correction step (see Eq. (15)), the pressure gradient is used to correct the velocity field, not the pressure itself. Hence, since IFE-CGP conserves the pressure gradient accuracy, the fidelity of the velocity field is also conserved.

Resolution P P
0 216576:216576 1.34241E6 6.00500E7 3.21871E12 1.81154E13
1 216576:54144 1.34241E6 6.00416E7 3.09665E12 1.90897E13
2 216576:13536 1.34241E6 6.00270E7 2.22832E12 2.43972E13
3 216576:3384 1.34241E6 6.00171E7 3.65721E12 5.78906E13
0 54144:54144 4.98048E6 2.38748E6 4.09860E11 3.48936E12
1 54144:13536 4.98048E6 2.38692E6 3.96108E11 3.81608E12
2 54144:3384 4.98048E6 2.38657E6 3.22830E11 6.43534E12
0 13536:13536 1.89334E5 9.42727E6 4.86969E10 7.21105E11
1 13536:3384 1.89334E5 9.42590E6 5.31351E10 8.42464E11
0 3384:3384 6.74979E5 3.66242E5 6.55823E9 1.66480E9
Table 2: The infinity and second error norms of pressure quantity and its gradient for the Taylor-Green vortex simulation at s. Resolution in the form of represents the grid resolution of the advection-diffusion solver, elements, and Poisson’s equation, elements.

The Cartesian coordinate of an element’s center at which the infinity error norm occurs for the velocity, pressure, and pressure gradient in the Taylor-Green vortex domain are tabulated in Table 3. For different levels of coarsening, the location changes (except for the pressure error, which remains unchanged). However, it is always near the boundaries of either the ring or the square. In particular, the maximum errors for the velocity and pressure gradient fields occur on the edge of the square for (). Thus, even using the prolongation operator of the IFE-CGP method, the maximum errors still occur around the boundaries with velocity Dirichlet conditions.

0 216576:216576 (0.576623, 0.808348) (0.482527, 0.110053) (0.340604, 0.695205)
1 216576:54144 (0.580109, 0.809028) (0.482527, 0.110053) (0.659396, 0.304795)
2 216576:13536 (0.601023, 0.813106) (0.482527, 0.110053) (0.306953, 0.664904)
3 216576:3384 (0.770613, 0.977284) (0.482527, 0.110053) (0.004544, 0.243629)
Table 3: Comparison of various locations of velocity magnitude, pressure, and gradient pressure infinity-norm errors for the Taylor-Green vortex problem at s. indicates the center of an element, which the errors occur in. Resolution in the form of demonstrates the grid resolution of the advection-diffusion solver, elements, and Poisson’s equation, elements.

From a general point of view, the spatial discretization of the advection-diffusion domain acts as a low-pass filter on the grid, and the Poisson solver also acts as a pre-filtering process [16]. The CGP procedure specifically uses the belief in order to increase saving in computational time without negatively affecting the properly-resolved advection-diffusion field. A visual demonstration of these effects is displayed in Figs. 4–5. If the difference between the exact solution and the numerical result is considered as a noise over the domain, Figure 4 and Figure 5 show the noise distribution, respectively, for the pressure and velocity variables. For example, by switching the grid resolution from 216576:216576 (see Fig. 4a and Fig. 5a) to 216576:54144 (see Fig. 4b and Fig. 5b), the maximum absolute value of the pressure noise roughly changes from to , implying a 7.14% noise increment, whereas there is no significant noise enhancement for the velocity field. Similar behavior is observed for further levels of the Poisson grid coarsening. For instance, as it can be seen in Fig. 4d and Fig. 5d, the noise of the pressure and velocity domains with three levels of coarsening increases, respectively, to 85.71% and 33.33%, demonstrating that the low-pass filtering feature of the discretization of the advection-diffusion part of the algorithm is why the large errors of the pressure field do not get transmitted to the velocity field in CGP methods.

Figure 4: Distribution of point wise pressure error for Taylor-Green vortex problem at s. Labels in the form of indicate the grid resolution of the advection-diffusion solver, elements, and the Poisson equation, elements. a Regular computation ; b IFE-CGP ; c IFE-CGP ; d IFE-CGP .
Figure 5: Distribution of point wise velocity magnitude error for Taylor-Green vortex problem at s. Labels in the form of indicate the grid resolution of the advection-diffusion solver, elements, and the Poisson equation, elements. a Regular computation ; b IFE-CGP ; c IFE-CGP ; d IFE-CGP .

Remark 3.1 (On the error order of the non incremental pressure-correction scheme). By calculating the slope of the discrete norms for the standard algorithm results, it is induced that the order of the overall accuracy is not fully a function of the spatial resolution. This is because the time error is dominant over the spatial step in the non-incremental pressure correction scheme with Dirichlet boundary conditions [1]. Numerically, in order to remove the time error, very small time intervals should be taken. For instance, this time error becomes invisible if the time step is chosen to be s for the simulation with a 3384:3384 spatial resolution. This severe condition is not odd because the fluid domain is surrounded by velocity Dirichlet boundaries from both the outside and inside. However, our computational resources do not allow us to choose such a tiny time increment for this test problem. Note that this discussion is not directly relevant to the topic of CGP methodology, but it can be concluded that a multigrid method (e.g., the CGP technique) provides optimal results in the presence of velocity Dirichlet boundary conditions if users take the incremental pressure-correction schemes introduced in Refs [1, 6, 27, 44, and 45].

3.2 Flow over a backward-facing step

To study the IFE-CGP algorithm efficiency in the presence of open boundary conditions, the flow over a backward facing step inside a channel is analyzed. Fig. 6 presents the problem geometry and imposed boundary conditions. Because an inlet channel upstream significantly affects the flow simulation at low Reynolds number [46], the inflow boundary is located at the step and is described by a parabolic profile:


The origin of the coordinate system is placed in the lower left corner of the step. Homogeneous natural Neumann conditions


are enforced at the exit. The Reynolds number is calculated as


where is the channel height and represents the space-averaged mean entrance flow velocity.

Figure 6: Schematic view of flow past a backward-facing step.

Although the stress free boundary condition is less restrictive than Dirichlet type boundary conditions at the outflow [27], it leads to the increased the possibility of a loss in spatial accuracy [6]. This unfavorable situation is mainly due to imposing artificial homogenous Dirichlet boundary conditions in the pressure Poisson solver [6]. On the other hand, because the IFE-CGP technique reduces the number of pressure unknowns at the outflow, it is one of the current research interests to check whether the IFE-CGP methodology provides a valid solution. For this purpose, the prediction of the reattachment length with respect to is plotted in Fig. 7. The obtained results for one () and two () levels of coarsening reveal good agreement with the numerical data of Kim and Moin [47], and Erturk [48]. At , the reattachment length just differs from that reported by Erturk [48] about 0.7% and 4.0%, respectively, for the IFE-CGP () and IFE-CGP () algorithms. However, the IFE-CGP () approach vastly overestimates the reattachment length after the Reynolds number has reached the value 300.

Figure 7: Normalized reattachment length as a function of Reynolds number.

To save space, detailed results are only presented for a Reynolds number of . The time step is chosen to be s. Based on our numerical experiments, the flow field reaches stationarity at s for the regular fine computations. The results for all other options are reported at this fluid flow time.

Figure 8 depicts the axial velocity contour maps of the flow simulated using both the normal and the IFE-CGP processes. Additionally, the efficiency and accuracy of the velocity field for the standard algorithm and the IFE-CGP approach are compared in Table 4. Flow velocity variables with one and two levels of coarsening agree well with the fine scale standard computations. Significantly, a 30-times reduction in computational cost of the IFE-CGP solution with the 102400:1600 resolution is obtained, as the corresponding velocity error norms are still in the acceptable range. As can be seen from Fig. 9a, interestingly, the advection-diffusion solver diverges for a simulation with the 1600:1600 resolution. This behavior will be discussed comprehensively in Sect. 3.4.

Resolution CPU (s) Speedup
0 102400:102400 - - 3447941 1.000
1 102400:25600 6.69447E6 7.88459E7 173052 19.924
2 102400:6400 1.43997E5 4.87393E6 120500 28.613
3 102400:1600 4.4854E5 2.23362E5 116745 29.534
Table 4: Comparison of relative norm errors and total CPU times between the standard and IFE-CGP algorithms for the flow past a backward-facing step at .
Figure 8: Horizontal velocity component contour plot of flow over backward-facing step at . Labels in the form of indicate the grid resolution of the advection-diffusion solver, elements, and Poisson’s equation, elements.
Figure 9: Demonstration of partial mesh refinement application of the CGP method, a comparison between axial velocity contours of a Regular computation, diverged, and b The IFE-CGP , converged, for flow past backward-facing step at . Labels in the form of indicate the grid resolution of the advection-diffusion solver, elements, and the Poisson equation, elements.

The CPU times devoted to the restriction/prolongation operators are tabulated in Table 5. By increasing the coarsening level , the prolongation operator becomes more expensive, whereas the time consumed by means of the restriction function decreases. To explain this fact, let’s consider, for instance, the data mapping procedure of the IFE-CGP strategy on the 102400:1600 grid resolution. The restriction operator directly injects the intermediate velocity field data from a grid with 102400 elements into the corresponding coarse grid with 1600 elements. From a programming point of view, this operation needs only 905 loops, which is the pressure node numbers of the coarse grid. The prolongation operator, in contrast, has to utilize two intermediate grids, associated with and , in order to extrapolate the pressure domain data. To handle this transformation, the prolongation operation requires 68659 loops, which is the summation of grid points belonged to the two intermediate as well as the finest meshes. As can be seen from Table 5, the IFE-CGP method slightly increases the time spent on the preprocessing block. Even though the and matrix assembling process is computationally cheaper in comparison with the standard algorithm, the mapping operator constructions ultimately overcome these savings. That is, the ratio is greater than 1.00 in Table 4. Appling AMG tools in order to establish the and matrices in the IFE-CGP technique is a way to optimize the preprocessing subroutine costs.

Resolution Restriction (s) Prolongation (s) Preprocessing (s) ()
0 102400:102400 - - 407.471 1.000
1 102400:25600 0.89 3.11 698.441 1.714
2 102400:6400 0.21 3.67 713.521 1.751
3 102400:1600 0.05 3.64 725.521 1.780
Table 5: CPU times consumed by restriction and prolongation operators, preprocessing segment, and its ratio () in comparison between the IFE-CGP and standard schemes for the backward-facing step flow simulation (see Sect. 2.4 for further details).

3.3 Flow past a circular cylinder

The concern of this section is a study of the effect of curvature on the IFE-CGP method. Although this capability of the method has been already investigated for decaying vortices in curved geometries in Sect. 3.1, an external unsteady flow over a cylinder is a more physically meaningful benchmark case [49]. San and Staples [16] have performed this fluid mechanics problem by means of the CGPRK3 solver, but for steady-state flows, at , and exclusively using one level of coarsening (). Here, the IFE-CGP methodology with three levels of coarsening is applied to this canonical flow problem at .

The computational field is considered as a rectangular domain . A circle with diameter represents the cylinder in two dimensions, and its center lies at the point (8, 16). A free stream velocity parallel to the horizontal axis is imposed at the inflow boundary. The circle is treated as a rigid boundary and no-slip conditions are enforced. The velocity at the top and bottom of is perfectly slipped in the horizontal direction. The outflow velocity is specified with a natural Neumann condition, Eq. (39). The Reynolds number is determined as


To set this dimensionless number to , the density, cylinder diameter, and free stream velocity are set to 1.00; and the viscosity is set to 0.01 in the International Unit System. The described geometry and boundary conditions are taken from the literature [50, 51] to satisfy far-field assumptions. A fixed time increment of s is selected and the numerical experiments are executed until time s. The grids utilized by the Poisson solver for , , , and are like those shown in Fig. 10, with 108352 nodes and 215680 elements, 27216 nodes and 53920 elements, 6868 nodes and 13480 elements and 1749 nodes and 3370 elements, respectively, with very fine grid spacing near the circle.

Figure 10: Computational mesh for the Poisson equation solution of the flow past a circular cylinder. a After two levels coarsening ; b After three levels coarsening . Details of the grids are reported.

A visual comparison between the obtained vorticity fields with and without the IFE-CGP method is made in Fig. 11 for several spatial resolutions at time s. For all levels of coarsening, the IFE-CGP field provides more detailed data compared to that modeled with a full coarse grid resolution. A comparison between the resulting vorticity fields with the resolutions of 215680:215680 and 13480:13480 demonstrate that the phases of periodic variation of these two fields are not equal to each other. Conversely, the fields computed by 215680:215680 and 215680:13480 oscillate with the same phase. However, there is a phase lag between the outcomes with 215680:215680 and 215680:53920 or 215680:3370 mesh resolutions. Our numerical experiments show that the phase lag between the standard and IFE-CGP results depends on the time step chosen. For instance, the velocity field phases, and consequently the vorticity ones, are the same in both the simulations performed by IFE-CGP () and IFE-CGP () tools when s.

Figure 11: Vorticity fields for the flow past a circular cylinder at s. Labels in the form of illustrate the spatial resolution of the advection-diffusion grid, elements, and the Poisson equation mesh, elements.

To more precisely analyze the IFE-CGP algorithm’s performance, velocity and vorticity distributions along the horizontal centerline, behind the cylinder and in the wake region, are shown in Fig. 12 at time s. By coarsening the advection-diffusion mesh at the exit of the fluid domain, the results of the 53920:53920 resolution includes spurious fluctuations. These fluctuations become stronger in the pure coarse case with the 13480:13480 resolution. However, they are successfully removed using the IFE-CGP approach. The computational times for the performed simulations are: 491068.2 s, 70809.2 s, 65492.7 s, and 56251.1 s, respectively for (215680:215680), (215680:53920), (215680:13480), and (215680:3370), leading to speedup factors between 6.935 and 8.730.

Figure 12: Comparison of various output variables in the wake of the flow over a cylinder at s for different values of the advection-diffusion and the Poisson grid resolutions. a Horizontal velocity for IFE-CGP (, 1, and 0); b Horizontal velocity for IFE-CGP (, 2, and 0); c Vertical velocity for IFE-CGP (, 1, and 0); d Vertical velocity for IFE-CGP (, 2, and 0); e Vorticity for IFE-CGP (, 1, and 0); f Vorticity for IFE-CGP (, 2, and 0). The grid resolution in the form of shows the element numbers of the advection-diffusion grid by , and the Poisson grid by .

Table 6 lists the calculated Strouhal number (), drag () and lift () coefficients compared with the experimental and numerical studies presented in [52–55]. The Strouhal number is based on the time evolution of the blunt body lift, and is formulated as:


where is the shedding frequency. The data presented in Table 6 demonstrates that as the coarsening level () increases, the drag and lift forces slightly decrease and increase, respectively. However, they still agree well with the values found in the literature. Because the IFE-CGP method computes the pressure variable on a coarsened mesh, there is a concern it might reduce the accuracy of the drag and lift coefficients. To check this issue, the time evolution of the viscous (), pressure (), and total lift coefficients are plotted separately in Fig. 13. Interestingly, the viscous and pressure lift diagrams for the full fine scale and the IFE-CGP () simulations become nearly identical after approximately time s; nevertheless, the IFE-CGP configuration provides the maximum magnitude of the lift force two periods of the flow cycle earlier. Even though a numerical computation performed on a pure coarse grid degenerates the viscous lift coefficient, choosing one level of coarsening for the IFE-CGP mechanism does not influence the viscous force integrity. For two further levels of the Poisson grid coarsening, spurious fluctuations can be observed at the beginning of the fluid flow simulation. For the IFE-CGP () results, when the simulation is using the 215680:3370 grid resolution, the lift coefficient reaches its final amplitude at approximately time s, roughly 60 s earlier than the standard full fine scale computations. Surprisingly, the lift force obtained by full coarse scale (3370:3370) computations oscillates around approximately 0.2 (instead of 0.0), showing the flow field is under-resolved. Contrarily, the IFE-CGP () lift force is more reliable since it fluctuates around 0.0.

Study C C
Braza et al. [52] 0.160 1.3640.015 0.25
Liu et al. [53] 0.165 1.3500.012 0.339
Hammache and Gharib [54] 0.158 - -
Rajani et al. [55] 0.156 - -
IFE-CGP () 0.156 1.2580.006 0.217
IFE-CGP () 0.152 1.2230.006 0.250
IFE-CGP () 0.149 1.1680.034 0.344
Table 6: Strouhal number, drag and lift coefficients for .
Figure 13: A Comparison between the time evolution of lift coefficients of the flow past a cylinder at for different combinations of the advection-diffusion and the Poisson mesh resolutions. a Viscous lift coefficient for IFE-CGP (, 1, and 0); b Viscous lift coefficient for IFE-CGP (, 3, and 0); c Pressure lift coefficient for IFE-CGP (, 1, and 0); d Pressure lift coefficient for IFE-CGP (, 3, and 0); e Lift coefficient for IFE-CGP (, 1, and 0); f Lift coefficient for IFE-CGP (, 3, and 0). The spatial resolution in the format of indicates the finite element numbers of the advection-diffusion grid by , and the Poisson grid by .

The magnitude of the centerline pressure and its gradient along the -axis in the wake region are shown in Fig. 14 for the full fine (215680:215680), IFE-CGP (, and 2) (215680:53920 and 215680:13480), and full coarse (53920:53920 and 13480:13480) simulations at time s. This figure reveals a key feature of the CGP methodology. The pressure magnitude estimated by the IFE-CGP (, and 2) method is far from that computed with the full-fine scale resolution. Comparably, the pressure magnitudes obtained using the standard algorithm with the full fine and coarse resolutions are relatively close to each other. Notwithstanding this, the pressure gradient magnitude of the full fine-scale and the IFE-CGP (, and 2) computations are indistinguishable for most of the domain (although with a phase lag in the case of one level coarsening). As discussed earlier, in contrast with the implicit pressure quantity, the pressure gradient is the relevant variable in the incompressible Navier-Stokes equations [27]. Hence, it does not matter if the IFE-CGP strategy does not retain absolute pressure values close to the full fine results. Importantly, it calculates the pressure gradients with a high level of accuracy and significantly better than those that are solely computed on a coarse grid.

Figure 14: The IFE-CGP methodology’s effect on the integrity of a The pressure variable and b Its gradient; The 215680:215680 label indicates the standard algorithm with full fine scale computations, The 215680:53920 and 215680:13480 labels show, respectively, the IFE-CGP outputs with one and two levels of coarsening, and 53920:53920 and 13480:13480 labels illustrate the standard algorithm with full coarse scale computations.

So far we have emphasized the fact that spurious oscillations of the velocity and vorticity fields are removed by IFE-CGP. Now by looking at Fig. 14, we see that this is also the case for the pressure fields. Although both the IFE-CGP (, and 2) and full coarse (53920:53920 and 13480:13480) simulations solve the pressure Poisson equation on the same coarse mesh, the artificial oscillations at the end of domain are removed only in the case of the IFE-CGP outputs. Thus, for the same number of degrees of freedom, when the pressure Poisson equation is fed with a smoother intermediate velocity field, the outcome pressure filed is also smoother. Lastly, there is a noticeable difference between the pressure gradient magnitude of IFE-CGP and full fine scale computations near m. This difference comes from homogenous artificial Dirichlet boundary condition for the pressure (). As can be seen from Fig. 14a and Fig. 14c, the pressure of IFE-CGP falls sharply near m to satisfy this boundary condition. It is conjectured that this issue would be eliminated by switching from non-incremental pressure projection methods to incremental ones.

3.4 CGP as partial mesh refinement

Let’s consider a condition that the standard numerical simulation diverged for the case due to a relative high Reynolds number or too coarse a mesh. Or the results obtained with a grid resolution are not sufficiently resolved and a fluid field with more detailed information is needed. The standard approach to resolving these common issues in projection methods is to refine both the advection-diffusion and the Poisson grids. In contrast with this approach, the CGP strategy suggests refining the advection-diffusion grid, without changing the resolution of the Poisson mesh. To be more precise from a terminology point of view, CGP does not propose a new mesh refinement method; however, it guides users to implement available mesh refinement techniques for the grids associated with the nonlinear equations. Here, we describe the concept by showing two simple examples.

Consider the simulation of the flow over a backward-facing step. Let’s assume that one is interested in the flow information at ; however, due to wall clock time or computational resource limitations, he is not able to run a simulation with the required pure fine 102400:102400 grid resolution. On the other hand, because a coarse 1600:1600 resolution is not high enough, the simulation diverges after 96 time steps, as depicted in Fig. 9a. The IFE-CGP framework with an intermediate resolution of 102400:1600 provides a converged solution as shown in Fig. 9b. The relative velocity error norms with reference to the full fine simulation are of order . Furthermore, the normalized reattachment length can be estimated around 14. Although this estimation is slightly different from that obtained by the standard computations, it is captured 30 times faster. Note that these results are achieved by only refining the advection-diffusion equation solver mesh. In this case, because the simulation on the coarse grid diverges, there is no real number for ; however, if a virtual considered, .

Let’s reconsider the flow over a circular cylinder computations described in Sect. 3.3. Coming back to the lift coefficient graphs, depicted in Fig. 13, and with an emphasis on the IFE-CGP () outputs, another interpretation of these results is discussed here. Let’s assume an exact measurement of the lift coefficient is needed for a specific engineering purpose. Using standard methods, this can be accomplished using either 215680:215680 or 53920:53920 grid resolutions. An implementation with the finer grid produces a more precise answer. It could be a user’s incentive to locally/globally refine the full coarse mesh. Obviously, this mesh refinement ends in an increase in CPU time for the processing part of the simulation. In this case, our numerical experiments show that the increase is equal to 339271.6 s (over 94 hr). As mentioned earlier, having a coarse mesh only degrades the level of the viscous lift accuracy. In fact, instead of refining the grids associated with both the nonlinear and linear equations, a mesh refinement of the nonlinear part is enough alone. Hence, in order to increase the precision of the lift force, one can refine the advection-diffusion grid and keep the resolution of the Poisson mesh unchanged. In this case, the IFE-CGP grid refinement cost factor is , whereas this factor for the regular mesh refinement is , illustrating a considerable saving of computational resources.

As a last point, obviously the types of two dimensional flow simulations described here are not challenging computation problems today. These problems are merely used as examples to explain one of the features of the IFE-CGP algorithm. Practical applications of the IFE-CGP mechanism as a mesh refinement tool are expected to be useful for three dimensional flow simulations on parallel machines.

3.5 Boundary conditions and data structure effects on CGP efficiency

Figures 15a–15d compare the CPU times consumed by various components of the processing segment for four test cases using different boundary conditions. The Taylor-Green vortex, flow over a backward-facing step, and flow past a cylinder are modeled by the IFE-CGP approach, while the double shear layer problem is simulated by the CGPRK3 technique [16]. Regarding the IFE-CGP strategy, by the coarsening level increment, the Poisson equation price lessens dramatically so that its portion becomes less than 5% after just one level of coarsening. For that reason, a considerable speedup cannot be achieved after . A similar trend occurs in application of the CGPRK3 method with the difference that the major reduction in the Poisson solver cost is obtained at . The significant difference between the IFE-CGP and CGPRK3 algorithms is associated with the computational expense of the data transfer between the advection-diffusion and the Poisson grids. According to Fig. 15d, the CGPRK3 method allocates more CPU resources to mapping the data for each subsequent level of coarsening, and finally the mapping process costs overcome the Poisson equation costs at . In contrast, this charge never exceeds 0.004% of the total algorithm computational cost using the IFE-CGP strategy.

Figure 15: The influence of boundary condition types and mapping algorithms on the achieved speedup by the CGP technique for a the Taylor-Green vortex problem with velocity Dirichlet boundary conditions, b the flow over a backward-facing step with stress free conditions, c the flow past a cylinder with open boundary conditions, and d the double shear layer problem [16] with periodic boundary conditions.

Concerning influences of the boundary condition, the maximum speedup is achieved when the outflow velocity in the backward-facing step problem is treated as a stress-free condition, and the lowest acceleration of the computations is observed when velocity Dirichlet boundary conditions are enforced for the Taylor-Green vortex simulation. The speedups when using periodic boundary conditions are in between these two extremes. This difference is expected because solving a Poisson equation with spurious homogeneous Dirichlet boundary conditions requires more computational effort in comparison with pressure Neumann conditions [56]. In terms of the accuracy level, the CGP class of methods retains the velocity field data close to that of full fine grid resolution simulations in the presence of less restrictive velocity boundary conditions, such as open and periodic ones. However, the CGP approach with pure velocity Dirichlet conditions results in a dampened velocity field, which has been also reported by Lentine et al. [15]. In the case of Euler’s equations, where the viscous term is neglected and the flow is allowed to slip on the solid surfaces, the CGP method acquires much less damped flows. For the same reason, the damping phenomenon disappears when the CGP solver is run at high Reynolds numbers.

4 Conclusions and future directions

The CGP method is a new multigrid technique applicable to pressure projection methods for solving the incompressible Navier-Stokes equations. In the CGP approach, the nonlinear momentum equation is evolved on a fine grid, and the linear pressure Poisson equation is solved on a corresponding coarsened grid. Mapping operators transfer the data between the grids. Hence, one can save a considerable amount of CPU time due to reducing the resolution of the pressure filed while maintaining excellent to reasonable accuracy, depending on the level of coarsening.

In this article, we proposed a semi-Implicit-time-integration Finite-element version of the CGP (IFE-CGP). The newly added semi-implicit time integration feature enabled CGP to run simulations with large time steps, and thus further accelerated the computations compared to the standard/previous CGP algorithms. The new data structure introduced resulted in nearly zero computational cost for the mapping procedures. Using the finite element discretization, CGP was adapted to be suitable for complex geometries and realistic boundary conditions. Moreover, the mapping functions were conveniently derived from the finite element shape functions.

In order to examine the efficiency of the IFE-CGP method, we solved three standard test cases: The Taylor-Green vortex in a non-trivial geometry, flow over a backward-facing step, and flow past a circular cylinder. The speedup factors ranged from 1.601 to 29.534. The minimum speedup belonged to the Taylor-Green vortex problem with velocity Dirichlet boundary conditions, while the maximum speedup was found for the flow over a backward-facing step with open boundary conditions. Generally, the outputs for one and two levels of the Poisson grid coarsening agreed well with those computed using full fine scale computations. For three levels of coarsening, however, only a reasonable level of accuracy was achieved.

The mesh refinement application of the CGP method was introduced for the first time in this work. Based on it, if for a given spatial resolution the numerical simulation diverged or the velocity outputs were not accurate enough, instead of refining both the advection-diffusion and the Poisson grids, the IFE-CGP mesh refinement suggests to only refine the advection-diffusion grid and keep the Poisson grid resolution unchanged. The application of the novel mesh refinement tool was shown in the cases of flow over a backward-facing step and flow past a cylinder. For the backward facing step flow, a three-level partial mesh refinement made a previously diverging computation numerically stable. For the flow past a cylinder, the error of the viscous lift force was reduced from 31.501% to 7.191% (with reference to the standard mesh refinement results) by the one-level partial mesh refinement technique.

We showed that the prolongation operator of IFE-CGP did not thicken the artificial layers that arose from the artificial Neumann boundary conditions. Additionally, we demonstrated that although CGP reduces the accuracy level of the pressure field, it conserves the accuracy of the pressure gradient, a key to the efficacy of the CGP method.

For future studies, the contribution of different incremental pressure projection schemes (such as standard [44], rotational [27], and vectorial forms [45]) to the CGP methodology will be analyzed. It is conjectured that because the above mentioned methods diminish the errors introduced by spurious Neumann conditions, applying a CGP technique to the incremental projection methods can lead to undamped flows in comparison with applying the CGP strategy to non-incremental pressure correction computations. Another objective of our future research is to perform a comparison between the CGP approach with one level of coarsening () and the standard finite element algorithm with Taylor-Hood mixed finite elements (P/P) (see e.g., Refs. [1, 6]). From a grid resolution point of view, for an assumed number of grid points of the velocity component, the Poisson solver utilizes a space with an equal pressure node numbers, discretized using either the IFE-CGP () method or Taylor-Hood elements. In this sense, a detailed investigation of the similarity/difference between these two concepts may introduce novel mapping functions for the IFE-CGP tool.

AK wishes to thank Dr. Michael Lentine, Dr. Peter Minev, Dr. Saad Ragab, and Dr. Omer San for helpful discussions.


  • (1) Guermond J, Minev P, Shen J (2006) An overview of projection methods for incompressible flows. Computer methods in applied mechanics and engineering 195 (44):6011-6045
  • (2) Chorin AJ (1968) Numerical solution of the Navier-Stokes equations. Mathematics of computation 22 (104):745-762
  • (3) Temam R (1969) Sur l’approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (II). Archive for Rational Mechanics and Analysis 33 (5):377-385
  • (4) Shen J (1992) On error estimates of projection methods for Navier-Stokes equations: first-order schemes. SIAM Journal on Numerical Analysis 29 (1):57-77
  • (5) Hugues S, Randriamampianina A (1998) An improved projection scheme applied to pseudospectral methods for the incompressible Navier–Stokes equations. International Journal for Numerical Methods in Fluids 28 (3):501-521
  • (6) Jobelin M, Lapuerta C, Latché J-C, Angot P, Piar B (2006) A finite element penalty-projection method for incompressible flows. Journal of Computational Physics 217 (2):502-518
  • (7) Fedkiw R, Stam J, Jensen HW Visual simulation of smoke. In: Proceedings of the 28th annual conference on Computer graphics and interactive techniques, 2001. ACM, pp 15-22
  • (8) Korczak KZ, Patera AT (1986) An isoparametric spectral element method for solution of the Navier-Stokes equations in complex geometry. Journal of Computational Physics 62 (2):361-382
  • (9) Guermond J-L, Minev P, Shen J (2005) Error analysis of pressure-correction schemes for the time-dependent Stokes equations with open boundary conditions. SIAM Journal on Numerical Analysis 43 (1):239-258
  • (10) Reusken A (1995) Fourier analysis of a robust multigrid method for convection-diffusion equations. Numerische Mathematik 71 (3):365-397
  • (11) Filelis-Papadopoulos CK, Gravvanis GA, Lipitakis EA (2014) On the numerical modeling of convection-diffusion problems by finite element multigrid preconditioning methods. Advances in Engineering Software 68:56-69
  • (12) Gupta MM, Kouatchou J, Zhang J (1997) A compact multigrid solver for convection-diffusion equations. Journal of Computational Physics 132 (1):123-129
  • (13) Gupta MM, Kouatchou J, Zhang J (1997) Comparison of second-and fourth-order discretizations for multigrid Poisson solvers. Journal of Computational Physics 132 (2):226-232
  • (14) Zhang J (1998) Fast and high accuracy multigrid solution of the three dimensional Poisson equation. Journal of Computational Physics 143 (2):449-461
  • (15) Lentine M, Zheng W, Fedkiw R A novel algorithm for incompressible flow using only a coarse grid projection. In: ACM Transactions on Graphics (TOG), 2010. vol 4. ACM, p 114
  • (16) San O, Staples AE (2013) A coarse-grid projection method for accelerating incompressible flow computations. Journal of Computational Physics 233:480-508
  • (17) San O, Staples AE (2013) An efficient coarse grid projection method for quasigeostrophic models of large-scale ocean circulation. International Journal for Multiscale Computational Engineering 11 (5)
  • (18) Jin M, Liu W, Chen Q (2014) Accelerating fast fluid dynamics with a coarse-grid projection scheme. HVAC&R Research 20 (8):932-943
  • (19) Losasso F, Gibou F, Fedkiw R Simulating water and smoke with an octree data structure. In: ACM Transactions on Graphics (TOG), 2004. vol 3. ACM, pp 457-462
  • (20) Moin P (2010) Fundamentals of engineering numerical analysis. Cambridge University Press
  • (21) Gottlieb S, Shu C-W (1998) Total variation diminishing Runge-Kutta schemes. Mathematics of computation of the American Mathematical Society 67 (221):73-85
  • (22) Gropp WD, Kaushik DK, Keyes DE, Smith BF (2001) High-performance parallel implicit CFD. Parallel Computing 27 (4):337-362
  • (23) Heys J, Manteuffel T, McCormick S, Olson L (2005) Algebraic multigrid for higher-order finite elements. Journal of Computational Physics 204 (2):520-532
  • (24) Becker R, Braack M (2000) Multigrid techniques for finite elements on locally refined meshes. Numerical linear algebra with applications 7 (6):363-379
  • (25) Reitzinger S, Schöberl J (2002) An algebraic multigrid method for finite element discretizations with edge elements. Numerical linear algebra with applications 9 (3):223-238
  • (26) Sampath RS, Biros G (2010) A parallel geometric multigrid method for finite elements on octree meshes. SIAM Journal on Scientific Computing 32 (3):1361-1392
  • (27) Timmermans L, Minev P, Van De Vosse F (1996) An approximate projection scheme for incompressible flow using spectral elements. International journal for numerical methods in fluids 22 (7):673-688
  • (28) Wanner G, Hairer E (1991) Solving ordinary differential equations II, vol 1. Springer-Verlag, Berlin
  • (29) Reddy JN (1993) An introduction to the finite element method, vol 2. vol 2.2. McGraw-Hill New York
  • (30) Jiang CB, Kawahara M (1993) A three-step finite element method for unsteady incompressible flows. Computational Mechanics 11 (5-6):355-370
  • (31) Brezzi F (1974) On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers. Revue française d’automatique, informatique, recherche opérationnelle Analyse numérique 8 (2):129-151
  • (32) Babuška I (1973) The finite element method with Lagrangian multipliers. Numerische Mathematik 20 (3):179-192
  • (33) Xu J, Chen L, Nochetto RH (2009) Optimal multilevel methods for H (grad), H (curl), and H (div) systems on graded and unstructured grids. In: Multiscale, nonlinear and adaptive approximation. Springer, pp 599-659
  • (34) Yserentant H (1986) On the multi-level splitting of finite element spaces. Numerische Mathematik 49 (4):379-412
  • (35) Bank RE, Xu J (1996) An algorithm for coarsening unstructured meshes. Numerische Mathematik 73 (1):1-36
  • (36) Hu J (2015) A robust prolongation operator for non-nested finite element methods. Computers & Mathematics with Applications 69 (3):235-246
  • (37) Besson J, Foerch R (1997) Large scale object-oriented finite element code design. Computer Methods in Applied Mechanics and Engineering 142 (1):165-187
  • (38) Bell N, Garland M (2008) Efficient sparse matrix-vector multiplication on CUDA. Nvidia Technical Report NVR-2008-004, Nvidia Corporation
  • (39) Saad Y, Schultz MH (1986) GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical computing 7 (3):856-869
  • (40) Van der Vorst HA (2003) Iterative Krylov methods for large linear systems, vol 13. Cambridge University Press
  • (41) Geuzaine C, Remacle JF (2009) Gmsh: A 3‐D finite element mesh generator with built‐in pre‐and post‐processing facilities. International Journal for Numerical Methods in Engineering 79 (11):1309-1331
  • (42) Taylor G, Green A (1937) Mechanism of the production of small eddies from large ones. Proceedings of the Royal Society of London Series A, Mathematical and Physical Sciences 158 (895):499-521
  • (43) Alam JM, Walsh RP, Alamgir Hossain M, Rose AM (2014) A computational methodology for two‐dimensional fluid flows. International Journal for Numerical Methods in Fluids 75 (12):835-859
  • (44) Goda K (1979) A multistep technique with implicit difference schemes for calculating two-or three-dimensional cavity flows. Journal of Computational Physics 30 (1):76-95
  • (45) Caltagirone J-P, Breil J (1999) Sur une méthode de projection vectorielle pour la résolution des équations de Navier-Stokes. Comptes Rendus de l’Académie des Sciences-Séries IIB-Mechanics-Physics-Astronomy 327 (11):1179-1184
  • (46) Barton I (1997) The entrance effect of laminar flow over a backward‐facing step geometry. International Journal for Numerical Methods in Fluids 25 (6):633-644
  • (47) Kim J, Moin P (1985) Application of a fractional-step method to incompressible Navier-Stokes equations. Journal of Computational Physics 59 (2):308-323
  • (48) Erturk E (2008) Numerical solutions of 2-D steady incompressible flow over a backward-facing step, Part I: High Reynolds number solutions. Computers & Fluids 37 (6):633-655
  • (49) Belov AA (1997) A new implicit multigrid-driven algorithm for unsteady incompressible flow calculations on parallel computers.
  • (50) Behr M, Hastreiter D, Mittal S, Tezduyar T (1995) Incompressible flow past a circular cylinder: dependence of the computed flow field on the location of the lateral boundaries. Computer Methods in Applied Mechanics and Engineering 123 (1):309-316
  • (51) Ding H, Shu C, Yeo K, Xu D (2004) Simulation of incompressible viscous flows past a circular cylinder by hybrid FD scheme and meshless least square-based finite difference method. Computer Methods in Applied Mechanics and Engineering 193 (9):727-744
  • (52) Braza M, Chassaing P, Minh HH (1986) Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder. Journal of fluid mechanics 165:79-130
  • (53) Liu C, Zheng X, Sung C (1998) Preconditioned multigrid methods for unsteady incompressible flows. Journal of Computational Physics 139 (1):35-57
  • (54) Hammache M, Gharib M (1989) A novel method to promote parallel vortex shedding in the wake of circular cylinders. Physics of Fluids A: Fluid Dynamics (1989-1993) 1 (10):1611-1614
  • (55) Rajani B, Kandasamy A, Majumdar S (2009) Numerical simulation of laminar flow past a circular cylinder. Applied Mathematical Modelling 33 (3):1228-1247
  • (56) Wang ZJ (1999) Efficient implementation of the exact numerical far field boundary condition for Poisson equation on an infinite domain. Journal of Computational Physics 153 (2):666-670
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description