Steady Rayleigh–Bénard convection

2D steady state solution of coupled momentum, continuity and energy equations. The fluid is driven by the Rayleigh–Bénard convection.

Governing equations:

Steady-state momentum in x and y, continuity, energy equation
The equations are non-dimensionalized with respect to Pr and Ra numbers. For the detailed description of the equations, please refer to "Ouertatani, et al, Numerical simulation of two-dimensional Rayleigh–Bénard convection in an enclosure. C. R. Mecanique 336 (2008)". The final results are validated against this paper as well.

Boundary conditions (2D square):

u=v=0 on all four sides
non-dimensional temperature = -0.5 on top and +0.5 on the bottom
zero temperature flux on left and right sides of the domain

Numerical method:

SIMPLE algorithm used for resolving velocity-pressure coupling. Please check out "Versteeg, Malalasekera: an introduction to computational fluid dynamics" text book. The velocity grids are staggered from the pressure grids. Temperature has the same grid as pressure. The velocity and temperature are updated using Jacobi method in every iteration and pressure correction equation is solved directly using a penta-diagonal matrix algorithm. Proper under-relaxation factors should be chosen for convergence. Jacobi method is the least efficient iterative method but it is parallelizable. You are encouraged to use other iterative methods such as, line by line TDMA, Guess-Seidel or SOR. The Power-law scheme has been used for upwinding momentum and energy equations. You may find the datils in "Patankar, Numerical Heat Transfer and Fluid Flow" text book.

I have written this code in Parallel using PETSc (C language). I am going to upload it on GitHub pretty soon. It's much faster than this version because:
1) It is parallel and not serial so you can use multiple processors.
2) The C version solves the momentum and energy equations by iterating a couple of times, unlike the MATLAB version that does only one Jacobi sweep.
3) Most importantly, PETSc uses Krylov subspace iterative method (KSP type:GMRES) and it also uses preconditioners (I found Additive Shwartz PC to be the best) to solve momentum, pressure correction and temperature, which is way more efficient than Jacobi method in terms of convergence rate.