Lubachevsky–Stillinger algorithm: Difference between revisions

Revision as of 12:53, 4 May 2013

Latest revision as of 13:31, 7 March 2024

(24 intermediate revisions by 21 users not shown)

Line 1:

{{Short description|Computational physics simulation algorithm}}

'''Lubachevsky-Stillinger (compression) algorithm''' (LS algorithm, LSA, or LS protocol) is a numerical procedure suggested by [[F. H. Stillinger]] and Boris D. Lubachevsky that simulates or imitates a physical process of compressing an assembly of hard particles.<ref name="StillingerLubachevskyJStat">{{cite journal|doi=10.1007/bf01025983|url=https://www.princeton.edu/~fhs/geodisk/geodisk.pdf|title=Geometric properties of random disk packings|journal=Journal of Statistical Physics|volume=60|issue=5–6|pages=561–583|year=1990|last1=Lubachevsky|first1=Boris D.|last2=Stillinger|first2=Frank H.|bibcode=1990JSP....60..561L|s2cid=15485746 }}</ref> As the LSA may need thousands of arithmetic operations even for a few particles, it is usually carried out on a [[digital computer]].[[File:1000 triangles packed in rectangle.png|thumb|Using a variant of Lubachevsky-Stillinger algorithm, 1000 congruent isosceles triangles are randomly packed by compression in a rectangle with periodic (wrap-around) boundary. The rectangle which is the period of pattern repetition in both directions is shown. Packing density is 0.8776]]

==Phenomenology (what is being simulated)==

A physical process of compression often involves a contracting hard boundary of the container, such as a piston pressing against the particles. The LSA is able to simulate such a scenario.<ref>F.{{cite Hjournal|doi=10. Stillinger and B. D. Lubachevsky, 1007/bf01054337|title=Crystalline-Amorphousamorphous Interfaceinterface Packingspackings for Disksdisks and Spheres,spheres|journal=Journal J.of Stat. Phys.Statistical Physics|volume=73,|issue=3–4|pages=497–514|year=1993|last1=Stillinger|first1=Frank 497-514H.|last2=Lubachevsky|first2=Boris (1993)D.|s2cid=59429012 }}</ref> However, the LSA was originally introduced in the setting without a hard boundary<ref name="StillingerLubachevskyJStat"/>B.<ref>{{cite Djournal|doi=10.1016/0021-9991(91)90222-7|title=How Lubachevskyto simulate billiards and F.similar Hsystems|journal=Journal of Computational Physics|volume=94|issue=2|pages=255–283|year=1991|last1=Lubachevsky|first1=Boris D.|bibcode=1991JCoPh..94..255L|arxiv=cond-mat/0503627|s2cid=6215418 Stillinger,}}</ref> Geometricwhere propertiesthe ofvirtual randomparticles diskwere packings"swelling" or expanding in a fixed, Jfinite virtual volume with [[periodic boundary conditions]]. StatisticalThe Physicsabsolute 60sizes (1990),of 561the particles were increasing but particle-583to-particle http://www.princeton.edu/~fhs/geodisk/geodisk.pdf</ref><ref>B.Drelative sizes remained constant. LubachevskyIn general, Howthe toLSA Simulatecan Billiardshandle an external compression and Similaran Systemsinternal particle expansion, Journalboth ofoccurring Computationalsimultaneously Physicsand Volumepossibly, 94but Issuenot 2necessarily, Maycombined 1991with http://arxiv.org/PS_cache/cond-mat/pdf/0503/0503627v2a hard boundary.pdf</ref> whereIn addition, the virtualboundary particlescan werebe "swelling" or expandingmobile.

in a fixed, finite virtual volume with [[periodic boundary conditions]]. The absolute sizes of the particles were increasing but particle-to-particle relative sizes remained constant. In general, the LSA can handle an external compression and an internal particle expansion, both occurring simultaneously and possibly, but not necessarily, combined with a hard boundary. In addition, the boundary can be mobile.

In a final, compressed, or "jammed" state, some particles are not jammed, they are able to move within "cages" formed by their immobile, jammed neighbors and the hard boundary, if any. These free-to-move particles are not an artifact, or pre-designed, or target feature of the LSA, but rather a real phenomenon. The simulation revealed this phenomenon, somewhat unexpectedly for the authors of the LSA. Frank H. Stillinger coined the term "rattlers" for the free-to-move particles, because if one physically shakes a compressed bunch of hard particles, the rattlers will be rattling.

In the "pre-jammed" mode when the density of the configuration is low and when the particles are mobile, the compression and expansion can be stopped, if so desired. Then the LSA, in effect, would be simulating a [[granular flow]]. Various dynamics of the instantaneous collisions can be simulated such as: with or without a full restitution, with or without tangential friction. Differences in masses of the particles can be taken into account. It is also easy and sometimes proves useful to "fluidize" a jammed configuration, by decreasing the sizes of all or some of the particles. Another possible extension of the LSA is replacing the hard collision [[force]] [[potential]] (zero outside the particle, infinity at or inside) with a piece-wise constant [[force]] [[potential]]. The LSA thus modified would approximately simulate [[molecular dynamics]] with continuous

short range particle-particle force interaction. External [[force field (physics)|force fields]], such as [[gravitation]], can be also introduced, as long as the inter-collision motion of each particle can be represented by a simple one-step calculation.

Differences in masses of the particles can be taken into account. It is also easy and sometimes proves useful to "fluidize" a jammed configuration, by decreasing the sizes of all or some of the particles. Another possible extension of the LSA is replacing the hard collision [[force]] [[potential]] (zero outside the particle, infinity at or inside) with a piece-wise constant [[force]] [[potential]]. The LSA thus modified would approximately simulate [[molecular dynamics]] with continuous

short range particle-particle force interaction. External [[force fields]], such as [[gravitation]], can be also introduced, as long as the inter-collision motion of each particle can be represented by a simple one-step calculation.

Using LSA for spherical particles of different sizes and/or for jamming in a non-commeasureable size container proved to be a useful technique for generating and studying micro-structures formed under conditions of a [[crystallographic defect]]<ref>F.{{cite Hjournal|doi=10. Stillinger and B. D. Lubachevsky. 1007/bf02183698|title=Patterns of Brokenbroken Symmetrysymmetry in the Impurityimpurity-Perturbedperturbed Rigidrigid-disk Diskcrystal|journal=Journal Crystal,of J.Statistical Stat.Physics|volume=78|issue=3–4|pages=1011–1026|year=1995|last1=Stillinger|first1=Frank PhysH.|last2=Lubachevsky|first2=Boris D.|bibcode=1995JSP....78,.1011S|s2cid=55943037 1011-1026 (1995)}}</ref> or a [[geometrical frustration]]<ref>Boris{{cite Djournal|doi=10. Lubachevsky and Frank H1103/physreve.70. Stillinger, 041604|pmid=15600418|title=Epitaxial frustration in deposited packings of rigid disks and spheres. |journal=Physical Review E |volume=70:44, 41604 (|issue=4|pages=041604|year=2004)|last1=Lubachevsky|first1=Boris http://arxivD.org/PS_cache/|last2=Stillinger|first2=Frank H.|bibcode=2004PhRvE..70d1604L|arxiv=cond-mat/pdf/0405/0405650v5.pdf0405650|s2cid=1309789 }}</ref><ref>{{cite journal|last1=Lubachevsky|first1=Boris D. Lubachevsky, |last2=Graham|first2=Ron L. Graham, and |last3=Stillinger|first3=Frank H. Stillinger, |title=Spontaneous Patterns in Disk Packings. |journal=Visual Mathematics, |year=1995. |url=http://vismath5.tripod.com/lub/}}</ref> It should be added that the original LS protocol was designed primarily for spheres of same or different sizes.<ref>A.R.{{cite Kansal, Sjournal|doi=10. Torquato, and F1063/1.H. Stillinger, 1511510|title=Computer Generationgeneration of Densedense Polydispersepolydisperse Spheresphere Packings,packings|journal=The J.Journal Chem.of Phys.Chemical Physics|volume=117, 8212-8218 (|issue=18|pages=8212–8218|year=2002)|last1=Kansal|first1=Anuraag http://catR.inist.fr/?aModele|last2=afficheN&cpsidtTorquato|first2=13990882Salvatore|last3=Stillinger|first3=Frank H.|bibcode=2002JChPh.117.8212K}}</ref>

Any deviation from the spherical (or circular in two dimensions) shape, even a simplest one, when spheres are replaced with ellipsoids (or ellipses in two dimensions),<ref>A.{{cite Donev, Fjournal|doi=10.H1103/physrevlett. Stillinger, P92.M. Chaikin, and S. Torquato, 255506|pmid=15245027|title=Unusually Dense Crystal Packings of Ellipsoids,|journal=Physical Phys. Rev.Review Letters |volume=92, |issue=25|pages=255506|year=2004|last1=Donev|first1=Aleksandar|last2=Stillinger|first2=Frank (2004)H.|last3=Chaikin|first3=P. M.|last4=Torquato|first4=Salvatore|bibcode=2004PhRvL..92y5506D|arxiv=cond-mat/0403286|s2cid=7982407 }}</ref> causes thus modified LSA to slow down substantially.

But as long as the shape is spherical, the LSA is able to handle particle assemblies in tens to hundreds of thousands

on today's (2011) standard [[personal computers]]. Only a very limited experience was reported<ref>M.{{cite Skoge, Ajournal|doi=10. Donev, F1103/physreve.H74. Stillinger, and S. Torquato, 041127|pmid=17155042|title=Packing Hypersphereshyperspheres in Highhigh-Dimensionaldimensional Euclidean Spaces,spaces|journal=Physical Phys. Rev.Review E |volume=74, |issue=4|pages=041127|year=2006|last1=Skoge|first1=Monica|last2=Donev|first2=Aleksandar|last3=Stillinger|first3=Frank (2006)H.|last4=Torquato|first4=Salvatore|bibcode=2006PhRvE..74d1127S|arxiv=cond-mat/0608362|s2cid=18639889 }}</ref>

in using the LSA in dimensions higher than 3.

==Implementation (how the calculations are performed)==

The state of particle jamming is achieved via simulating a [[granular flow]]. The flow is rendered as a [[discrete event simulation]], the events being particle-particle or particle-boundary collisions. Ideally, the calculations should have been

performed with the infinite precision. Then the jamming would have occurred [[ad infinitum]]. In practice, the precision is finite as is the available resolution of representing the real numbers in the [[computer memory]], for example, a [[double-precision]] resolution. The real calculations are stopped when inter-collision runs of the non-rattler particles become

smaller than an explicitly or implicitly specified small threshold. For example, it is useless to continue the calculations when inter-collision runs are smaller than the roundoff error.

The LSA is efficient in the sense that the events are processed essentially in an [[Event-driven programming|event-driven]] fashion, rather than in a time-driven fashion. This means almost no calculation is wasted on computing or maintaining the positions and velocities of the particles between the collisions. Among the [[Event-driven programming|event-driven]] algorithms intended for the same task of simulating [[granular flow]], like, for example, the algorithm of D.C. Rapaport,<ref>{{cite journal|doi=10.1016/0021-9991(80)90104-7|title=The event scheduling problem in molecular dynamic simulation|journal=Journal of Computational Physics|volume=34|issue=2|pages=184–201|year=1980|last1=Rapaport|first1=D.C|bibcode=1980JCoPh..34..184R}}</ref> the LSA is distinguished by a simpler [[data structure]] and data handling.

The LSA is efficient in the sense that the events are processed essentially in an [[event-driven]] fashion, rather than in a

time-driven fashion. This means almost no calculation is wasted on computing or maintaining the positions and velocities

of the particles between the collisions. Among the [[event-driven]] algorithms intended for the same task of simulating [[granular flow]], like, for example, the algorithm of D.C. Rapaport,<ref>D.C. Rapaport, The Event Scheduling Problem in Molecular Dynamic Simulation, Journal of Computational Physics Volume 34 Issue 2, 1980</ref> the LSA is distinguished by a simpler [[data structure]] and data handling.

For any particle at any stage of calculations the LSA keeps record of only two events: an old, already processed committed event, which comprises the committed event [[time stamp]], the particle state (including position and velocity), and, perhaps, the "partner" which could be another particle or boundary identification, the one with which the particle collided in the past, and a new event proposed for a future processing with a similar set of parameters. The new event is not committed. The maximum of the committed old event times must never exceed the minimum of the non-committed new event times.

and a new event proposed for a future processing with a similar set of parameters. The new event is not committed. The maximum of the committed old event times must never exceed the minimum of the non-committed new event times.

Next particle to be examined by the algorithm has the current minimum of new event times. At examining the chosen particle,

Line 35 ⟶ 31:

As the calculations of the LSA progress, the collision rates of particles may and usually do increase. Still the LSA successfully approaches the jamming state as long as those rates remain comparable among all the particles, except for the rattlers. (Rattlers experience consistently low collision rates. This property allows one to detect rattlers.) However,

it is possible for a few particles, even just for a single particle, to experience a very high collision rate along the approach to a certain simulated time. The rate will be increasing without a bound in proportion to the rates of collisions in the rest of the particle ensemble. If this happens, then the simulation will be stuck in time, it won't be able to progress toward the state of jamming.

The stuck-in-time failure can also occur when simulating a granular flow, without the particle compression or expansion. This failure mode was recognized by the practitioners of granular flow simulations as an "inelastic collapse" <ref>S.{{cite McNamara, and Wjournal|doi=10.R1103/physreve. Young, 50.r28|pmid=9962022|title=Inelastic collapse in two dimensions, |journal=Physical Review E |volume=50: pp|issue=1|pages=R28–R31|year=1994|last1=McNamara|first1=Sean|last2=Young|first2=W. R28-R31, 1994R.|bibcode=1994PhRvE..50...28M}}</ref> because it often occurs in such simulations when the [[restitution coefficient]] atin collisions is low (and hence the collisions arei.e. inelastic). The failure is not specific to only the LSA algorithm. Techniques to avoid the failure have been proposed.<ref>{{cite thesis|last=Drozd|first=John J. Drozd, |title=Computer Simulation of Granular Matter: A Study of An Industrial Grinding Mill, Thesis, |publisher=Univ. Western Ontario, |___location=Canada, |year=2004|url=http://imaging. robarts.ca/~jdrozd/thesisjd.pdf|access-date=2011-05-25|url-status=dead|archive-url=https://web.archive.org/web/20110818102135/http://imaging.robarts.ca/~jdrozd/thesisjd.pdf|archive-date=2011-08-18 }}</ref>

== History ==

The LSA was a by-product of an attempt to find a fair measure of [[speedup]] in [[Parallel algorithm|parallel]] [[simulations]]. The [[Time Warp (algorithm)|Time Warp]] parallel simulation algorithm by David Jefferson was advanced as a method to simulate asynchronous spacialspatial interactions of fighting units in combat models on a [[parallel computer]].<ref>F. Wieland, and D. Jefferson, Case studies in serial and parallel simulations, Proc. 1989 Int'l Conf. Parallel Processing, Vol.III, F. Ris, and M. Kogge, Eds., pp. 255-258. </ref> Colliding particles models<ref>P. Hontales, B. Beckman, et al., Performance of the colliding pucks simulation of the Time Warp operating systems, Proc. 1989 SCS Multiconference, Simulation Series, SCS, Vol. 21, No. 2, pp. 3-7. </ref> offered similar simulation tasks with spacialspatial interactions of particles but clear of the details that are non-essential for exposing the simulation techniques. The [[speedup]] was presented as the ratio of the execution time on a [[uniprocessor]] over that on a [[multiprocessor]], when executing the same parallel [[Time Warp]] algorithm. Boris D. Lubachevsky noticed that such a speedup assessment might be faulty because executing a [[parallel algorithm]] for a task on a uniprocessor is not necessarily the fastest way to perform the task on such a machine. The LSA was created in an attempt to produce a faster uniprocessor simulation and hence to have a more fair assessment of the [[parallel speedup]]. Later on, a parallel simulation algorithm, different from the Time Warp, was also proposed, that, when run on a uniprocessor, reduces to the LSA.<ref>{{cite journal|last=Lubachevsky|first=B.D.|title=Simulating Billiards: Serially and in Parallel|journal=International Journal in Computer Simulation|volume=2|year=1992|pages=373–411}}</ref>

different from the Time Warp, was also proposed, that, when run on a uniprocessor, reduces to the LSA.<ref>B.D. Lubachevsky, Simulating Billiards: Serially and in Parallel, Int.J. in Computer Simulation, Vol. 2 (1992), pp. 373-411.</ref>

== References ==

Line 49 ⟶ 44:

==External links==

* [http://cims.nyu.edu/~donev/Packing/Animations.html LSA in action. A collection of animations by Alexander Donev]

* [httphttps://cherrypitcims.princetonnyu.edu/~donev/Packing/C++/ Source C++ codes of a version of the LSA in arbitrary dimensions]

* [http://people.physics.anu.edu.au/~tas110/Pubblicazioni/ReggioC_10a.pdf Volume fluctuation distribution in granular packs studied using the LSA]

* [httphttps://wwwweb.archive.org/web/20110203002410/http://pack-any-shape.com/ LSA generalized for particles of arbitrary shape]

* [http://onlinelibrary.wiley.com/doi/10.1002/pamm.200610180/pdf LSA used for production of representative volumes of microscale failures in packed granular materials]

<!--- Categories --->

{{DEFAULTSORT:Lubachevsky-Stillinger Algorithmalgorithm}}

[[Category:Computational physics]]

[[Category:Scientific modelingSimulation]]

[[Category:GranularGranularity of materials]]

[[Category:Articles created via the Article Wizard]]