## 3 Performance Evaluation
### Triaxial Test
The loading velocity remained at 1 cm/s. The stress-strain curves obtained with 31,360 and 250,880 particles present better convergence and stability (Figure 5b). The average runtime of the CPU sequential simulation in 10,000 incremental steps was investigated for a comparison with the CPU and GPU parallel simulations. For the CPU sequential calculations, the runtime is nearly linearly proportional to the number of particles (Table 1). The CPU parallel simulations were also performed using the OpenMP with up to 10 cores. The speedups of the CPU parallel simulations increase with the number of CPU cores (see Table 1 and Figure 7a). The maximum speedup with 10 threads is 5.7. With more CPU threads mobilised, the acceleration efficiency wanes, which is due to the overhead on the scheduling manipulations. The significant reduction in runtimes by the GPU parallel strategy is clear. The speedup increases to less than 31,360 particles; if the computational scale is enlarged further, the GPU runs with a full load, and the improvement presents a good scaling behaviour and gradually stabilises to an average speedup of 84 (Figure 7b). The total runtime for the whole simulation process within an axial strain of 20% is proportional to the square of \(D/d_{50}\), e.g., 8.3 h for \(D/d_{50}=\) 50 (31,360 particles) and 18.75 h for \(D/d_{50}=\) 75 (100,880 particles), and that for \(D/d_{50}=\) 158 (1,003,520 particles) is estimated as 75 h. However, for the real Fujian sand with \(d_{10}=\) 0.32 mm, \(d_{50}=\) 0.43 mm and \(d_{90}=\) 0.67 mm, the runtime is further increased by the leap of an incremental step. Therefore, the total runtimes for a standard sample with \(H=\) 8 cm and \(D=\) 4 cm accommodating 1,650,000 real sand grains is extrapolated to be more than 92 days, which remains unacceptable for quantitative analysis of the real sand. Figure 4: Schematic for discrete element model of the triaxial test. (**a**) specimen setup (**b**) particle size distribution. Therefore, compromising techniques, such as decreasing the particles' stiffness, expanding the particle sizes and mass scaling, may remain necessary. The performance of the GPU parallel computing varies between the two computational procedures, which are listed in Table 2. The efficiency of the procedure \"Update kinematic state of particles\" is enhanced tremendously with speedups of more than 150. In contrast, the acceleration effect of the \"Calculate contact forces\" procedure is much less, mainly due to the low arithmetic intensity in the operations. The overhead on updating the contact pair list accounts for about 15% of that on the whole \"Calculate contact forces\" procedure. The acceleration effect of the GPU parallelisation can be influenced by various factors, e.g., the size of the coarse cells and the number of empty cells. To accommodate all the potential connecting particles, the coarse mesh needs to be larger than \(2R^{\text{max}}\). Given that a larger coarse mesh size is used, the cell number reduces and the particle number in each cell accordingly increases (Figure 8). The speedup increases from 82 to 123 with the coarse mesh size varying from \(2R^{\text{max}}\) (0.005 m) to \(4R^{\text{max}}\) (0.01 m), and then its value decreases with a coarse mesh larger than \(4R^{\text{max}}\). For consistency purposes, the coarse mesh size remained as \(2R^{\text{max}}\) in the following simulations. The speedups can be undermined by the existence of empty cells, which means some GPU threads are idle in the calculations. The influence of the empty cells can be roughly estimated by the ratio of the number of empty cells to engaged ones. Figure 5: Deviator stress changes versus axial strain during triaxial tests with confining stress of 300 kPa. (**a**) with different tools (**b**) influence of model size. Figure 6: Deformation and mechanical profiles of sample at axial strain of 6.2% (3920 particles). (**a**) transient velocity (**b**) accumulative displacement.