Re: [Question #676451]: optimize camputational time for vibrated granular media

```Question #676451 on Yade changed:

Bruno Chareyre posted a new comment:
I see.
Well, as I just reduced the radius a little to make them all fit in I can give  you some results with -j4  (also reduced verletDist but the speedup is small, script below):

200 spheres: 0.52 sec for 10k iterations
2000 spheres: 5.5 sec for 10k iterations

That's nearly linear.
Unless the the cost in LAMMPS is sub-linear (how is that possible?!) I can't imagine how yade becomes 100x slower. Do you have approximately similar numbers?

Could it be that sphere insertions makes it *very* different?

Bruno

import numpy as np
import math
import sys
import random

if(len(O.bodies)<Nball+nBodyCont):
sp=pack.SpherePack()
sp.makeCloud(minCorner=(-rcont/1.4,-rcont/1.4,hcone),maxCorner=(rcont/1.4,rcont/1.4,hcone+hcil),rMean=rball*0.8,num=Nball+60-len(O.bodies))
sp.toSimulation()

#Dimensioni utili
rball=0.002
hcone=6.37*rball
rcont=22.5*rball
rconelow=4*rball
hcil=17.13*rball

Nball=int(sys.argv[1]) #Number Of Ball
A=0.00025 #amplitude shaker
fr=200 #freq shaker

#Steel
densSteel=8000
ySteel=21e7 #originale e9
poisSteel=0.293
shearModStell=ySteel/(2*(1+poisSteel))

#Plaexiglass
densPMMA=1190
yPMMA=33e6
poissPMMA=0.37

#plexiglass as lammps
O.materials.append(FrictMat(young=yPMMA, poisson=poissPMMA,
frictionAngle=frictAngle,density=densPMMA, label='lmpPMMA'))

#Steel as lammps
O.materials.append(FrictMat(young=ySteel, poisson=poisSteel,
frictionAngle=frictAngle,density=densSteel, label='lmpSteel'))

#rayleigh time
tRay=math.pi*rball*(densSteel/shearModStell)**(0.5)/(0.1631*poisSteel+0.8766)

#container
contenitore=coneId+cyliId

nBodyCont=len(O.bodies) #60 facets

track=[]
for k in range(nBodyCont,len(O.bodies)):
track.append(O.bodies[k])

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),],label='collider',verletDist=rball*0.3),
InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),
], [Ip2_FrictMat_FrictMat_MindlinPhys(betan=0.4, betas=0.4),], [Law2_ScGeom_MindlinPhys_Mindlin(),],),
NewtonIntegrator(gravity=(0,0,-9.8),damping=0.0),
]

O.engines = O.engines + [HarmonicMotionEngine(ids = contenitore, A = (0,0,A), f = (0,0,fr), label='shaker')]
O.engines = O.engines + [PyRunner(command = "print O.iter, O.time,O.speed, O.realtime,kineticEnergy()",iterPeriod=10000)]+[PyRunner(command =

O.dt=0.2*tRay
O.timingEnabled=False
O.run(20000,1)
timing.reset()
O.timingEnabled=True
O.run(10000,1)
timing.stats()

--