Luding, Stefan

Ogarko, V.

English

We investigate the influence of the packing fraction of highly polydisperse particle systems on the performance of a high-performance multilevel contact detection algorithm as applied for molecular dynamics type simulations. For best performance, this algorithm requires two or more hierarchy levels in order to cope with the strongly different particle size classes. In order to predict the optimal number of levels, an empirical parameter corresponding to the "overhead" is identified. For homogeneous systems, the density is not much affecting the performance of the algorithm, however, the optimal number of levels slightly increases with density as well as the speed-up as compared to a single-level method

University of Twente Research Information

II International Conference on Particle-based Methods - Fundamentals and ApplicationsPARTICLES 2011E. Oñate and D.R.J. Owen (Eds)A STUDY ON THE INFLUENCE OF THE PARTICLEPACKING FRACTION ON THE PERFORMANCE OF AMULTILEVEL CONTACT DETECTION ALGORITHMV. OGARKO AND S. LUDINGMulti Scale Mechanics (MSM), CTW, UTwente,PO Box 217, 7500 AE Enschede, Netherlands;v.ogarko@utwente.nl, s.luding@utwente.nlKey words: contact detection, discrete element, polydisperse, particle size distribution,packing fractionAbstract. We investigate the influence of the packing fraction of highly polydisperseparticle systems on the performance of a high-performance multilevel contact detectionalgorithm as applied for molecular dynamics type simulations. For best performance,this algorithm requires two or more hierarchy levels in order to cope with the stronglydifferent particle size classes. In order to predict the optimal number of levels, an empiricalparameter corresponding to the “overhead” is identified. For homogeneous systems, thedensity is not much affecting the performance of the algorithm, however, the optimalnumber of levels slightly increases with density as well as the speed-up as compared to asingle-level method.1 INTRODUCTIONContact detection is a fundamental computational problem arising in computer sim-ulations of systems consisting of many discrete objects such as particles or atoms. Theparticle based modeling methods like the Discrete Element Method (DEM) [1] play an im-portant role for physics-based simulations in many diverse fields. The performance of thecomputation relies on several factors, which include the physical model, on the one hand,and the contact detection method used, on the other. The collision detection of short-range pairwise interactions between particles is usually one of the most time-consumingtasks in calculations [2].The most commonly used method for contact detection of nearly monosized particleswith short-ranged forces is the Linked-Cell method [3]. Nevertheless, the Linked-Cellmethod is unable to efficiently deal with particles of greatly varying sizes [4]. This caneffectively be addressed by the use of methods based on hierarchical grids [4, 5, 6].1V. Ogarko and S. LudingThe multilevel contact detection algorithm presented in Ref. [6] contains the numberof hierarchy levels and the cell sizes at each level as adjustable parameters. A methodto find the optimal parameters for an arbitrary polydisperse size distribution of objectsis also shown in [6]. As an input parameter, this method requires the “overhead” of thedata structure. In addition to the detailed study of different size distributions in Ref.[6], we examine the influence of the density of the particle system on the performanceof the algorithm [6] and on the “overhead” of the data structure by means of moleculardynamics simulations.2 CONTACT DETECTION ALGORITHMIn this section we briefly describe the main ideas of the multilevel contact detectionalgorithm and the method to find the optimal parameters, which are necessary for under-standing the results of this study. For more details see Ref. [6].2.1 AlgorithmThe algorithm is designed to determine all the pairs in a set of N spherical particles ina d-dimensional Euclidean space that are in close neighborhood and possibly do overlap.A hierarchial grid structure used in the algorithm is a set of L regular grids with differentcell sizes. Every regular grid is associated with a hierarchy level h ∈ [1, L], where L is theinteger number of hierarchy levels. Each level h has a different cell size sh ∈ R, where thecells are d-dimensional cubes. Grids are ordered with increasing cell size so that h = 1corresponds to the grid with smallest cell size, i.e. sh < sh+1. For a given number of levelsand cell sizes, the hierarchical grid cells are defined by the following spatial mapping ofpoints ~x ∈ Rd to a cell at specified level h:(~x, h) 7→ (bx1/shc , ..., bxd/shc , h), (1)where bxc denotes the floor function 1.The algorithm is made up of two phases. In the first “mapping phase” all the particlesare mapped into hierarchical grid cells based on their location and size as follows. Thediameter of the smallest particle within a given level of hierarchy should be larger thanthe cell size of the previous level of hierarchy, and the diameter of the largest particleat some hierarchical level should be smaller or equal to the cell size of this level. Thecell-sizes and the distribution of particles among levels are detailed below.In the second “contact detection phase” for every particle in the system the potentialcontact partners are determined, and the geometrical intersection tests with them aremade. The contact detection is split into two steps, and the search is done by loopingover all particles and performing the first and second steps consecutively for each particle.The first step is the contact search at the level of insertion of a given particle, usingthe classical Linked-Cell method [3]. The second step is the cross-level search. For a1The largest integer not greater than x2V. Ogarko and S. Ludinggiven particle, one searches for potential contacts at every level h lower than the level ofinsertion. For detailed description of these steps and further optimizations see Ref. [6].2.2 Optimal parametersIt is shown in Ref. [6] that in the case of polydisperse systems, i.e., systems whereall the particles’ sizes are different, the performance of the multilevel contact detectionalgorithm strongly depends on the selected parameters (the number of levels L and cellsizes sh). Furthermore, a hypothesis on the optimal distribution of particles by levels wasformulated in [6], which states:Hypothesis 1 Let mh be the average number of particles per cell at level h, that is,mh = Nh/Nch, where Nh is the number of particles at level h, and Nch is the number ofcells at this level. Then the optimal distribution of particles by levels satisfies the followingcondition:m := mi = mj, for all i, j ∈ [1, L]. (2)If the particles are mapped to the levels, such that Eq. (2) is approximated, the CPUtime spent for contact detection, TCD, scales asTCD ∼ NL(m + K), (3)where K is a constant corresponding to the “overhead” of the algorithm, i.e., the timespent to access cells to be tested, and m = m(L, PSD) is the number of particles percell (see Ref. [6]). To compute m for a given L and for the particle system at hand, oneneeds to choose cell sizes sh so that Eq. (2) is approximately satisfied. How to do thisis explained in Ref. [6]. In the next section we investigate the influence of the volumefraction of the particle system, i.e., the ratio between the volume of the particles and thevolume of the system, on the value of the “overhead” K.3 SIMULATION RESULTSWe use homogeneous and isotropic disordered systems of colliding elastic sphericalparticles in a cubical box with hard walls. The motion of particles is governed by Newton’ssecond law with a linear elastic contact force during overlap. Every particle undergoesonly translational motion (without rotation) and gravity is set to zero.3.1 Experimental setupThe systems are prepared in two stages. Starting from a random uniform distributionof points in a cubical box, the radius of the particles grows linearly with time. We usenon-overlapping spheres [7], with an Event-Driven code [8], whose growing rate conformsto and conserves a prescribed size distribution. Initial velocities are set randomly in orderto keep the system dynamic and random, for details see Ref. [8]. When the target volumefraction is reached, the growth process is stopped.3V. Ogarko and S. Luding(a) (b)Figure 1: Particle systems with N = 125001 with (a) volume fraction ν = 0.1, (b) volume fractionν = 0.4, with uniform volume distribution. The size ratio ω is 50. Colour is by relative size.With the final configuration from the first step, the simulation switches to the relaxationstage with “soft” particles, i.e., particles move according to interparticle forces [1]. Thelinear elastic normal contact force model is used [9], which leads to a certain contactduration. At the beginning of this stage the velocities of the particles are scaled in a waythat a collision between two of the smallest particles would reach an average maximumoverlap of one percent of their radius. We let the simulations run for a few collisions perparticle for equilibration before making the measurement of the performance, in order tohave contacts (overlaps) between particles in the system.We use uniform volume size distribution, i.e., the probability distribution of the volumesof the particles is constant. This type of size distribution implies that the particle systemcontains rather many small particles. We denote the size ratio between the radius of thelargest and the smallest particle in the system as ω.3.2 Numerical performance resultsWe measure the total CPU time for varying number of hierarchy levels L for particlesystems with different volume fractions. For this, we calculate the value of m for eachL ∈ [1, 50] using the method given in Ref. [6], and run simulations where the multilevel gridis used with L levels, and the cell sizes sh are computed in accordance with hypothesis(1). To present the total CPU time, we use the slowdown factor SF, that is the totalCPU time divided by the smallest CPU time for a given system. In Fig. 2 the results ofthis numerical experiment are shown for systems with uniform volume distribution, withN = 125001 and ω = 50. The analytical prediction (3) is also plotted with K = 0.2,4V. Ogarko and S. Ludingscaled in such a way that SF = 1 corresponds to the minimum of the right hand side ofEq. (3). Note that even though the prediction (3) is for CPU time spent only for contactdetection, TCD, the total CPU time for fixed number of particles also scales as TCD. Thisis because the CPU time spent in the force calculation and integration does not depend onthe grid parameters used. From the numerical experimental results shown in Fig. 2 it canbe seen that the prediction of Eq. (3) used with K = 0.2 can be used for determining theoptimum number of levels, and the value of K is not dependent on the packing fraction.We also confirm that the dependence on density is very weak for rather large packingfractions of 0.62 (data from Ref. [6] – not shown here). Furthermore, the optimal numberof levels, L∗, slightly increases with the volume fraction from L∗ = 5 for ν = 0.1 to L∗ = 8for ν = 0.62, see Fig. 3(a). Finally, the best speed-up over the Linked-Cell method(L = 1), S, increases nearly linearly with the volume fraction, as it can be seen from Fig.3(b). Note that even for relatively low volume fraction ν = 0.1, the speed-up is high, i.e.,S ≈ 150, so the use of the multilevel algorithm is highly advantageous for this type ofparticle size distributions (many small particles) for the whole range of volume fractions. 1 5 10 50 100 500 1  2  5  10  15  20  30  50SFLSimulations, ν = 0.1Prediction(a) 1 5 10 50 100 500 1  2  5  10  15  20  30  50SFLSimulations, ν = 0.4Prediction(b)Figure 2: Slowdown factor for different numbers of levels for systems with N = 125001 and ω = 50 withuniform volume distribution; (a) volume fraction ν = 0.1, (b) volume fraction ν = 0.4. The prediction ofEq. (3) is used with K = 0.2 for plotting solid lines. Note that the data can be obtained only for integervalues of L.5V. Ogarko and S. Luding 5 6 7 8 0.1  0.2  0.3  0.4  0.5  0.6  0.7L*ν(a) 140 150 160 170 180 190 200 210 220 230 0.1  0.2  0.3  0.4  0.5  0.6  0.7Sν(b)Figure 3: (a) The optimal number of levels for different densities. (b) The speed-up relative to thesingle-level contact search for different packing fractions. The data for ν = 0.62 is taken from Ref. [6].4 SUMMARY AND CONCLUSIONSWe have shown that the multilevel algorithm presented in Ref. [6] can be used tosignificantly improve the performance of simulations of particle systems with rather manysmall particles at very different volume fractions. The speed-up relative to the single-levelLinked-Cell method increases nearly linearly with the volume fraction for the homogeneoussystems tested. Moreover, the value of the “overhead” parameter, which is necessary toanalytically find the optimal number of hierarchy levels is not dependent on the packingfraction for these systems. This leads to an increasing number of levels with the densityfor optimal performance of the contact detection.ACKNOWLEDGMENTSThe authors would like to thank S. Gonzalez, A. R. Thornton and O. Imole for helpfuldiscussions. This research is supported by the Dutch Technology Foundation STW, whichis the applied science division of NWO, and the Technology Programme of the Ministryof Economic Affairs (STW MuST project 10120).REFERENCES[1] P. A. Cundall and O. D. L. Strack. A discrete numerical model for granular assem-blies. Géotechnique, 29 (1), 47-65 (1979).[2] J. R. Williams and R. O’Connor. Discrete element simulation and the contact prob-lem. Archives of Computational Methods in Engineering, 6 (4), 279-304 (1999).6V. Ogarko and S. Luding[3] M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids. Clarendon Press,1989.[4] T. Iwai, C. W. Hong, and P. Greil. Fast particle pair detection algorithms for particlesimulations. International Journal of Modern Physics C, 10 (5), 823-837 (1999).[5] B. Mirtich, Efficient Algorithms for Two-Phase Collision Detection in Practical Mo-tion Planning in Robotics: Current Approaches and Future Directions, Ed. K. Guptaand A. P. del Pobil, Wiley, New York, 1998.[6] V. Ogarko and S. Luding, A fast multilevel algorithm for contact detection of arbi-trarily polydisperse objects, submitted 2011.http://www2.msm.ctw.utwente.nl/sluding/PAPERS/HGrid.pdf[7] B. D. Lubachevsky and F. H. Stillinger. Geometric Properties of Random Disk Pack-ings. J. Stat. Phys., 60, 561-583 (1990).[8] M. Skoge, A. Donev, F. H. Stillinger and S. Torquato. Packing hyperspheres in high-dimensional Euclidean spaces. Phys. Rev. E, 74 (4), 041127 (2006).[9] S. Luding. Cohesive, frictional powders: contact models for tension. Granular Matter,10 (4), 235-246 (2008).7

A study on the influence of the particle packing fraction on the performance of a multilevel contact detection algorithm

https://ris.utwente.nl/ws/files/5384752/Ogarko11study.pdf

A study on the influence of the particle packing fraction on the performance of a multilevel contact detection algorithm

Abstract

Similar works

Full text

Available Versions

University of Twente Research Information