A Galerkin/POD Reduced-Order Model from Eigenfunctions of Non-Converged Time Evolution Solutions in a Convection Problem

: A Galerkin/POD reduced-order model from eigenfunctions of non-converged time evolution transitory states in a problem of Rayleigh–Bénard is presented. The problem is modeled in a rectangular box with the incompressible momentum equations coupled with an energy equation depending on the Rayleigh number R as a bifurcation parameter. From the numerical solution and stability analysis of the system for a single value of the bifurcation parameter, the whole bifurcation diagram in an interval of values of R is obtained. Three different bifurcation points and four types of solutions are obtained with small errors. The computing time is drastically reduced with this methodology.

The Rayleigh-Bénard convection is a bifurcation problem where a two-dimensional horizontal fluid layer is uniformly heated from the lower plate. To model this problem, we use the continuity equation, the two-dimensional Navier-Stokes equations, and the heat equation, where the Boussinesq approximation is considered. The conductive quiescent state becomes unstable for a critical vertical temperature gradient, measured by the Rayleigh number R [32]. Beyond this threshold, a convective motion sets in. This state suffers transitions as the value of R is further increased. Some industrial setups concerning the use of Rayleigh-Bénard convection equations can be found in [33][34][35][36][37][38][39][40]. Reduced-order models based on POD usually operate in two phases: (1) an off-line phase, where proper bases for the problem unknowns are computed from snapshots, and (2) an on-line phase, where the original partial differential equations are projected over the aforementioned where Ω = {(x, z) ∈ R 2 : 0 < x < Γ, 0 < z < 1}, Γ = L/d, e z is the upwards vertical unit vector, Pr = ν/k is the Prandtl number, R = d 3 αg∆T/(νk) is the Rayleigh number, α is the thermal expansion coefficient, and g is the acceleration of the gravity. The boundary conditions (bc) are: v = 0, θ − 1 = 0 at z = 0, θ = ∂ z u = w = 0 at z = 1, ∂ x θ = ∂ x w = u = 0 at x = 0 and at x = Γ.

Stationary Solutions and Linear Stability Analysis
There exist several stationary solutions for this problem [54]. We denote with superscript b a stationary solution (v b , θ b , P b ). The linear stability of this solution is studied through the perturbation: where the tilde refers to the perturbation fields [50,55]. To quantify v := ( u, w), θ, P, and σ, Expression (11) is introduced into (6)-(8), together with (9) and (10). As a result, the following equations are obtained: where the second-order term e 2σt v · ∇ θ in the third equation will be neglected because we are interested in linear stability. Noticing that (v b , θ b , P b ) is a stationary solution, and removing the tildes to simplify notation, the perturbation fields and σ satisfy the following eigenvalue problem: with bc (9) and (10). The eigenvalues and their corresponding eigenfunctions depend on the stationary solution and, therefore, on the Rayleigh number R. We denote σ 1 as the largest real part of all the eigenvalues σ. To determine the stability of the stationary solution, we study the sign of σ 1 . If σ 1 < 0, the stationary solution is stable. If σ 1 > 0, the solution is unstable. In this context, we first name the eigenfunction of the stationary solution the eigenfunction related to σ 1 .

The POD Reduced-Order Method
One of the solutions of Problem (6)-(10) is the trivial one, with the fields equal to zero, named the conductive solution, denoted by Φ 0 . This solution is stationary and exists for any value of R. We fix a value of the Rayleigh number R 0 . We solve the linear stability analysis Problem (12)- (14) with bc (9) and (10) for the conductive numerical solution Φ 0 with a standard numerical method. A linear combination of the two eigenfunctions associated with the eigenvalues with a maximum real part is the initial condition to obtain a different non-trivial steady solution, solving Problem (6)-(10) with a standard time evolution discretization. This procedure generates some transitory numerical states before convergence φ 0 , φ 1 , φ 2 , . . . . We are taking the transitory states φ j such that ||φ j − φ j−1 || ∞ > . The linear stability analysis for those transitory states is numerically performed, and the first eigenfunctions for these states are the snapshots for the POD analysis. Thus, we obtain where i denotes the time step. This procedure is summarized in Algorithm 1 below.
Algorithm 1: Calculation of the snapshots.
Step 1: Fix a value of the Rayleigh number R 0 ; Step 2: Solve the linear stability analysis Problem (12)- (14) with bc (9) and (10) for the conductive numerical solution, Φ 0 , with a standard numerical method; Step 3: Solve Problem (6)-(10) with the standard time evolution numerical scheme using a linear combination of the two eigenfunctions related to the eigenvalues with a maximum real part, obtained in the previous step, as the initial condition to obtain a different non-trivial steady solution; Step 4: Take K transitory numerical states, obtained with the previous step, before Step 5: Solve, numerically, the linear stability analysis for those transitory states. The first eigenfunctions for these states are the snapshots for the POD analysis, These snapshots construct, in columns, the corresponding snapshot matrices, S, for the first and the second components of the velocity field, the temperature, and the pressure fields, S u , S w , S θ , and S P , whose columns are the velocity, temperature, and pressure snapshots, respectively. Once these matrices are obtained, we proceed to calculate the modified snapshot matrices M u , M w , M θ , and M P in terms of the snapshots, as is defined in Algorithm 2 below. These matrices are not the usual correlation matrices used in other works, and it has been proved they are more efficient [46]. We obtain the singular values λ i , i = 1, 2, . . . , K and the corresponding eigenvectors E i , i = 1, 2, . . . , K of the modified snapshot matrices. If the singular values decay quickly to zero, this means that the model reduction approach will be worth implementing.
Algorithm 2 describes, for the thermal snapshots S θ , the numerical method applied to compute λ i , i = 1, 2, . . . , K and their corresponding POD modes. The procedure is similar for the hydrodynamic fields.
Algorithm 2: Computation of λ i and POD modes.
In practice, matrices S P and M P are not considered and Algorithm 2 is not applied to the pressure field because the pressure field disappears from the variational formulation used to solve the problem; see the next section.

The POD/Galerkin Projection Procedure for the Stationary Problem
A Galerkin method is developed to solve the stationary Problem (6)-(10), where the temporal derivative is disregarded. The unknown fields are expanded in the orthonormalized POD bases calculated following Algorithm 2.
The variational form of the stationary problem can be written as follows [51]: The non-linearity is solved with a Newton procedure, and the velocity and temperature fields are expanded into the POD modes: v I = (u I , w I ) = ∑ I j=1 α j Q j v and θ J = ∑ J j=1 β j Q j θ . We note that the pressure field disappears from the variational formulation of Equation (7) because the POD modes are incompressible and the normal components of the velocity POD modes vanish at the boundaries [16].
Each linear problem of the Newton iteration of the stationary problem in its variational form (15) and (16) becomes a (I + J) × (I + J) algebraic system of equations.

Linear Stability Analysis of POD/Galerkin Solutions
Let us explain the POD/Galerkin method developed to solve the eigenvalue Problem (12)- (14) together with its corresponding boundary conditions. It is based on the POD modes B POD The variational form of the eigenvalue problem can be written as follows [56]: where the velocity and temperature fields that we seek are expanded into the POD modes: As in Section 3.1, the pressure field disappears from the variational formulation of Equation (13) because the POD modes are incompressible and the normal components of the velocity POD modes vanish at the boundaries [16]. B POD v,I and B POD θ,J are the orthonormalized POD bases for the velocity and temperature fields, where I and J are the number of unsaturated POD modes for each field.
The eigenvalue problem in its variational form (17) and (18) becomes a (I + J) × (I + J) algebraic eigenvalue problem as follows: where A 1 ∈ M I×I , A 2 ∈ M J×I , B 1 ∈ M I×J , B 2 ∈ M J×J , α = (α 1 , . . . , α I ); β = (β 1 , . . . , β J ) and I J is the identity matrix of size J. The matrix M is (I + J) × (I + J). Problem (19) has eigenvalues σ i , i = 1, . . . , I + J. As explained in Section 2, we are interested in studying the sign of the real part of the eigenvalue with largest real part, σ 1 . If σ 1 > 0, the solution is unstable; otherwise, the solution is stable.

Numerical Results
A Legendre collocation method is used as a space discretization. Two different Legendre grids are considered, one with 18 nodes in the x-direction and 14 nodes in the z-direction, and another with 36 nodes in the x-direction and 14 nodes in the z-direction. Thus, the maximum number of grid points is N = 36 × 14 = 504. The integrals in Equations (15)- (18) and the usual L 2 scalar product are performed with the Legendre-Gauss-Lobatto quadrature formulas [57]. Then, G is the diagonal matrix whose elements are the Legendre-Gauss-Lobatto weights. These expansions have good convergence properties [55]. The time evolution scheme is a finite-differences Euler implicit. We consider an aspect ratio Γ = 3.495, as in Refs. [50,51], and the parameter R takes values in the interval [1000; 2000].

First Bifurcation
The exhibited reduced-order model is applied to find a first bifurcation point. The following numerical results were obtained with a grid of 18 × 14 Legendre-Gauss-Lobatto nodes. First, we follow Algorithm 1 for the Rayleigh number R 1 = 1250. We solve the linear stability analysis Problem (12)- (14) with bc (9) and (10) for the conductive solution Φ 0 with a Legendre collocation method. A linear combination of the two eigenfunctions associated with the eigenvalues with a maximum real part is the initial condition to obtain the solution Φ 1 (R 1 ), solving Problem (6)-(10) with the time evolution discretization. This procedure generates some transitory numerical states before convergence φ 0 , φ 1 , φ 2 , . . . . We are initially taking the transitory states φ j , such that ||φ j − φ j−1 || ∞ > 10 −2 . A total of 66 transitory states meet this condition; however, we only consider the first 21, so K = 21. In Figure 2a, the horizontal component of the velocity in the middle of the cell u((m + 1)/2, (n + 1)/2), depending on the transitory time step, is plotted. We observe that the transitory states converge to the solution Φ 1 (R 1 ). In Figure 2b, the infinity norm of the difference between two consecutive transitory states of the horizontal component of the velocity field is displayed. As expected, these differences tend to zero. Isotherms of the transitory states at temporal steps t = 0.4, 0.8, 1.2, and 1.6 can be seen in Figure 3. These transitory states are far from the converged solution.
The linear stability analysis for the transitory states is numerically performed, and the first eigenfunctions for these states are the snapshots for the POD analysis. Therefore, a set of numerical snapshots, Figure 4 shows isotherms of the first eigenfunctions of the transitory states in Figure 3. It can be observed that some eigenfunctions are very similar because the transitory states are close in time.
Algorithm 2 is followed to calculate the modified snapshot matrices, M θ , M u , M w , the singular values and eigenvectors for those matrices, and the resulting POD bases. The singular values for the modified snapshot matrix M θ (R 1 ) are plotted in Figure 5. The number of unsaturated modes for temperature is J = 12, and for velocity is I = 11. We observe that the singular values rapidly decay to zero, so the POD analysis is reliable. The thermal POD modes, computed at R 1 = 1250, corresponding with the first four singular values λ i , i = 1, 2, 3, and 4, are represented in Figure 6. The complexity of the modes increases as the singular values decrease.

Validation and Results
Now, we solve Problem (6)- (10) for values of the Rayleigh number in the interval [1000; 2000] with the POD/Galerkin method and with Legendre collocation. For instance, the isotherms and the velocity field of the POD solution Φ 1 (R) for R = 1500 can be seen in Figure 7a. The same solution obtained with Legendre collocation is shown in Figure 7b, as in Ref. [50]. The L 2 norm of the difference between both solutions is order O(10 −3 ). Now, we solve the Galerkin/POD linear stability Problem (17) and (18) for Φ 1 (R). The first eigenfunction of the solution Φ 1 (1500), computed with the POD method, can be seen in Figure 8a. The norm of the horizontal component of the velocity field for the solutions Φ 1 (R) in the interval can be seen in Figure 9a; the continuous line has been obtained with Legendre collocation and the dashed line with POD. The difference between both solutions is order O(10 −1 ). Finally, a linear stability analysis for those stationary solutions has been performed with the POD method. A plot of the real part of the eigenvalue with the largest real part, σ 1 , can be seen in Figure 9b; the continuous line has been obtained with Legendre collocation and the dashed line with POD. The relative difference between both values of σ 1 is plotted in Figure 10; it is order O(10 −2 ), except at the bifurcation point R c1 = 1100, where we divide by a value near zero. The POD stability analysis in Figure 9b shows that σ 1 is negative in the interval, so solutions Φ 1 (R) are stable for any value of R in the interval The POD method is applied to perform the linear stability analysis of the conductive solutions Φ 0 (R). A plot of σ 1 can be seen in Figure 11; the continuous line has been obtained with Legendre collocation and the dashed line with POD. The relative difference between both values of σ 1 is plotted in Figure 12; it is order O(10 −2 ), except at the bifurcation point R c1 = 1100, where we divide by a value near zero. The POD stability analysis in Figure 11 shows that σ 1 is negative for R < R c1 and positive for R > R c1 . Then, the conductive solutions Φ 0 (R) are stable for R < R c1 and unstable for R > R c1 . We observe that, at R = R c1 = 1100, an exchange of stability takes place; a stable branch of solutions, Φ 1 (R), emerges while the conductive branch, Φ 0 (R), becomes unstable. These results coincide with those from Ref. [50], obtained with Legendre collocation.
Note that, with the information obtained in just a single value of the Rayleigh number, R 1 = 1250, the bifurcation diagram for Φ 0 (R) and Φ 1 (R) solutions, R ∈ [1000; 2000], and its linear stability has been computed. Indeed, only 21 transitory states and their corresponding eigenfunctions have been considered at R 1 . The number of thermal POD modes is J = 12, and the number of hydrodynamic POD modes is I = 11.

Second Bifurcation
Now, the presented reduced-order method is applied to calculate a second bifurcation point. The following numerical results were obtained with a grid of 18 × 14 Legendre-Gauss-Lobatto nodes. First, we follow Algorithm 1 for the value of the Rayleigh number R 2 = 1300. We solve the linear stability analysis Problem (12)- (14) with bc (9) and (10) for the conductive solution, Φ 0 , with a Legendre collocation method. A linear combination of the two eigenfunctions corresponding to the eigenvalues with the largest real part is the initial condition to obtain the steady solution Φ 2 (R 2 ) with the time evolution discretization. This procedure generates some transitory numerical states before convergence φ 0 , φ 1 , φ 2 , . . . . We are taking the transitory states φ j such that ||φ j − φ j−1 || ∞ > 10 −2 . A total of 64 transitory states meet this condition. The first 13 are disregarded. Half of the remaining states have been taken, namely, the even ones. Then, K = 26. In Figure 13a, the horizontal component of the velocity in the middle of the cell u((m + 1)/2, (n + 1)/2), depending on the transitory time step, is plotted. We observe that the transitory states converge to the solution Φ 2 (R 2 ). In Figure 13b, the infinity norm of the difference between two consecutive transitory states of the horizontal component of the velocity field is displayed. As expected, these differences tend to zero. Isotherms of the transitory states at temporal steps t = 0.4, 0.8, 1.2, and 1.6 can be seen in Figure 14. The exhibited transitory states are far from the converged solution.
The linear stability analysis for the transitory states is numerically performed, and the first eigenfunctions for these states are the snapshots for the POD analysis: Ψ i (R 2 ) ≡ (V i (R 2 ), Θ i (R 2 ), P i (R 2 )), i = 1, . . . , K. Figure 15 shows isotherms of the first eigenfunctions of the transitory states in Figure 14. Some eigenfunctions are similar because the transitory states are close in time.
As in Section 4.1, Algorithm 2 is followed to calculate the modified snapshot matrices M θ , M u , M w , the singular values and eigenvectors for those matrices, and the resulting POD bases. The singular values for the modified snapshot matrix M θ (R 2 ) are plotted in Figure 16. The number of unsaturated modes for temperature is J = 18, and for velocity is I = 10. We observe that the singular values rapidly decay to zero, so the POD analysis is reliable. The thermal POD modes, computed at R 2 = 1300, corresponding with the first four singular values λ i , i = 1, 2, 3, and 4, are represented in Figure 17. The complexity of the modes increases as the singular values decrease.

Validation and Results
Now, we solve Problem (6)- (10) for values of the Rayleigh number in the interval [1000; 2000] with the POD/Galerkin method and with Legendre collocation. For instance, the isotherms and the velocity field of the POD solution Φ 2 (R) for R = 1500 can be seen in Figure 18a. The same solution obtained with Legendre collocation is shown in Figure 18b, as in Ref. [50]. The L 2 norm of the difference between both solutions is order O(10 −3 ). Then, we solve the Galerkin/POD linear stability Problem (17) and (18) for Φ 2 (R). The first eigenfunction of the solution Φ 2 (1500), computed with the POD method, can be seen in Figure 8b. The norm of the horizontal component of the velocity field for these solutions Φ 2 (R) in the interval is displayed in Figure 19a; the continuous line has been obtained with Legendre collocation and the dashed line with POD. The difference between both solutions is order O(10 −1 ). Finally, a linear stability analysis for those stationary solutions is performed with the POD method. A plot of σ 1 can be seen in Figure 19b; the continuous line has been obtained with Legendre collocation and the dashed line with POD. The relative difference between both values of σ 1 is plotted in Figure 20; it is order O(10 −1 ), except at the bifurcation points R c1 = 1100 and R c2 = 1558, where we divide by a value near zero. In this case, the computed eigenvalue shown in Figure 19b, σ 1 , corresponds to the stability of Φ 0 (R) in the interval [1000; 1252], and to the stability of Φ 2 (R) in the interval [1252; 2000]. As we already calculated, Φ 0 (R) is stable for values of R < R c1 and unstable for R > R c1 . At R = R c2 = 1252, a new solution Φ 2 (R) emerges from Φ 0 (R). This solution is unstable in the interval [1252; 1558] because σ 1 is positive in this interval for this solution; this solution is stable for R > R c3 = 1558 because σ 1 is negative there. These results coincide with those in Ref. [50], obtained with Legendre collocation.
As in Section 4.1, we remark that, with the information obtained in just a single value of the Rayleigh number, R 2 = 1300, the bifurcation diagram for Φ 0 (R) and Φ 2 (R) solutions, R ∈ [1000; 2000], and its linear stability has been computed. Only 26 transitory states and their corresponding eigenfunctions have been considered at R 2 . The number of thermal POD modes is J = 18, and the number of hydrodynamic POD modes is I = 10.

Third Bifurcation
Finally, the presented reduced-order method is applied to calculate a third bifurcation point. The following numerical results were obtained with a grid of 36 × 14 Legendre-Gauss-Lobatto nodes. This time, a finer numerical discretization was required to capture the bifurcation point. We follow Algorithm 1 for the value of the Rayleigh number R 3 = 1300. We solve the linear stability analysis Problem (12)- (14) with bc (9) and (10) for the conductive solution, Φ 0 , with a Legendre collocation method. A linear combination of the two eigenfunctions corresponding to the eigenvalues with the largest real part is the initial condition to obtain the steady solution Φ 3 (R 3 ) with the time evolution discretization. This procedure generates some transitory numerical states before convergence φ 0 , φ 1 , φ 2 , . . . . We are taking the transitory states φ j such that ||φ j − φ j−1 || ∞ > 2 · 10 −2 . A total of 41 transitory states meet this condition. Only the last 21 are considered, so K = 21. In Figure 21a, the horizontal component of the velocity in the middle of the cell u((m + 1)/2, (n + 1)/2), depending on the transitory time step, is plotted. We observe that the transitory states converge to the steady solution Φ 3 (R 3 ). In Figure 21b, the infinity norm of the difference between two consecutive transitory states of the horizontal component of the velocity field is displayed. Again, these differences tend to zero. Isotherms of the transitory states at temporal steps t = 0.4, 0.8, 1.2, and 1.6 can be seen in Figure 22. The exhibited transitory states are far from the converged solution.     The linear stability analysis for the transitory states is numerically performed, and the first eigenfunctions for these states are the snapshots for the POD analysis: where i is the temporal step. Figure 23 shows isotherms of the first eigenfunctions of the transitory states in Figure 22. Some eigenfunctions are very similar because the transitory states are close in time. As in previous subsections, Algorithm 2 is followed to calculate the modified snapshot matrices: M θ , M u , M w , the singular values and eigenvectors for those matrices, and the resulting POD bases. The singular values for the modified snapshot matrix M θ (R 3 ) are plotted in Figure 24. The number of unsaturated modes for temperature is J = 15 and for velocity is I = 7. We observe that the singular values rapidly decay to zero, so the POD analysis is reliable. The thermal POD modes, computed at R 3 = 1300, corresponding with the first four singular values λ i , i = 1, 2, 3, and 4, are represented in Figure 25. The complexity of the modes increases as the singular values decrease.

Validation and Results
Now, we solve Problem (6)- (10) for values of the Rayleigh number in the interval [1000; 2000] with the POD/Galerkin method and with Legendre collocation. For instance, the isotherms and the velocity field of the POD solution Φ 3 (R) for R = 1900 can be seen in Figure 26a. The same solution, obtained with Legendre collocation, is shown in Figure 26b, as in Ref. [50]. The L 2 norm of the difference between both solutions is order O(10 −2 ). Then, we solve the Galerkin/POD linear stability Problem (17) and (18) for Φ 3 (R). The first eigenfunction of the solution Φ 3 (1900), computed with the POD method, can be seen in Figure 8c. The norm of the horizontal component of the velocity field for this solution Φ 3 (R) in the interval can be seen in Figure 27a; the continuous line has been obtained with Legendre collocation and the dashed line with POD. The difference between both solutions is order O(10 −1 ). Finally, the POD method is applied to perform the linear stability analysis. The corresponding plot of σ 1 can be seen in Figure 27b; the continuous line has been obtained with Legendre collocation and the dashed line with POD. The relative difference between both values of σ 1 is plotted in Figure 28; it is order O(10 −2 ), except at the bifurcation points R c1 = 1100 and R c3 = 1558, where we divide by a value near zero. The computed eigenvalue shown in Figure 27b, σ 1 , corresponds to the stability of Φ 0 (R) in the interval [1000; 1252], the stability of Φ 2 (R) in the interval [1252; 1558], and the stability of Φ 3 (R) in the interval [1558; 2000]. As we already calculated, Φ 0 (R) is stable for values of R < R c1 and unstable for R > R c1 . At R = R c2 = 1252, solution Φ 2 (R) appears from Φ 0 (R). These solutions are unstable in the interval [1252; 1558] because σ 1 is positive in this interval for these solutions and it is stable for R > R c3 = 1558, as can be seen in Figure 19b. Solution Φ 3 (R) emerges from Φ 2 (R) at R c3 . These solutions are unstable for R > R c3 because σ 1 is positive there. These results coincide with those in Ref. [50].
As in previous subsections, it should be noted that, with the information obtained in just a single value of the Rayleigh number, R 3 = 1300, the bifurcation diagram for Φ 0 (R), Φ 2 (R), and Φ 3 (R) solutions, R ∈ [1000; 2000], and its linear stability has been computed. Only 21 transitory states and their corresponding eigenfunctions have been considered at R 3 . The number of thermal POD modes is J = 15, and the number of hydrodynamic POD modes is I = 7.

Bifurcation Diagram
The different branches of solutions and bifurcations are summarized in the bifurcation diagram ( Figure 29) obtained with the POD method. The bifurcation digram represents the norm of the component u of the velocity field as function of the parameter R (horizontal axis). The horizontal line represents the conductive solution Φ 0 ; the velocity field for the conductive solution is zero (in particular u). The branches of stable stationary solutions are marked with solid lines and the unstable solutions with dashed lines. The conductive stationary solution constitutes the conductive branch of solutions, Φ 0 (R), which exhibits two pitchfork bifurcations as R increases in [1000; 2000]. The first takes place at R c1 = 1100, where Φ 0 (R) loses its stability and a new branch of solutions, Φ 1 (R), emerges. These solutions Φ 1 (R) are stable for any value of R in the interval [1000; 2000]. The second is a pitchfork bifurcation at R c2 = 1252 that produces a branch of solutions Φ 2 (R) emerging from Φ 0 (R) at this critical threshold R c2 . At R c3 = 1558, a subcritical pitchfork bifurcation occurs. Solutions Φ 3 (R) emerge from Φ 2 (R). Solutions Φ 2 (R) are unstable for R < R c3 = 1558, but at higher values, they become stable. Solutions Φ 3 (R) are unstable for R > R c3 . These results coincide with those in Refs. [46,50].

Computational Cost
Now, we approximate the number of operations considering only products at each step of the algorithms. We name N * = (m + 1) × (n + 1), N = 4N * , N t : number of time steps computed by the standard time evolution solver, N R : number of values of R for which the problem is solved, and N N : number of Newton iterations for solving the variational problem. We consider first the number Breakdown of the computational cost of Algorithm 1:

•
Step 2: Solve the linear stability analysis Problem (12)- (14) with bc (9) and (10) with a standard method: O(N 3 ); • Step 3: Solve Problem (6)-(10) with the standard time evolution scheme: N t · O(N 2 ); • Step 5: Solve, numerically, the linear stability analysis for K transitory states: Breakdown of the computational cost of Algorithm 2: • Step 3: Construct the matrices M θ , M u , and M w from thermal and hydrodynamic snapshots: 3KN * ; • Step 4: Apply a SVD decomposition to M θ , M u , and M w to obtain their eigenvectors and singular values: 3 · O(N * 3 ); • Step 5 (15) and (16): Solving the eigenvalue Problem (17) and (18): The cost of the off-line calculation is O((K + 1)N 3 ). The cost of the on-line calculation is O(N N N R N * 2 ). The calculation of a branch of the bifurcation diagram with the standard method is N R O(N 3 ) and with the POD method is O(N N N R N * 2 ). In our case, the reduction goes from O(10 12 ) to O(10 8 ). However, taking into account the on-line and the off-line, the POD computational cost is O((K + 1)N 3 ), in our case O(10 11 ). The Legendre collocation method calculates 201 solutions in the interval of R, [1000; 2000], by steps of five with a computational cost of 0.0279 h in a 3.6-GHz Intel Core i9 microprocessor. The calculation of the eigenvalues for those solutions takes 1.862 h. Therefore, the total computational cost for the Legendre collocation method is 1.8899 h. Regarding the POD method, the off-line computational cost, including the snapshots, the singular values decomposition, the POD bases, and the stability calculation is 0.2008 h. The on-line computational cost of the POD solutions and the corresponding eigenvalue problem with the same interval and steps is 0.4344 s, or 1.21 · 10 −4 h. Then, the total computational cost for the POD method is 0.2009 h.
We obtain a reduction in the CPU time of a factor of 9.4, i.e., O(10), as previously estimated in the computational cost. Taking into account only the on-line computational cost, the factor of reduction is 15619, i.e., O(10 4 ), as was seen in the computational cost.

Conclusions
A Galerkin/POD reduced-order model from eigenfunctions of non-converged time evolution states in a Rayleigh-Bénard problem is presented. This convection problem is modeled in a rectangular domain with a heat equation coupled with the incompressible Navier-Stokes equations, depending on the Rayleigh number R as a bifurcation parameter. This problem is solved with a time evolution Legendre collocation method for a value of the Rayleigh number. The eigenfunctions of the linear stability analysis of the non-converged states are the snapshots of the POD analysis. From a single value of the bifurcation parameter, the whole bifurcation diagram in an interval of values of R is obtained. Three different bifurcation points and four types of solutions are obtained with small errors. The calculations use less than 30 modes, so the matrices for the eigenvalue problem solved with the POD method are very small, and the computational cost is drastically reduced, by a factor of 10 4 .