Question #341304 on Yade changed:
https://answers.launchpad.net/yade/+question/341304
Status: Open => Answered
Jan Stránský proposed the following answer:
Hi Jonathan,
not studied your script very much in detail, but you may try:
- decrease number of cores (contrary to what one expects :-). Sometimes
using more cores actually make the simulation take longer time, see [1],
page 470.
- not using polyhedrons, as from its nature it takes much longer time to
evaluate them than spheres
- increase time step (instead of 0.5 use higher factor)
- increase particle size (increasing time step)
- decrease number of particle number (using some tricks with boundary
conditions etc.)
- decrease stiffness or increase mass to increase time step (but this
changes also the physics), perhaps increase mass of just a few smallest
particles if you use some PSD (cannot be determined from your script)
- use timing stats to check if some of PyRunner functions does not use too
much time (not likely since you call them not very often)
Anyway, none of the options I gave decrease the simulation time
drastically..
cheers
Jan
2016-08-12 18:17 GMT+02:00 Jonathan Pergoli <
question341...@answers.launchpad.net>:
> New question #341304 on Yade:
> https://answers.launchpad.net/yade/+question/341304
>
> Hi guys,
>
> I am performing simulation in which a satellite lands on the soil. The
> simulation deal with more than 70.000 particles and take a week to make 10
> seconds of simulations using 10 cores. I'd like to know if it is possible
> to speed it up.
> This is the code:
>
> from yade import utils,ymport,export,plot
> import math as m
>
>
> # Material
> E1=1e+8
> E2=5e+7
> mli=FrictMat(density=643,frictionAngle=0.1489,label="MLI",young=E2)
> MLI=O.materials.append(mli)
> gravel=FrictMat(density=1700,frictionAngle=0.7853,label="gravel",young=E1)
> GRAVEL=O.materials.append(gravel)
>
> # Ground
> s=ymport.textExt('1cm_0g2.txt',format='x_y_z_r')
> sphere=O.bodies.append(s)
>
>
> # Create a vector o spheres to eliminate those that have COG above the
> container
> for i in O.bodies:
> if isinstance(i.shape,Sphere):
> if i.state.pos[2]>.30:
> O.bodies.erase(i.id)
> print i.state.mass
> for i in O.bodies:
> if isinstance(i.shape,Sphere):
> if i.state.pos[2]<0:
> O.bodies.erase(i.id)
> for i in O.bodies:
> if isinstance(i.shape,Sphere):
> x=i.state.pos[0]
> y=i.state.pos[1]
> r=m.sqrt(x**2+y**2)
> if r>.75:
> O.bodies.erase(i.id)
>
> aa=[]
> for i in O.bodies:
> if isinstance(i.shape,Sphere):
> aa.append(i.id)
>
> print len(aa)
>
> # Cylinder
> hc=.30
> c=geom.facetCylinder((0,0,.15),radius=.75,height=hc,
> segmentsNumber=100,wallMask=6,material="gravel")
> O.bodies.append(c)
>
> # SAT
> a=.2774
> b=.2922
> c=.1973
> aa=a/2
> bb=b/2
> cc=c/2
> h=.65
> dist=0
> theta=0
> thetav=0
> v1=(aa,bb,c)
> v2=(aa,-bb,c)
> v3=(-aa,-bb,c)
> v4=(-aa,bb,c)
> v5=(aa,bb,0)
> v6=(aa,-bb,0)
> v7=(-aa,-bb,0)
> v8=(-aa,bb,0)
> V=[v1,v2,v3,v4,v5,v6,v7,v8]
> vz=.19
> R=[[m.cos(theta),0,m.sin(theta)],[0,1,0],[-m.sin(theta),0,m.cos(theta)]]
> v1=(R[0][0]*V[0][0]+R[0][1]*V[0][1]+R[0][2]*V[0][2],R[1][0]*
> V[0][0]+R[1][1]*V[0][1]+R[1][2]*V[0][2],R[2][0]*V[0][0]+R[
> 2][1]*V[0][1]+R[2][2]*V[0][2])
> v2=(R[0][0]*V[1][0]+R[0][1]*V[1][1]+R[0][2]*V[1][2],R[1][0]*
> V[1][0]+R[1][1]*V[1][1]+R[1][2]*V[1][2],R[2][0]*V[1][0]+R[
> 2][1]*V[1][1]+R[2][2]*V[1][2])
> v3=(R[0][0]*V[2][0]+R[0][1]*V[2][1]+R[0][2]*V[2][2],R[1][0]*
> V[2][0]+R[1][1]*V[2][1]+R[1][2]*V[2][2],R[2][0]*V[2][0]+R[
> 2][1]*V[2][1]+R[2][2]*V[2][2])
> v4=(R[0][0]*V[3][0]+R[0][1]*V[3][1]+R[0][2]*V[3][2],R[1][0]*
> V[3][0]+R[1][1]*V[3][1]+R[1][2]*V[3][2],R[2][0]*V[3][0]+R[
> 2][1]*V[3][1]+R[2][2]*V[3][2])
> v5=(R[0][0]*V[4][0]+R[0][1]*V[4][1]+R[0][2]*V[4][2],R[1][0]*
> V[4][0]+R[1][1]*V[4][1]+R[1][2]*V[4][2],R[2][0]*V[4][0]+R[
> 2][1]*V[4][1]+R[2][2]*V[4][2])
> v6=(R[0][0]*V[5][0]+R[0][1]*V[5][1]+R[0][2]*V[5][2],R[1][0]*
> V[5][0]+R[1][1]*V[5][1]+R[1][2]*V[5][2],R[2][0]*V[5][0]+R[
> 2][1]*V[5][1]+R[2][2]*V[5][2])
> v7=(R[0][0]*V[6][0]+R[0][1]*V[6][1]+R[0][2]*V[6][2],R[1][0]*
> V[6][0]+R[1][1]*V[6][1]+R[1][2]*V[6][2],R[2][0]*V[6][0]+R[
> 2][1]*V[6][1]+R[2][2]*V[6][2])
> v8=(R[0][0]*V[7][0]+R[0][1]*V[7][1]+R[0][2]*V[7][2],R[1][0]*
> V[7][0]+R[1][1]*V[7][1]+R[1][2]*V[7][2],R[2][0]*V[7][0]+R[
> 2][1]*V[7][1]+R[2][2]*V[7][2])
> p=utils.polyhedron((v1,v2,v3,v4,v5,v6,v7,v8),fixed=False,
> color=(.6,.45,0),material="MLI",wire=False)
> SAT=O.bodies.append(p)
> p.state.vel=(vz*m.sin(thetav),0,-vz*m.cos(thetav))
> p.state.ori=((0,-1,0),theta)
> p.state.pos=(-dist,0,h)
> Ixx=0.081026
> Iyy=0.10031
> Izz=0.12116
> p.state.inertia=(Ixx,Iyy,Izz)
> M=p.id
> r=m.sqrt(aa**2+bb**2)
> Rj=m.sqrt(r**2+cc**2)
> Ri=0.05
> Rr=Rj*Ri/(Rj+Ri)
> mu_rM=0.016
> mu_rG=2.05
> KN=E1*2*Ri*E1*2*Ri/(E1*2*Ri+E1*2*Ri)#6.5e+4
> KR=3*Ri**2*mu_rG**2*KN/4
> print "SAT's mass = ",p.state.mass
> print "SAT's position = ",p.state.pos
> print "SAT's orientation = ",p.state.ori
> print "SAT's inertia = ",p.state.inertia
> print "Timestep = ",O.dt
>
> # Functions
>
> def forces():
> # rotation of axis
> q1=p.state.ori[0]
> q2=p.state.ori[1]
> q3=p.state.ori[2]
> q4=p.state.ori[3]
> RR=[[q1**2-q2**2-q3**2+q4**2,2*(q1*q2+q3*q4),2*(q1*q3-q2*
> q4)],[2*(q1*q2-q3*q4),-q1**2+q2**2-q3**2+q4**2,2*(q2*q3+q1*
> q4)],[2*(q1*q3+q2*q4),2*(q2*q3-q1*q4),-q1**2-q2**2+q3**2+q4**2]]
> e1=(RR[0][0],RR[0][1],RR[0][2])
> e2=(RR[1][0],RR[1][1],RR[1][2])
> e3=(RR[2][0],RR[2][1],RR[2][2])
> massa1=0
> massa2=0
> massa3=0
> for i in O.bodies:
> if isinstance(i.shape,Sphere):
> if i.state.vel[2]>.001:
> massa1+=i.state.mass
> if i.state.pos[2]>.4:
> massa2+=i.state.mass
> if i.state.vel[2]>.01:
> massa3+=i.state.mass
> # forces
> #Fx=utils.sumForces([MASCOT],e1)
> #Fy=utils.sumForces([MASCOT],e2)
> #Fz=utils.sumForces([MASCOT],e3)
> #Tx=utils.sumTorques([MASCOT],axis=e1,axisPt=(p.state.pos[0]
> -aa,p.state.pos[1],p.state.pos[2]))
> #Ty=utils.sumTorques([MASCOT],axis=e2,axisPt=(p.state.pos[0]
> ,p.state.pos[1]-bb,p.state.pos[2]))
> #Tz=utils.sumTorques([MASCOT],axis=e3,axisPt=(p.state.pos[0]
> ,p.state.pos[1],p.state.pos[2]-cc))
> # energy
> vx=p.state.vel[0]
> vy=p.state.vel[1]
> vz=p.state.vel[2]
> wx=p.state.angVel[0]
> wy=p.state.angVel[1]
> wz=p.state.angVel[2]
> K_lin=.5*p.state.mass*(vx**2+vy**2+vz**2)
> K_rot=.5*(Ixx*wx**2+Iyy*wy**2+Izz*wz**2)
> # plot
> plot.addData(i=O.iter,t=O.time,vz=p.state.vel[2],depth=
> p.state.pos[2]-.4,E=O.energy.total(),Ek=utils.kineticEnergy(),wx=p.state.
> angVel[0],wy=p.state.angVel[1],wz=p.state.angVel[2],x=p.
> state.pos[0],y=p.state.pos[1],vx=p.state.vel[0],vy=p.state.
> vel[1],massa1=massa1,massa2=massa2,massa3=massa3,q1=q1,q2=
> q2,q3=q3,q4=q4,Ek_lin=K_lin,Ek_rot=K_rot)
> plot.saveDataTxt('SATflatflat.txt',vars=('i','t','vz','
> depth','E','Ek','wx','wy','wz','x','y','vx','vy','massa1','
> massa2','massa3','q1','q2','q3','q4','Ek_lin','Ek_rot'))
> export.textExt('particles.txt',format='x_y_z_r',comment='Final
> position of the spheres')
>
>
> def checktime():
> if O.time>20:
> O.pause()
> print ("Simulation time = %f" % O.time)
> def itercount():
> print ("Simulation iter 0deg = %f" % O.iter)
>
> # Engines
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()
> ,Bo1_Polyhedra_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_
> Sphere_Polyhedra_ScGeom()],
> [Ip2_FrictMat_FrictMat_MindlinPhys(en=.55,es=.55,
> krot=KR,frictAngle=.7853)],#MatchMaker(matches=((GRAVEL,
> MLI,.1489),(GRAVEL,GRAVEL,.7853))))],
> [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True)]
> ),
> NewtonIntegrator(gravity=(0,0,-2.5e-4),damping=0),
> PyRunner(command='checktime()',realPeriod=3600),
> PyRunner(command='forces()',realPeriod=180),
> PyRunner(command='vtkExporter.exportPolyhedra()',iterPeriod=
> 15000),
> PyRunner(command='vtkExporter.exportSpheres(ids="all",what=[
> (velocities,s_velocities)])',iterPeriod=15000),
> PyRunner(command='vtkExporter.exportInteractions(ids="all",
> what=[(n,NF)])',iterPeriod=15000),
> #PyRunner(command='vtkExporter.exportInteractions(ids="all",
> what=[(cr,MB)])',iterPeriod=15000)
> ]
>
> O.trackEnergy=True
> O.save('SATflatflat')
> Prova="SATflatflat"
> vtkExporter = export.VTKExporter(Prova)
> vtkExporter.exportFacets(ids="all",what=[('pos','b.state.pos')])
> #vtkExporter.exportSpheres(ids="all",what=[('vel','b.state.vel')])
> velocities='vel'
> s_velocities='b.state.vel'
> n='normalstress'
> NF='i.phys.normalForce'
>
> # Final
> O.dt=.05*PWaveTimeStep()
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : yade-users@lists.launchpad.net
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
