<p>Accurate simulation of geomaterials requires 3D meso-mechanical random fields that couple rock particle morphology with internal heterogeneity—a goal hindered by methodological decoupling and high computational costs. To address this, we present an integrated, GPU‑accelerated framework for the closed‑loop generation of such coupled random fields for rock particles. The framework unifies 3D scanning, spherical harmonics, and random discrete element theory into a seamless workflow from physical digitization to random field realization. Its key innovations are (1) a scheme coupling stochastic morphology (via an anisotropic spherical harmonic representation) with meso-mechanical heterogeneity; (2) a grayscale image-based method for determining the scale of fluctuation (SOF); and (3) optimized GPU modules enabling remarkable speedups of 1714.2 × in key steps and 113.4 × overall. Statistical characterization reveals multi-scale controls: the coefficient of variation (COV) and SOF govern strength dispersion and spatial correlation, while macroscopic rock size and morphology modulate higher-order statistics and heterogeneous distribution. The framework is rigorously validated by its unique capability to replicate experimental particle crushing behavior, complex fragmentation patterns, and the Weibull-type statistical distribution of particle strength, thereby establishing a direct, quantitative link between meso-scale heterogeneity and macro-mechanical response. Critically, the framework reveals that material homogeneity enhances stiffness and peak strength at the cost of localized brittle collapse, whereas spatial heterogeneity promotes lower stiffness and peak strength but enables ductile, distributed failure—a fundamental insight linking meso‑scale randomness to macro‑scale ductility. This work provides a transformative numerical tool for the probabilistic analysis of granular materials, linking particle-scale stochasticity to engineering-scale performance in geotechnics and granular physics.</p>