← Back to team overview

yade-users team mailing list archive

  1. yade-users team
  2. Mailing list archive
  3. Message #13255

[Question #341304]: Landing of a satellite

  Thread Previous • Date Previous • • Thread Next 
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.
  Thread Previous • Date Previous • • Thread Next