The reliability of fluid–structure interaction (FSI) computational models depends on the accurate solution of both the fluid (Navier–Stokes) and solid governing equations, as well as their mutual coupling. The resulting numerical systems typically involve millions, if not billions, of degrees of freedom and require substantial computational resources. Even on modern supercomputing platforms, such simulations often demand long execution times. In recent years, graphics processing units (GPUs) have been widely adopted to accelerate high-performance computing applications; however, their use in FSI simulations remains limited, due to the multiphysics nature of the problem and the challenges associated with efficient parallelization. In this work, we present the open-source release of IBflows, our multi-GPU CUDA Fortran code designed to address large-scale FSI problems [1]. The solver integrates the incompressible Navier–Stokes equations using a fractional-step method with second-order central finite differences on a staggered Cartesian grid. No-slip conditions on solid boundaries are enforced through several immersed boundary (IB) techniques, including a moving least-squares (MLS) and related improvements [2]. Both rigid and deformable solids are supported. For the latter, a spring-network formulation based on the Fedosov interaction potential [3] is used. The accuracy, stability, and parallel scalability of the framework are demonstrated through representative FSI benchmarks, with particular attention to challenging added-mass regimes [4]. [1] Viola, F., Verzicco, R, et al. (2022). Computer physics communications, 273, 108248. [2] Vagnoli, G., Viola, F., et al. (2025). Computer Physics Communications, 109741. [3] Fedosov, D. A. (2010). Multiscale modeling of blood flow and soft matter. Brown University. [4] Sotiropoulos, F., \& Yang, X. (2014). Progress in Aerospace Sciences, 65, 1-21.