280 likes | 807 Views
An introduction to scripting DEM simulations. ESyS-Particle:. D. Weatherley Earth Systems Science Computational Centre University of Queensland. Outline. ESyS-Particle features Why Python? Scripting a DEM simulation Constructing a particle assembly Specifying interactions
E N D
An introduction to scripting DEM simulations ESyS-Particle: D. Weatherley Earth Systems Science Computational Centre University of Queensland
Outline • ESyS-Particle features • Why Python? • Scripting a DEM simulation • Constructing a particle assembly • Specifying interactions • Simulation initialisation & execution • Real-time processing • Visualisation
Features Software: • 3D C++ DEM engine • Verlet list neighbour search algorithm • Domain decomposition parallelism via MPI • Python scripting interface • Open source license
Features • Particle Physics: • Spheres & aggregates • Contact models: • Linear elastic • Friction • Rotational bonds • Thermal diffusion • Linear and tri-mesh walls • Gravity, bulk viscosity
Why Python? http://www/python.org/ • An interpretive scripting language • Nice syntax and easy to learn • Object-oriented design • Reusable packages and modules • Extensive standard and 3rd party modules • Email, web, XML, I/O, math etc. • Numarray, SciPy, matplotlib, ... • Existing C/C++/Fortran libraries can be “boosted” into Python easily http://www.boost.org/libs/python/doc/
Scripting a simulation 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
Constructing particle assemblies 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() )
Particle assemblies (2) 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() 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
Particle Assemblies (3) • Other assemblies: • Spherical aggregate grains • Cylinders via trimming of rectangular blocks • Limitations: • Lacks some of the tricks of PFC • Particle expansion • Packing with specified porosity • Random packing with specified radius size-distribution
Specifying Interactions • Interaction Groups specified via IGprms objects • Types of IGprms: • Body forces • Particle-particle • Particle-wall my_gravity = GravityPrms( name=”earth-gravity”, acceleration=Vec3(0,-9.81,0) ) wall_interaction = NRotElasticWallPrms( name = “SpringyWall”, wallName = “bottom”, normalK = 10000.0 # N/m ) particle_interactions = NRotElasticPrms( name = “SpringySpheres”, normalK = 10000.0 # N/m )
Bonded particle-wall interactions utilise particle “tags” to denote those particles to be bonded to walls Specifying Interactions (2) • Bonded particle interactions require a list of particle-pairs • A NeighbourSearcher class provides a mechanism to compute particle-pair lists bondGrp = \ sim.createInteractionGroup( NRotBondPrms( name = “SphereBonds”, normalK = 10000.0, # N/m breakDistance = 10*r ) ) neighFinder = SimpleSphereNeighbours(maxDist=0.01*r) idPairs = neighFinder.getNeighbours(particles) bondGrp.createInteractions(idPairs)
To initialise a simulation, we create an instance of the LsmMpi class and add the various features we need Simulation Initialisation • The LsmMpi class is a “container” for particle assemblies, walls, interactions etc. 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) sim.createParticles(my_particle_collection) sim.createInteractionGroup(my_gravity) sim.createInteractionGroup(particle_int) ...
Initialisation (2) & Execution sim.createWall( name = “bottom”, posn=Vec3(0,-20,0), normal=Vec3(0, 1, 0) ) sim.createInteractionGroup( NRotElasticWallPrms( name = “SpringyWall”, wallName = “bottom”, normalK = 10000.0 # newtons per meter ) ) ... sim.setTimeStepSize(0.000025) sim.setNumTimeSteps(600000) sim.run()
Particle/Wall FieldSavers provide mechanism to store useful scalar/vector data at regular intervals Particle savers include: Kinetic Energy Displ./Veloc. Field Wall savers include: Force on wall Wall displacement Real-time processing • A Runnable class permits users to script tasks to perform before/after each timestep • Runnables: • control wall displacement/force (e.g. Servos) • Output data to disk • Move non-inertial particles/walls etc.
Real-time processing (2) 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)
Visualisation 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( Vec3(10,15,5) ) camera.setPosn( Vec3(10,100,5) ) scene.render( offScreen=True, fileName="cubicBlock.png", size=[800,600] ) • Three Visualisation modules provide batch-mode visualisation for large/long runs • Visualisation of: • Particle displacements • Velocity fields • Force chains • Other ... • Wrapping visualisation tasks as Runnables results in a library of re-usable visualisation tools
Resources ESyS-Particle Twiki pages: https://shake200.esscc.uq.edu.au/twiki/bin/view/ESSCC/ParticleSimulation EsyS-Particle API Documentation: http://iservo.edu.au/esys/esys_particle_python_doc/esys_particle-1.x.x/pythonapi/html/ HELP! d.weatherley@uq.edu.au