<p>The efficiency of the discrete element method (DEM) of non-convex particles is restricted by the time-consuming multipoint-contact resolution due to complex surface morphologies. In this work, we first analyze the runtime occupancy of different phases in DEM computational loops by our previously developed Graphics Processing Unit (GPU)-accelerated framework, which adopts the grid method, the bounding sphere method, and the level set algorithm in the broad, middle and narrow phases of contact detection respectively. Next, to enhance computational efficiency, we propose a hierarchical optimization protocol that integrates two acceleration data structures into the middle and narrow phases of contact detection in a DEM computational loop. At the first level, an oriented bounding box (OBB) is employed in the middle phase to efficiently eliminate numerous non-overlapping neighbor pairs during the preliminary overlap detection. At the second level, a particle surface bounding volume hierarchy (PSBVH) is applied during the advanced contact resolution (narrow phase) to rapidly localize contact regions thereby avoiding exhaustive point-by-point traversal in the traditional level set algorithm. Subsequently, we apply the protocol to accelerate DEM simulations of gravitational packing and column collapse of non-convex particles and the optimization results are discussed in detail. The results show a maximum runtime reduction of 64.40% among numerical examples in this work, demonstrating substantial efficiency gains enabled by our two-level optimization strategy compared with the previously developed GPU framework. This efficiency improvement facilitates large-scale, high-resolution DEM simulations of actual non-convex granular materials, expanding their applicability in fields such as civil, hydraulic and geotechnical engineering.</p>