Laminar-Turbulent Bifurcation Scenario in 3 D Rayleigh-Benard Convection Problem

We are considering two initial-boundary value problems for Rayleigh-Benard convection in Oberbeck-Boussinesq approximation for incompressible fluid in 3D-rectangular domain with 4:4:1 geometric ratio with periodicity in two directions and cubic domain with 1:1:1 ratio and zero velocity and temperature gradient boundary conditions. For this purpose, we use two numerical method: one is a Pseudo-Spectral-Galerkin method with trigonometric-Chebyshev polynomials and the other is finite element/volume method with WENO interpolation for advection term. Numerical methods are presented shortly and are benchmarked against known DNS data and against one another (for quasi-periodic domain problem). Then we perform stability analysis using analytical expression for main stationary solutions and eigenvalue numerical analysis by applying Implicitly Restarted Arnoldi (IRA) method. The IRA is used to perform linear stability analysis, find bifurcations of stationary points and analyze eigenvalues of monodromy matrices. Thus characteristic exponents of the system for time periodic solutions (limited cycles of various periods and resonance invariant tori) are computed. We show, numerically, the existence of multistable rotes to chaos through chaotic fractal attractors, full Feigenbaum-Sharkovski cascades and multidimensional torus attractors (Landau-Hopf scenario). The existence of these attractors is shown through analysis of phase subspaces projections, Poincare sections and eigenvalue analysis of numerically computed DNS data. These attractors burst into chaos with the increase of Rayleigh number either through resonance and phase-locking or through emergence of singular chaotic attractors.


Introduction
The bifurcation analysis of Rayleigh-Benard convection was inspired by the 0-th modal approximation-the Lorenz system.The latter is known to show chaotic behavior and is a classic example of such ODEs (formulated as the 14th Smale problem).For the rigorous proof of the Smale's 14th Problem, you can see [1].The first analysis of the phenomenon was conducted by Lord Rayleigh in [2].Details about Rayleigh-Benard convection in general are described in [3], where reader can find information about linear analysis, secondary flows, experiments and other useful information.The bifurcation analysis of the full system of Navier-Stokes equations was formulated in some papers later, see [3] [4] [5] [6].Note that most of these papers are dedicated to 2D convection [6] or low mode problems [5].Good review on the problem in general is presented in the Paul Manneville's chapter, see [7].However, a large amount of review is dedicated to intermittency with almost no focus on Landau-Hopf scenario and Feigenbaum-Sharkovskii scenario.On the other hand, we were able to obtain some results in previous papers, see [8].We show that the problem branches itself with the transition to chaos either through bifurcations of limited cycles (so called Feigenbaum-Sharkovskii-Magnitskii scenario, see [8]) or through Landau-Hopf scenario with the formation of high dimensional tori in the phase space.The present work is a revision of these results for cubic wall bounded domain and presentation of new results obtained for cuboid periodic domain.In this paper, we also apply analysis of Monodromy matrix eigenvalues to confirm some bifurcations and transition mechanisms.
The paper is laid out as follows.Firstly, the Initial-Boundary value problem is posed.Then, the analytical data concerning linear stability are presented in order to perform benchmark of numerical methods.The next section includes numerical methods: Pseudo-Spectral-Galerkin method, Finite Element/Volume method and the IRA matrix-free eigenvalue solver.We outline some properties of these methods.Then we show some benchmark results for Rayleigh-Benard convection problem.We compare eigenvalues with some known data and analytical expressions; we also compare DNS results for moderate Rayleigh numbers with known statistical results.In the last section, we show result for bifurcation analysis in the considered two domains.The final sections are discussion and conclusion.

Initial-Boundary Value Problem
We are considering Oberbeck-Boussinesq approximation for incompressible Navier-Stokes equations, i.e. the temperature dependence of the fluid parameters is nulled for all but the density in the buoyancy term.This term varies with temperature linearly.Let equipped with initial-boundary conditions:  is a solution on a 2D torus, i.e. periodic, in appropriate direction, l n -unit outward normal to the l x boundary, p is the pressure; Pr µ χ = , ( ) ( ) ( ) are Prandtl and Rayleigh numbers, where g is the magnitude of gravity; β is the fluid thermal expansion coefficient; h is the length in the direction of gravitational force ( z L in our case); χ is the fluid thermal con- ductivity; µ is the fluid kinematic viscosity; h θ is the temperature on the hot plane of the domain and c θ is the temperature on the cold plane of the domain, ( ) . The nondenominational form is derived if the time scale is chosen as characteristic time for momentum transfer by viscosity through the layer of height Pressure is not explicitly defined and is treated differently for every numerical method that we use.

Stability of the Main Solution
We are following [3] to show the analysis of stability for the main stationary solution.
Considering a layer of fluid that is defined on 1 Ω with 0 1, 2 z z L = − = and boundary conditions of temperature set to 1 h θ = , 0 c θ = .The boundary conditions for velocity can be either set up as non-slip, i.e.Ω we get ( ) the stability is defined for the linearly varying temperature profile with no dependence on , x y .So we will consider the stability in the form of normal modes with fixed wave vector { } , ,0 where λ is the increment and ( ) , w x y is the planar waveform function being a solution to the Helmholtz equation ( ) Applying normal mode analysis to (1) gives the following boundary conditions for (4): The case for the free slip condition can be easily solved by choosing eigenfunctions in the form of ( ) that satisfy boundary conditions.The solution of the quadratic equation for λ is: for the case of stability loss 0 n λ = one can find that ( ) ( ) and also see that the stability of the main solution is not dependent on Pr .However coming up with eigenfunctions for the non-slip case is much harder so we use numerical approach, since we are only interested in the value for cr Ra itself.We seek solution for the critical point 0 λ = in the form e qz f = .Plugging conditions into (4), one gets: that has the following six roots: So the solution for f is presented in the form: The solution for the constants can be organized in the homogenoius system of linear equations using boundary conditions (5): 0, = Ab (10) where is the column vector of the unknown constants A, B,•••, F and matrix A depends on the boundary conditions.We are not interested in finding b explicitly, but rather determine the condition for the homogenoius liner system to have a solution and from that condition derive critical values of Ra and k , called cr Ra and cr k , respectively.For the free-slip boundary the expression for the matrix is: e e e e e e q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q wise.For non-slip boundary conditions the matrix is: q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q k q q k q q k q q k q q k q q k q q k q q k q q k q q k q q wise.The determinant of A must be zero in order for the system (10) to have a nontrivial solution.So for a fixed value of k we must solve the nonlinear equation ( ) It is done by using Newton method as: here the Jacobi matrix is found using analytical differentiation technique.Iterations stop if On the other hand, knowing from the literature [3] the critical values of Rayleigh number, one can find a critical wavenumber k using the same method for fixed cr Ra .It should be noticed, however, that the solution should be pure real despite that the roots are complex.In Figure 1 we show neutral curves for

Numerical Methods
In this section we give information about numerical methods that are used to solve the problem.

Pseudo-Spectral-Galerkin Method
The Pseudo-Spectral-Galerkin Fourier-Chebyshev method is applied to the domain 1 Ω (Fourier-Galerkin method or FGM for short).All vector-functions in (1) are expanded by the divergence free basis functions ( )

, , jkl x y z w
, constructed analogues to [5]: We can check that the following bases functions are analyticaly divergence free, i.
, , 0 matter the expressions for them.For the case of 1 Ω we use the following set of scalar bases functions: , , , hence, d .
In order to meet homogeneous Dirichlet conditions we must have: With all this together we get the following scalar basis functions in z direction: It is a straightforward way to check that (16) and, hence, (12) are complied with (14) and (15).Such functions form a basis in From here we turn our attention to the truncated series for basis functions, so that the problem can be solved on the computer.We assume that there are , , x y z N N N number of polynomial modes, i.e. x y N N is the number of Fourier modes and z N is the number of Chebyshev polynomials in z direction.Please note that the number of degrees of freedom is lesser since all complex coefficients for Fourier modes are subject to reality condition, i.e. ( ) ( ) ˆ, , . So the total number of degrees of freedom is ( ) . In this case the basis scalar components (13) are rewritten as: , , Since the basis is divergence free and the mean flow through periodic boundaries is zero, this implies that the integral of velocity over these boundaries is zero.In this case it is straightforward to show that the pressure is eliminated from (1) by projecting pressure gradient into the divergence free functional subspace that is formed by the span of divergence free basis functions.
We use the following scalar basis functions for the temperature expand: , , e , where ( ) l z χ are constructed such that boundary conditions for temperature are satisfied: With the correction term ( ) z − we satisfy boundary conditions for temperature on 1 z = ± .The cost of the full Bubnov-Galerkin method being applied to Equation (1) is of ( ) ( ) due to the non-linear term.These multiplication terms become convolution in the functional space thus causing multiplications of tensors rank 3.Such computation complexity is very limiting.In order to reduce the computational cost of calculations we use two stage transfer from physical space to functional space and use Fourier collocations with pseudo-spectral approach in xy direction.
First we span the functions in (1) using Discrete Fourier Transfer (DFT) in regular grid points, forming the following system (taking into account divergence-free nature of basis): with B being a convolution term.All coefficients For every , j k point we apply Bubnov-Galerkin projection in z direction: where we denote: as projections of DFT basis corrected coefficients to polynomial space and − and minus sign in diffusion matrices is due to the in- tegration by parts and zero boundary conditions.Since these matrices are independent of , j k indexes, we can apply them for each point , j k in z -direction.The mass matrices are full since the basis we use is not orthogonal but all matrices are not singular.It can be proved by the fact that we are using linear combinations of polynomials that form basis in 2  .So the vectors are linear intendant and so the matrices are positive-definite, hence, invertible and the system can be solved for unknown coefficients ( ) jkl t u , ( ) jkl t θ . In order to perform integration we use exact symbolic integrals for diffusion and inverse mass matrices that are calculated in Wolfram Mathematica and stored for further use in the program.We use Gauss-Chebyshev quadrature in z direction.This quadrature is exact for polynomials of degree 2 1 N − and can be efficiently used.We use it to perform projection from domain to image and back for (22), (23) and we find DFT divergence-free coefficients at points p z as where p z points are Gauss-Chebyshev quadrature points and where [ ] ( ) ( )   so boundary points belong to the boundary.
The nonlinear (multiplication) term is calculated using pseudo-spectral approach.
We calculate derivatives in polynomial space, then return to the physical space at specific points and calculate multiplication in physical space with computational difficulty ( ) . Having known coefficients ( ) at time 0 t we perform the following steps to calculate ( ) (the advection term in scalar energy equation is calculated analogously).
1) Calculate derivatives in functional spaces: where D1 is the differentiation matrix for our basis functions in z direction.
2) Increase the size of arrays for u and , , x y z ∂ u in x and y direction by 1 3 and fill added elements with zeros, so we have arrays with sizes of 3 2 3 2 3) Return from functional space to physical space step by step to get ( ) 0 , , , x y z t u and ( ) , , 0 , , , 4) Calculate multiplication in physical space at every point to get: ) 6) Apply ( 22) using Gauss-Chebyshev quadrature to get This leads to no aliasing of frequencies.
The following approach requires less operations then the exact Bubnov-Galerkin method of ( ) ( ) , and for most practical use (say ) the limiting factor would be ( ) ( ) N N N  so we assume that our method requires ( ) ( ) operations for every time step.
Explicit Runge-Kutta 4 (RK-4) method that is used to integrate the semi-discrete system (21) in time.In order to satisfy stability condition we analyze spectra of the linear operator and bounds for the nonlinear part and derive stability condition.A necessary condition for stability by the linear change of mapping (by the change of the time step) is the location of all discrete spectra of the spacial operator inside the RK-4 stability region on a complex plane.The discrete Fourier spectra has a standard estimate [9] and spectral norm for Chebyshev matrices is used for the estimate.This part is beyond the scope of the paper.

Finite Element/Volume Method
We consider another method that is applied in 1 Ω and 2 Ω domains.It uses nodal Finite Element approach to reconstruct pressure and Finite Volume/Difference method for the reset part of the equations.Finite element method uses values of pressure and velocity in vertexes of elements to form matrix equations.Finite volume/difference method uses values of velocities in centers of elements to approximate integral/ differential operators.So this is a combination of finite element method and finite volume method that we call "FEM" for short.
We start with discretization of Stokes operator: , 0, in bounded domain s Ω with appropriate initial-boundary conditions.We use some discretization method that we discuss later and projection method (see [10]) to translate (30) into: , 0.
Here M is a mass matrix; Q is a gradient matrix; A is a diffusion matrix and G is a divergence matrix.We introduce time slices with ( ) l being l -the time slice with time-step t ∆ and derive the following system: p + is unknown on l -th time slice, we split the system as: N Q , and we get correction equation for the potential function φ : where matrix After the solution of (34) we correct velocity and pressure functions in such way, that where 0 1 β ≤ ≤ is a parameter.This projection method is another way to write Helmholtz-Hodge decomposition for specific type of equations.The whole step of solution consists of using three steps (33), (34), (35).In general this results in first order of approximation for (30) in time.The stability and consistence of the method depend on the discretization.It is known that the discretization must obey Ladijenskaya-Babuska-Brezzi (LBB) inf -sup condition in order to be stable and consistent.In case of the provided approximation it means that sponds to sup condition) and that ( ) , where κ is a condition number and h is a discretization parameter (e.g.grid size).In general one usually takes T = G Q , chooses mixed finite element ap- proximation (for finite element methods) or staggered grid (for finite difference approximation) and changes matrix N to N that can be easily inverted, e.g.
= N E .Here we use different strategy that is dealing with different approximations for different operators.

Let us introduce rectangular cuboids jkl
W that from a 3D tessellation of a rectangular domain , where j is a multi index with j being a center of W j .We introduce another set of tessellation U k that is constructed from swapping central nodes and vertexes, thus each vertex of W j becomes a center for U k and vice versa.
We define basis functions in an element Q k , where Q can be W or U , as follows: Now we use the following expansion in this element space for a scalar function ( ) Now we consider set of points j formed by the centers of W j or vertexes of U k .
Such elements can be considered as finite volumes, i.e. with : Q W = and the second operation is an inverse of (37) with : Q U = .In this case the scheme can be written as follows: In this work we use compact finite difference scheme of 4-th order to approximate A and G using method of alternating directions.The approximation of L and Q is done using Bubnov-Galerkin projection (integration over Ω ), e.g. for the second equation in (38): Inserting (37) into (39) and doing integration by parts: It is a straightforward way to check the BBL condition (we don't consider this in the papaer).Trivial kernel of G and Q is proved by considering space of finite elements and 4-th order compact differences schemes.The condition number of finite element approximation can be estimated from the space of finite elements and can be shown that it has a marginal bound.Now it is a straightforward way to return to Navier-Stokes equations by applying some approximation for the nonlinear term in (38) (defined bellow as B ).We use 7-th order WENO scheme that has good spectral properties and guaranties TVB behavior of the solution (on each WENO stage we use Runge-Kutta 3rd order SSP method [11]).In order to increase the temporal accuracy we also use Runge-Kutta 3rd order explicit method [11] for which the projection step (38) is applied on every stage.The stability of the method is deduced from CFL condition since diffusion is considered in implicit way.The usage of this type of finite elements gives one more positive result.There is no artificial boundary layer near wall boundary.
Equation on φ requires zero Neumann boundary conditions on the wall, but in this FEM setup it just means skipping values of φ on zero Neumann boundary.The same is true for the pressure gradient approximation.Then the total scheme for (1) is given as: ( ) , . 3 3

Solution of the Eigenvalue Problem
We briefly give description of the matrix free eigenvalue solver that is used for the problem.More information about this approach can be found in many papers, for example [12] [13] [14] [15].Considering discrete systems ( 21) and (41).Let us consider another discrete systems for small temporal perturbations s u and s θ formed in the vector: Inserting those into discrete systems and linearizing one can gets the following set of equations: where ( ) where , , , , , , , , .
for FEM system.Here  stands for linearization of the operator (40).The size of the perturbation vector is 4 and is used in the Implicitly Restarted Arnoldi (IRA) method.Note that if we are using Fourier method, this vector includes real and imaginary parts of Fourier modes that are treated as real values and the size changes to ( ) . The system (43) is used as an operator A that maps the vector (42) as Av .These perturbations are automatically divergence-free since we are using divergence-free basis for Fourier-Galerkin method.For FEM method these perturbations are guaranteed to be divergence-free since we are applying projection algorithm with pressure initialized during the projection.There's no need to introduce pressure perturbations in both cases.The algorithm for IRA to find eigenvalues of A is based on [16] and goes as follows: 1) Initialization.Initialize vector v with different random values.Then normalize the vector, so 2) Arnoldi Step.We form the Krylov subspace as . We use the following process: , , , This precess generates the following decomposition: of H are the eigenvectors (eigenvalues) of A .If the process is not stopped (which is usually true), then continue to 2. Find eigenvalues of Hessenberg matrix 0 1 1 , , , k m λ λ λ + − , sort them in an appropriate order (either maximum real part or maximum magnitude in our applications).Perform QR algorithm with shifts for H using polynomials with number of shifts q m = : 1, , Set 1 s k m q = + − − .Note that at this point q m ≠ since some of j λ ∈  .Apply shifts: , At this point we have matrices V , H with additional value In order to find eigenvalues of Jacobi matrix we consider the system (43) with zero temporal derivative.So the system (43) is called on every stage of Arnoldi process by applying the linearized equations to the vector in Arnoldi step.If we are computing eigenvalues of Monodromy matrix then the system is applied as follows.Let the system (43) have a periodic solution with period  .Then the system is formed as where ( ) M  is a state transition operator or a Monodromy matrix.To find M in a matrix-free way we integrate the linearized system in time for a period  : Now we apply (48) for the input vector in Arnoldi step of IRA, while integral is evaluated by applying the selected time stepper.

Implementation Details
Our goal is to perform direct numerical simulation (DNS) of problems and trace transition from initial stationary point in the phase space (main solution) to the chaos through cascades of bifurcations.In order to detect and clarify bifurcations we use the analysis of eigenvalues of Jacobi and Monodromy matrices and phase portraits with Poincare sections.The calculation of the whole IRA algorithm and achievement of the statistically quasi-periodic solution regimes for relatively high Ra numbers require significant computational power.We use the same idea as in [8] to adapt number of harmonics or number of elements for the problem during the analysis of bifurcations.
The control of the accuracy is done by the analysis of the energy spectrum for the whole simulation time.We define the correlation tensor and its Fourier transfer For the discrete problem we are using FFT to calculate (51).In order to get the energy spectra we integrate over the spherical shell: that becomes a summation for the discrete problem.Then we check, that for [ ] with N being the size of FFT discretization.The last relation (53) in physical interpretation means a track of the solution to be in a deep dissipation regime.It is an overdiscretization from a standard DNS point of view, however it is essential to obtain complex bifurcations in near chaos region.For all calculations we are using  39) is carried out using geometric multigrid approach [17].
The IRA algorithm is using dot product, matrix-vector operations and matrix matrix operations from MAGMA library across multiGPUs, wheres QR routine for the upper Hessenberg matrix is taken from LAPACK.The visualization is done using GMSH [18], N. M. Evstigneev Gnuplot and LibreOffice, simple calculations are done in MATLAB.

Benchmarks
At first we perform the verification of our methods vs. known results.The first benchmark is to obtain the neutral curve by applying different cr k as it is done in Section 3.For this purpose we consider domain 1 Ω with the following perturbations of temperature: , , 2 sin , 50 where cr k is a given critical wave number, and we take Pr 1 = .The results for numerical methods are brought in Figure 1 with comparison to the linear stability analysis.We can see that maximum deviation for FEM ( 64 64 30 × × elements) method is about 11.2% and for FGM is about 6.5% ( 32 32 20 × × modes).Please note that further increase of degrease of freedom did not improve the results much.We are able to trace exact points of transition with the application of the IRA solver.The results for leading eigenvalues are brought into Table 1.Please note that for both cases we have leading eigenvalues of multiplicity two.So the supercritical pitchfork bifurcation is observed, that complies with well known data about Rayleigh-Benard convection.In physical space we can observe the formation of rolls that run parallel to one of the axis x or y , depending on the direction of perturbation vector in (54).Velocity vectors for these rolls are shown in Figure 2. One can observe small amplitude of velocity but the amplitude increases with the increase of Ra , with the formation of a mushroom type distribution for temperature.
Velocity vectors and temperature distribution for , , 1, 2, 3, here R is an averaged Raynolds number taken analogous to [20] for the considered problems, ( ) Instantaneous snapshots of temperature distributions are shown in Figure 4 and

Bifurcations and Route to Chaos
Some bifurcations for the problem in 2 Ω were presented in [8] as a survey of the results with the main emphasis on the Feigenbaum-Sharkovskiy sequence of cycles.In here we give new results for both setups.

Bifurcations in Domain with XY Periodicity
After the first bifurcation that is shown in Benchmarks section the flow stays steady that corresponds to the point in the phase space that remains in this point for  for FEM with the formation of limited cycle in the whole phase space.The following process for FGM is depicted in Figure 7.The limited cycle increases its amplitude and for 5100 = Ra is shown in Figure 8 alongside with the eigenvalues of Monodromy matrix.The corresponding velocity and temperature distributions are shown in Figure 9 and trajectories in physical space with the magnitude of leading eigenvector are shown in Figure 10.It is clear that one eigenvalue is located at the unit cycle at the point ( ) 1, 0 + that corresponds to the limited cycle.It is true for both numerical methods, although the convergence of IRA is faster for FGM method.Other eigenvalues are situated inside the unit circle and are stable.One can see the difference in eigenvalues for FEM and FGM methods, however both methods give correct leading eigenvalue.As one can see from the eigenvector, the flow is formed by the bending of the rolls with the subsequent oscillation.Although the attractor dimension is just one (limited cycle), the flow in physical space is already complicated.It can be traced with the stream lines obtained by Lagrangian particle tracers (Figure 10) that are forming a complected path in the physical space.
At this point one may observe the effect of multistability.If perturbations of magnitude With the further increase of Ra number we see the increase of amplitude of the limited cycle and with the sequential formation of the invariant two dimensional torus.
It can be expected since complex-conjugate eigenvalues of the Monodromy matrix in for FEM.This leads to the formation of the limited torus in the phase space, whose projection into three-dimensional subspace is shown in Figure 12 and physical space functions are shown in Figure 13.Please note that in many papers the bifurcation that leads to the formation of invariant torus is called Neimark-Sacker bifurcation.This is a    Such exact value was chosen in order to perform eigenvalue analysis (was found using quasi-Newton method).The period during integration  in (48) was defined by the return map in the Poincare section.Please note that the calculation of these eigenvalues took about a month on a 5GPU cluster.Corresponding eigenvalues are presented in Figure 15.Two pure real eigenvalues have magnitude close to unity (0.999 and 0.997 and we assume that these eigenvalues are of magnitude one) and are located at the ( ) 1, 0 + point on the unit circle.This is another way of saying that the system has two zero Lyapunov exponents and all other exponents are negative (inside the unit circle on the complex plane).This corresponds to the phase-lock of three frequencies: two are connected through period doubling and another was irrational to them both.It is interesting that there are two more eigenvalues are closing the point ( ) mean a possibility of period doubling bifurcation.However with the increase of Ra , after the phase-lock, the system continuous to maintain the same attractor (2D two period torus with increasing magnitude) and suffers anther phase-lock (with another frequency of the period doubling) at about 8350 Ra = (see Figure 15).Again, two pure real eigenvalues are of unit magnitude (1.019 and 0.998) correspond to zero Lyapunov exponents.However there are more eigenvalues closing to the unit circle in Please note that we present data only from FGM here since FEM was not used for this branch.Evolution of return maps is shown in Figure 16.One can see that the onset of chaos emerges very fast in parameter space through local hyperbolicity.With this the attractor dimension remains bounded between 2 and 3.This can be stated since the  ODEs are constructed in [22].However we may assume another scenario of fast and multiple period doubling bifurcations with the formation of singular Feigenbaum attractor in the Poincare section.This can be implied through the analysis of some leading eigenvalues that are located closer to the ( ) 1, 0 − point for resonance case, see Figure 15.But for this case (close to the resonance) the IRA algorithm filed to converge so we cannot justify this scenario.We are unable to reveal any localized structures for 8382 Ra > in the phase space for this branch.
Another branch is observed if small perturbations (of  22, the idea of its construction is taken from [6].We will not stop on detail analysis of every solution and only point out that further increase of Ra leads to the   ).Example of such intermittency is presented in Figure 23 along with cycle period three from Sharkovsky cascade and its Monodromy matrix eigenvalues.Finally the system emerges into chaos through singular attractor and from 9550 Ra > we are unable to detect any regular structures in the phase space.

Bifurcations in Bounded Cubic Domain
We are only discussing results here that were not mentioned in [8].This subsection is focused on eigenvalues since we were unable to perform eigenvalue analysis in our previous papers due to computational limitations.All results are obtained with FEM.
There were three series of experiments conducted in [8] for various Pr numbers.Each series resulted in different scenarios of transition to chaos.As a matter of fact, this was not only due to Pr change, but also due to the symmetries in the system.At first, all scenarios have a common initial stage-supercritical pitchfork bifurcation.The first one was found in Section 5.The Rayleigh number of all bifurcations is higher compared to periodic domain due to the wall stabilization effect.If the initial perturbations are given as discussed in Benchmark Section 5 (dominated along one of the axis), then the flow develops a symmetrical solution in one plane, see Figure 24.
Another set of initial conditions may lead to other symmetries.For example, corner structures are formed if an initial condition are taken with constant perturbation of magnitude  in the form: ( )       intermittency, as suggested in [7].The dependence on Pr was not investigated, however it is known, see [7], that the route to turbulence depends on it and intermittency emerges more for high Pr (spacial instabilities are more elaborated in this case).The emergence of a 4D torus in Landau-Hopf scenario is still a question.It can be confirmed by the author's new idea (to be published soon) of constructing an ε -net of splines over the attractor and tracking its evolution on the attractor.If the structure in question is a torus and its mapping in phase space is diffeomorphic, then the ε -net will converge.But it requires enormous amount of computational power and time.
Another question is the automatization of the process of bifurcation detection and eigenvalue analysis as it was suggested in [12] [23].Another question is the emergence of traveling waves for periodic problem and explicit study of it's influence on bifurcation scenarios.All these questions are topics for further research.

Conclusions
In this paper, we present results for laminar-turbulent transition in Rayleigh-Benard We show that there are many routes to turbulence in 3D Rayleigh-Benard convection problem.If the flow is bounded only from one direction, then the onset of turbulence emerges for small Ra numbers.For all bounded domain, the typical values of Ra number are 30 -50 times higher.In all cases, the initial stage is analogous to 2D Rayleigh-Benard convection with the emergence of pitchfork bifurcation followed by Andronov-Hopf bifurcation, see [6].However, the scenario is different from here on.
There is more symmetry in 3D problem and symmetry groups may become generators for Hopf bifurcations.It was studied in [24] on an ABC flow.If the symmetry is preserved exactly (that can only be achieved numerically by the application of high order methods of quasi-spectral accuracy and detailed discretization), then the system undergoes bifurcations of limited cycles.In this case, we confirm the existence of multiple Figenbaum-Sharkovsky sequences of cycle periods with direct and inverse directions, see bifurcation diagram in Figure 22.We also observed multistability and existence of intermittency.Multistability can be explained by the neutral curve in Figure 1, since for high Rayleigh numbers there are many possible attracting sets.
Intermittent solutions are only confirmed for limited cycles.If the symmetry is broken, then an invariant torus emerges through the secondary Hopf bifurcation.It was confirmed by the analysis of Monodromy matrix eigenvalues with the emergence of complex conjugate eigenvalue with magnitude greater than unity.At this point, the solution may vary and be continuous either through Landau-Hopf scenario of a N-d tori cascade or through the formation of resonant torus with further phase-locking and possible local hyperbolicity.We also show how the mapping to itself of the Poincare section in the singular torus works using Lameray diagram.All these "tori" routes lead to chaos much faster than that is for cycle route.It was noted in [6], that there are chaotic attractors and stationary points coexisting for 2D case after the development of chaotic solution for high Ra numbers.We are unable to confirm this yet for the considered 3D cases.However, it seems likely in 3D as well, due to multistability.But in 3D case, the basis of attraction of these point attractors is smaller and, hence, it is more difficult to find this kind of system behavior.

2 Ω
with almost everywhere Lipschitz continuous boundaries s ∂Ω .We are interested in the solution of the following problem: Problem 1.For given values of Pr, , Ra τ + ∈ » , 1, 2 s = , find fluid velocity vector- function [ ]  and scalar function of fluid temperature [ ] where * () is a complex conjugate.Solving (3) for ( ) f z one gets the following homogeneous equation[3]:

Figure 1 .
Figure 1.Neutral curves for free-slip and non-slip boundary conditions.Dot plots are numerically calculated points of first bifurcation for given cr k and non-slip boundary conditions for 1 z = ± .Asterisks for Fourier-Galerkin method and circles for Finite Element/Volume Method.
in[5] for 2D and 3D Dirichlet box cases.We use the relation for Chebyshev polynomials of the first and the second kind to get: that we use to approximate the solution of the initial-boundary value problem (2) for (1).Note, that the basis (12) is not orthogonal in z direction but orthogonal in , x y directions due to the use of Fourier basis functions.It can be shown by the construction of the Mass matrix using an 2  projection.

D
are mass and diffusion matrices for velocities and temperature, respectively for every point , j k .Matrix sizes are z z N N × and they are formed as: ].The most operation-hungry part is the nonlinear term calculation using 3/2 padding.Now we assume that FFT can be used for DFT with ( )2 log N N operation, for example one can use fftw or cufft on GPU.Calculation of derivatives is done using return to the image from the domain of the mapping requires the same difficulty as for the image to domain transfer.So we get maximum difficulty as and introduce velocity correction vector C u and scalar potential function φ : We now define the following differential operators: Q and S in space of nodal finite elements W , A and G in space of finite volumes W , we use Laplace operator L for the approximation of S and identity matrix for N and M .Define projection of nodal operator to central finite volume point as  and projection of central finite difference/volume operator to nodes as  .The first operation is performed using (37)

1 g
j are coefficients of expansion for Dirichlet and Neumann boundary conditions and ( ) j G u are coefficients of divergence operator projection into the space of finite elements.Other operators are derived analogously.

.
Select number of eigenvalues k that are desired and number of additional vectors m for the implicit procedure, so the dimension of Krylov subspace is k m + accelerate computations we are using multiple Graphic Processor Units (multiGPU).The GPUs used are k40 NVIDIA GPUs, all programs are implemented on C++ with CUDA C. The application of DFTs is done by the CUFFT library on 2 or 4 GPUs.The matrix vector products are conducted using MAGMA library.The solution of the Poisson Equation (

2 Ω are shown in Figure 3
using FEM.The critical pertur-

Figure 5 .
Figure 5.Note that for the latter distribution the boundary layer is thinner.Comparison of statistical characteristics with available data from [19] is shown in Figure 6.It is clear that the statistical data is close to the provided simulation data.Please note, that the current results are obtained with higher order methods and, thus can be more exact.But in general the distribution is similar and so one can conclude that the proposed numerical methods for relatively hight Raylight numbers are correct.The leading direction of the flow is determained by the initial perturbations.In order to compare results with [19] we used x -dominated perturbations.
leading eigenvalue with multiplicity two crosses the imaginary axis.This results in another supercritical pitchfork bifurcation, forming solution that is symmetrical in another plane.Velocity vectors and temperature distribution are shown in Figure7.We can see that the solution is now rotated π 4 with the formation of similar rolls.

Figure 7 . 1 Ω
Figure 7. Velocity vectors and temperature distribution in

5 1 10 Figure 11 .Figure 8 .
Figure 11.It is clear that a formation of distorted square structures is presented.If we trance this solution back by decreasing Ra we will get the square tile structures.For

Figure 12 are
Figure 12 are closing to the unit circle on the complex plane.In order to obtain these results we had to run IRA with the same randomly initialised vector for all values of Ra .Secondary Hopf bifurcation occurs at 5893.4 Ra = for FGM and 5943.3Ra =

Figure 10 .Figure 11
Figure 10.Lagrangian particle movement in 1 Ω and magnitude of leading eigenvector that corresponds to maximum magnitude eigenvalue of Monodromy matrix, located at 1,0 + on the complex plane for FGM.(a) Lagrangian particle movement.(b) Modulus of the leading eigenvector for velocities.

Figure 12 .
Figure 12.Leading eigenvalues of Monodromy matrix near second Hopf bifurcation as functions of Ra and 2D invariant torus in phase subspace.(a) Evolution of Monodromy matrix leading eigenvalues.(b) Projection of invariant 2D torus into three-dimensional phase subspace.

Figure 13 .Figure 14 .
Figure 13.Temperature and velocity distributions in 1 Ω for 6000 Ra = that coresponds to the invariant torus in the phase space.(a) Velocity vectors.(b) Temperature isosurfaces.

Figure 14 .
Figure 14.Invariant 2D torus period 2 and Poincare section for it with close resonant torus.(a) Projection to the three dimensional phase subset, 7300 Ra = .(b) Poincare section using plane

Figure 15 .
Figure 15.Leading eigenvalues for resonant 2D torus period two in two cases of phase locking.(a) Leading eigenvalues of the resonant torus, 7421.5623Ra = .(b) Leading eigenvalues of the resonant torus, 8350.0712Ra = .

..
Figure 22, the idea of its construction is taken from[6].We will not stop on detail

Figure 18 .
Figure 18.Transition to chaos through the Landau-Hopf scenario.(a) Phase subspace and Poincare sections for8350 Ra = .(b) Poincare section for

Figure 19 .
Figure 19.Transition to chaos through the Landau-Hopf scenario-formation of 4D torus.(a) Second Poincare sections for 8590 Ra = that corresponds to stable 3D invariant torus.(b) Second Poincare section for 8618 Ra = and 8618.5 Ra = .

Figure 25 .Figure 21 .
Figure 25.These solutions are symmetric relative to planes: , and 1 z x y x y ∀ = =− .The direction of symmetry is selected only by initial conditions.Further increase of Ra number leads to the increase of amplitude of velocities.Finally, the solution can either go through supercritical pitchfork bifurcation or through Andronov-Hopf bifur-

Figure 25 .
Figure 25.Symmetry solutions for modified velocity initial conditions at 130000 Ra = .(a) Temperature distribution.(b) Sections of temperature distribution.(c) Sections of velocity magnitude distribution.

Figure 26 .
Figure 26.Eigenvalues of Jacobi matrix for pitchfork and Hopf bifurcations and Monodromy matrix for 150250 Ra = .(a) Hopf bifurcation at 150208.58Ra = .(b) Pitchfork bifurcation at

Figure 27 .
Figure 27.Singular attractor phase subspace projection and Poincare sections for 376370 Ra = .(a) Cycle projection into 2D phase subspace.(b) Poincare section using velocity.(c) Poincare section using velocity and curl.

Figure 28 .
Figure 28.Cycle period 3 and singular attractor phase subspace projections, Poincare section of the singular attractor and eigenvalues of Monodromy matrix.(a) Cycle period 3 projection into 2D phase subspace for 308565 Ra = .(b) Singular attractor projection into 2D phase subspace for 308672 Ra = .(c) Poincare section of the singular attractor.(d) Evolution of Monodromy matrix eigenvalues.
convection from the nonlinear dynamics point of view.In order to analyze simple bifurcations, we use Implicitly Restarted Arnoldi eigenvalue solver implemented on top of Navier-Stokes solvers.All methods are extended to MultiGPU architecture for acceleration. d , . ,

Table 1 .
Leading eigenvalues for the first bifurcation for four critical wave numbers.