The axial flow in rod bundles has been the object of several investigations, since they are relevant for various industrial applications (e.g., heat exchangers, nuclear reactor cores). For tight configurations (pitch-to-diameter ratio smaller than 1.1) the large difference in velocity within the cross section creates the possibility of a Kelvin-Helmholtz instability and the generation of a vortex street in the gap ("gap instability"). While the presence of spacing devices may perturb the flow, simple grid spacers-grids of the type usually encountered in sodium fast reactors-do not prevent this instability. This work seeks to improve our understanding of the gap instability in large rod bundles by using high-fidelity largeeddy simulation (LES). The application of LES methods has been historically limited to single-pin calculations because of the large cost of full-bundle calculations. However, using Nek5000, a massively parallel spectral-element computational fluid dynamics code, we are able to carry out well-resolved LES simulations for this class of geometries. Results are compared with experimental data. For the largest bundle case, over 8 billion collocation points at a polynomial order of 20 are judged necessary to achieve overall excellent accuracy. More than 500,000 processors are used to carry out the simulations (costing approximately 50,000,000 CPU-hours), and 3,000 snapshots of the flow field have been collected to apply several coherent structure recognition techniques, including one of the largest proper orthogonal decompositions carried out to date. Also investigated here is the effect of geometry, including the (presence of spacers, changing pitch-to-diameter ratio, and canister wall. In particular, global linear stability analysis is applied to a series of simplified geometries in order to gain insight into the physics of the gap instability.