Running code with 4 threads

Asked by Carine Tanissa

How can i edit this code to run it in parallel using 4 threads:

from yade import pack, plot
import math

# 7-15. The code block is trying to create a directory named "vtk" if it doesn't exist already. This directory is used for storing VTK files.

import os
import errno
try:
   os.mkdir('./vtk/')
except OSError as exc:
   if exc.errno != errno.EEXIST:
      raise
   pass

#Enable the storage of potential particles in the simulation.

Gl1_PotentialParticle().store=True

O.engines=[
 ForceResetter(),
 InsertionSortCollider([PotentialParticle2AABB()],verletDist=0.01, avoidSelfInteractionMask=2),
 InteractionLoop(
  [Ig2_PP_PP_ScGeom(twoDimension=False, unitWidth2D=1.0, calContactArea=True, areaStep=5)],
  [Ip2_FrictMat_FrictMat_KnKsPhys(kn_i=1e8, ks_i=1e7, Knormal = 1e8, Kshear = 1e7, useFaceProperties=False, viscousDamping=0.05)],
  [Law2_SCG_KnKsPhys_KnKsLaw(label='law',neverErase=False)]
 ),
 NewtonIntegrator(damping=0.1,exactAsphericalRot=True,gravity=[0,-9.81,0]),
    # call the checkUnbalanced function (defined below) every 2 seconds
    PyRunner(command='checkUnbalanced()', realPeriod=2),
    # call the addPlotData function every 200 steps
    PyRunner(command='addPlotData()', iterPeriod=100),
    PotentialParticleVTKRecorder(fileName='./vtk/cubePPscaled',label='vtkRecorder',twoDimension=False,iterPeriod=1000,sampleX=50,sampleY=50,sampleZ=50,maxDimension=0.2)
]

#Define the density of the powder material.
powderDensity = 2000
#Define the distance to the center for particles.
distanceToCentre= 0.5
# Define the mean size of particles.
meanSize = 1.
#Calculate the wall thickness based on the mean size.
wallThickness = 0.5*meanSize
#this material represents frictionless particles and sets some parameters such as density.
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(0.0),density=powderDensity,label='frictionless')) #The normal and shear stifness values are determined in the IPhys functor, thus the young, poisson parameters of the FrictMat are not used.
lengthOfBase = 9.0*meanSize
heightOfBase = 14.0*meanSize
# create empty sphere packing
# sphere packing is not equivalent to particles in simulation, it contains only the pure geometry
sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),7.0*2*heightOfBase,0.5*(lengthOfBase-wallThickness))
R=sqrt(3.0)*distanceToCentre
sp.makeCloud(mn,mx,R,0,700,False)

#initialize the counter variable
count= 0
#Calculate the radius of particles based on the mean size.
r=0.05*meanSize

for s in sp:
 b=Body()
 b.mask=1
 b.aspherical=True
 wire=False
 color=Vector3(random.random(),random.random(),random.random())
 highlight=False
 b.shape=PotentialParticle(k=0.2, r=r, R=R, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r], isBoundary=False, color=color, wire=wire, highlight=highlight, minAabb=sqrt(3)*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), maxAabb=sqrt(3)*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), maxAabbRotated=1.02*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), minAabbRotated=1.02*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), AabbMinMax=True, id=count)
 V=(2*distanceToCentre)**3 # (approximate) Volume of cuboid
 geomInert=(1./6.)*V*(2*distanceToCentre)**2 # (approximate) Principal inertia of cuboid to its centroid
 utils._commonBodySetup(b, V, Vector3(geomInert,geomInert,geomInert), material='frictionless', pos=[0,0,0], fixed=False)
 b.state.pos = s[0] #s[0] stores center
 b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
 O.bodies.append(b)
 count=count+1
# for s in sp:
     #creates a new instance of the Body class to represent a particle.
 # b=Body()
     #Assigns a mask value of 1 to the body.
 # b.mask=1
     # Sets the body as aspherical, meaning it has a non-spherical shape.
 # b.aspherical=True
     #Specifies whether the body is displayed as a wireframe or not
 # wire=False
     # Generates a random RGB color for the body.
 # color=Vector3(random.random(),random.random(),random.random())
     #highlight = False: Specifies whether the body should be highlighted or not
 # highlight=False
     # Defines the shape of the body using a PotentialParticle shape. It specifies various parameters such as radius (r), outer radius (R), shape coefficients (a, b, c, d), color, and AABB (Axis-Aligned Bounding Box) dimensions.
     #k=0.2: The spring constant (k) determines the stiffness of the interactions between particles. Higher values result in stronger interactions.
     #r=r: The r parameter represents the radius of the sphere.
     #R=R: The R parameter represents the outer radius of the sphere.
     #a=[1, -1, 0, 0, 0, 0]: The a parameter defines the shape coefficients for the particle. These coefficients control the particle's non-sphericity by modifying its shape in terms of the Legendre polynomials.
     #b=[0, 0, 1, -1, 0, 0]: The b parameter represents additional shape coefficients, similar to a, but affecting different Legendre polynomials.
     #c=[0, 0, 0, 0, 1, -1]: The c parameter represents yet another set of shape coefficients affecting different Legendre polynomials.
     #d=[distanceToCentre - r, distanceToCentre - r, distanceToCentre - r, distanceToCentre - r, distanceToCentre - r, distanceToCentre - r]: The d parameter specifies the distance between the particle center and the surface of the particle in different directions. It is based on the distanceToCentre variable.
     #isBoundary=False: This parameter determines if the particle is considered a boundary particle or not.
     #color=color: The color parameter specifies the RGB color of the particle.
     #wire=wire: Specifies whether the particle should be displayed as a wireframe or not.
     #highlight=highlight: Determines whether the particle should be highlighted or not.
     #minAabb=sqrt(3) * Vector3(distanceToCentre, distanceToCentre, distanceToCentre): The minAabb parameter represents the minimum axis-aligned bounding box (AABB) coordinates for the particle. It is based on the distanceToCentre variable.
     #maxAabb=sqrt(3) * Vector3(distanceToCentre, distanceToCentre, distanceToCentre): The maxAabb parameter represents the maximum AABB coordinates for the particle. It is also based on the distanceToCentre variable.# Calculates the approximate volume of the cuboid based on distanceToCentre
 # b.shape=PotentialParticle(k=0.2, r=r, R=R, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r], isBoundary=False, color=color, wire=wire, highlight=highlight, minAabb=sqrt(3)*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), maxAabb=sqrt(3)*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), maxAabbRotated=1.02*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), minAabbRotated=1.02*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), AabbMinMax=True, id=count)
  # V=(2*distanceToCentre)**3 # (approximate) Volume of cuboid
  #geomInert=(1./6.)*V*(2*distanceToCentre)**2 # (approximate) Principal inertia of cuboid to its centroid
     #Sets up common properties of the body, such as volume (V), geometric inertia (geomInert), material, position, and fixity
  #utils._commonBodySetup(b, V, Vector3(geomInert,geomInert,geomInert), material='frictionless', pos=[0,0,0], fixed=False)
     # This line assigns the position of the body (b) to the value stored in s[0]. Here, s[0] represents the center position of the sphere.
  #b.state.pos = s[0] #s[0] stores center
     # This line sets the orientation of the body (b) using a Quaternion. The quaternion is initialized with random values for the quaternion components and a random rotation angle.
  #b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
     # Here, the body (b) is added to the O.bodies list, which represents the collection of bodies in the simulation.
  #O.bodies.append(b)
     # This line increments the count variable by 1, which is used as an identifier for the bodies
  #count=count+1

# create rectangular box from facets
# O.bodies.append() adds particle to the simulation; returns id of the particle() added
# create rectangular box from facets
# geom.facetBox(center, extents, orientation=Quaternion((1, 0, 0), 0), wallMask=63, **kw)
# wallMask determines which wall sides will be created
box=geom.facetBox((0, 49.25, 0), (20, 49.25, 20), wallMask=31)
O.bodies.append(box)

#Bottom faces of the box
r=0.1*wallThickness
bbb=Body()
bbb.mask=3
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialParticle(k=0.1, r=0.1*wallThickness, R=0.2*lengthOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[lengthOfBase/2.0-r,lengthOfBase/2.0-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/2.0-r,lengthOfBase/2.0-r], id=count,isBoundary=True,color=color,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(lengthOfBase/2.0,0.5*wallThickness,lengthOfBase/2.0),maxAabb=1.02*Vector3(lengthOfBase/2.0,0.5*wallThickness,lengthOfBase/2.0),maxAabbRotated=1.02*Vector3(lengthOfBase/2.0,0.5*wallThickness,lengthOfBase/2.0),minAabbRotated=1.02*Vector3(lengthOfBase/2.0,0.5*wallThickness,lengthOfBase/2.0),fixedNormal=False)
length=lengthOfBase/1.
V=length*length*wallThickness
geomInertX=(1./12.)*V*(length**2+wallThickness**2)
geomInertY=(1./12.)*V*(length**2+length**2)
geomInertZ=(1./12.)*V*(length**2+wallThickness**2)
utils._commonBodySetup(bbb, V, Vector3(geomInertX,geomInertY,geomInertZ), material='frictionless', pos=[0,0,0], fixed=True)
bbb.state.pos = [0.0,0,0]
lidID = O.bodies.append(bbb)
count=count+1

#Vertical faces A-B-C-D of the box
bA=Body()
bA.mask=3
wire=False
color=[0,0.5,1]
highlight=False
bA.shape=PotentialParticle(k=0.1, r=0.1*wallThickness, R=0.5*heightOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], id=count,isBoundary=True,color=color,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),maxAabb=1.02*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),maxAabbRotated=1.02*Vector3(0.5*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),minAabbRotated=1.02*Vector3(0.5*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),fixedNormal=False)
#length=lengthOfBase
V=lengthOfBase*heightOfBase*wallThickness
geomInertX=(1./12.)*V*(heightOfBase**2 + lengthOfBase**2)
geomInertY=(1./12.)*V*(wallThickness**2 + lengthOfBase**2)
geomInertZ=(1./12.)*V*(wallThickness**2 + heightOfBase**2)
utils._commonBodySetup(bA, V, Vector3(geomInertX,geomInertY,geomInertZ), material='frictionless', pos=[0,0,0], fixed=True)
bA.state.pos = [0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bA)
count=count+1

bB=Body()
bB.mask=3
wire=False
color=[0,0.5,1]
highlight=False
bB.shape=PotentialParticle(k=0.1, r=0.1*wallThickness, R=0.5*heightOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], id=count,isBoundary=True,color=color,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),maxAabb=1.02*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),maxAabbRotated=1.02*Vector3(0.5*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),minAabbRotated=1.02*Vector3(0.5*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),fixedNormal=False)
#length=lengthOfBase
V=lengthOfBase*heightOfBase*wallThickness
geomInertX=(1./12.)*V*(heightOfBase**2 + lengthOfBase**2)
geomInertY=(1./12.)*V*(wallThickness**2 + lengthOfBase**2)
geomInertZ=(1./12.)*V*(wallThickness**2 + heightOfBase**2)
utils._commonBodySetup(bB, V, Vector3(geomInertX,geomInertY,geomInertZ), material='frictionless', pos=[0,0,0], fixed=True)
bB.state.pos = [-0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bB)
count=count+1

bC=Body()
bC.mask=3
wire=False
color=[0,0.5,1]
highlight=False
bC.shape=PotentialParticle(k=0.1, r=0.1*wallThickness, R=0.5*heightOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*lengthOfBase-r,0.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], id=count,isBoundary=True,color=color,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness),maxAabb=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness),maxAabbRotated=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.5*wallThickness),minAabbRotated=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.5*wallThickness),fixedNormal=False)
#length=lengthOfBase
V=lengthOfBase*heightOfBase*wallThickness
geomInertX=(1./12.)*V*(heightOfBase**2 + lengthOfBase**2)
geomInertY=(1./12.)*V*(wallThickness**2 + lengthOfBase**2)
geomInertZ=(1./12.)*V*(wallThickness**2 + heightOfBase**2)
utils._commonBodySetup(bC, V, Vector3(geomInertX,geomInertY,geomInertZ), material='frictionless', pos=[0,0,0], fixed=True)
bC.state.pos = [0,0.5*heightOfBase,0.5*lengthOfBase]
O.bodies.append(bC)
count=count+1

bD=Body()
bD.mask=3
wire=False
color=[0,0.5,1]
highlight=False
bD.shape=PotentialParticle(k=0.1, r=0.1*wallThickness, R=0.5*heightOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*lengthOfBase-r,0.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], id=count,isBoundary=True,color=color,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness),maxAabb=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness),maxAabbRotated=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.5*wallThickness),minAabbRotated=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.5*wallThickness),fixedNormal=False)
#length=lengthOfBase
V=lengthOfBase*heightOfBase*wallThickness
geomInertX=(1./12.)*V*(heightOfBase**2 + lengthOfBase**2)
geomInertY=(1./12.)*V*(wallThickness**2 + lengthOfBase**2)
geomInertZ=(1./12.)*V*(wallThickness**2 + heightOfBase**2)
utils._commonBodySetup(bD, V, Vector3(geomInertX,geomInertY,geomInertZ), material='frictionless', pos=[0,0,0], fixed=True)
bD.state.pos = [0,0.5*heightOfBase,-0.5*lengthOfBase]
ThisID=O.bodies.append(bD)

O.dt = 0.2*sqrt(0.3*O.bodies[0].state.mass/1.0e8)

escapeNo=0
def myAddPlotData():
 global escapeNo
 global wallThickness
 global meanSize
 uf=utils.unbalancedForce()
 if O.iter>25000:
  removeThisID()

from yade import plot
plot.plots={'timeStep1':('unbalancedForce'),'timeStep2':('kineticEn'),'time':('outsideNo')}
#plot.plot() #Uncomment to view plots
O.engines=O.engines+[PyRunner(iterPeriod=100,command='myAddPlotData()')]

def removeThisID():
 global ThisID
 if (O.bodies[ThisID]):
  O.bodies.erase(ThisID)
# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently

O.trackEnergy = True

# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
 if unbalancedForce() < .05:
  O.pause()
  plot.saveDataTxt('bbb.txt.bz2')
  # plot.saveGnuplot('bbb') is also possible

# collect history of data which will be plotted
def addPlotData():
 # each item is given a names, by which it can be the unsed in plot.plots
 # the **O.energy converts dictionary-like O.energy to plot.addData arguments
 plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

##Control the rendering quality
#Gl1_PotentialParticle.aabbEnlargeFactor=1.1
#Gl1_PotentialParticle.sizeX=30
#Gl1_PotentialParticle.sizeY=30
#Gl1_PotentialParticle.sizeZ=30

from yade import qt
qt.Controller()
v=qt.View()

# define how to plot data: 'i' (step number) on the x-axis, unbalanced force
# on the left y-axis, all energies on the right y-axis
# (O.energy.keys is function which will be called to get all defined energies)
# None separates left and right y-axis
plot.plots = {'i': ('unbalanced', None, O.energy.keys)}

# show the plot on the screen, and update while the simulation runs
plot.plot()
O.saveTmp()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Carine Tanissa
Solved:
Last query:
Last reply:
Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#1

Hi Carine,

I could not run your example, because I don't have 'PotentialParticleVTKRecorder', but I am pretty sure that you don't have to modify the code. Just run the script with '-j' parameter. For example

yade -j4 yourscript.py

Cheers,
Karol

Revision history for this message
Carine Tanissa (carinatanissa) said :
#2

Thank you!