{'fraction_non_alphanumeric': 0.032386402112819124, 'fraction_numerical': 0.015575269181282589, 'mean_word_length': 4.435130658814869, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finite volume adaptive solutions using SIMPLE as smoother

This paper describes a new multilevel procedure that can solve the discrete Navier-Stokes system arising from finite volume discretizations on composite grids, which may consist of more than one level. SIMPLE is used and tested as the smoother, but the multilevel procedure is such that it does not exclude the use of other smoothers. Local refinement is guided by a criterion based on an estimate of the truncation error. The numerical experiments presented test not only the behaviour of the multilevel algebraic solver, but also the efficiency of local refinement based on this particular criterion.


INTRODUCTION
The SIMPLE algorithm ( [1]) and its variants (see [2] for descriptions and references) are popular methods for the solution of discrete non-linear systems of equations arising from the discretisation of the incompressible Navier -Stokes equations. Like most single-grid solution methods they have the disadvantage that the larger the grid (in terms of number of nodes), the larger the number of iterations required to attain a specified level of convergence.
Multigrid ( [3], [4], [5]) is a class of methods which overcome this disadvantage by using a series of progressively coarser grids in addition to the finest one. Any solver for systems of algebraic equations may be used in the context of a multigrid method, and in this case it is called a smoother rather than a solver. In a multigrid method, the error of a current approximation compared to the exact solution of the system is considered to be composed of Fourier components, and each grid is responsible for the reduction of those components which may be represented on it and which oscillate rapidly compared to the grid spacing. The efficiency of the method depends on the ability of the iterative solver to reduce the rapidly oscillating components of the error faster than the smooth components on each grid. After a number of iterations of such a solver the error appears smoother than before, and therefore in the context of multigrid such solvers are called smoothers. By transferring the solution procedure between grids of various densities, every component of the error will appear oscillatory on some grid, where it will be efficiently reduced by the smoother. Therefore, ideally, the number of iterations (or cycles) required in multigrid methods to attain a specified level of convergence is small and independent of the density of the finest grid.
Originally Gauss-Seidel and line relaxation smoothers were used, but beginning in the late 1980s several studies have shown that SIMPLE-like procedures also have smoothing properties. In [6] and [7] this was shown for staggered grids, while in [8] an algorithm was derived for colocated grids. In [9] the method was tried on curvilinear grids, demonstrating high efficiency, and in [10] and [11] the method was tried on three-dimensional cases. In [12], [13] and [14] it was shown that the method can handle turbulence models successfully, although the gains are not in general as high as in the laminar case. In [15] the effect of the discretisation scheme for the convection terms on the overall efficiency was studied. In [16] and [17] it was shown that the performance of SIMPLE as a smoother is comparable to that of SCGS (Gauss-Seidel -like smoother) and SCAL (line-relaxation -type smoother). Also, in [18] SIMPLE is compared with other methods of the same family and it is shown that their performance as smoothers is not significantly different. A more specialised use of the procedure is presented in [19], where SIMPLE is the smoother of a multigrid procedure which is used as a preconditioner for a Newton-Krylov method.
In all these studies, except maybe [13], each of the grids used by the multigrid procedure covers the entire computational domain. There are cases though where it would be convenient that each grid is allowed to cover a subset of the domain. Such grids are more appropriately termed levels, and the corresponding solution techniques are called multilevel rather than multigrid techniques. Multilevel techniques are convenient in the case of locally refined grids, where finer levels are added on top of an original grid only in regions of the flow field where enhanced spatial resolution is required. Multilevel techniques date back to the work of Brandt [4], who defined the technique known as MLAT (MultiLevel Adaptive Technique) for finite difference discretisations. In the case of finite volume discretisations some caution is needed in constructing the multigrid equations for the control volumes which are located at interfaces between different levels. In [5] a method is proposed for composite grids (i.e. grids composed of different levels) whose nodes identify with the centres of the control volumes. In the present work a method will be proposed for grids whose nodes identify with the vertices of the control volumes.
The classical SIMPLE/multigrid formulation as described in the studies mentioned above cannot readily be applied in the case of locally refined (composite) grids. In [13] results of using SIMPLE as smoother on composite grids are presented, but the method described is precisely the classic one, and the difficulties which arise at the interfaces between different levels are not addressed. These will be addressed in the present work. A different approach is adopted in [20], where the composite grid is partitioned into blocks, each assigned to a different processor, and each block is covered by a single level. The classic SIMPLE/multigrid method is used within each block, and information is exchanged at the interfaces between blocks, between successive multigrid cycles. However, in the present work the composite grid will not be partitioned into blocks but treated as a whole.
In this paper the emphasis is on how to extend the multigrid procedure to be used on composite grids, but one also has to choose an appropriate refinement criterion on which the construction of the composite grid will be based. In some studies, e.g. [20], [21], the grid is refined before the calculations are performed, at regions where the flow field is expected to exhibit strong variation e.g. boundary layers. A more useful approach is to solve the equations on a given grid and then refine those regions of the grid which are indicated by a criterion which assesses the quality of the solution. This can be done by examining the variation of the solution in relation with the grid spacing. For example, a second order accurate discretisation, like the one used in this paper, is constructed by assuming that the control volumes are small enough such that the flow variables vary almost linearly within the neighbourhood of each control volume. Therefore, after a solution has been obtained on a given grid, the refinement criterion can refine those control volumes where this assumption does not hold within a selected tolerance. Such refinement criteria can be based on estimates of the truncation error (e.g. [22], [5]), the finite element residual (e.g. [23], [24]) or some other similar quantity (e.g. [13]). These quantities have the same dimensions as the flux of the transported quantity (e.g. mass or momentum) but some criteria (e.g. [24], [13]) normalise them by some factor (usually the main diagonal coefficient of the linearised algebraic system) to convert them into a quantity which has the same dimensions as the flow variables. In this paper the refinement criterion is based on the truncation error estimate presented in [25], which is not normalised but the goal is to reduce the integral of the truncation error in the domain below some selected value.

DISCRETISATION OF THE DOMAIN INTO FINITE VOLUMES
The continuous space of the flow field under study is divided into quadrilateral control volumes (CVs), which are organised into levels (see Figure 1). The grid is constructed as follows: First, level 1 is constructed (say, by a structured grid generator) which covers the entire computational domain and is relatively coarse. Then a second level is added, whose CVs are formed by subdivision of some (or all) of the CVs of level 1. The subdivision procedure produces four (4) "children" CVs out of each chosen CV of the previous level, by 3 joining the parent CV centre with the midpoints of its 4 sides. Level 3 is similarly constructed by subdividing some (or all) of the CVs of level 2 and so on. A restriction is posed on this procedure: at an interface between two different levels, say k and l, the two levels which meet must be successive, i.e. |k-l|=1.
Once a CV is subdivided into children it is not destroyed but it is retained in the data structure because it will be needed by the multilevel procedure. Therefore, each level consists of its global part, i.e. the CVs which do not have children, and its local part, i.e. the CVs which have children. The sum of the global parts of all levels forms the composite grid, which is the actual grid used for the discretisation of the Navier-Stokes equations.
A composite grid will be denoted by a letter such as h which will also be abstractly interpreted as the distribution of the sizes of the CVs of the given composite grid in the computational domain. A multiple of h, say a•h, will denote a composite grid whose CV at any point of the domain has size (length) a times the size of, and is geometrically similar to, the corresponding CV of grid h at the same location. Further, it is required that the CV relates to its neighbours in a similar way that the corresponding CV of grid h relates to its own neighbours. For an arbitrary grid h it is generally not possible to find a grid which fulfils the requirements for being a•h, but often one may find grids which approximately fulfil these requirements. An important special case is the following: Some composite grids have an underlying grid which is formed by removing all global CVs from the original grid. Therefore, the global CVs of the underlying grid are precisely those local CVs of the original grid which do not have "grandchildren". For an underlying grid to exist, all local CVs of the original grid must have only local siblings (siblings are the CVs which have a common parent). For example, the grid of Figure 1 does not have an underlying grid. This grid is reproduced in Figure 2(a), where the global siblings of local CVs are highlighted in grey. Once these CVs are refined, as in Figure 2(b), the underlying grid shown in Figure 2(c) becomes available. One can ensure that an underlying grid exists at all stages of the finite volume solution by requiring that whenever a CV is marked for refinement, say due to high truncation error, all it's siblings are marked as well. In the present method, the underlying grid is used to estimate the truncation error.
In the present method the CVs are considered to be logically polygonal, an approach suggested for example in [13], [2] and [26]. According to this approach the surface which forms the boundary of each CV is considered to consist of a number of faces, each of which separates the given CV from another single distinct CV. Therefore, although the CVs of the present method are quadrilaterals, those CVs which are located at interfaces between different levels may have more than 4 faces. For example, CV P of Figure 3 has 6 faces and 6 neighbours. The faces which separate two CVs of different level will be called exterior, while the rest will be called interior, except for boundary faces which coincide with the domain boundary.
If  is a function defined on the computational domain and the computational domain is discretised by a grid h, then the grid function h is the discrete function which assumes the values of  at the centres of the CVs of grid h. The P-th component of h, which equals the value of  at the centre of CV P, is denoted as h,P or (h)P.

THE EQUATIONS AND THEIR DISCRETISATION
In Cartesian coordinates, the two-dimensional stationary incompressible Navier -Stokes equations and the continuity equation integrated over a control volume P may be written as: P N (3.8) where:  (3.7). The definition (3.8) of the discrete mass flux is taken from [25], and it contains a pressure term to avoid the appearance of pressure oscillations in the discrete solution. The real factor ami is included for better control of the pressure term. Usually one can use ami=1, except for some rare cases of very coarse grids where ami < 1 is necessary to avoid convergence problems (see [25] or [27] for more details).
In the experiments of Section 6, ami=1 is used unless otherwise specified.
With the above definitions, the exact operators N x , N y , N c (3.1) 12) where ΔΩP is the volume of CV P, and similarly for N y , y h N , y h  . If (3.12) were substituted into (3.1) and the system was solved one would obtain the exact values uh, vh, ph at the CV centres. However, this is not possible since the truncation errors are not known. According to the finite volume methodology, the truncation errors are dropped from (3.12) under the assumption that they have small magnitude, and instead of (3.1) the following algebraic system is solved: The solution h u  , h v  , h p  of (3.13) differs from the exact solution uh, vh, ph of (3.1) by the For arbitrary grids, it is shown in [27] that the contributions of the convection and pressure terms of (3.7) to the truncation error have a magnitude of O(h), while the contributions of the viscous terms have a magnitude of O(1). However on structured grids which come from the discretisation of smooth fields of dimensionless curvilinear coordinates, where the grid skewness and expansion vary smoothly between opposite faces of a CV and tend to zero with grid refinement, parts of these contributions cancel out between opposite faces of a CV and what remains is that . Therefore many numerical experiments, e.g. in [25], [27], show that on smooth and therefore the truncation error contributions of (3.7) and (3.8) of the face which is opposite to the boundary face are not cancelled out. However, as numerical experiments have shown, it appears that these truncation errors of the boundary CVs do not affect the overall O(h 2 ) convergence rate of the discretisation errors. Theoretical explanations for this are proposed in [28] and [29] h . This is because skewness and expansion of exterior faces are large and do not reduce with grid refinement (see Figure 3). However, experiments in Section 6 of the present work suggest that again this does not affect the O(h 2 ) convergence rate of the discretisation error. This may be due to cancellation between parts of the truncation errors of the coarse CV and the pair of fine CVs which meet at a level interface (e.g. CVs P, N3, N4 of Figure 3). Of course these truncation errors occur at different locations in space, but as the grid spacing h reduces such triads of CVs tend to identify with single points in space.

TRUNCATION ERROR ESTIMATE AND LOCAL GRID REFINEMENT CRITERION
The refinement criterion which will be used in the context of this work is based on an estimate of the truncation error. The truncation error is estimated with the method proposed in [25]. The truncation error estimate x h   is given by: [ΔΩ] -1 is a diagonal matrix whose P-th diagonal element equals 1/ΔΩP, the reciprocal of the volume of CV P of grid 2h. The linear operator 2h h I is called a restriction operator and it transfers a grid function from grid h to grid 2h. Likewise the prolongation operator 2 h h I transfers a grid function from grid 2h to grid h. It was found that the particular choice of the restriction and prolongation operators does not have a strong impact on the overall efficiency of the local refinement method. The 3 rd order accurate restriction operator and the prolongation operator used in [25] are also adopted in the present work. Finally, p is the order of accuracy of the discretisation, which is assumed to be p = 2.
The local refinement criterion used here has as an ultimate goal to reduce the integral of the truncation error below some limit, while simultaneously distributing it roughly equally among all CVs of the grid. The limit can be set as follows: Given reference momentum and mass fluxes Qmon and Qmas respectively, which are of the order of the momentum and mass fluxes in the computational domain, and an appropriate real number rτ << 1, the goal is to reduce the integrals of the momentum and continuity truncation errors below rτQmom and rτQmas respectively. The (approximate) integrals of the truncation errors are defined as P|τh,P|ΔΩP for each truncation error, where summation is over all CVs P of grid h, and ΔΩP is the volume of CV P. This goal may be achieved by marking for refinement every CV P of grid h which does not fulfil the following condition: where #Vh is the number of CVs in grid h. When (4.2) is satisfied by all CVs of the grid then P|τh,P|ΔΩP  rτQ as can be easily seen by summing (4.2) of all CVs of the grid (#Vh in number). Of course since the exact truncation errors are not known one has to use the estimates (4.1). After refinement of all CVs which do not satisfy (4.2) a new grid arises on which the equations are solved. Then again refinement takes place according to the criterion (4.2) and so on, until a grid is reached such that all its CVs satisfy (4.2). The criterion (4.2) is dynamic in the sense that a CV of a grid h which satisfies (4.2) may not satisfy it in a subsequent grid h', and thus be refined, because the number of CVs has increased (#Vh' > #Vh). In this way the procedure is aimed at reducing the integrals of the truncation errors below rτQ by refining the CVs with the largest contribution to these integrals. The discrete operators (3.11) of any particular CV, and therefore also the truncation errors, depend not only on the size of P itself but also on the sizes of its neighbours, so if P does not satisfy (4.2) then its neighbours will be refined as well.
A complication arises near the interfaces between different levels of the composite grid where, as has been mentioned before, Therefore, as far as the truncation error is concerned the assumption that p=2 in (4.1) is not accurate there. Furthermore in these regions the truncation error distribution is discontinuous and so 2 h h I in (4.1), which is based on an assumption of smoothness, is not appropriate. Despite these problems, numerical experiments show that the estimate (4.1) predicts the zeroorder accuracy of the discretisation, i.e. that the truncation error estimate is large and does not reduce with refinement. This poses a problem to the local refinement procedure: According to the criterion (4.2) the CVs of the set V e h  V ne h must be refined, but this offers no real benefit as the truncation error will not reduce (τh causing refinement to take place and so on, resulting in perpetual refinement with no real benefit. This behaviour is indeed observed in practice. To overcome this problem in the present work it has been chosen not to allow refinement at these regions. Since the truncation error is effectively estimated on the underlying grid 2h, refinement must not be allowed at the CVs of grid h which cover 2 V e h  2 V ne h . That is, refinement is not allowed within a depth of 4 CVs from the level interface. A similar situation occurs at domain boundaries, where also τh  O(1) as has been mentioned. However, in this case because it is difficult to determine a priori the required grid fineness at the boundary, local refinement is allowed and the number of refinements is limited by setting a maximum number of grids to be constructed.

SIMPLE / MULTILEVEL PROCEDURE FOR LOCALLY REFINED GRIDS
In this section the solution of the discrete system (3.13) on a fixed composite grid, denoted by the capital letter H, is considered. Lowercase letters such as h will now refer to a particular level of the composite grid H. Also, since in this section the interest is not in the exact differential solution u, v, p of (3.1), let now H û , H v , H p be used instead of h u  , h v  , h p  to denote the exact solution of the algebraic system (3.13), and let uH, vH, pH denote an estimate of this solution in the iterative solution procedure. With this notation the algebraic system (3.13) to be solved is written as: Suppose for simplicity that the composite grid H consists of only 3 levels with global parts, h, 2h and 4h, as depicted in Figure 5. During the multilevel cycle, the iterative procedure visits one by one all the levels of the grid and it solves on both the global and local part of each level. Therefore a grid function h defined on a level h will be written as h = (actually the left hand sides of (5.2) also involve variables of some of the CVs of level 2h). If there were a still finer level then level h would also have a local part, and apart from the subset of (5.1) which correspond to level h, (5.2) would also include equations for the local CVs, which would also be satisfied up to a residual. The form of these local equations will be described shortly. Now, the iterative procedure moves to level 2h. First the current solution is restricted from level h to level 2h (including the current solution at the local part of level h) and stored as 2h h h I p  (defined only on the local part of 2h), where 2h h I is a restriction operator which is different from that of Section 4. The experience of the authors has shown that operator (5.3) is suitable also for high Reynolds number flows: where CP is the set of 4 children (of level h) of CV P (of level 2h).
On level 2h the equations solved are such that their solution, say u2h, v2h, p2h, is such that: On the global part of 2h, the new solution 2  which when prolonged back to level h give an estimate which satisfies the equations of level h more closely than the previous estimate The equations of level 2h are the following: Global CVs with only global neighbours (CVs indicated as • and ○ in Figure 5): At these CVs the equations solved are precisely those of the composite grid, i.e. the corresponding subset of (5.1). For CVs ○, the variables at their coarse neighbours + are not treated as unknowns but as fixed Dirichlet boundary values ( 4 . The discretisation stencils of CVs ○ also involve the gradients at CVs +, which in turn involve the variables at CVs ○; therefore in our procedure we have chosen to update the values of the gradients at CVs + after each iteration on level 2h. Concerning the gradients, it should also be mentioned that because the equations of some of these global CVs also make use of the gradients at global CVs of type □, the gradients at CVs □ must be calculated grid-wise (like in Figure 3) and not levelwise, i.e. assuming that the neighbours of the CVs □ are the fine CVs  and not the local CVs ■. For the calculation of these gradients in the present work the values at CVs  were assumed to be fixed at g h u  , g h v  , g h p  rather than updating them after each iteration according to the corrections produced at CVs ■.
Local CVs (CVs indicated as ■ in Figure 5): First we consider the momentum equations. Suppose a local CV P of level 2h, and CP the set of its children for which equations of the type (5.2) hold. According to the multigrid methodology (see e.g. [5] or [27]) if the iteration errors and the residuals per unit volume which correspond to the estimate h u  , h v  , h p  are smooth enough compared to the spacing of level h, then the subset of 4 x-momentum equations of (5.2) which correspond to the children of P may be approximated by a single xmomentum equation of CV P as: The y-momentum equations are treated in the same way. Now we turn to the continuity equation. Normally the continuity equation of local CVs should be defined as (5.4) with superscript x replaced by c. However, following the usual SIMPLE/multigrid algorithm, we depart slightly from this rule and define the local continuity equation similarly to the global continuity equation, as: where: In the above, fP is the set of 4 local faces of the local CV P, and cf is the set of 2 child faces of the local face f. In (5.5), F2h,f are defined by (3.8). If f separates P from a global CV then for g h p  will also be used. In (5.6), Fh,c is to be replaced by Fh,c if the child c is itself local. In other words, the mass flux 2h f F  , is defined as the sum of the mass fluxes through the children of f at the time of the restriction from level h to level 2h; and the mass flux F2h,f is defined as equal to the restricted mass flux 2h f F  , plus a correction (the term in square brackets in (5.5)) which is due to the improvement of the flow field estimate. A slight improvement of this definition of the local mass fluxes will be proposed in Section 6. Summarising, the local continuity equation is defined as F2h,f = 0, like its global counterpart, but the mass fluxes though local faces, F2h, are defined differently than F2h, using (5.5) instead of (3.8). This is necessary so that the local continuity equations will produce corrections which are driven by the residuals of the finer level. In fact, by summing the continuity equations of the children of P, it is not hard to see from definition (5.6) that: Using (5.7), equation (5.5) becomes: which is precisely (5.4) with superscript x replaced by c. So in fact it is not the continuity operator which is defined differently on local and global CVs, but the momentum operator because 2 x h N uses the mass fluxes F2h while 2 x h N uses the mass fluxes F2h.
Global CVs with local neighbours (CVs indicated as □ in Figure 5): This is the most difficult case. The equations must involve unknowns of the same level, and therefore the unknowns at CVs ■. However, at the overall convergence of the multigrid procedure the composite grid equations must be satisfied, which involve the global unknowns at CVs  of level h and not the local unknowns at CVs ■ of level 2h. This is easy to achieve as far as the continuity equation is concerned. Indeed, if P is a CV of type □ and g P, sP are the sets of global and local (sub-exterior) faces of P (regarding P in a level-wise manner, bounded by 4 faces) then its continuity equation has the natural form: At overall convergence of the multigrid cycles the "correction" term in square brackets in (5.5) of the local mass fluxes becomes zero and F2h = 2h F  . Substituting this into (5.9) and using (5.6), (5.9) becomes: (Fh,c also involve values from u2h, v2h, p2h) which is precisely the composite grid equation, i.e. the P-th continuity equation of the system (5.1), as is required. A similar reasoning lies behind the construction of the momentum equations. Let the net x-momentum fluxes (3.7) be written as  Then, with reasoning similar to that for the continuity equation, the x-momentum equation for a CV P of type □ is defined as: That is, the momentum flux through each local face s of P is defined as equal to the sum of This completes the definition of the equations at all CVs of level 2h. After the required number of iterations have been performed at level 2h, iteration may pass to a still coarser level 4h etc. At some point of the multilevel cycle, iteration will return to level 2h and then it will move up to level h. At this point prolongation of the corrections produced at the local part of 2h will occur, to give an improved estimate of the solution of the equations of level h: where P is the parent of C, and P, C are their centres. The initial guess for the solution of the equations of a particular level, say 2h, is û , H v , H p of (5.1), then the initial guess at any level satisfies the equations of that particular level (both local and global) and zero corrections are produced at the local parts. In other words, the multilevel cycle does not alter the exact solution.
Since the adopted discretisation uses the central difference scheme (3.5) (CDS), the solution of the system of equations of any particular level may exhibit oscillations, as is discussed in [25]. The problem of pressure oscillations is taken care of by the use of momentum interpolation for the mass fluxes in both local and global equations of each levelsee [25] (however, see Section 6 for further discussion). Some care is needed to avoid velocity oscillations which may appear especially at coarse levels and high Reynolds numbers. For the multilevel procedure, the problem is twofold: First, during the restriction phase, say from level h to level 2h, any velocity oscillations in the field p   which causes the multilevel procedure to fail to converge. This has indeed been observed in practice at high Reynolds numbers when other restriction operators are used, but it seems that the operator (5.3) overcomes this problem. Conversely, during the prolongation of the corrections from level 2h to level h, if the solution u2h, v2h contains oscillations then the corrections will also be oscillatory. These oscillations will appear to have a greater wavelength compared to the mesh spacing on level h than on level 2h. In fact, since at any level the equations are not solved exactly but only a few iterations are performed, these oscillations will survive until prolongation to the finest levels. On these finest levels the wavelength of the oscillations will appear to be very large compared to the mesh spacing, and therefore the smoother will be unable to reduce them. The impact on the overall efficiency of the multilevel procedure will be detrimental. One known solution to this problem (see [2]) is to define the local momentum fluxes using a blend of the CDS with the 1 st order upwind scheme (UDS) for the convection terms. That is, replace Fh,fuh,c by u , is the value of u at the centre of the face as given by the 1 st -order UDS scheme (it is equal to the value at the centre of the adjacent CV from which fluid flows towards the face). The global fluxes remain unaltered so that the composite grid equations (5.1) still use pure CDS. This change results in a smoother solution within the local part of each level, but because the local operator N2h departs from Nh, the corrections produced are less effective and the convergence rate of the multilevel procedure drops. This will be shown in Section 6. Alternatively, it has been found that it is more efficient to leave the discretisation of the local momentum fluxes as it is and to smooth the corrections prior to prolongation using a smoothing operator S2h, so instead of (5.13) prolongation takes place as follows: The corrections of v are similarly smoothed, while the corrections of p need not be smoothed. The smoothing operator was inspired by the fact that, as discussed in [25], in the absence of momentum interpolation an oscillatory pressure field results, but which results in near-correct values of pressure at face centres. Therefore, Sh is defined as:  (5.16) where fP is the set of 4 faces of CV P (in level-wise treatment every CV has exactly 4 faces) and the overbar denotes linear interpolation at point c' of face f from the values at adjacent CV centres. The efficiency of this technique will be demonstrated in Section 6.
Any smoother can be used in the context of the multilevel procedure as described so far. Now some particular issues which concern the use of SIMPLE will be addressed. First of all, a problem arises in the presence of outlet boundaries, when the outlet boundary condition is implemented in the usual way (e.g. as in the CAFFA code provided with [2] -see also [27]). In this implementation, after solution of the velocity linear systems within each SIMPLE outer iteration the outlet mass fluxes are scaled to satisfy overall mass conservation through the boundaries of the domain. This actually interferes with the discretisation of the equations and produces a final solution which does not exactly fulfil the zero-gradient condition at the outlet -see [27] for more details. The problem with the multilevel procedure is that since a particular level may not include all global outlet faces, it is difficult to find a suitable way to update the outlet mass fluxes of the particular level while simultaneously ensuring global mass conservation. To overcome this problem, it was chosen in this work not to alter the outlet mass fluxes during SIMPLE smoothing sweeps on any particular level, but to perform global SIMPLE sweeps on the whole composite grid between multilevel cycles. These will be termed composite-grid smoothing sweeps in contrast to the level smoothing sweeps which occur within each multilevel cycle.
Another complication arises at the exterior CVs of each level. Because the coarse CVs + ( Figure 5) do not contribute unknowns to the system of equations of level 2h, it so happens (see [27]) that the exterior faces of 2h may have negative contributions to the main diagonal of the matrix of the velocity linear systems of each SIMPLE sweep on level 2h. Thus for the solution of these systems we used unpreconditioned GMRES which does not require strict diagonal dominance. It was observed that usually 4-5 inner iterations suffice to achieve the full rate of convergence of outer iterations. Also, for the construction of the pressure correction system of SIMPLE, the mass flux correction through a face assumes a form which is a function of the main diagonal coefficients of the velocity linear systems of the CVs which lie on either side of the face. For exterior faces, during level sweeps, the coefficient of the fine CV is unsuitable as mentioned above, and the coefficient of the coarse CV does not exist. However, since composite-grid smoothing sweeps occur between multilevel cycles, we store the contributions of exterior faces to the matrix of coefficients of the pressure correction system during each composite-grid sweep and use them within the immediately following cycle. Since at each level sweep the exterior face mass fluxes are also corrected, the pressure correction system becomes diagonally dominant and always has a single solution (except if the particular level does not have exterior faces). We use ILU(0)-preconditioned conjugate gradients for the pressure correction system.
Finally, the matter of the type of multilevel cycle must be briefly addressed. For structured grids it has been found that for medium and high Reynolds numbers W cycles are usually more efficient than V cycles. However, for composite grids this may not be the case if the finest levels have very few CVs. In such cases it is often more efficient to use a cycle which resembles a W cycle at the lower levels and a V cycle at the higher levels, like the one shown in Figure 6. In the following, the notation W k /V(ν1,ν2)-s means a multilevel procedure where on levels up to the k-th 2 cycles are used to solve the problem of the immediately finer level, while 1 cycle is used on levels > k. Also ν1, ν2 are the number of pre-smoothing and post-smoothing iterations, and s is the number of composite-grid sweeps between multilevel cycles.

Square Lid-Driven Cavity
The square lid-driven cavity problem is probably the most often used case to test new methods for incompressible flows. It is used in many of the papers mentioned in Section 1 ([6], [7], [12], [13], [15], [16], [18], [19]). Usually methods are validated by comparing their results against those of [30], but these results are not very accurate. For the particular case of Re=1000 the very accurate results of [31] are available. In the following, the cavity has side L = 1 m while the lid velocity and fluid density are fixed at V = 1 m/s (in the positive xdirection) and ρ = 1 kg/m 3 . The viscosity μ is varied according to the required Reynolds number Re = ρVL/μ. The quantities Qmom and Qmas are defined as Qmom = ρLV 2 = 1 kgm/s 2 and Qmas = ρLV = 1 kg/s. For comparison, r , /ΔΩP is the residual per unit volume after cycle k, and m, m0 are large enough so that the rate of convergence has stabilised but q is not influenced by the first m0 cycles where the rate of convergence may be irregular. The maximum value over the momentum and continuity equations is displayed (usually the reduction factors are nearly equal for the three sets of equations). Unless otherwise indicated Full Multigrid (FMG) is used with parameters au=0.8, ap=0.2 (SIMPLE underrelaxation factors), amg=1, ami=1, and the coarsest level (CL) is level 1 (88 CVs). For low Re numbers V cycles are more efficient, while W cycles are more efficient for medium and high Re numbers. For Re  5000 velocity oscillations make it impossible to attain convergence unless either the velocity corrections are smoothed according to (5.15) or a blend of CDS/UDS is used for the local convective momentum fluxes. We have found that it is always more efficient to smooth the corrections, as demonstrated in Figure 7, and it is this technique which was used to obtain the results for Re=5000, 10000 in the Tables. Table 2 shows the same data for the new SIMPLE/multilevel method described in Section 5. For low and medium Re the results are similar, but the method becomes very inefficient at high Re numbers. In fact, it is not possible to attain convergence unless ami is dropped to 0.5 for Re=2000 (not shown) and 0.1 for Re=5000, resulting in great inefficiency. The larger the distance between the finest and the coarsest level, the smaller ami must be. For example, Table 2 shows that by letting the coarsest level CL=2 (1616 CVs) it is possible to increase ami to 0.2, increasing the efficiency. The problem seems associated with momentum interpolation, and an outline of a possible explanation is that it is caused by the fact that local mass fluxes (5.5) contain not one but two separate pressure terms. The fixed pressure term of p   ) may be relatively large, especially at coarse levels because it is calculated from the restricted pressure field, and it must be balanced by the "dynamic" pressure term of F2h,f(u2h,v2h,p2h) and by the mass flux change which is due to the velocity correction. This spoils the corrections to be prolonged back to the finer levels. On the other hand, in the classic SIMPLE/multigrid method the local mass fluxes contain a single pressure term like global mass fluxes, and its magnitude should in general be negligible compared to the total flux (see [25]).
then the sum of the two pressure terms would equal the single pressure term of the classical SIMPLE/multigrid method involving the pressure correction, and the two methods would become equivalent. However this would require that F2h,f( 2h u   , 2h v   , 2h p   ) be recalculated at every smoothing sweep (because the value of Af(u2h,v2h) would change) and to avoid this it was chosen instead to replace Af(u2h,v2h) by Af( 2h u   , 2h v   ) in the definition (3.8) of F2h,f(u2h,v2h,p2h) of (5.5). This again makes the sum of the two pressure terms equivalent to a single pressure term involving the pressure correction, but it is not equivalent to the pressure term of the classic method. The results shown in Table 3 were obtained with this method and it is clear that the situation has improved, with ami=1 used throughout. In the case Re=10000 amg=0.8 had to be used for the velocity, and efficiency is somewhat lower than that of the classic method, but still it is considered adequate. This modification will be used in all subsequent experiments. It should be noted here that the difference between Af(u2h,v2h) and becomes greater as Re increases as seen from (3.10) (because their convective parts are different but their viscous parts are the same). This may explain why the problem occurred at high Re. The table also shows the CPU time required by the FMG procedure on a 1.4 GHz Pentium 4 processor. The time required per cycle for the results of tables 1 and 2 is nearly the same as for table 3.
The next step is to allow local refinement. Starting from level 4 (6464) local refinement levels are added according to the criterion (4.2). At the upper corners of the cavity the flow field is discontinuous and this causes the truncation error to actually increase near the corners as the grid is refined. Therefore, the refinement criterion cannot be satisfied and the process of refinement must be ended by setting a maximum allowable number of grids. For Re = 100 and 5000, rτ = 0.01 was used and two additional grids were allowed (with level 6 being the finest). The underlying grids of the corresponding final grids are shown in Figure 8. Local refinement occurs at the regions where the estimate (4.1) predicts a high truncation error, that is near the top lid and corners for Re = 100 and at the circumference of the main vortex for Re = 5000 (see [27]). The number of cycles required to reduce the maximum residuals per unit volume below 10 -8 and the corresponding reduction factors are shown in Table 4, where grid 1 comes from grid 6464 with local refinement, and grid 2 (whose underlying grid is shown in Figure 8) comes from grid 1. For each grid the finest level is written in parentheses. The FMG procedure is used in each case starting from grid 88. For Re = 100 V(2,2)-1 cycles were used and comparison with Table 3 shows that the number of cycles and the reduction factors are the same as for the non-composite grids. For Re = 5000 W k /V(2,2)-1 cycles were used with different values of k as shown in the Table. In this case the highest possible value k=5 proves to be the most efficient, also in terms of CPU time. A slight degradation of the convergence rate is observed for grid 2 compared to grid 256256.
For Re = 1000 rτ = 0.001 was used and four additional grids were allowed (level 8 being the finest) to study the benefits of local refinement comparing with the accurate results of [31]. The underlying grids of the resulting series of 4 grids are shown in Figure 9. As expected, because the refinement criterion is dynamic, not only is a new level added on each subsequent grid but also existing levels are extended. The truncation error is highest at the lid and at the circumference of the main vortex, especially at its right and lower part. Figure 10 shows the u-velocity discretisation error | u h  | (h again denotes the whole grid, not a particular level) along the vertical centreline on various composite and non-composite grids, calculated at the points and from the uvalues given in [31], which are regarded as "exact". The distance between the distributions of the non-composite grids is in accordance with the 2 nd order accuracy of the method. Composite grid 1 and grid 128128 have about the same number of CVs and offer nearly identical accuracy. The accuracy offered by composite grid 2 is very close to that of grid 256256 although it has half as many CVs. Composite grid 3 has about the same number of CVs as grid 256256 but offers a clearly more accurate solution. The accuracy of composite grid 4 is comparable to that of the much larger (in terms of number of CVs) grid 512512, except near the centre of the cavity where however the error of the 512512 solution is already very small and so this difference is not so significant. The addition of a relatively large number of CVs to get grid 4 from grid 3 does not seem to reflect in an equivalent reduction in error. To interpret this one must keep in mind that the new CVs are mostly located near the lid corners (see Figure 9), where the increase of the grid density is not so efficient in improving the accuracy due to the singularities. Also, the introduction of many refinement levels increases the area of the domain which is near level interfaces, where there are large O(1) truncation errors as mentioned in Section 3. Already this area is comparable to the area occupied by the whole of the finest level of grid 4. Finally, according to Tables 3 and 4 the convergence rates of SIMPLE/multilevel are similar to those on the noncomposite grids but convergence is achieved in less cycles because successive composite grids differ less than successive non-composite grids (since the finest levels cover only a small proportion of the domain) and so the initial guess due to FMG has a smaller error on the composite grids. This time using cycles W k /V with an increased number k does not seem to pay off, since the finest levels cover less space, and it is more efficient to use k = 3 or 4.
Next the procedure was tested on non-uniform grids. Figure 11 (left) shows a nonuniform 3232 grid, where the ratio of the heights or the widths of consecutive CVs is r32 = 1.1561. Again for Re = 1000, using rτ = 0.05 and allowing 3 additional grids produces Grid 1 shown in Figure 11 (right). Then Grid 2 is produced by refining every CV of Grid 1, and Grid 3 is produced by refining every CV of Grid 2. The number of global CVs of each level of Grid 1 and the percentage of the domain that they occupy are displayed in Table 5. Also, for comparison a series of non-composite grids are constructed, up to 256256, which are such that r64 = r32, r128 = r64 etc. The grid lines of any given grid are also grid lines of the immediately finer grid. Figure 12 again displays | u h  | along the vertical centreline on various grids. Again the distance between the distributions of successive non-composite grids is in accordance with the 2 nd -order accuracy of the method. However, the distance between the distributions of the composite grids also suggests 2 nd -order accuracy. This is a very interesting observation in view of the fact that the truncation error has magnitude of O(1) near level interfaces, as noted in Section 3. Unfortunately in [31] only results along the centrelines are given, so it can not be strictly verified that convergence is 2 nd -order throughout the domain. Of course, 2 nd order accuracy with respect to refinement of the whole grid does not mean that the addition of any single level may not even cause an increase of the discretisation error due to the truncation error increase at the new level interface. Comparing the error distributions of Figure 12, there does not appear to be any clear benefit in using the composite rather than the non-composite grids. A possible explanation is that the non-composite grids are already nearly optimal. This can be seen also from the fact that the two finest levels cover very small percentages of the domain. These finest levels are located near the top corners where the flow field is singular which limits their contribution towards the increase of accuracy. Also, despite the fact that they cover only a very small percentage of the domain, they consist of a relatively large number of CVs due to the non-uniformity of the grid. Also, in consecutive noncomposite grids the maximum width to height ratio of the CVs increases (from 8.8 in 3232 to 10 in 256256) while this ratio remains constant at 8.8 in the composite grids because the child CVs inherit it from their parent during refinement. This means that the spacing of the non-composite grids near the walls is smaller than that of the composite grids. Figure 13 shows . Similarly, the estimate shows high truncation error in a region near the boundaries which is twice as wide as the actual one. Table 6 displays the number of cycles required to drop the maximum residual per unit volume below 10 -8 (FMG is used) and the residual reduction factors for various cases. The reduction factors are somewhat worse than for the uniform grids, and this is due to the aspect ratio of the CVs (it is a well known problem that multigrid smoothers loose their efficiency when the grid spacing in one direction is much smaller than in the other -see [5]). For the composite grids only the results for the most efficient type of cycle (in terms of k) are displayed, in terms of CPU time (level 3232 is level 3). It was observed that in general increasing the number of k causes the number of cycles and the reduction factor to drop, but since it also causes the cycle to become more expensive it does not pay off to increase k as far as possible.

Skew lid-driven cavities
Next the procedure is tested on a problem which is often used for testing on nonorthogonal grids, the flow in a cavity whose side walls are inclined at an angle of β=45 or β=30 to the horizontal level. This problem is proposed in [32]. All sides of the cavity have a length of L = 1 m and the flow and solution parameters are the same as in Section 6.1 unless otherwise stated. Starting from a uniform 6464 grid (level 4) two additional grids are allowed according to a criterion rτ = 0.001. Simulations are performed for Re = 100 and Re = 1000. The underlying grids of the resulting final grids are shown in Figure 14. Again for Re = 100 the truncation error is higher near the lid and top corners, while for Re = 1000 this also occurs at the circumference of the main vortex which is smaller in size than that of the square cavity and is located near the top right corner. Figure 15 shows the u discretisation error distributions along the centrelines of the cavities which are parallel to the side walls. To calculate the discretisation error the exact u values were estimated with Richardson extrapolation (see [2] or [22]), using the solutions of the 128128 and 256256 grids and assuming 2 nd order accuracy. This was done because it was observed that the results presented in [32] are not significantly more accurate than the 256256 solution. For Re = 100, the benefits of using local refinement are not evident, as the solution on the final composite grid is of comparable accuracy as on grid 128128 and the two grids also have a comparable number of CVs. On the other hand it seems advantageous to use composite grids for Re = 1000: in regions of high discretisation error the solutions on composite grids 1 and 2 have nearly the same accuracy as the solutions on grids 128128 and 256256 respectively, although they have significantly less CVs. It is not surprising that local refinement works better with higher Reynolds numbers: Convection tends to transport discretisation errors generated at regions of high truncation error to distant locations without alteration, while diffusion (viscous forces) tends to transport discretisation errors with an ever decreasing magnitude as the distance from the source increases. Therefore, it pays off more to use refinement to reduce locally high truncation errors in the presence of strong convection.
The convergence of the SIMPLE/multilevel procedure, again for a criterion of 10 -8 , is shown in Tables 7 and 8, for the non-composite and composite grids respectively. Level 1 (88) is used as the coarsest level in the FMG solution procedure. Each outer SIMPLE iteration includes a 2 nd pressure correction step to account for grid non-orthogonality as suggested in [2]. For β=45 convergence is similar as for the square cavity, but for β = 30 some additional difficulties arise. For Re = 100 V(2,2) cycles no longer converge and either W(2,2) or V (3,3) have to be used. Also on composite grids for β = 30 it pays to increase the number of composite-grid SIMPLE sweeps to 4. Except for the case β=30 Re=100, which does not exhibit typical multigrid convergence, in the other cases the convergence rates on the composite grids are about the same as on the non-composite grids. The gains of the local refinement procedure in terms of CPU time are not as high as in terms of number of CVs. This may be due to the compiler used. The difference is that our code allocates memory for the CVs of the non-composite grids all at once (in the form of arrays), while the CVs of local refinement levels are allocated one at a time to allow for unrefinement (a feature not used in the present work).

Backward facing step
Finally, a brief description of the solution of the backward facing step problem (see e.g. [22], [33], [34]) is provided to test the method on problems with outlet boundaries. The domain and level 1 are shown in Figure 16. Like in the aforementioned studies, the narrow channel before the step is not included as part of the computational domain, but a fully developed velocity profile is assumed at the height of the step (x=0). More specifically, the top and bottom boundaries are solid walls, the right boundary is the outlet, while the left boundary (step) is a solid wall from y = -0.5 to y = 0 and an inlet boundary from y = 0 to y = 0.5 with a parabolic velocity profile u = 24y(0.5-y)U, U = 1 m/s. The grid is uniform with square CVs up to the middle of the domain (x = 15) and afterwards it becomes stretched under a mild constant stretching factor. Simulations were performed for Reynolds numbers Re = ρUH/μ of 133 (low), 400 (medium) and 800 (high) (transition starts at Re  1150). Again ρ = 1 kg/m 3 and μ is varied according to the Reynolds number.
Tables 9 and 10 contain convergence data for non-composite and composite grids respectively. Again FMG is used and the algebraic convergence criterion is 10 -8 . Level 1 is the coarsest level of the multilevel procedure. The composite grids are constructed starting from the 44024 grid (level 3) and using Qmom = 1 kgm/s 2 , Qmas = 1 kg/s, rτ = 0.01, allowing two additional grids. This time it is useful to perform smoothing of corrections according to (5.15) for the whole range of Reynolds numbers, and the results of Tables 9 and 10 include such smoothing. It is observed that, in general, increasing the grid fineness causes the reduction factor to increase, contrary to the ideal multigrid properties. The situation improves significantly if the number of composite-grid smoothing sweeps is increased. This makes likely the following explanation: Since outlet mass fluxes and velocities are not updated during each multilevel cycle, this means that the composite grid sweeps are totally responsible for the reduction of the errors of these outlet mass fluxes and velocities. These errors also contain smooth components, which would normally be reduced on coarse levels in a multilevel procedure, but now they have to be reduced by the composite grid SIMPLE sweeps. Like most single-grid solvers, SIMPLE becomes less efficient as the grid density increases which reflects in increased overall reduction factors in Tables 9 and 10. Fortunately, it seems that a significant improvement of the convergence rate results with only a small increase in the number of composite-grid sweeps.

CONCLUSIONS
A multilevel algorithm has been proposed for locally refined grids, and tested using SIMPLE as the smoother. The tests show that the existence of local refinement levels does not adversely affect the convergence rates, which are similar to those of the classic SIMPLE/multigrid method. Contrary to the algorithm of [20], the present algorithm does not require the partitioning of the domain into blocks, and solution takes place on the whole composite grid simultaneously and not block-by-block. Therefore, one expects that the present algorithm is more efficient than that of [20]. However, the algorithm of [20] is more easily parallelizable. Of course, the two approaches can be combined by partitioning the domain into blocks each of which is itself a composite grid.
Although the main focus of this work has been on the multilevel solver, several issues related to the local refinement criterion itself have risen: Local refinement was driven by a criterion which tries to minimise the integral of the truncation error by refining the CVs which contribute the greatest to this integral. The tests have shown rather moderate gains from the use of this local refinement technique, especially at low Reynolds numbers. The method suggested in [23] appears to be more efficient, although a direct comparison cannot be made since different test cases are treated. A factor which may limit the usefulness of local refinement is the zero-order accuracy of usual discretisation schemes at interfaces between different levels. This should be further investigated: The present tests suggest that this does not affect the overall 2 nd order of accuracy with respect to refinement of the whole grid, but one suspects that the high truncation errors at level interfaces are indeed sources of additional discretisation error, thus limiting the gains from local refinement, especially at low Reynolds numbers because viscous terms contribute O(1) to the truncation error while convection terms contribute O(h). A possible way around this problem would be to construct more accurate discretisation schemes for the CVs which are close to the level interfaces. Such schemes would also be more expensive, but they would not result in significant overhead because the number of such CVs is very small compared to the total number of CVs of the domain (the dimensionality of the region covered by these CVs is one less than the dimensionality of the computational domain).
An important issue related to the efficiency of the local refinement technique is how to choose the local refinement parameter rτ in (4.2). Normally rτ should somehow be related to the required level of discretisation error, but this issue was not investigated in the present work. Alternatively, one may normalize the truncation error estimate by the main diagonal coefficients of the linearised discrete system as is done in [13], [24], to obtain a quantity which has the same dimensions as the discretisation error, and which is an indication of the contribution of each CV to the discretisation error. A difficulty with this approach is how to obtain an indication of the discretisation error of pressure since the continuity equation does not contain pressure terms.
The O(1) magnitude of the truncation error near level interfaces prevents the use of the local refinement criterion there, but in the present work this problem was overcome by simply not allowing local refinement near level interfaces. The same problem arises at boundary CVs, and there the problem is more serious because a suitable way must be found to identify regions where refinement is actually needed. Otherwise, the use of the local refinement criterion may result in the pile-up of local refinement levels near the boundaries, unnecessarily increasing the number of CVs and causing a large increase of the truncation error near the boundaries. A possible remedy is to use a refinement criterion which is based on the finite element residual (e.g. [23]) instead of the truncation error estimate. However, the truncation error estimate may also be used to obtain a more accurate solution (see [25]). An alternative remedy would be again to construct more accurate discretisation schemes for boundary CVs.  Table 1: Square lid-driven cavity, uniform grids: For each case, the left column shows the number of cycles of the classic SIMPLE/multigrid procedure ( [8]) required to reduce the maximum residual per unit volume below 10 -8 , and the right column (in italic) shows the corresponding reduction factor (6.1).
Full multigrid (FMG) is used in each case starting from the coarsest level 88.