Computational fluid dynamics (CFD) provides the mathematical and computational tools for the numerical simulations of fluid flows, an essential component in a number of applications ranging from geophysics to industrial applications. It can be used to forecast the areas subject to potential inundation by natural flows (tsunamis, lava flows, landslides, pyroclastic flows), to evaluate the impact and cost of strategies aimed at reducing long- and short-term risk from natural and man-made hazards, or to test mechanical designs and verify their efficiency before the costly implementation processes. In many applications, multiple fluids (e.g. oil, water, air) may be present in the same problem, and the very different rheological properties of these fluids pose a challenge for some CFD methods. The work described in this report is focused on the implementation and validation of two state-of-the-art multi-fluid models for the Smoothed Particles Hydrodynamics (SPH) numerical method. The SPH method, in its weakly-compressible formulation, is highly parallelizable, and is therefore an excellent candidate for impementation on parallel computing hardware, such as Graphic Processing Units (GPUs). While GPUs have been primarily developed for high-performance three-dimensional graphic rendering and video-games, since 2007 the two major manufacturers have also exposed their massive parallel computing capabilities to more generic applications. Our work is integrated into GPUSPH [Hérault et al., 2010], the first implementation of SPH to use CUDA-enabled GPUs as high-performance parallel computing devices. The report begins with an introduction to the SPH method and its application to the Navier-Stokes equations for fluid dynamics, highlighting the issues that the classical SPH formulation encounters in the case of multiple fluids. This is followed by an overview of GPUSPH, the implementation of SPH on CUDA GPUs which is the basis for our implementation of the multi-fluid SPH models: the Hu & Adams model [Hu and Adams, 2006], and the more recent one developed by Grenier [Grenier et al., 2008; 2009; Grenier, 2009]. A description of these models, highlighting the differences over the classic SPH formulation, is included in the subsequent chapter. Finally, validation tests are shown for problems in both two and three dimensions, covering the classic Rayleigh-Taylor instability, the lock exchange gravity current problem, and finally a rising bubble example. To the best of our knowledge, the work described here represent the first fully three-dimensional implementation of these multi-fluid SPH models on the high-performance computing platform provided by CUDA-enabled GPUs.