Due to significant multi-scale heterogeneity, understanding sub-grid structures is critical to effective continuum-based description of gassolid flow. However, it is challenging for both physical measurements and numerical simulations. In this article, with the macro-scale pseudo-particle method (MaPPM) implemented on a GPU-based HPC system, up to 30,000 fluidized solids are simulated using the NS equation directly. The destabilization of uniform suspensions and the formation of solids clusters are reproduced in two-dimensional suspensions. Distinct scale-dependence of the statistical properties in the systems at moderate solid/gas density ratio is observed. Obvious cluster formation and its effect on drag coefficient are shown in a system at high solid/gas density ratio. On the computational side, about 19 folds speedup is obtained on one GT200 GPU, as compared to a mainstream CPU core. The necessity for investigating even larger systems is prospected. (C) 2010 Elsevier Ltd. All rights reserved.