In the present study, we developed a CPU-based discrete element method (DEM) to tackle the challenges associated with modelling powder compaction, in particular large scale systems with wide size distributions. In the model, a multi-grid searching method specifically designed within the GPU architecture was proposed for particle neighbour searching. A memory layout was designed to ensure coalesced memory access for neighbour list and associated contact history. The proposed GPU implementation was able to achieve a three-level parallelism, from single GPU to CPUs within a computing node and to GPUs across nodes. The model was applied to powder compaction and the simulation results showed significant gain in computational efficiency and reliable prediction of the compaction behaviour. (C) 2017 Elsevier B.V. All rights reserved.