9.4 Multigrid 多重网格[cfd-9-4]

Before the geometric multigrid scheme can be applied, the coarser grids have to be generated. The standard way is to coarsen the grid evenly in all coordinate directions. However, Mulder [49] proposed an approach called semicoarsening. Here, the grid is coarsened only in one direction, which is changed from one coarse level to another. The semicoarsening methodology is especially suited for flow problems, where the governing equations are stiff in one spatial direction. An example is the direction normal to the wall in boundary layers. Applications of semicoarsening to the Navier-Stokes equations were reported, e.g., in Refs. [26], [50], [51].

The discretised governing equations on the fine grid read in accordance with the relationship (6.1)

In the following, the finest grid will be denoted by subscript \(h\) in reference to the spacing of the grid lines (characteristic dimension of the control volume on unstructured grids). Starting from a known solution \(\vec{W}^{n}_{h}\), a new solution \(\vec{W}^{n+1}_{h}\) is obtained after one time step with some suitable iterative scheme. A new residual \(\vec{R}^{n+1}_{h}\) is evaluated with this solution. In order to improve the solution \(\vec{W}^{n+1}_{h}\) using a coarse grid, the following three steps are carried out:

The solution is transferred to the coarse grid by means of the interpolation

where the subscript \(2h\) denotes the coarse grid¹ and \(\hat{I}^{2h}_{h}\) is the interpolation operator. The residuals have to be transferred to the coarse grid as well, so that their low-frequency error components can be smoothed. A conservative transfer operator is employed for this purpose. This means that when the control volume size increases, the value of the residual must increase by the same amount. The residuals of the fine grid are also required in order to retain the solution accuracy of the fine grid on the coarse grid. For this purpose, a source term, the so-called forcing function [19], [22], is formed as the difference between the residual transferred from the fine grid and the residual computed using the initial solution \(\vec{W}^{(0)}_{2h}\) (Eq. (9.17)) on the coarse grid, i.e,

Here, \(I^{2h}_{h}\) represents the restriction operator which transfers residuals from the fine to the coarse grid. This type of multigrid scheme is known as the Full Approximation Storage (FAS) method [19]. The FAS method is particularly suited for non-linear equations because the nonlinearities in the system are carried down to the coarse levels through the re-discretisation.

The solution is evaluated on the coarse grid in the same way as on the fine grid. The forcing function in Eq. (9.18) is added to the residual of the coarse grid, i.e.,

Hence, the time-stepping scheme can be written in the form

In the case of the explicit multistage scheme (Subsection 6.1.1), this results in

in accordance with Eq. (6.5). It has to be noted that during the first iteration (stage in Eq. (9.21)), \((\vec{R}_F)_{2h}\) is identical to the residual transferred from the fine grid (i.e., \((\vec{R}_F)_{2h} = I^{2h}_{h}\vec{R}^{n+1}_{h}\) from Eq. (9.18)). This guarantees that the solution on the coarse grid depends on the residual of the fine grid and thus retains the accuracy of the fine grid.

An important question is the accuracy of the spatial discretisation scheme on coarse grids. Since the coarse grids do not influence the accuracy of the fine-grid solution, first-order schemes are sufficient. The advantages of first-order accurate discretisation on coarse grids are the increased robustness, better damping properties, and lower numerical effort in comparison to higher-order schemes.

After one or several time steps (iterations) were carried out on the coarse grid, the correction with respect to the initial - interpolated - solution (Eq. (9.17)) is computed. This so-called coarse grid correction is given by

The coarse-grid correction is interpolated to the fine grid in order to improve the solution there. Hence, the new solution on the fine grid reads

where \(I^{h}_{2h}\) is denoted as the prolongation operator.

The basic multigrid scheme described above consists of one coarse grid only. If multiple coarse grids are present, steps 1 and 2 are repeated until the coarsest grid is reached. It is important to realize that the forcing function on the coarse grids is formed from the restricted corrected residual of Eq. (9.19)). For example, on the coarse grid \(4h\), the forcing function is obtained from

In this way, the residual of the finest grid controls the accuracy of the solution on all coarse grids. After a given number of time steps on the coarsest grid, step 3 can be successively repeated until the finest grid is reached again. This procedure is known as a saw-tooth or V-cycle (see Fig. 9.4a). However, it is also possible to conduct more cycles on the coarse grids. This strategy, termed the W-cycle, is displayed in Fig. 9.4b. It is employed particularly frequently for transonic flows. In the case of supersonic and hypersonic flows, the V-cycle proved to be more efficient.

The optimum number of time steps before the restriction and after the prolongation depends on the type of the time-stepping scheme. In the case of the explicit multistage scheme (Subsection 6.1.1 or 6.1.2), it is common to carry out only one time step before the restriction of residuals and no time step after the prolongation. However, the robustness of the multigrid scheme can be improved by smoothing the coarse grid corrections (Eq. (9.22)) before adding them to the fine grid solution \(\vec{W}^{n+1}_{h}\) (Eq. (9.23)). The same central implicit smoothing (with constant coefficients) as described in Section 9.3 is utilised.

The other popular time-stepping method, the implicit LU-SGS scheme (see Subsection 6.2.4), requires two iterations before the restriction for the best multigrid efficiency [34]. The number of time steps after the prolongation depends on the spatial discretisation. In the case of the central scheme (Subsection 4.3.1), no time step is necessary [34]-[36], but the solution correction can be smoothed. On the contrary, one time step should carried out after the prolongation if an upwind spatial discretisation is used. This (2,1)-strategy proved to be an optimum with respect to robustness and computing time for various flow conditions [35], [36].

It should be pointed out that in practice the multigrid scheme is not started directly from the finest grid. Instead, several multigrid cycles are executed from one of the coarse grids. The approximate solution is interpolated to the next finer grid (using the same operator as for the prolongation), few more cycles are performed, the solution is again interpolated to the next finer grid and so on, until the finest grid is reached. In this way, a good starting solution is obtained on the finest grid with only a moderate numerical effort. This very efficient procedure is termed the Full Multigrid (FMG) method [19].

The restriction (9.18) and the prolongation (9.23) operator must fulfil certain accuracy requirements, namely [52]

where \(m_R\) and \(m_P\) denote the degree plus 1 of the polynomial, which is exactly interpolated by the restriction and the prolongation operator, respectively. For example, \(m_R\) or \(m_P\) are equal to two in the case of linear interpolation. Furthermore, \(m_E\) represents the order of the governing equations. Thus, \(m_E = 1\) for the Euler equations, and \(m_E = 2\) in the case of the Navier-Stokes equations. If the condition (9.25) is violated, the additional errors introduced by the restriction and/or prolongation will disturb the fine-grid solution. Hence, such multigrid scheme will converge only slowly or it will even diverge.

The implementation of multigrid on structured grids is straightforward since the coarse grids can be easily generated by deleting every second grid line in the respective coordinate direction. The spacing of the grid lines is therefore \(2h\), \(4h\), etc. This guarantees that the numerical effort on the coarse grids stays low as compared to the finest grid. Several representative examples are provided in Refs. [11]-[13], [20]-[22], [26], [32]-[36].

As we can conclude from Fig. 9.5, the operators for the solution interpolation, the restriction of residuals and the prolongation of corrections have to be defined differently for cell-centred and node-centred (cell-vertex) schemes. For example, two successive grids have every second grid point in common. On the contrary, the cell centres are always at different locations. Therefore, we shall discuss the standard forms of the transfer operators separately for the cell-centred and node-centred (identical to cell-vertex) finite-volume schemes.

Apart from the symmetrical, purely geometrically defined restriction and prolongation operators, which will be presented next, upwind-biased forms were suggested in [16], Chapter 4. Upwind restriction and prolongation, which account both for the characteristics of the flow equations, improve the robustness of the multigrid scheme for hypersonic flows. In order to save space, we shall discuss the implementation of an upwind prolongation operator for the node-centred spatial discretisation only.

The solution is transfered from the fine to the coarse grid by using a volume weighted interpolation. In 2D, Eq. (9.17) becomes (see Fig. 9.6a)

Similar transfer operator is employed in 3D, where the summation is over the eight fine-grid control volumes, which form one coarse-grid cell.

The restriction operator is defined as a sum of the residuals from all cells which are contained in one coarse-grid control volume. Hence, in 2D we have (cf. Fig. 9.6a)

and likewise in 3D. Those residuals \(\vec{R}^{n+1}_{h}\) in Eq. (9.27), which are located outside the physical domain, have to be set to zero.

The prolongation of the coarse-grid correction Eq. (9.23) can be conducted in two different ways. First, a zeroth-order prolongation operator results, if \(\delta\vec{W}_{2h}\) is equally distributed to all surrounding cell centres as indicated in Fig. 9.6b. Thus, for instance

The second possibility, which leads to a faster convergence of the multigrid scheme, consists of two steps. In a first step, \(\delta\vec{W}_{2h}\) is interpolated to the grid nodes like for the cell-vertex scheme (see below). In a second step, the nodal values are averaged to obtain the value in the centre of the fine-grid cell. Referring to Fig. 9.6c, the following final relationship can be derived

The corresponding expression in 3D can be found by a similar procedure. It reads

Because of the common nodes between the fine and coarse grid, the solution can be transfered simply by injection, i.e.,

The standard central restriction operator represents a linear interpolation from the nodes of all four (eight in 3D) fine-grid cells, which resemble one coarse-grid cell. According to Fig. 9.7, the restricted residual is computed in 2D as

In 3D, the fine-grid residuals are collected as follows

with the factors

In the above Eq. (9.34), only those indices are shown which are different from \(i, j, k\). It should be noted that the residuals \(\vec{R}^{n+1}_{h}\) must be set to zero at all physical boundary and dummy points before the restriction.

The prolongation of the coarse-grid correction can be implemented as a loop over the points of the coarse grid. Within the loop, the values \((\delta\vec{W}_{2h})_{i,j,k}\) are distributed to the fine-grid points using the same weights as for the restriction (cf. Fig. 9.7b). The particular contributions are summed up in order to obtain the complete transfered correction at each point of the fine grid.

In principle, upwind prolongation can be formulated either in characteristic or in conservative variables. A particularly efficient implementation in conservative variables was proposed in [16]. The methodology employs upwind-biased interpolation of the corrections in Eq. (9.23) according to the Mach number and the velocity direction. The numerical effort is very low, nevertheless the robustness of the multigrid scheme can be significantly improved for high Mach-number flows [16].

The interpolation of the solution corrections \(\delta\vec{W}_{2h}\) to the finer grid is accomplished in two steps. In the first step, the corrections at the points \(A\), \(B\), \(C\), and \(D\) (see Fig. 9.8), which are common to both grids, are transferred directly to the finer grid. In the second step, the corrections are interpolated to the points \(e\), \(f\), \(e'\), \(f'\) and \(g\), which are contained only on the finer grid. The interpolation depends on the sign and the absolute value of the Mach number at the corresponding point. The Mach number in this case is calculated using

is employed. The normal vector \(\vec{n}\) in Eq. (9.35) can be obtained either by averaging the face vectors of the control volume (in the direction \(A-B\)), or by normalising the vector from point \(A\) to point \(B\). Then, the correction is transferred to point \(e\) according to the rule

The same procedure applies also to the points \(f\) to \(f'\). The upwinding helps to match the information exchange between the grids better to the real physics. The interpolation to the point \(g\) is more difficult. In Ref. [16], the values at the surrounding points \(A\) to \(D\) were simply averaged

However, some sort of upwind weighted interpolation would be more appropriate. The upwind prolongation can be implemented in similar way also in 3D. Despite the simplification, encouraging results were obtained in a number of test cases [16].

As compared to the structured grids, the construction of the coarse grids is much more involved in the case of unstructured grids. The problem is how to construct an uniformly coarsened grid from a set of elements (grid cells) which have no particular ordering. Additionally, the ratio of the cell volumes of the coarse to the fine grid has to stay within a certain margins (about 4 in 2D and 8 in 3D). One possibility how to solve this problem is to apply the AMG methodology [41]-[48], which we briefly discussed at the beginning of Section 9.4. However, the geometric multigrid is still more widely used. Therefore, we shall concentrate here on this approach.

Three main methods for the generation of coarse grids can be identified:

• nonnested-grids approach,
• topological methods, and
• agglomeration of control volumes.

Reviews of the above methods were presented in Refs. [53] and [54].

The standard restriction and prolongation operators are based on purely geometrically defined interpolation. Leclerq and Stoufflet [55] suggested upwind transfer operators, which are particularly promising for flows with strong shocks. Their upwind restriction/prolongation is based on the transformation of the residuals/corrections into the characteristic variables. After an upwind-biased interpolation, the restricted/prolongated values are transformed back into the physical variables. Numerically much less expensive upwind multigrid method was presented in Ref. [16] (see also Subsection 9.4.3, Eq. (9.36)).

The most obvious idea is to generate a sequence of completely independent, increasingly coarser grids [56]-[61]. It is not necessary that the grids contain any common nodes. Therefore, we speak of nonnested grids. However, it is important that the main geometrical features (leading and trailing edge, fuselage nose, etc.) are retained on all coarse grids. This is not easy to accomplish, particularly in the case of a geometrically complex configuration. Multigrid based on nonnested grids is hardly used today.

One particular approach applies graph-based algorithms in order to remove certain nodes from the fine grid. The remaining nodes are then re-triangulated [62], [63]. On the contrary to the nonnested-grids approach, the interpolation between the grids becomes easier, since the successive grids contain common nodes. However, the method inherits the drawback of the nonnested grids with respect to geometry conformance.

A further topological method employs grid refinement [64], [29], [65]. The technique starts from a coarse grid and generates finer grids by element division. The methodology can be applied either over the whole physical domain or only locally (e.g., at boundary layers). The disadvantage of this approach is that the quality of the finest grid strongly depends on the initial coarse grid and the refinement procedure. The problem can be partially cured by edge swapping [66], [67].

Another idea for the generation of coarse grids is based on edge collapsing [68]. It was initially developed for inviscid flows on tetrahedral grids. The edge-collapsing method was further extended to viscous flows on mixed-element grids in [69].

A very efficient methodology for unstructured grids is the so-called agglomeration multigrid. It was first presented by Lallemand [70], Lallemand et al. [71] and by Koobus et al. [72]. Later on, the agglomeration multigrid was adopted by various authors [73]-[77], [31]. The method generates a coarse grid by fusing the control volumes of the finer grid with their neighbours. The resulting coarse grids consist of successively larger, irregularly shaped polyhedral cells. This is depicted in Fig. 9.9. As we can see, the agglomeration technique retains the full discretisation of the boundary surfaces. This represents a significant advantage over all previously discussed unstructured multigrid methods. However, it should be mentioned that up to now, implementations of the agglomeration multigrid were based mostly on the median-dual cell-vertex scheme (Subsection 5.2.2). The application of the agglomeration multigrid to a cell-centred scheme (Subsection 5.2.1) was described in [73].

The volume agglomeration for a node-centred scheme proceeds in the following steps:

• 1. build a list of the so-called seed points. Seed points are grid points selected to agglomerate the surrounding control volumes. The list of seed points can contain either those points which form an approximate maximal independent set [75], or simply all points of the current grid level.
• 2. Loop over all seed points.
• 3. If the seed point is unagglomerated, agglomerate all its nearest neighbours (connected by an edge), which were not already agglomerated.
• 4. Check the coarsening ratio (i.e., how many fine-grid control volumes are contained within a coarse-grid volume). If the ratio is less than four (eight in 3D), the neighbours of the already agglomerated nearest neighbours are added (if not associated with another seed point), until the optimum coarsening ratio is achieved. In Ref. [27], it was proposed to agglomerate those distance-two neighbours first, which are connected to at least two (three in 3D) agglomerated nearest neighbours.
• 5. If there are still seed points in the list, goto step 2.
• 6. Eliminate singletons. These are single control volumes which could not be agglomerated, because there were no unagglomerated neighbours. A singleton can be eliminated by agglomeration with such neighbouring control volume, which has the smallest coarsening ratio. This leads to coarse-grid levels with a more regular distribution of control-volume areas [31].

The above procedure is repeated until all coarse grids are generated.

The volume agglomeration has to start from the boundary in order to preserve grid isotropy. In 3D, user intervention may be required to prescribe the agglomeration direction depending on the shape of the boundary. To overcome this difficulty, Okamoto et al. [77] proposed another algorithm denoted as global coarsening. The method employs a global partitioning scheme, which is based on edge colouring. The partitioning scheme is used to generate an independent set of edges. In a second step, all control volumes, which share an edge of the independent set, are agglomerated. The procedure is repeated until the prescribed coarsening ratio is achieved. The method does not require the specification of an initial seed point or agglomeration direction. It can also treat any type of grid cells.

The implementation of agglomeration multigrid presents little difficulty for inviscid flows. The Euler equations are in general discretised as fluxes over individual control-volume faces. In this respect, it does not matter how complex is the shape of the control volume (in fact, averaged face vectors are used). Furthermore, a first-order accurate spatial scheme requires the knowledge of flow quantities in the neighbouring control volumes only. This information is readily available on the coarse grids.

In the case of viscous flows, the discretisation of the diffusive fluxes on arbitrary shaped control volumes is no longer straightforward. The problem is the evaluation of gradients at face midpoints (see the discussion in [31]). Another, and even more serious, difficulty is related to the required accuracy of the prolongation operator. The inequality in Eq. (9.25) suggests \(m_P = 2\), since the common restriction operator (sum of residuals) leads to \(m_R = 1\) only. However, the construction of a linear interpolation is not easy on the coarse grids.

In order to circumvent the construction of a first-order accurate prolongation operator, Mavriplis [53] proposed to use constant prolongation (i.e., all points of the fine grid contained within an agglomerated coarse-grid volume get the same solution correction) and a scaling of the viscous fluxes. However, this approach does not lead to optimal multigrid efficiency.

Haselbacher [31] suggested to retain the fine-grid discretisation of the viscous fluxes also on the coarse grids and to enforce the boundary conditions like on the finest grid. Moreover, he proposed a piecewise linear prolongation operator. For a scalar quantity \(U\), it can be written as

The gradient \((\nabla\delta U_{2h})_i\) in Eq. (9.38) is calculated by using the linear least-squares reconstruction described in Subsection 5.3.4, Eq. (5.55). The values of the limiter function \(\Psi_i\) are evaluated according to the Barth-Jespersen limiter function presented in Subsection 5.3.5.