Buoyancy-driven flows are widespread in diverse engineering applications. Such flows have been studied in great detail theoretically, experimentally, and numerically. However, the fluid-dynamic instabilities and flow reversals of thermosiphon are still actively investigated. The presence of such instabilities limits the effectiveness of such devices for decay heat removal. Traditionally the stability analysis of natural convection loops has been confined to one-dimensional calculations, on the argument that the flow would be mono-dimensional when the ratio between the radius of the loop and the radius of the pipe is much larger than 1. Nevertheless, accurate velocity measurements of the flow in toroidal loops have shown that the flow presents three-dimensional effects. Previous works of the authors have shown that these structures can be seen in thermosiphons. In this paper, we aim to use modern CFD methods to investigate the three-dimensional flow in thermosiphons. This paper focuses on rectangular thermosiphons. In particular, we perform a series of high-fidelity simulations using the spectral element code Nek5000 to investigate the stability behavior of the flow in a rectangular thermosiphon. We compare the results with available existing experimental data from the L2 facility in Genoa. We examine in detail the flow structures generated. Moreover, in the past various authors have demonstrated that the overall behavior of the thermosiphon depends strongly on the boundary conditions (BCs). The simulation campaign was carried out with different BCs to investigate and confirm this effect. In particular, simulations with Dirichlet, Neumann and Robin BCs for heater and sink were performed.