2D Newton Schwarz Legendre Collocation Method for a Convection Problem

: In this work, an alternate Schwarz domain decomposition method is proposed to solve a Rayleigh–Bénard problem. The problem is modeled with the incompressible Navier–Stokes equations coupled with a heat equation in a rectangular domain. The Boussinesq approximation is considered. The nonlinearity is solved with Newton’s method. Each iteration of Newton’s method is discretized with an alternating Schwarz scheme, and each Schwarz problem is solved with a Legendre collocation method. The original domain is divided into several subdomains in both directions of the plane. Legendre collocation meshes are coarse, so the problem in each subdomain is well conditioned, and the size of the total mesh can grow by increasing the number of subdomains. In this way, the ill conditioning of Legendre collocation is overcome. The present work achieves an efﬁcient alternating Schwarz algorithm such that the number of subdomains can be increased indeﬁnitely in both directions of the plane. The method has been validated with a benchmark with numerical solutions obtained with other methods and with real experiments. Thanks to this domain decomposition method, the aspect ratio and Rayleigh number can be increased considerably by adding subdomains. Rayleigh values near to the turbulent regime can be reached. Namely, the great advantage of this method is that we obtain solutions close to turbulence, or in domains with large aspect ratios, by solving systems of linear equations with well-conditioned matrices of maximum size one thousand. This is an advantage over other methods that require solving systems with huge matrices of the order of several million, usually with conditioning problems. The computational cost is comparable to other methods, and the code is parallelizable.


Introduction
Different scales, complex domains, and boundary layers are computational challenges for laminar and turbulent flow models [1][2][3][4][5][6][7][8][9]. One of these flows is Rayleigh-Bénard convection, in which a layer of fluid is heated from below. The first to systematically study this phenomenon was H. Bénard [10], who observed a pattern of hexagons when heating exceeded a threshold. Lord Rayleigh [11] explained this phenomenon by identifying the destabilizing mechanism as buoyancy. However, this explanation does not correspond to what happens in Bénard's experiments. Pearson [12] discovered the second destabilizing mechanism, which is changes in surface tension. In this case, the problems are called Benard-Marangoni. This mechanism is relevant in crystal growth processes [13,14] or ferrofluids in a magnetic field [15]. If the heating is increased further, different states occur: stationary, periodic, or chaotic states, up to turbulence. Some of these phenomena have been studied in REfs. [16][17][18][19][20][21], and the monograph [22] provides detailed information about the understanding of this problem. If the mechanism is buoyancy, the problem is called Rayleigh-Bénard, which is the one we address in this work and is also addressed in Refs. [23][24][25][26][27]. The model equations are the incompressible Navier-Stokes equation and the heat equation under the Boussinesq approximation in a bounded domain [28]. Three parameters appear in the model governing the convective motion: the Rayleigh number R (ratio between the buoyancy and viscous force), the Prandtl number Pr (ratio between the kinematic viscosity and thermal diffusivity), and the aspect ratio Γ (ratio between the length and the depth of the fluid layer). In geophysical contexts studying the interior of the stars and planets, Pr is large and the convective flow is dominated by plumes [29,30]. This is context of this work.
In the case of order one aspect ratio, this problem has been numerically solved with different numerical techniques, mainly finite elements, spectral elements, and collocation methods [26,30]. The collocation methods are high-order methods that allow one to find precise numerical solutions for this problem with coarse meshes [31][32][33][34]. The resulting algebraic systems for these methods are ill conditioned when the mesh size [32,34] is increased. The condition number of the matrices resulting from the discretization is large. The maximum size of the mesh is restricted, and it is not possible to obtain solutions with many oscillations that appear in domains with a large aspect ratio or solutions with different scales that arise when the Rayleigh number increases until reaching turbulent states. The same problem can be found in other numerical techniques, i.e., finite differences, finite elements, or spectral elements for fine meshes. The linear systems that appear in the procedure are huge and usually poorly conditioned. One way to avoid this drawback is to use a domain decomposition strategy; examples of this procedure applied to some partial differential equations can be found in Refs. [2,[35][36][37][38][39][40]. The domain is divided into several subdomains, and the problem is numerically solved in each subdomain with a coarse mesh. In this way, ill conditioning is avoided because the numerical method is used on a coarse mesh. Domain decomposition is widely used for fluid dynamics problems. Alternating Schwarz is the most common domain decomposition technique used to solve incompressible Navier-Stokes. The use of this technique for incompressible Navier-Stokes with different numerical methods can be found in Refs. [41][42][43][44][45][46][47][48][49][50]. Schwarz domain decomposition methods have been broadly used as a parallelization strategy of fluid flows [4,5,[51][52][53]. The model equations are nonlinear. A way to deal with the nonlinearity is the use of a Newton method. Newton's method is a powerful nonlinear solver used successfully in many engineering communities. Recent research has introduced improvements and variants on Newton's method. In Ref. [54], inexact Newton-Krylov and quasi-Newton algorithms for nonlinear elasticity equations are presented. In Ref. [55], a spatial additive Schwarz preconditioned exact Newton method as nonlinear preconditioner for Newton's method is applied to multiphase flows. In Ref. [56], the advanced Jacobian-Free Newton-Krylov method in a steam generator is studied. Refs. [57,58] modify the Newton increment to ensure global convergence. Newton's method considered in this work is the standard method that also provides good results for this type of problems (see Refs. [59][60][61][62][63]).
A partial alternating Schwarz method for Navier-Stokes-Boussineq with aspect ratio one, low Prandtl number, and 4 × 4 subdomains solved with the method of approximate particular solutions can be found in Ref. [64]. In Ref. [30], this problem is solved with a spectral element method for the 3D and the 2D problem for different values of the parameters, showing that the 3D Rayleigh-Bénard problem is quasi two-dimensional for an infinite Prandtl number.
In this work we, achieve an efficient algorithm for the alternating Schwarz method to solve each step of the Newton scheme used to deal with the nonlinearity solved with a Legendre collocation method of a 2D Rayleigh-Bénard problem with infinite Prandtl number [63]. In Ref. [3], the convergence of this methodology was theoretically proven in an infinite domain when an overlap is considered for a linear approximation to the problem, and some numerical tests were included. The present work achieves an efficient algorithm to solve the fully nonlinear problem with an alternating Schwarz-Newton-Legendre collocation method. The number of subdomains can be indefinitely increased in both directions of the plane. The numerical solutions have been validated by comparing them with solutions obtained with other methods and with physical experiments [29]. Thanks to this domain decomposition, the aspect ratio and the Rayleigh number can be increased considerably by adding domains in any direction. Values of the Rayleigh number for near turbulent regimes can be reached and domains with large aspect ratio can be considered. The great advantage of this method is that we get solutions close to turbulence, or in domains with large aspect ratios, by solving systems of linear equations with wellconditioned matrices of maximum size one thousand. This is the main advantage over other methods that require solving systems with huge matrices of the order of several million, usually badly conditioned that require the use of preconditioners. The method is parallelizable, and, even without parallelizing, the computational cost is affordable.
The article is organized in the following manner. Section 2 provides the mathematical formulation of the Rayleigh-Bénard problem. The numerical method that includes the domain decomposition method is explained in Section 3. The numerical results and a discussion appear in Section 4. Finally, in Section 5, the conclusions are presented.

Formulation of the Problem
A two-dimensional (2D) fluid layer of depth d (z coordinate) lies between two parallel plates of length L (x coordinate). The upper plate is at temperature T 0 and the bottom plate is at T 1 such that T 1 > T 0 . A usual way to non-dimensionalize in Rayleigh-Bénard problems is the following: and θ = (θ − T 0 )/∆T, where the primes correspond to the quantities with their dimensions, x , z space coordinates, t time, p pressure, u velocity, and θ temperature, κ is the thermal diffusivity, ρ 0 is the mean density, ν is the viscosity, and ∆T = T 1 − T 0 . The Rayleigh and Prandtl numbers are defined as R = d 3 δg∆T/(νκ) and Pr = ν/κ, where g is the gravity constant and δ is the thermal expansion coefficient. The domain becomes The dimensionless forms of the momentum, mass, and energy balance equations that model the problem are [63]: where e z is the upwards vertical unit vector. The boundary conditions come from mantle convection related models [65,66]. They correspond to a rigid bottom wall, and free-slip, non-deformable surfaces at the top and lateral walls, the environment temperature is fixed at the upper plate, and the lateral walls are insulated, and at the bottom, a constant temperature is imposed, The boundary conditions for pressure are obtained using the normal component of the momentum equations on x = 0, Γ and z = 1, and the continuity equation at z = 0 (see [67,68]). The positivity of the inf-sup constant is guarantueed because these conditions eliminate the parasite modes for pressure [31].

Numerical Method
Problem (6)-(9) is nonlinear. The nonlinearity is solved with a Newton iteration [67]. Each step of the Newton method is a linear problem in which a Schwarz alternating domain decomposition method is introduced as explained below, and each Schwarz iteration at each subdomain is discretized with a Legendre collocation method [70].
This is the linear problem to be solved at each step of the Newton method. The numerical solution fields are then u n = u 0 + ∑ n i=1ũ i , p n = p 0 + ∑ n i=1p i , θ n = θ 0 + ∑ n i=1θ i .

Alternating Schwarz Domain Decomposition
In this section, we introduce an alternating Schwarz domain decomposition method to solve a step of the Newton method (10)- (13). The computational domain Ω is split into np × mp subdomains of the same size, where mp is the number of vertical partitions and np the number of horizontal partitions; Ω is a rectangle, and each subdomain is a rectangle of the same size. In Ref. [3], it is proven that two adjacent domains have to share an overlapping area for the algorithm to converge. Partitions marking the position of each subdomain are performed with double points to handle the overlap. A partition with 2np points on the horizontal interval of the domain, [0, Γ], and another with 2mp points on the vertical interval of the domain, [0, 1], are calculated. These partitions are denoted as . . , γ z 2(mp−1) , γ z 2mp = 1}, respectively. Note that neither γ x 0 , nor γ z 0 , nor γ x 2np−1 , nor γ z 2mp−1 are considered due to future calculations. These values depend not only on mp and np, but also on the Legendre-Gauss-Lobatto (LGL) collocation points. Therefore, they will be calculated in the next section. Then, each subdomain Ω q,p is defined as Ω q,p = (γ x 2p−3 , γ x 2p ) × (γ z 2q−3 , γ z 2q ), p = 1, . . . , np, q = 1, . . . , mp. Each subdomain has length Γ x and height Γ z . An example of these partitions is shown in Figure 1 with mp = 2 and np = 2. Now, the problem has to be adapted to each subdomain Ω q,p . The equations inside each subdomain are kept the same. The boundary conditions at the edges of the subdomain, which are also boundaries of the original domain, (∂Ω q,p ∩ ∂Ω), remain the same. The interfaces are the boundaries of the subdomains, which are inside the original domain (∂Ω q,p ∩ Ω). The new interface conditions are given by the continuity of the fields u n , p n , and θ n . The previous values in the exact same point of the adjacent rectangle are considered. The following system of equations is solved at each subdomain and each step of the iterative Schwarz method: −∆u s q,p + ∇p s q,p − Rθ s q,p e z = −∇p n−1 q,p + ∆u n−1 q,p + Rθ n−1 q,p e z in Ω q,p , ∇ · u s q,p = −∇ · u n−1 q,p in Ω q,p , −∆θ s q,p + u s q,p · ∇θ n−1 q,p + u n−1 q,p · ∇θ s q,p = ∆θ n−1 q,p − u n−1 q,p · ∇θ n−1 q,p in Ω q,p , B u s q,p , θ s q,p , p s q,p = −B u n−1 q,p , θ n−1 q,p , p n−1 q,p IB u s q,p , θ s q,p , q s q,p = h on ∂Ω q,p ∩ Ω.
The interface boundary conditions IB and h on ∂Ω q,p ∩ Ω are: where q = 1, . . . , mp, p = 1, . . . , np, (u s q,p , p s q,p , θ s q,p ) are the unknown fields at each subdomain Ω q,p , (q , p ) are the subscripts of the adjacent subdomains, which contain the interface (p = p − 1, p or p + 1 and q = q − 1, q or q + 1), s denotes the iteration index in the Schwarz algorithm, and n means the Newton iteration. Note that the continuity of the Newton iterations is considered at the interface, i.e., u n q,p = u n q ,p , where u n q,p = u n−1 q,p + u s q,p and u n q ,p = u n−1 q ,p + u s−1 q ,p , and similarly for the rest of the fields. These continuity conditions improve iteration errors.

Legendre Collocation
At each subdomain, the problem is numerically solved with a Legendre collocation method. The Legendre polynomials are defined on the interval [−1, 1], therefore, for computational convenience, each subdomain Ω q,p is transformed into (−1, 1) × (−1, 1) through the following change of variables . These changes of coordinates introduce the corresponding scaling factors in the equations and boundary conditions as follows: Then, the unknown fields in each subdomain Ω q,p are expanded in a truncated series of Legendre polynomials The corresponding expansions for the four different fields at the different subdomains are introduced into the equations, the boundary, and the interface conditions, and they are evaluated at the Legendre-Gauss-Lobatto collocation points (x l p , z k q ) with l = 1, 2, . . . , N and k = 1, 2, . . . , M. The Schwarz method requires an overlap between adjacent subdomains. The values of N, M, and the number of overlap points are known a priori. With all these numbers we can calculate the partitions for the Schwarz method of the (14)- (17) is defined for each subdomain independently. The interface condition requires evaluating the fields at the exact same points in two different subdomains.
The symmetry of the LGL collocation points is very useful for this purpose. If we couple two sets of N LGL collocation points in such a way that the first point of one of them matches the (N − overlap)th point of the other, it will occur that the (1 + overlap)th point of the first set matches the last point of the other set, as shown in Figure 2. This process allows building subdomains with overlap shared points. With this information, the partition values can be calculated as follows: as it should be. The other partition values can be calculated recursively knowing these four values, as γ The solution strategy is mixed with the approximation method.

Algorithm
A pseudocode for the algorithm can be summarized as follows: 1.
Define input parameters: R; Γ; number of subdomains mp × np; size of the mesh at each subdomain N (x direction), M (z direction); number of points of overlap, o;
Introduce the initial solution in the whole domain: u x0 , u z0 , t 0 , p 0 ; 4.
Calculate the part of the solution at each subdomain; 5.
Introduce the corresponding Legendre expansions (22) into equations and boundary conditions (14)- (21), and evaluate at the Legendre-Gauss-Lobatto points, the result is a system of linear algebraic equations; B.
Solve the linear system of equations to get the Newton perturbations at each subdomain (q, p): U x {q, p}, U z {q, p}, T{q, p}, P{q, p}; iii. End; iv.
Calculate the Schwarz errors, errorS: (e) Update the solution at each subdomain (q, p): Calculate the Newton error:
Save the results.

Numerical Results
Different values of the overlap can be considered with similar results, therefore overlap is fixed to 4 points. Different mesh sizes, N and M, between 10 and 18, are considered. Different values of the maximum number of iterations can be taken for the Newton and Schwarz methods. A balance between computational cost and accuracy is obtained by considering a maximum of 40 iterations for Newton's method and a maximum of 5 iterations for Schwarz's method. The iterative initial value for Newton method is a numerical solution calculated with Legendre collocation with expansions 18 × 14 in a domain with aspect ratio Γ = 3.495 and R = 1300 adapted to the new multi-domain.

Validation
Legendre collocation is an ill conditioned method, and the maximum mesh size that can be used in a double precision calculation for parameter values in a laminar regime with a small aspect ratio is N × M = 50 × 50. Some of these values of the parameters are R = 1300 and Γ = 3.495. For this case, we have compared the numerical solutions obtained with Legendre collocation in a single domain with the maximum size of the expansion, N × M = 50 × 50, with numerical solutions obtained with the Newton alternating Schwarz-Legendre collocation method with different numbers of subdomains and different sizes of the Legendre expansions. Table 1 shows the relative L 2 norm of the difference between the horizontal component of the velocity solution, u x , obtained with Legendre collocation with a mesh 50 × 50 in a single domain, and the solution obtained with the Newton-Schwarz method with different numbers of subdomains on the horizontal axis, np, a single domain on the vertical axis, mp = 1, and different Legendre collocation meshes, n. The optimal result, taking into account a balance between the coarsest mesh and the smallest error, is obtained with two subdomains on the horizontal axis, one subdomain on the vertical axis, and a Legendre mesh of size 10 × 10. In this case, the total mesh contains only 200 nodes. Increasing the number of subdomains and/or the size of the mesh at each subdomain makes the solution worse. The same can be seen in Table 2 for two subdomains on the vertical axis, i.e., the relative L 2 norm of the difference between the horizontal component of the velocity solution, u x , obtained with Legendre collocation with a mesh 50 × 50 in a single domain, and the solution obtained with the Newton-Schwarz method with different numbers of subdomains on the horizontal axis, np, two domains on the vertical axis, mp = 2, and different Legendre collocation meshes n. Cases with a single domain on the horizontal axis and two subdomains on the vertical do not converge. All other errors are worse than for a single subdomain on the vertical axis. Table 1 shows the optimal result taking into account the size of the total mesh, i.e., two subdomains on the horizontal axis, one subdomain on the vertical axis, and a Legendre mesh of size 10 × 10. Numerical solutions for larger values of the aspect ratio or larger values of the Rayleigh number cannot be obtained with a single domain, even with the finest mesh. For those cases, the results of the domain decomposition method have been compared with the numerical solutions obtained by means of a finite elements method with a fine mesh with COMSOL. The solutions obtained with both methods can be seen in Figures 3-10 for different values of the parameters. Figure 3 shows solutions obtained with finite elements for R = 1300 and Γ = 3.495. In Figure 4, solutions for the same values of the parameters obtained with alternating Schwarz domain decomposition are presented. This solution has been calculated with two subdomains on the x axis, one subdomain on the z axis, and a collocation mesh of size N × M = 10 × 10. Three rolls are observed in both cases. The difference of the maximum value of the temperature field between both solutions is 0.035. Figure 5 shows solutions obtained with finite elements method for R = 1300 and Γ = 27.96. In Figure 6, solutions for the same values of the parameters obtained with the Schwarz can be seen. This solution has been calculated with 14 subdomains on the horizontal direction, two subdomains on the vertical direction and a collocation mesh of size N × M = 16 × 10. Twenty-five rolls are observed in both solutions. The difference of the maximum value for the temperature field between both solutions is 0.3.    In Figure 7, solutions obtained with finite elements for R = 5000 and Γ = 3.495 are shown. In Figure 8, solutions for the same values of the parameters obtained with the Schwarz method can be seen. This solution has been calculated with eight subdomains on the x axis, two subdomains on the z axis, and a collocation mesh of size N × M = 16 × 10. Three rolls are obtained, and the isotherms are mushroom-shaped in both cases. The difference of the maximum value for the temperature field between both solutions is 0.7.
In Figure 9, solutions obtained with finite elements for R = 5000 and Γ = 27.96 can be seen. Figure 10 shows solutions for the same values of the parameters obtained with the Schwarz method. This solution has been calculated with 13 subdomains on the x direction, three subdomains on the z direction, and a collocation mesh of size N × M = 18 × 12. Twenty-nine rolls are obtained in this case, and the isotherms are mushroom-shaped. The difference of the maximum value for the temperature field between both solutions is 0.6.    The same solutions are obtained with the two different numerical methods. There is an optimal value of the number of subdomains that coincides with the minimum number of subdomains.

Numerical Convergence
The error in the Newton iterations for each field is defined as the Newton perturbations, similarly to the Newton error in Section 3.4, for the first domain: where F refers to any of the fields. The error in the Schwarz iterations for each field is defined similarly to the Schwarz error in Section 3.4 for the first domain and adjacent ones: We determine the convergence rate of the Schwarz alternating method as the limit of the ratio

Legendre Collocation
In Tables 3 and 4, we have varied the mesh sizes for two different values of the parameters and a number of subdomains, and we observe that the errors decrease to a minimum and increase as the mesh size increases. The optimum in the case of R = 5000 and Γ = 3.495 with eight subdomains on the horizontal axis and 2 = two subdomains on the vertical axis, is N = 14 and M = 12, which is of the same order as the one used in the figures N × M = 16×10 with the same number of subdomains. The optimum in the case of R = 1300 and Γ = 27.96 with 14 subdomains on the horizontal axis and two subdomains on the vertical axis, is N = 14 and M = 10, which is of the same order as the one used in the figures N × M = 16×10 for the same number of subdomains. The optimal mesh size of the expansions is a coarse mesh of order 10 × 10 in each subdomain. Legendre collocation is a high-order method on a single domain, but the global multidomain method with Schwarz and Newton iterations loses this property. Table 3. Newton errors of the last iteration for u x obtained with Legendre collocation with different meshes (N denotes the number of horizontal collocation points and M denotes the number of vertical collocation points); Ra = 5000, Γ = 3.495, np = 8, and mp = 2.

Schwarz
In this section, we analyze the behavior of the Schwarz parts of the algorithm, the errors, convergence rate, number of subdomains, and influence of the parameters.
First, we look at the minimum number of subdomains necessary to obtain convergence depending on the parameters present in the problem. Two examples are shown in Tables 5 and 6. In Table 6, the last error in the Newton iteration for the field u x , for R = 1300 and Γ = 27.96, increases in the number of subdomains in the vertical and horizontal directions can be seen. A minimum number of subdomains is required to achieve convergence. The minimum number of subdomains to get convergence is mp = 2 and np = 3, or mp = 3 and np = 3. A single domain in the vertical direction or one or two subdomains in the horizontal direction never converges. In the case of R = 5000 and Γ = 3, 495, in Table 5, the minimum number of subdomains corresponds to mp = 1 and np = 8, or mp = 2 and np = 6, or mp = 3 and np = 2. The number of Schwarz subdomains required increases with R and Γ, i.e., one domain is enough for R = 10 3 , and around one hundred subdomains are required for R = 10 5 or Γ = 100 with the same size of the mesh at each subdomain.  Second, we study the influence of the parameters in the Schwarz errors. First, we vary the Rayleigh number. In Figure 11a, the decimal logarithm of the errors in the velocity, temperature, and pressure fields for the Schwarz iterations in the case R = 1300, Γ = 3.495, np = 2, mp = 1, and N × M = 10 × 10, can be seen. In that figure, we can check the convergence of Schwarz's method for each Newton problem. The errors of the Schwarz solver reach values between 10 −8 and 10 −9 for all fields as Newton's problems progress. Figure 11b shows the decimal logarithm of the errors in the four fields for the Schwarz iterations in the case R = 5000, Γ = 3.495, np = 8, mp = 2, and N × M = 16 × 10. The errors for the Schwarz iterations reach values between 10 −5 and 10 −6 for all the fields. In that figure, the convergence of Schwarz for each Newton problem can be seen. We have increased further the Rayleigh number for a fixed value of the aspect ratio, Γ = 3.495. Figure 12 shows a solution for R = 10 4 obtained with 25 subdomains on the horizontal axis and two subdomains on the vertical axis for a collocation mesh of 12 × 10. The three rolls with the mushroom shaped isotherms still remain. In Figure 13, a solution for R = 10 5 obtained with 25 subdomains on the horizontal axis and five subdomains on the vertical axis for a collocation mesh of 16 × 12 can be seen. Eleven rolls with plume shaped isotherms is a new solution in the near turbulent regime. Figure 14a, shows the decimal logarithm of the errors in the velocity, temperature, and pressure fields for the Schwarz iterations in the case of R = 10 4 and Γ = 3.495. These errors take values between 10 −4 and 10 −3 . Figure 14b shows the decimal logarithm of the errors in the velocity, temperature, and pressure fields for the Schwarz iterations in the case of R = 10 5 and Γ = 3.495. In that figure, the convergence of Schwarz for each Newton problem is tested. The errors for the Schwarz solver reach values between 10 −5 and 10 −4 . The increase in R needs an increase in the number of subdomains to reach convergence. There is a loss of accuracy of the method with R or with the increase in the number of subdomains. Now we increase the values of the aspect ratio and the Rayleigh number. The decimal logarithm of the errors in the velocity, temperature, and pressure fields for the Schwarz iterations in the case of R = 1300, Γ = 27.96, np = 14, mp = 2, and N × M = 16 × 10, can be seen in Figure 15a. The convergence of the Schwarz method for each Newton problem is shown in this figure. The errors for the Schwarz solver reach values between 10 −5 and 10 −4 for all the fields. Figure 15b shows the decimal logarithm of the errors in the four fields for the Schwarz iterations in the case of R = 5000, Γ = 27.96, np = 13, mp = 3, and N × M = 18 × 12. These errors take values between 10 −4 and 10 −3 . If we compare with the previous results in the case of Γ = 3.495, the increase in Γ needs an increase in the number of subdomains to reach convergence. There is a loss of accuracy of the method with Γ or with the increase in the number of subdomains.
We have increased further the aspect ratio for a fixed value of the Rayleigh number, R = 1300. Figure 16 shows a solution for Γ = 50 obtained with 20 subdomains on the horizontal axis and two subdomains on the vertical axis for a collocation mesh of 16 × 10. The number of rolls increases. Figure 17 shows a solution for Γ = 100 obtained with 35 subdomains on the horizontal axis and two subdomains on the vertical axis for a collocation mesh of 16 × 10. The number of rolls doubles. The number of necessary subdomains increases with Γ in a nearly logarithmic way. A noteworthy aspect is that the values of the aspect ratio can be increased as much as required.
The convergence rate for the Schwarz iterations appears in Figure 18a,b. The convergence rate are 0.7 in the first case and 0.9 in the second case. There is a dependence of the convergence rate with the parameters; it increases with R and Γ, although it always remains less than 1.

Newton
In this section, we analyze the behavior of the Newton method, the errors, convergence rate, number of subdomains, and influence of the parameters.
First, we investigate the influence of the Raylegh number for a fixed value of the aspect ratio Γ = 3.495. Figure 19a shows the decimal logarithm of the errors in the velocity, temperature, and pressure fields for the Newton iterations for R = 1300, Γ = 3.495, np = 2, and mp = 1. Newton's errors are 10 −8 for the temperature field and 10 −7 for both components of the velocity field and the pressure field. This figure shows the convergence of the nonlinear method. Figure 19b shows the decimal logarithm of the errors in the four fields for the Newton iterations for R = 5000, Γ = 3.495, np = 8, and mp = 2. The Newton errors are around 10 −3 for all the fields. Figure 20a shows the decimal logarithm of the errors for the Newton iterations for R = 10 4 and Γ = 3.495. These errors are between 10 −2 and 10 −1 . The convergence of the nonlinear problem is shown in this figure. Figure 20b shows the decimal logarithm of the errors for the Newton iterations for R = 10 5 and Γ = 3.495. The Newton errors are between 10 −2 and 10 −1 . The method is convergent. The increase in R needs an increase in the number of subdomains to reach convergence. There is a loss of accuracy of the method with R or with the increase in the number of subdomains.
Now we increase the Rayleigh number and the aspect ratio. The decimal logarithm of the errors in the velocity, temperature, and pressure fields for the Newton iterations for R = 1300, Γ = 27.96, np = 14, and mp = 2 are shown in Figure 21a. The Newton errors are around 10 −2 . Convergence of the nonlinear problem is shown in this figure. In Figure 21b, the decimal logarithm of the errors in the fields for the Newton iterations for R = 5000, Γ = 27.95, np = 13, and mp = 3 can be seen. These errors are around 10 −2 . The method is convergent. If we compare with the previous results in the case of Γ = 3.495, the increase in Γ needs an increase in the number of subdomains to reach convergence. There is a loss of accuracy of the method with Γ or with the increase in the number of subdomains.
Newton's method has a second-order local rate of convergence in the single-domain case. This has been shown in the case of the Legendre collocation in a single domain in Figure 22, where the numerical sequence |u n+2 x − u n+1 x | |u n+1 x − u n x | 2 , n = 1, 2, . . . for R = 1300, Γ = 3.495, np = 1, mp = 1, N × M = 50 × 50, has been plotted, and the tendency to 1 can be seen. For a single domain, the iterations are just Newton iterations. However, in the multi-domain case, Schwarz iterations are introduced at each Newton iteration, and the errors due to those iterations become relevant. In the multi-domain case with Schwarz iterations, Newton's method becomes first order, as can be seen in Figure 23, where the numerical sequence |u n+2 , n = 1, 2, . . . for the same case and two horizontal subdomains has been plotted, the tendency to 0.17 can be seen.

Computational Cost
Next, we study the computational cost. We approximate the number of operations considering only products. We name n = N × M, N t = 4n, N S : number of iterations of the Schwarz solver, N N : number of iterations of the Newton procedure. Each Newton and Schwarz iteration solves a linear system of size N t with a computational cost O(N 2 t ); then, to calculate the total computational cost, we multiply by the number of subdomains and the number of Newton and Schwarz iterations, then the total number of operations is No = np · mp · N S · N N · N 2 t . In the smallest case considered in this work, R = 1300, Γ = 3.495, N t = 4 × 10 × 10, N S = 5, N N = 12, np = 2, mp = 1, No = 2 × 10 7 . In the case of large aspect ratio, Γ = 50, np = 25, mp = 2, N S = 5, N N = 40, then No = 4 × 10 9 . In the case of large aspect ratio, Γ = 100, np = 50, mp = 2, N S = 5, N N = 40, then No = 6 × 10 9 . The case of large Rayleigh number R = 10 5 , np = 25, mp = 5, N S = 5, N N = 40, then No = 10 10 . In a parallel computation, if we do not take into account the computational cost penalty due to parallelization, the calculation would run at the same time in each subdomain, then the number of subdomains would be saved, and the number of operations would be reduced to Nop = N S · N N · N 2 t . Regarding computing time, the calculations have been performed with a 4 GHz Intel Core i7 microprocessor. For instance, in the case of R = 1300, Γ = 100, with 50 subdomains on the x axis and two subdomains on the z axis, with a Legendre mesh of N t = 4 × 16 × 10, the computing time is 612 s. This time could be reduced the order of magnitude associated with the total number of subdomains in an ideal parallel computation. The order of reduction is np · mp, and the time computing is N S × N N the time it takes to solve a linear system of size N t .

Discussion
This method includes several elements, validation, Legendre collocation method, Schwarz iterations, Newton iterations, turbulence, and computational cost.
From the point of view of validation, the equations have been solved with finite elements and with the present method. The same solutions are obtained with both methods. The same states of rolls and thermal plumes are observed in the real physical experiments [29].
We look now at the Legendre collocation method. The optimal mesh size of the expansions is a coarse mesh of order 10 × 10 in each subdomain. Legendre collocation is a high-order method in a single domain, but the global multi-domain method with the Schwarz and Newton iterations loses this property.
With respect to the Schwarz method, the Schwarz parts of the method converge, the errors are on the order of 10 −4 or less, and the convergence rate is less than 1 for all the fields and values of the parameters present in the problem. There is a dependence of the convergence rate with the parameters; it increases with R and Γ, although it always remains less than 1. The minimum distribution of necessary subdomains to reach convergence depends on R and Γ; it increases with both parameters with a logarithmic dependence.
In relation with the Newton method, it converges, and the final errors increase with the parameters R and Γ, varying from 10 −8 for the lower values R = 1300 and Γ = 3.495, to 10 −2 for the larger values R = 5000 and Γ = 3.495. Newton's method has a secondorder local rate of convergence in the single-domain case, but, in the multi-domain case that includes Schwarz iterations, Newton's method becomes first order. There is a loss of accuracy with the number of subdomains that can be related with some lack of regularity at the interfaces for Schwarz.
In relation to turbulence, the model can be near turbulence for the convenient values of R, but it cannot be turbulent because it is stationary. In Ref. [29], the authors consider the maximum value R = 10 5 and plumes are observed in the experiments, but the flow is not yet turbulent. In other references, such as [71], a power-law relation between the Reynolds number and the Rayleigh number, Re = R 1/2 , was found. A flow starts to be turbulent at approximately Re = 2000, following this power law; in our case, it corresponds to R = 4 × 10 6 , higher than our highest value R = 10 5 . Turbulence can be modeled with direct numerical simulation (DNS), which would be our model including time dependency. As we do not include the time dependency in this work, we only look for stationary solutions. Now, the highest value of R is 10 5 which is not yet turbulent. The time dependent model for turbulence can be easily implemented, and it will be addressed in future work. Because of this domain decomposition, the aspect ratio and the Rayleigh number can be increased considerably by adding subdomains in any direction. Then, values of the Rayleigh number for near turbulent regimes can be reached, and domains with large aspect ratio can be considered. Namely, the great advantage of this method is that we get solutions close to turbulence, or in domains with large aspect ratios, by solving systems of linear equations with well-conditioned matrices of maximum size one thousand. The main benefit of the developed method compared with finite differences, finite elements, spectral elements, continuous Galerkin, and discontinuous Galerkin is the small well-conditioned matrices that appear in this method. The maximum size of the matrix in the systems that are solved with this method is 10 3 in the case of a mesh grid 16 × 12 for a near turbulent flow with R = 10 5 . This is the main advantage over other methods that require solving systems with huge matrices of the order of several million, usually badly conditioned. In this case, ill conditioning is avoided and there is no need for preconditioners; even domain decomposition is used as preconditioner for those methods, as, for example, in Refs. [55,72].
Regarding the computational cost, the size of the collocation mesh is around 100 points, the problem has four unknowns, then the size of the matrices is 400 × 400; the linear solver needs this size squared operations, i.e., O(10 5 ), the Schwarz and Newton iterations are fixed to 5 and 40, respectively. Then, the computational cost is around the number of subdomains times 10 7 . The computational cost is comparable to other methods. The method is parallelizable; the number of subdomains can be eliminated from the count in this case, but then the computational cost penalty due to parallelization should be included.

Conclusions
In this work, an alternating Schwarz domain decomposition method to solve a Rayleigh-Bénard problem is presented. The problem is modeled with the stationary version of incompressible Navier-Stokes equations and the heat equation in a rectangular domain under the Boussinesq approximation. The system of partial differential equations is nonlinear. The nonlinearity has been solved numerically with a Newton's method. Each iteration of Newton's method has been discretized with an alternating Schwarz scheme, and each Schwarz iteration is solved with a collocation Legendre method. The Legendre collocation method is ill conditioned; for this reason it is not possible in a single domain to increase the mesh size appropriately to obtain solutions with many oscillations or solutions with structure at different scales. These solutions appear for large domains or large values of the Rayleigh number in the turbulent regime. This difficulty has been overcome with the use of an alternating Schwarz domain decomposition method. The original domain has been divided into several subdomains in both directions of the plane. Each step of the Newton and Schwarz iteration problems in each subdomain has been solved by Legendre collocation with a coarse mesh, then the problem in each subdomain is well conditioned. The present work achieves an efficient alternating Schwarz algorithm such that the number of subdomains can be increased indefinitely in both directions of the plane. The method has been validated with a benchmark with numerical solutions obtained with other methods and with real experiments. Because of this domain decomposition method, the aspect ratio and Rayleigh number can be increased considerably by adding subdomains. Then Rayleigh values close to the turbulent regime can be reached. The great advantage of this method is that we obtain solutions close to turbulence, or in domains with large aspect ratios, by solving systems of linear equations with well-conditioned matrices of maximum size one thousand. This is an advantage over other methods that require solving systems with huge matrices on the order of several million, usually with conditioning problems. The computational cost is comparable to other methods, and the code is parallelizable.

Conflicts of Interest:
The authors declare no conflict of interest.