Scripting Parallel Discrete Element Simulations with ESyS_Particle

1. Shane Latham Australian Computational Earth Systems Simulator Earth Systems Science Computational Centre/QUAKES University of Queensland Brisbane, Australia Scripting Parallel Discrete Element Simulations with ESyS_Particle

2. ESyS_Particle discrete element simulation: Particles and interactions Parallelisation Python Scripting of Simulations Initial particle configuration Model parameterisation Running simulation Visualisation Data output Overview

3. C++ DEM implementation with Python scripting interface Neighbour lists and Verlet-distance for interaction/contact detection Three dimensional Particles/elements are currently spheres Parallelisation uses Message Passing Interface (MPI) ESyS_Particle Simulation (LSM)

4. Spherical particles Selection of contact physics: Non-rotational and rotational dynamics Friction interactions Linear elastic interactions Bonded interactions Particles and Interactions

5. Model elastic continuum as linear elastically bonded particles Non-rotational bonds only impart compression/tension forces Rotational bonds impart additional shear, torsion and bending forces Bonds can break – simulate fracturing Bonded Particle Interactions

6. Decompose rectangular domain into rectangular subdomains Parallel Implementation For 12 processes/CPUs divide domain into 3x4 grid (or 2x6 or 1x12)

7. Particles and interactions divided among separate MPI processes Parallel Implementation Distribute particles in each subdomain to individual MPI processes

8. Identify particles on boundary of the subdomains. Parallel Implementation

9. Duplicate boundary particles on neighbouring MPI processes Parallel Implementation

10. Why Python (http://www.python.org)? Python is “...a tool engineered with good balance, the right combination of rigor and flexibility...” A simple but effective approach to object-oriented programming High-level concepts such as collections and iterators High-level encapsulation facilities (packages and modules) to support the design of re-usable libraries Exception-handling for effective management of error conditions Extensive standard library (serialization, XML parsing, process management, file i/o, sets,...) Extensive range of third party libraries/applications with Python interfaces (VTK, MayaVi, numarray, SciPy,...) Interpreter is easily extended with new functions and data types implemented in C or C++ using Boost.Python (http://www.boost.org/libs/python/doc/index.html) Python Scripting ESyS_Particle Simulations

11. The esys.lsm Python package can be used in conjunction with existing Python modules for: Model initialisation geometrical configuration of particles configuration of inter-particle bonds Running simulations extract model data during simulation dynamic setting of model parameters, eg. creation of particles and bonds Visualisation Interactive visualisation for model debugging Batch option for large data sets and offline simulations Python Scripting ESyS_Particle Simulations

12. Create initial configuration of particles Setup boundary conditions and particle interactions Run simulation and generate output Perform online/offline visualisation Typical Simulation Steps Use basic example scripts to illustrate these steps.

13. Initial particle configuration from esys.lsm.geometry import SimpleSphere from esys.lsm.util import Vec3 r = 2.0 ssList = [] sphereId = 0 for i in xrange(0, 10): for j in xrange(0, 7): for k in xrange(0, 5): ssList.append( SimpleSphere( id = sphereId, centre = Vec3(i*2*r, j*2*r, k*2*r), radius = r ) ) sphereId += 1 # Save the particle config in a custom format: f = file(“RegularGrid.txt”, “w”) for ss in ssList: f.write( “%d %f %f %f %f\n” % ((ss.getId(),) + tuple(ss.getCentre()) + (ss.getRadius(),)) ) f.close() SimpleSphere objects will be used to initialise particles within the LSM model. Integer which is used to uniquely identify a particle

14. Initial particle configuration from esys.lsm.geometry import CubicBlock, HexagBlock from esys.lsm.util import Vec3 # Create a Face Centred Cubic (FCC) # rectangular block fcc = CubicBlock([20,15,10], radius=1.0) # Translate all particles in block fcc.translate(Vec3(10, 5, 10)) # Create a Hexagonal Centred Cubic (HCC) # rectangular block hcc = HexagBlock([40,30,20], radius=0.5) # Rotate all particles in block about centre # of block from math import pi hcc.rotate( Vec3(pi/2, pi/4, pi/8), hcc.getParticleBBox().getCentre() )

15. Initial particle configuration from esys.lsm.geometry import RandomBoxPacker from esys.lsm.util import Vec3, BoundingBox # Create randomly-sized spheres inside # a rectangular block rndPacker = \ RandomBoxPacker( minRadius = 0.1, maxRadius = 1.0, cubicPackRadius = 1.0, bBox = BoundingBox(Vec3(0,0,0), Vec3(20,15,10)), maxInsertFails = 10000 ) rndPacker.generate() rnd = rndPacker.getSimpleSphereCollection() # Use pickling to save the configuration # to file import pickle pickle.dump(rnd, file(“RndPacking.pkl”, “w”)) # Load configuration back in at some later stage rnd = pickle.load(file(“RndPacking.pkl”, “r”)) Parameters for creating a random packing of SimpleSphere objects inside a box Generate the packing, uses a random-location with fit-to-neighbours algorithm Assign the collection of SimpleSphere objects

16. (2-3) Set up and run simulation from esys.lsm import * from esys.lsm.geometry import CubicBlock from esys.lsm.util import Vec3, BoundingBox sim = LsmMpi(numWorkerProcesses=4, mpiDimList=[2,2,1]) sim.initVerletModel( particleType = “RotSphere”, gridSpacing = 2.5, verletDist = 0.5 ) domain = BoundingBox(Vec3(-20,-20,-20), Vec3(20, 20, 20)) sim.setSpatialDomain(domain) particles = CubicBlock([20,20,20], radius=0.5) sim.createParticles(particles) # Create gravity acting on particles. sim.createInteractionGroup( GravityPrms(name=”earth-gravity”, acceleration=Vec3(0,-9.81,0)) ) sim.setTimeStepSize(0.000025) sim.setNumTimeSteps(600000) sim.run() Create the (empty) model object, including MPI worker processes. Set the domain in which particles may exist. Create rotational spheres in the model from collection of SimpleSphere objects Create gravity body force in the model which acts on all particles. user@lam$ mpiexec -machinefile=mf.txt -np 1 /opt/lsm/bin/mpipython grav.py user@mpt$ mpirun -up 5 -np 1 /opt/lsm/bin/mpipython grav.py GravSim.mpg

17. (2-3) Set up and run simulation from esys.lsm import * from esys.lsm.geometry import CubicBlock from esys.lsm.util import Vec3, BoundingBox . . . sim.createInteractionGroup( GravityPrms(name=”earth-gravity”, acceleration=Vec3(0,-9.8,0)) ) # Create a “wall”, that is, an infinite plane sim.createWall( name = “bottom”, posn=Vec3(0,-20,0), normal=Vec3(0, 1, 0) ) # Create interactions between wall and particles sim.createInteractionGroup( NRotElasticWallPrms( name = “SpringyWall”, wallName = “bottom”, normalK = 10000.0 # newtons per meter ) ) . . . sim.run() Name is used to identify/reference the wall Name for the group of interactions. Specify that interactions occur between particles and “bottom” wall. GravSimBotWall.mpg

18. (2-3) Set up and run simulation from esys.lsm import * from esys.lsm.geometry import CubicBlock from esys.lsm.util import Vec3, BoundingBox . . . sim.createInteractionGroup( NRotElasticWallPrms( name = “SpringyWall”, wallName = “bottom”, normalK = 10000.0 # newtons per meter ) ) # Create interactions between particles sim.createInteractionGroup( NRotElasticPrms( name = “SpringySpheres”, normalK = 10000.0 # newtons per meter ) ) . . . sim.run() Create linear elastic interactions between particles which come into contact. GravSimSpringySpheres.mpg

19. (2-3) Set up and run simulation from esys.lsm import * from esys.lsm.geometry import CubicBlock from esys.lsm.util import Vec3, BoundingBox . . . # Create bonds between particles bondGrp = \ sim.createInteractionGroup( NRotBondPrms( name = “SphereBondage”, normalK = 10000.0, # N/m breakDistance = 10*r ) ) neighFinder = SimpleSphereNeighbours(maxDist=0.01*r) idPairs = neighFinder.getNeighbours(particles) bondGrp.createInteractions(idPairs) . . . sim.run() Linear elastic bonds can be created between specified particles. Determine pairs of particles (pairs of particle id's) which are less than a specified distance appart. Create the bonds between the specified pairs of particles. GravSimBondedSpheres.mpg

20. (3) Run and Save Particle Data to File from esys.lsm import * import cPickle as pickle class Saver(Runnable): def __init__(self, sim, interval): Runnable.__init__(self) self.sim = sim self.interval = interval self.count = 0 def run(self): if ((self.sim.getTimeStep() % self.interval) == 0): pList = self.sim.getParticleList() fName = “particleList%04d.pkl” % self.count pickle.dump(pList, file(fName, “w”)) self.count += 1 . . . saver = Saver(sim=sim, interval=100) sim.addPostTimeStepRunnable(saver) . . . pList = pickle.load(file(“particleList0000.pkl”,”r”)) for particle in pList: # do something with particle Create a class which inherits from esys.lsm.Runnable. Override the run method to perform an action during each time-step of a simulation.

21. (4) Visualisation from esys.lsm.vis import vtk, povray from esys.lsm.util import Vec3 from esys.lsm.geometry import CubicBlock #pkg = vtk pkg = povray sphereBlock = CubicBlock([20,30,10], radius=0.25) scene = pkg.Scene() for ss in sphereBlock: ss.add( pkg.Sphere(ss.getCentre(), ss.getRadius()) ) camera = scene.getCamera() camera.setLookAt(sphereBlock.getParticleBBox().getCentre()) camera.setPosn(sphereBlock.getParticleBBox().getCentre() + Vec3(0,0,8)) scene.render(offScreen=False, fileName="cubicBlock.png", size=[800,600]) Select rendering package: (1) VTK: Allows user interaction with data (2) Povray: High quality ray-traced rendering Add sphere objects to the scene. Setup the camera position and point it at the centre of the block of spheres.

22. Python interface for executing LSM simulations Minimal user exposure to parallelisation – only specify number of processes and use the mpirun command Separation of initial particle configuration from running simulation Flexibity of choosing between different boundary conditions and particle contact physics Basic visualisation support Conclusion