Introduction

Team: Bernard Kleynhans, Shaan Desai, Yaniv Toledano and Sébastien Lemieux-Codère

Background

Problem Description (Heavily Inspired by Sondak et al., 2015):

Thermal convection is a phenomenon found in a wide variety of physical processes. Examples include, but are not limited to, the motion of the Earth’s molten core, atmospheric and oceanic dynamics and solar phenomena (Sondak 2015).

A common approximation in studies of convection is that fluid parameters are constant except for the fluid density in the buoyancy force (Drazin & Reid 2004; Lappa 2009; Chandrasekhar 2013). This approximation, known as Oberbeck–Boussinesq, is valid across many physical situations and allows systems to be described by the incompressible Navier–Stokes equations subject to a buoyancy force and an advection–diffusion equation for the temperature field (Sondak 2015).

Specifically, Rayleigh–Bénard convection focuses on flow fields driven by buoyancy and confined between two infinite horizontal plates. A key question in thermal convection is how much heat is transported for a given temperature difference. In context of Rayleigh–Bénard convection, this is asking how the non-dimensional heat flux Nu depends on the Rayleigh number Ra (Sondak 2015).

Goal of Analysis:

The goal of this project is to parallelize an existing Fortran code (by Sondak et al.) using a combination of MPI and OpenMP in order to enable simulations of Rayleigh–Bénard convection at high Rayleigh numbers not feasible (or highly impractical) using the serial code base.

Need for High Performance Computing

Steady state flows that optimize heat transport have been obtained using a serial implementation for Rayleigh number up to $Ra \approx 10^9$. For higher Rayleigh numbers, serial execution can take in the order of days to months of compute time depending on granularity and problem size.

Specifically, the implemented time integrator uses an accurate (and stable) Implicit-Explicit Runge-Kutta algorithm which frequently calls Fast Fourier Transforms. As a result, the time complexity of the overall algorithm scales worse than $O(n \,log \,n)$ (given FFT’s scale as such) and likely in the order of approximately $O(n^2 log \,n)$.

Based on these observations, it should be evident execution time of the serial code base is a bottleneck for large-scale problems.

Profiling Current Implementation

Profiling of the existing serial Fortran code has illustrated that a significant majority of the compute time is spent on the time integrator as illustrated in the table below:

With this observation in mind, the main focus of the parallel implementation will deal with this module primarily.

Theory

The below theoretical background aims to set out the implementation used in the existing serial Fortran code and as such heavily relies on details provided in (Sondak 2015).

Governing Equations:

We wish to determine the temperature distribution $T(x,y,t)$ and the velocity field $\mathbf{u}(x,y,t) = (u,v)$ in the channel. The $x$-direction is periodic and the $y$-direction is bounded by two horizontal boundaries. The temperature boundary conditions at the top and bottom wall are given by

The boundary conditions for the velocity field will be described after introducing the governing equations. The fields evolve according to the Oberbeck-Boussinesq equations, which just means incompressible Navier-Stokes with a buoyancy force due to gravity and temperature changes. The temperature evolves according to an advection-diffusion equation. Note that all fields are dimensionless.

The variable $\phi$ has been introduced to help enforce the no-slip boundary conditions on the vertical velocity field, which are given by

The appendix of (Sondak 2015) contains details on how $\phi$ is used to enforce the boundary conditions.

Solution Approach:

The above PDEs are discretized in space and time separately. The temporal discretization uses the aforementioned Implicit-Explicit Runge-Kutta algorithm. In terms of the spatial discretization, recall the geometry of the Rayleigh–Bénard problem is periodic in the $x$-direction and bounded at $y=\pm 1$. As such, the implemented spatial discretization in $x$ is accomplished via Fourier transforms. Contrarily, in the wall-normal direction (i.e., $y$-direction) a non-uniform (using on a cosine grid) finite difference stencil based on Lagrange polynomials is used.

References

Ascher, U.M., S.J. Ruuth, and B.T.R Wetton. 1995. “Implicit-explicit methods for time-dependent partial differential equations”. SIAM Journal on Numerical Analysis 32 (3): 797–823

Chandrasekhar, S. 2013 Hydrodynamic and Hydromagnetic Stability. Dover

Drazin, P. G. & Reid, W. H. 2004 Hydrodynamic Stability. Cambridge University Press

Lappa, M. 2009 Thermal Convection: Patterns, Evolution and Stability. John Wiley & Sons

Sondak, David, Leslie M Smith, and Fabian Waleffe. 2015. “Optimal heat transport solutions for Rayleigh–Bénard convection”. Journal of Fluid Mechanics 784:565–595