Scalable Fully Implicit Runge-Kutta Solvers for Time-Dependent Rayleigh-Bénard Convection

We consider the efficient numerical solution of time-dependent Rayleigh-Bénard convection governed by the incompressible Boussinesq equations. The coupled velocity-pressure-temperature system presents significant computational challenges arising from incompressibility, buoyancy coupling, convection, and the simultaneous resolution of multiple physical time scales. We present a systematic study of fully implicit Runge-Kutta methods for time-dependent Rayleigh-Bénard convection, together with scalable finite element solvers for the resulting stage-coupled systems. We employ high-order, L-stable Radau IIA methods for temporal integration and investigate Taylor-Hood and Alfeld-Sorokina finite elements in two dimensions, along with Taylor-Hood and enriched macroelement discretizations in three dimensions. The macroelement approaches offer strong local inf-sup stability and patch structures well suited to Vanka-type relaxation. In two dimensions, monolithic multigrid achieves nearly stage-independent outer Krylov iteration counts and maintains robust convergence under mesh refinement. On the same mesh with incremental pressure correction, selected Radau configurations reduce elapsed solve time by 18-32% and yield smaller sampled transient differences relative to a numerical reference solution. Three-dimensional experiments demonstrate favorable parallel scalability with Taylor-Hood elements and substantially improved discrete mass conservation with the enriched macroelement discretization. These results identify regimes in which coupled multistage fully implicit integration is computationally competitive when equipped with appropriate multigrid preconditioners.