Boris D. Lubachevsky

Geometric properties of random disk packings

Random packings ofN⩽2000 rigid disks in the plane, subject to periodic boundary conditions on a square primitive cell, have been generated by a concurrent construction which treats all disks on an equal footing, as opposed to previously investigated sequential constructions. The particles start with random positions and velocities, and as they move about they grow uniformly in size, from points to jammed disks. The collection of packings displays several striking geometric features. These include (for largeN) typically polycrystalline textures with irregular grain boundaries and linear shear fractures. The packings occasionally contain monovacancies and trapped but unjammed “rattler” disks. The latter appear to be confined to the grain boundaries. The linear shear fractures preserve bond orientational order, but disrupt translational order, within the crystalline grains. A new efficient event-driven simulation algorithm is employed to generate the histories of colliding and jamming disks. On a computer which can process one million floating-point instructions per second the algorithm processes more than one million pairwise collisions per hour.

How to simulate billiards and similar systems

Abstract An N -component continuous-time dynamic system is considered whose components evolve independently all the time except for discrete asynchronous instances of pairwise interactions. Examples include colliding billiard balls and combat models. A new efficient serial event-driven algorithm is described for simulating such systems. Rather than maintaining and updating the global state of the system, the algorithm tries to examine only essential events, i.e., component interactions. The events are processed in a non-decreasing order of time; new interactions are scheduled on the basis of the examined interactions using preintegrated equations of evolutions of the components. If the components are distributed uniformly enough in the evolution space, so that this space can be subdivided into small sectors such that only O (1) sectors and O (1) components are in the neighborhood of a sector, then the algorithm spends time O (log N ) for processing an event which is the asymptotic minimum. The algorithm uses a simple strategy for handling data: only two states are maintained for each simulated component. Fast data access in this strategy assures the practical efficiency of the algorithm. It works noticeably faster than other algorithms proposed for this model.

Dense packings of congruent circles in a circle

Abstract The problem of finding packings of congruent circles in a circle, or, equivalently, of spreading points in a circle, is considered. Two packing algorithms are discussed, and the best packings found of up to 65 circles are presented.

Disks vs. Spheres: Contrasting Properties of Random Packings

Collections of random packings of rigid disks and spheres have been generated by computer using a previously described concurrent algorithm. Particles begin as infinitesimal moving points, grow in size at a uniform rate, undergo energy-onconserving collisions, and eventually jam up. Periodic boundary conditions apply, and various numbers of particles have been considered (N⩽2000 for disks,N⩽8000 for spheres). The irregular disk packings thus formed are clearly polycrystalline with mean grain size dependent upon particle growth rate. By contrast, the sphere packings show a homogeneously amorphous texture substantially devoid of crystalline grains. This distinction strongly influences the respective results for packing pair correlation functions and for the distributions of particles by contact number. Rapidly grown disk packings display occasional vacancies within the crystalline grains; no comparable voids of such distinctive size have been found in the random sphere packings. “Rattler” particles free to move locally but imprisoned by jammed neighbors occur in both the disk and sphere packings.

An analysis of rollback-based simulation

We present and analyze a general model of rollback in parallel processing, The analysis points out three possible modes where rollback may become excessive; we provide an example of each type. We identify the parameters that determme a stability, or efficiency region for the simulation. Our analysis suggests the possibility of a dangerous “phase-transition” from stabil ity to instability y in the parameter space. In particular, a rollback algorlthm may work efficiently for a small system but become inefficient for a large system. Moreover, for a given system, it may work quickly for a while and then suddenly slow down On the positive side, we give a tunable algorlthm, Filtered Rollback, that is designed to avoid the failure modes, Under appropriate assumptions, we provide a rigorous mathematical proof that Faltered Rollback m efficient, if implemented on a reasonably efficient multiprocessor. In particular, we show that the average time r to complete the simulation of a system with N nodes and R events on a p-processor PRAM satisfies

Algorithms for unboundedly parallel simulations

New methods are presented for parallel simulation of discrete event systems that, when applicable, can usefully employ a number of processors much larger than the number of objects in the system being simulated, Abandoning the distributed event list approach, the simulation problem is posed using recurrence relations. We bring three algorithmic ideas to bear on parallel simulation: parallel prefix computation, parallel merging, and iterative folding. Efficient parallel simulations are given for (in turn) the G/G/l queue, a variety of queueing networks having a global first come first served structure (e.g., a series of queues with finite buffers), acyclic networks of queues, and networks of queues with feedbacks and cycles. In particular, the problem of simulating the arrival and departure times for the first N jobs to a single G/G/l queue is solved in time proportional to N/P + log P using P processors.

A chaotic asynchronous algorithm for computing the fixed point of a nonnegative matrix of unit spectral radius

Given a nonnegative, irreducible matrix P of spectral radius unity, there exists a positive vector π such that π = πP. If P also happens to be stochastic, then π gives the stationary distribution of the Markov chain that has state-transition probabilities given by the elements of P. This paper gives an algorithm for computing π that is particularly well suited for parallel processing. The main attraction of our algorithm is that the timing and sequencing restrictions on individual processors are almost entirely eliminated and, consequently, the necessary coordination between processors is negligible and the enforced idle time is also negligible. Under certain mild and easily satisfied restrictions on P and on the implementation of the algorithm, x(.), the vectors of computed values are proved to converge to within a positive, finite constant of proportionality of π. It is also proved that a natural measure of the projective distance of x(.) from π vanishes geometrically fast, and at a rate for which a lower bound is given. We have conducted extensive experiments on random matrices P, and the results show that the improvement over the parallel implementation of the synchronous version of the algorithm is substantial, sometimes exceeding the synchronization penalty to which the latter is always subject.

Efficient parallel simulations of dynamic Ising spin systems

Abstract An event-driven method for parallel simulation of a class dynamic Monte Carlo models is presented. The method can be applied to several models studied in the computational physics such as Ising spin simulations by the method of Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller, continuous time Ising spin simulation by Glauber, and the dynamic binary alloy simulation. Unlike previously known parallel multi-spin algorithms, the proposed algorithms do not change the simulated model. For example, the asynchrony and randomness of update time arrivals which are present in Glaubers formulation, are not disturbed here and the simulated update history is precisely the same as it is in the serial algorithm. The theoretical efficiency evaluation is encouraging: for 768 × 768 spins using a parallel processor with 256 processing elements, the estimated efficiency is not lower than 71%. This means a parallel speed-up of 180 in the computations which were previously believed inherently serial. The algorithm by Bortz, Kalos, and Lebowitz can be incorporated, further contributing to speed-up.

Synchronization barrier and related tools for shared memory parallel programming

The synchronization barrier is a point in the program where the processing elements (PEs) wait until all the PEs have arrived at this point. In a reduction computation, given a commutative and associative binary operationop, one needs to reduce valuesa0,...,aN-1, stored in PEs 0,...,N-1 to a single valuea*=a0op a, op...op aN-1 and then to broadcast the resulta* to all PEs. This computation is often followed by a synchronization barrier. Routines to perform these functions are frequently required in parallel programs. Simple and efficient, workingC-language routines for the parallel barrier synchronization and reduction computations are presented. The codes are appropriate for a CREW (concurrent-read-exclusive-write) or EREW parallel random access shared memory MIMD computer. They require only shared memory read and write; no locks, semaphores etc. are needed. The running time of each of these routines isO(logN). The amount of shared memory required and the number of shared memory accesses generated are botO(N). These are the asymptotically minimum values for the three parameters. The algorithms employ the obvious computational scheme involving a binary tree. Examples of applications for these routines and results of performance testing on the Sequent Balance 21000 computer are presented.

Synchronous relaxation for parallel simulations with applications to circuit-switched networks

Synchronous relaxation, a new, general-purpose, efficient method for parallel simulation, is proposed. The method is applied to obtain a new parallel algorithm for simulating large circuit-switched communication networks. To show that synchronous-relaxation method is efficient, we present the results of circuit-switched network simulation experiments, and analytic approximations derived from a mathematical model of the simulation method.