Many cardiovascular flow systems involve tight interaction between blood flow and soft tissue mechanics (e.g. heart valves, pulse propagation). The computational modelling of these systems is challenging as they often comprise complex flow patterns (e.g. flow separation, instabilities, turbulence), complex constituive laws (e.g. anisotropy, hyperelasticity), and fluid-structure instabilities (e.g. valve leaflet fluttering). Heart valves present also a topological problem as valve closure splits the flow domain in two separate volumes. Here, we present a high-performance high-fidelity solver for fluid-structure interaction (FSI) which was designed to study heart valve mechanics. It comprises a Navier–Stokes solver with sixth-order finite differences on a staggered grid and a third-order Runge–Kutta time integration scheme. The structure solver uses a finite-element discretization of the full elasto-dynamics equations for the soft tissue, a Newmark scheme for time stepping, and hyperelastic material models. A motility formulation of the immersed boundary method (IBM) is used to couple the flow and structure solvers: Flow velocities are passed to the structure solver as Dirichlet boundary conditions, and reaction forces of the deformed structure are returned to the flow solver as a force density. This is in contrast to many other FSI solvers which pass fluid pressures to the structure solver, and structure velocities back to the flow solver. Benefits of the motility formulation for simulating heart valves will be critically discussed. Furthermore, a variational formulation is used to transfer velocities and forces between fluid and structure meshes in a weak sense. This variational transfer has benefits over pointwise interpolation schemes of the classical IBM because it conserves integrals over the overlapping regions. We will discuss the effects of this approach on the well-known problem of leaky structures in the IBM. Finally, we will present examples of heart valve simulations including FSI instabilities and laminar-turbulent transition and we will demonstrate high-performance computing capabilities of a GPU-accelerated version of the solver.