Collocation and Least Residuals Method and Its Applications

The collocation and least residuals (CLR) method combines the methods of collocations (CM) and least residuals. Unlike the CM, in the CLR method an approximate solution of the problem is found from an overdetermined system of linear algebraic equations (SLAE). The solution of this system is sought under the requirement of minimizing a functional involving the residuals of all its equations. On the one hand, this added complication of the numerical algorithm expands the capabilities of the CM for solving boundary value problems with singularities. On the other hand, the CLR method inherits to a considerable extent some convenient features of the CM. In the present paper, the CLR capabilities are illustrated on benchmark problems for 2D and 3D Navier– Stokes equations, the modeling of the laser welding of metal plates of similar and different metals, problems investigating strength of loaded parts made of composite materials, boundary-value problems for hyperbolic equations.


Introduction
The purpose of the present communication is to demonstrate the capabilities of the method combining collocation and least residuals (CLR) at the numerical solution of boundary-value problems for partial differential equations (PDE) systems.In the collocation method (CM), the method of collocation with least squares (CLS), and in the CLR methods, the original boundary-value problem for a system of differential equations is projected onto a finite-dimensional linear functional space serving to the construction of approximate solutions.More often than not, the approximate solution is sought as a linear combination with indeterminate coefficients of the elements of a basis defined in this space.In the CM, the solution follows from the collocation equations -the requirements for the approximate satisfaction of the equations of the original problem at separate points of the computational regionfor finding the solution.The collocation equations are obtained by substituting the representation of the sought approximate solution into the differential equations.In doing so the nonlinear equations of the problem are linearized.As a result, the collocation equations are linear algebraic equations for the coefficients of the linear combination sought for.The substitution of the approximate representation of the solution into the linearized boundary conditions also yields linear equations with respect to its coefficients.This makes it possible to approximate the differential problem by a discretized problem the solution of which is found from the corresponding system of linear algebraic equations (SLAE).This provides approximants to the solution of the original differential problem.Usually, the approximants are constructed from polynomial bases defined on finite-dimensional spaces in view of the a e-mail: vshapeev@ngs.rusimplicity of their application.The approximate solution may be defined in terms of a polynomial of a sufficiently high degree which is uniquely defined over the entire region of the problem solution (the p-variant of the method).Alternatively, the solution region may be partitioned into subregions inside which one searches for "local" polynomials approximating the solution of the differential problem (the h-and hp-variants of the method).In this case, the global solution is got as a piecewise polynomial.The subregions may consist either of single cells or of unions of several cells of the discretization grid which covers the global solution region.
The theoretical prerequisites of the method were first developed in the work of Vainniko [1] who proved the possibility of the exponential convergence of polynomial approximants with the increasing degrees of the polynomials of the basis.
In any CM version, the number of equations written down in the SLAE coincides with the number of coefficients employed in the expression for the discrete problem solution.The CM has a simple implementation and gives good results at the solution of boundary-value problems for ordinary differential equations, both linear and nonlinear [2][3][4] as well as at the solution of the problems for PDEs without singularities.But at the solution of problems with singularities, resulting in large gradients for instance, their discretization by the collocation method yields an approximate solution characterized by a high error and ill-conditioned SLAE.
In the CLS and CLR methods, the coefficients of the sought approximate solution are got from an overdetermined SLAE, that is the number of equations exceeds the number of unknowns.Besides the collocation equations and the boundary conditions, matching conditions for the solutions found inside neighboring subregions, are to be written down at suitable chosen points on the common boundaries.The requirements of continuity of the values of the solution and/or its derivatives at the chosen points can be taken for such matching conditions.They can be implemented in different forms.The only condition is to be consistent with the definition of the differential problem solution.After the substitution of the tentative solution into them, they also yield equations for the coefficients of the solution sought for.
In the CLR method, the solution of the overdetermined SLAE is sought from the requirement of the minimization of the residual functional of its equations on the solution.The fulfillment of this requirement minimizes the maxima of the numerical solution error.The staightforward consequence of this fact is the better conditioning of the SLAE resulting within the CLR method as compared to that resulting within the CM method for the discretization of a same problem on a same grid.As a matter of fact, the error of the SLAE solution, its condition number, and the magnitude of the residual functional are interrelated by direct dependencies and correlate with one another.This circumstance is similar to that encountered in the problem of constructing an approximant of a one-dimensional function from a discrete sampling involving data affected by considerable errors.In this case, the error of the approximation spanned by an interpolating Lagrange polynomial is frequently higher than that based on the least squares method of Gauss-Legendre.In the latter case, the requirement of the minimization of the functional of the problem solution leads to the reduction of the maxima of its error.It is distributed more uniformly over the solution region.
The idea of combining the CM and the finite element method (FEM) with the least-squares method for the numerical solution of the mathematical physics problems was proposed independently and simultaneously in [5] and [6].Similar to the CM, in the CLR method one should distinguish among the versions indexed by the letters h, hp, and p.In the p-version, one can assume that the entire computational region represents a single grid cell.In doing so one tries to reach the needed accuracy of the numerical solution at the expense of a high approximation accuracy of the differential problem solution.To this end, one takes an increased dimension of the basis of the employed functional space onto which the differential problem is projected.The main merit of this version is the much smaller size of the computer storage as compared to the other versions.Its shortcoming lies in the fact that the SLAE matrix of solving the discrete problem is densely filled.In more flexible h-and p-versions, the computational region is covered by a grid which, similar to the difference methods, may have different shapes: it can be either regular or irregular, it may consist of rectangular, triangular, or other shape cells.These versions are well parallelized.The necessary accuracy increase in the h-version is obtained at the expense of a substantial refinement of the grid cells.A merit of this approach to the numerical solution stems from the occurrence of a sparse coefficient matrix of the corresponding SLAE.This enables the application of efficient algorithms for the SLAE solution.The main features of the hp-version stay inbetween those of the p-and h-versions.The necessary accuracy increase of the approximation can be obtained on fairly crude grids.An example of such a version, which makes use of sixth-degree polynomials, is discussed below in section 3.
A number of works of different authors have been recently published in which the minimization of the residual functional of the overdetermined systems of equations of the discrete problem is used in conventional versions of the FEM (see, e.g., [7]) and Galerkin techniques.Their discussion goes beyond the scope of the given communication.

P-versions of the CLR method
In view of the limitations for the paper size we illustrate the presented concepts on simple problems.More complex problems will be discussed briefly with some comments on published results.
Consider the Dirichlet problem for the Poisson equation for in the two-dimensional region Ω with boundary ∂Ω.For simplicity, in the numerical experiments mentioned in the following, Ω = [0, 1] × [0, 1] is a unit square with its left corner at the coordinate origin.In the p-version of the CLR method, we search for the solution in the space of polynomials P k (Ω) of degree k defined over Ω. Bases of power bivariate monomials, products of two orthogonal or non-orthogonal polynomials, one of which depends on the variable x 1 and the other on x 2 , were tested in numerical experiments.In particular, in the case of power monomials, the solution is sought in the form The coefficients b i+ j( j+1)/2 are determined from the collocation equations and the boundary conditions on ∂Ω.Several choices of the point sets for writing the collocation equations were taken in Ω and the boundary conditions on ∂Ω: uniformly distributed and at the zeros of Chebyshev and Legendre orthogonal polynomials.It is to be noted that, unlike the interpolation problem, the accuracy of the obtained results were only slightly different from each other in the above-mentioned cases.For problems with smooth test solutions, e.g., u(x 1 , x 2 ) = e x 1 + e x 2 , the solution accuracy grows rapidly with the increasing powers of the polynomials of the basis.Double precision solution accuracy in the C norm (the uniform norm) was reached at k = 10, 11, and 12.With increasing k, the number of coefficients in (2), the number of collocation points, and the number of equations for boundary conditions is increasing.In the implemented versions, the number of the equations in SLAE was by a factor of about 1.5 higher than the number of coefficients b i in (2).When the solution error is near the machine zero, the attainable accuracy depends substantially on the method of solving the SLAE.Under the use of orthogonal bases on a computer with double precision arithmetic, one can reach the accuracy of about 10 −15 , while under the use of the least-squares method, the occurence of significant roundoff errors prevents reaching an accuracy higher than 10 −13 .
Computer experiments were also conducted on test problems with large solution gradients taken from available literature sources.In some of them, the double precision was reached only at k near 50.At high k values, the rates of convergence of the problem solution and the accumulation of the rounding errors depend substantially on the specific expressions used for the spanning polynomials and the order of the arithmetic operations entering the computation of their values.
The capabilities of the different h-and hp-versions of the CLR method were also tested on the case study problem (1).For example, the approximation with an accuracy of the order of 10 −7 of the test solution u(x 1 , x 2 ) = e x 1 +x 2 by the h-version of the CLR method in the space of second-degree polynomials asks for a grid of 160 × 160 cells.The CPU time needed to reach this accuracy exceeds by far the CPU time asked by the p-version to reach the double precision accuracy for the computation of such a smooth solution.

Versions of the CLR method for solving the Navier-Stokes equations
The 2D and 3D problems of viscous flows in a cavity with the top cover moving with a constant speed provide benchmarks for the capabilities of different numerical methods in the solution of Stokes and Navier-Stokes equations.Various versions of the CLS and CLR methods have been tested against such benchmarks [9][10][11][12][13][14].For simplicity, the following discussion is confined to a square cavity Ω = Let the Navier-Stokes equations be satisfied in Ω (Figs. 1, 2), and the velocity v = g be set on its boundary ∂Ω: g = (1, 0) on the top side of the square, while g = (0, 0) is asked on the remaining sides to fulfil the smoothness requirement.
In Ω, a grid consisting of N rectangular cells Ω 1 , . . ., Ω N is defined.Let 2h 1i and 2h 2i be the sizes of the cell Ω i along the x 1 and x 2 directions, respectively, i = 1, . . ., N. The h-and hp-variants of the CLS method are applied in combination with a domain decomposition method.Specifically, the derivation of the solution of the problem implies the construction of a converging iterative process at each step of which auxiliary problems are sequentially solved in all the decomposition subdomains Ω 1 , . . ., Ω N entering the problem domain Ω.These problems are called local.The derivation of the solution of the problem starts with the choice of an initial guess of the solution in Ω.Let a subdomain be an individual cell Ω i .The problem in the cell Ω i is formulated for the Navier-Stokes equations linearized after Newton where j = 1, 2 and (x 1 , x 2 ) ∈ Ω i .Here, v k+1 1 , v k+1 2 , and p k+1 are sought quantities at the (k  in Ω i are obtained from the requirement of continuity of suitable combinations of the components of the velocity and the pressure on the sides separating Ω i from the neighboring cells.If the side ∂Ω i of the cell Ω i lies on the boundary of Ω, then the following boundary conditions of the original problem are satisfied on it: EPJ Web of Conferences Here , and τ denote respectively the outward normal and tangent unit vectors to the boundary ∂Ω i , ∂x 2 , j = 1, 2; vn , vτ , p are the approximations of the solution components in the neighboring cells of the reference cell Ω i within which the (k + 1)th iteration of the solution is constructed, h i = (h 1i h 2i ) 1/2 .The conditions ( 5) and ( 6) are analogous to those proposed in [9].The pressure is uniquely determined provided the following relationship is satisfied, It is convenient to define inside the cell Ω i the new variables where are the local coordinates, and (x 1ci , x 2ci ) is the cell center.In the versions of the method which are considered here, the velocity (u 1 , , where Ωi is the image of Ω i in the local coordinate system, P m v ( Ωi ) is the space consisting of all polynomials of degree not higher than m v defined in Ωi .As basis functions for the velocity (u 1 , u 2 ) we use the elements (ξ α 2 2 , 0), (0, They and any of their linear combinations identically satisfy the continuity equation div u = 0. The pressure q is sought in the space P m p ( Ωi ) spanned by the basis elements Denote the basis elements of the V m v ( Ωi ) space by ϕ v 1 , . . ., ϕ v d v and the basis elements of the P m p ( Ωi ) space by ϕ p 1 , . . ., ϕ p d p .An approximate solution in the cell Ωi is represented as a linear combination of the basis functions: C p ik ϕ p k , where C v i1 , . . ., C v id v and C p i1 , . . ., C p id p are undefined coefficients.They are found from the overdetermined system of linear algebraic equations (SLAE), which is obtained after the substitution of the desired approximate solution in equations ( 3), ( 5) - (7), and (8).The collocation points, the matching and the boundary condition points are chosen to be uniformly distributed inside each cell and on its boundary, respectively.Their number was chosen such that the number of the equations in the SLAE for determining the expansion coefficients outnumbered the number of the unknowns.
The SLAE was solved in the present work by the orthogonal method of linear algebra rather than the least squares method.If the elimination of the under diagonal elements of the matrix does not ask for permutations of the equations, then the orthogonal method yields the same manifold of solutions (pseudosolutions) as the method of the least squares provided the round off errors are absent.It is to be noticed that, in contradiction to the least squares method, the orthogonal method does not worsen the coefficient matrix conditioning.As a consequence, it is more stable with respect to the round off errors and allows the derivation of numerical solutions of SLAE of larger dimensions and to build up more accurate solutions of the initial problem.Due to the decomposition method applied, the computations can be parallelized.
The class of driven cavity flows at various Re is large and diverse.The solution has singularities in the upper corners of the cavity.There is a large number of published numerical results in the form of tables of characteristic values of the solution components which have been obtained for this problem Mathematical Modeling and Computational Physics 2015 01009-p.5 by various methods, see, e.g., [15][16][17][18].In these equations, the quantity 1/Re occurs as a coefficient of the highest partial derivatives with respect to the velocities.Due to the fact that the velocity component u 1 undergoes a discontinuity at the boundary ∂Ω the solution has singularities: in the neighborhood of corners A(0, 1) an D(1, 1) (Fig. 1), the pressure and the first derivatives of the velocity components are unbounded.To increase the computation accuracy, the elimination of the solution singularities in the upper corners was done following with minor differences [16].There is in the driven cavity flow a central vortex (denoted PE as in [15]) as well as chains of Moffatt vortices in the left (BL 1 , BL 2 , . . . ) and in the right (BR 1 , BR 2 , . . . ) lower corners.At high values of the Reynolds number, one or two eddies (T L 1 and T L 2 ) can appear in the left upper corner.They are indexed in order of decreasing size: the larger the index, the smaller the eddy.The results presented here were obtained on a sequence of three grids, M 1 , M 2 , and M 3 .The grid M 1 with the maximum and minimum cell side lengths equal to 1/16 and 1/256, respectively, is shown in Fig. 1 left.The grid M 2 is obtained from M 1 by dividing each of its cells into four equal parts.The grid M 3 is constructed from M 2 in a similar manner.The driven cavity flow was computed using the CLR version with m v = 6, m p = 5.To compare our results with previously published eddy characteristics, the stream function was computed using the formula In each cell, the velocity components of the solution were polynomials.The integrals of them were calculated analytically.Thus, no additional error was associated with their computation.Figure 1 right plots streamlines of the cavity flow at Re = 1000 for the approximate solution obtained on M 3 .The central eddy PE and the corner eddies BL1 and BR1 are apparent.Other eddies are poorly distinguishable in Fig. 1 because of their smallness.For illustrative purposes, the left part of Fig. 2 depicts an enlarged fragment of the flow pattern computed in the lower left corner of the cavity.The right part of Fig. 2 depicts an enlarged fragment of the flow pattern in the lower right corner of the cavity.Here, the vortices BL 2 , BL 3 , BR 2 , and BR 3 are well resolved.Table 1 provides the characteristics of the eddy PE: the coordinates of its center (x 1 , x 2 ) and the value of the stream function ψ at the eddy center.The results reported in [16][17][18] are also given for comparison.The comparison of the outputs got on the grids M 1 , M 2 , and M 3 points to a good convergence rate.Moreover, exhaustive comparison of various flow parameters produced by the CLR method with the results presented in [16][17][18] besides the sampling provided in Table 1 points to a good agreement of the outputs.For example, the stream function value at the center of PE as calculated on M 3 coincides with that presented in [18] to within 2•10 −8 .Here, we used the hp-version −0.118938 0.5300 0.5650 [18] −0.1189366 0.5307901 0.5652406 of the CLR method with 56 coefficients of the approximate solution in a single grid cell.This allowed the accurate description of small details of the flow.Thus, although the eddies BL 3 and BR 3 are small and weak, their relative error on M 3 is sufficiently small.The driven cavity flow in this paper and in [16][17][18] were computed using different high-order accurate methods.A comparison of the solutions produced on the sequence of grids suggests that they converge to the same solution.This result can be regarded as proof of the existence of a unique solution to the problem by numerical experiment.
In this study, we also implemented versions of the CLR method for solving 3D Navier-Stokes equations and conducted the flow simulation in a 3D cavity.The results were compared with various published results.The closest match was with those of [19].It was found out that the difference between the values of the velocity components obtained here and in [19] at Re=1000 did not exceed 10 −3 .The obtained accuracy of the numerical solution on a grid of moderate dimension points to the good performance of the present CLR method for the numerical simulation of the flow of the viscous fluids.

Applications
The problem of the 3D modeling of laser welding of metal plates consisting either of homogeneous or heterogeneous metals has been solved by a conservative version of the CLR method in different physical formulations [23,24].Cases have also been computed, which are important for technology, when there is a thin insert of a third metal between the plates of heterogeneous metals to be welded.Figure 3 shows a schematic longitudinal cut of the welding zone.The laser beam is directed perpendicularly to Mathematical Modeling and Computational Physics 2015 01009-p.7 the upper surface of the plates and moves along their joint at a constant velocity in a direction opposite to the x-axis.Under the assumption of the quasi-stationarity of the process, without considering the secondary factors, the computations of the temperature field, the shape of the welding joint, and other important technological parameters of welding, were done.The computations were based on a thermophysical model in which the welding was simulated by solving a 3D heat equation.The thermal interaction of the laser with the welding zone and of the remaining subregions with one another and with the ambient medium was modelled by different boundary conditions on the external and internal boundaries of the computational region, including the nonlinear ones, and on the curvilinear surface of the vapor channel.The involved physical processes are characterized by large temperature gradients, several bodies with different thermophysical properties are present, there is a complex interaction of the laser beam with the melt in the welding pool, a large peripheral zone responsible for the heat removal from the welding zone into the ambient medium.It is necessary to take into account in the process description many physical parameters affecting its quantitative characteristics.To deal with all these intricacies, nonuniform grids were defined.These resulted in an acceptable accuracy of the computation on a desktop computer.In the implemented versions, special attention was paid to the control of the heat conservation law in the computational region.An accuracy of better than 1% for the check of this law was reached.The melt dynamics in the welding pool was modelled by solving the Navier-Stokes equations.A satisfactory agreement of the results of numerical modeling of various characteristics of the process with the data of physical experiments was obtained.
Various problems of the theory of the strength of loaded plates of homogeneous materials [25] and composite materials [26] were modelled by using the h-and p-versions of the CLR method and different known models.In [26], the problem of the bend of rectangular layered plates consisting of an arbitrary number of orthotropic layers has been considered.The arbitrary orthotropy directions of each layer in the plate plane lead to the anisotropy of the overall design.The problem statements have been considered within the framework of the classical Kirchhoff-Love theory, the Timoshenko theory, and the broken line theory.The comparisons were done with the spatial elasticity theory.The p-version of the CLR method was applied in [26] in a unit square with the polynomial basis in the Lagrange form.Collocation points were chosen at the roots of Chebyshev polynomials along each coordinate line passing through the collocation nodes lying on the coordinate lines transverse to them.The points at writing the boundary conditions on ∂Ω i were located in a similar manner.The control of the accuracy of the solutions obtained by the CLR method was done preliminarily on the solution of known test problems computed with a high accuracy by other methods.To check the reliability of EPJ Web of Conferences 01009-p.8 the solution of the above specific problems, comparisons were done with the solutions obtained for some of them on sufficiently fine grids both by using the above theories and the model of the threedimensional elasticity theory.An exponential convergence was observed in the numerical experiments with increasing number of the collocation points.Thus, the comparison with the numerical results presented in [26], yielded already a relative error of the order of 10 −5 for the value of the maximum bend of plates at the number of collocation points equal to 16 × 16 and, respectively, 64 points for writing down the boundary conditions.Similarly, a high convergence rate and high accuracy were observed for the stress quantities in the solution of this and other problems.The merits of the pversions of the method have manifested themselves well.

On other capabilities of the CLR methods
It is appropriate to note here also some other capabilities of the CLR method which follow from its two independent properties: 1) the presence of a piecewise polynomial solution and 2) minimization of the residual functional of the equations of the discrete problem on their solution.Both of these properties give positive effects on the solution of boundary-value problems for PDEs.
The fact that the approximate solution of the global problem consists of local analytical pieces enables a relatively simple construction of the method versions on irregular and adaptive grids [5,8,10,11].
The CLR method also makes it possible to construct relatively simply the versions of the higher accuracy method in regions with curvilinear boundaries [8].Within the framework of methodic developments, several conservative and non-conservative versions of the CLR method were tested over regions with curved boundaries for the solution of two-dimensional heat equations, including the Poisson and convection-diffusion equations over a region containing a curved internal boundary at which the equation coefficients have a jump.In all cases, we have succeeded in constructing versions of the method which have the second order of convergence.The simplest approach among the tested ones was that in which the computational region was covered by a regular grid with rectangular cells.Irregular (boundary) cells arose only at the region boundary at the intersection of the regular grid cells with the boundary.Then for writing the collocation equations at such boundary cells, their external part was used, which was cut off by the boundary.The boundary conditions were written exactly on the boundary arc, which cut off the external cell part.The neighborhood of small-size cells with the large ones affects negatively the SLAE condition number of the discrete problem.A criterion for the smallness of the irregular cell was introduced according to which a small irregular cell was attached to the neighboring cell with which they had a common side of the largest length.It was assumed further that the local solution of this neighboring cell is continued into the attached part.
The use of a regular grid inside the region greatly simplifies the algorithm and the computer code in comparison with the case of the use of irregular grids.In the first case, at the determination of neighboring cells there is no need to use the graph algorithms which complicate the code.
Figure 4 shows a fragment of a region, in which the Dirichlet problem for the Poisson equation with a smooth solution was solved, with the indication of the collocation points, the matching points, and the points for writing the boundary conditions.The approximate solution constructed by the CLR method with a basis of second-degree polynomials converged with the second order on the given sequence of grids under the use of the parts of boundary cells which were outside the region.An error of 1.630 • 10 −8 on a test solution was reached on the grid of 512 × 512 cells.The application of the Krylov's algorithm has reduced the CPU time by a factor of five.A more complex approach which does not use the parts of the boundary cells lying outside the region boundary was reported in [8].In this case, all the collocation points were placed approximately uniformly in the part of the irregular Mathematical Modeling and Computational Physics 2015 01009-p.9cell which belonged to the computational region.Figure 5 left shows a region in which the Dirichlet problem for a two-dimensional convection-diffusion equation with a small multiplicative parameter at the higher derivative (with the value of 10 −4 ) was solved.The equation coefficients of the first derivatives were chosen in such a way that a narrow boundary layer in the equation solution is located inside the region and is curved at a right angle, see figure 5.This has been described in more detail in [8].The definition of a grid adapted to large solution gradients and to the curved boundary resulted in highly accurate solutions.Subsequently it was found that the order of accuracy of the solution equates the order of accuracy of solution of the same problem over a square region under the use of a difference scheme of the fourth-order approximation on a grid with the number of cells approximately equating their number in the first case.In Fig. 5 (right) the plot of a fragment of the surface of the case study solution is shown.

01009-p.10
In the implemented variants of the method described above, the alternating Schwarz algorithm was used together with global Gauss-Seidel iterations.At each global iteration, all the cells of the discretization have been successively visited.The SLAE defining the solution over the cell discretization has been solved by the orthogonal method.
The implementation of some variants of the CLR method for solving, in particular, the 2D and 3D Navier-Stokes equations, used the efficient convergence acceleration of the iterative process based on the multigrid Fedorenko approach [21].At the first stage, the problem was solved on a coarse grid which consisted of one, nine or one hundred cells, until the convergence to a given accuracy was reached.After that, each cell was partitioned into four equal parts.The availability of the approximate solution in analytic form within a reference cell makes it possible to rewrite the solution formula in local variables inside each of the four finer cells by an accurate substitution of variables.In doing so, there occurs no loss of the solution accuracy reached on the foregoing grid.This positive property of the CLR method distinguishes it from some other methods which are asking for approximate numerical operations of interpolation or averaging under the use of multigrid complexes.Within such methods, accuracy is lost at the intermediate stages of the passage in the complex from one grid to another, while the convergence rate of the iterations decreases.
For different versions of the CLR method for solving the 2D and 3D Navier-Stokes equations, the acceleration algorithm of the iteration processes based on Krylov subspaces [10,11,14,20] was also used with the aim at reducing the CPU time.Its application on a grid of the 256 × 256 size in problems with the Reynolds number Re = 1000 sometimes resulted in an acceleration of the CPU time by a factor of up to 25. Numerous experiments were devoted to the application of the Krylov method on each grid of the multigrid complex [21] with the use of only the prolongation operation on it.In some versions, at Re = 1000, the reduction of the CPU time up to factors exceeding 300 was reached on grids of size 160 × 160 as compared to the CPU time asked for the straighforward derivation of the numerical solution over the finest grid of the complex under the use of the same criteria for the termination of the iterations.Such an effect may be explained by the fact that under the combination of these two algorithms for the acceleration of the iteration processes the minimization of the residual functional in the CLR method acts on a substantially improved process and introduces a significant contribution to the convergence acceleration of the solution.
We finally note one more positive property of the CLR method which has manifested itself at the solution of boundary-value problems for hyperbolic equations.Nonstationary problems for the nonlinear one-dimensional Burgers and Korteweg-de Vries equations were solved in [22].Fairly simple versions of the method were used with the bases defined in the space of third-degree polynomials and without adding any artificial terms to these equations.No manipulations with the solution of the type of the averaging over the cell and other similar operations were used.It was noticed that in the computations by the CLR method of the generalized solutions there were no oscillations of the solution, and the discontinuity was smeared over no more than 1.5 grid cells within grids consisting of up to 200 cells over the spatial variable.
The above-mentioned properties of the CLR method point to its capability to provide accurate modeling of large classes of physical processes.The development of dedicated program packages is feasible and highly desirable.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences, 201 and p k are the solution components obtained at the kth iteration step; ( a, b) is the scalar product of the vectors a and b; and Re is the Reynolds number.The boundary conditions in the local problem for v k+1 1 and v k+1

Figure 1 .
Figure 1.The adaptive grid and the picture of the streamlines in the cavity

Figure 2 .
Figure 2. Enlarged pieces of the flow pattern in the left lower and right corners of the cavity