particlesSilo.py
You can view and download this file on Github: particlesSilo.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: test with parallel computation and particles
5#
6# Author: Johannes Gerstmayr
7# Date: 2021-11-01
8#
9# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13import exudyn as exu
14from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
15import exudyn.graphics as graphics
16
17import numpy as np
18
19
20SC = exu.SystemContainer()
21mbs = SC.AddSystem()
22
23nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
24
25np.random.seed(1) #always get same results
26
27useGraphics = True
28
29useRigidBody = True #needed for friction
30staticTriangles = True #True speeds up in case of many static objects (silo)
31
32L = 1
33n = 8000 #fast
34# n = 500000
35row = 8*2
36a = L*0.5*0.5*0.75
37stepSize= 0.0001 #Velocity Verlet also works with 2e-4 with 8000 particles
38m = 0.05
39ss=16*3 #the more cells, the more time the search tree needs to build, but contact search is more efficient
40holeRad = 3*a
41SHsilo = 4*L #height of container; adapt to fit all particles
42
43if n >= 4000*8:
44 a*=0.4
45 m *=0.2
46 row=40
47 ss = 16*4
48 holeRad *= 1.4 #better 1.4 !
49 stepSize *= 0.5
50
51if n >= 4000*64:
52 #a*=0.75
53 SHsilo = 8*L
54 row=int(1.3*row)
55 stepSize *= 0.5
56 ss=100
57
58radius = 0.5*a
59t = 0.5*a
60k = 16e4
61d = 0.0004*k
62
63frictionCoeff = 0.5
64if not useRigidBody:
65 frictionCoeff = 0
66
67markerList = []
68radiusList = []
69gDataList = []
70
71
72gContact = mbs.AddGeneralContact()
73#gContact.verboseMode = 1
74gContact.SetFrictionPairings(frictionCoeff*np.eye(1))
75gContact.SetSearchTreeCellSize(numberOfCells=[ss,ss,ss])
76# gContact.computeExactStaticTriangleBins = False #default = True; speeds up a lot for static triangles
77
78
79#%% ground
80LL=6*L
81p0 = np.array([0,0,-0.5*t])
82color4wall = [0.6,0.6,0.6,0.5]
83addNormals = False
84hw=10*a
85gFloor = graphics.Brick(p0,[LL,LL,t],graphics.color.steelblue,addNormals)
86gFloorAdd = graphics.Brick(p0+[-0.5*LL,0,0.5*hw],[t,LL,hw],color4wall,addNormals)
87gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
88gFloorAdd = graphics.Brick(p0+[ 0.5*LL,0,0.5*hw],[t,LL,hw],color4wall,addNormals)
89gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
90gFloorAdd = graphics.Brick(p0+[0,-0.5*LL,0.5*hw],[LL,t,hw],color4wall,addNormals)
91gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
92gFloorAdd = graphics.Brick(p0+[0, 0.5*LL,0.5*hw],[LL,t,hw],color4wall,addNormals)
93gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
94
95gDataList = [gFloor]
96
97
98nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0] ))
99mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
100
101[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gFloor)
102#[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
103gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
104 pointList=meshPoints, triangleList=meshTrigs, staticTriangles=staticTriangles)
105
106if True: #looses color
107 gFloor = graphics.FromPointsAndTrigs(meshPoints, meshTrigs, color=color4wall) #show refined mesh
108 gDataList = [gFloor]
109
110
111color4node = graphics.color.blue
112print("start create: number of masses =",n)
113for i in range(n):
114 kk = int(i/int(n/8))
115 color4node = graphics.colorList[min(kk%9,9)]
116
117 if (i%20000 == 0 and i>0): print("create mass",i)
118 offy = 0
119
120 iz = int(i/(row*row))
121 ix = i%row
122 iy = int(i/row)%row
123
124 if iz % 2 == 1:
125 ix+=0.5
126 iy+=0.5
127
128 offz = 5*L+0.5*a+iz*a*0.74 #0.70x is limit value!
129 offx = -0.6*a-row*0.5*a + (ix+1)*a
130 offy = -0.6*a-row*0.5*a + (iy+1)*a
131
132 valueRand = np.random.random(1)[0]
133 rFact = 0.2 #random part
134 gRad = radius*(1-rFact+rFact*valueRand)
135 v0 = [0,0,-2]
136 pRef = [offx,offy,offz]
137
138 if not useRigidBody:
139 nMass = mbs.AddNode(NodePoint(referenceCoordinates=pRef,
140 initialVelocities=v0,
141 visualization=VNodePoint(show=True,drawSize=2*gRad, color=color4node)))
142
143 oMass = mbs.AddObject(MassPoint(mass=m, nodeNumber=nMass,
144 #visualization=VMassPoint(graphicsData=[gSphere,gSphere2])
145 # visualization=VMassPoint(graphicsData=gData)
146 ))
147 mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
148 else:
149 RBinertia = InertiaSphere(m, radius)
150 RBinertia._nTilesGraphics = 4 #reduce from 16, if drawn
151 dictMass = mbs.CreateRigidBody(
152 inertia=RBinertia,
153 nodeType=exu.NodeType.RotationRotationVector,
154 referencePosition=pRef,
155 initialVelocity=v0,
156 returnDict=True,
157 show=False, #nodes are drawn!
158 )
159 [nMass, oMass] = [dictMass['nodeNumber'], dictMass['bodyNumber']]
160
161 mbs.SetNodeParameter(nMass, 'VdrawSize', 2*gRad)
162 mbs.SetNodeParameter(nMass, 'Vcolor', color4node)
163 mbs.SetNodeParameter(nMass, 'Vshow', True)
164 mThis = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
165
166 mbs.AddLoad(Force(markerNumber=mThis, loadVector= [0,0,-m*9.81]))
167
168 gContact.AddSphereWithMarker(mThis, radius=gRad, contactStiffness=k, contactDamping=d,
169 frictionMaterialIndex=0)
170
171
172
173if True: #add Silo
174 SR = 3.1*L
175 SH2 = 1*L #hole
176 SR2 = holeRad #hole
177 ST = 0.25*L
178 #contour=8*np.array([[0,0.2],[0.3,0.2],[0.5,0.3],[0.7,0.4],[1,0.4],[1,0.]])
179 contour=np.array([[0,SR2],[0,SR2+ST],[SH2-ST,SR2+ST],[2*SH2-ST,SR+ST],[2*SH2+SHsilo,SR+ST],
180 [2*SH2+SHsilo,SR],[2*SH2,SR],[SH2,SR2],[0,SR2]])
181 contour = list(contour)
182 # contour.reverse()
183 gSilo = graphics.SolidOfRevolution(pAxis=[0,0,3*L], vAxis=[0,0,1],
184 contour=contour, color=[0.8,0.1,0.1,0.5], nTiles=32)
185
186 [meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gSilo)
187 gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
188 pointList=meshPoints, triangleList=meshTrigs, staticTriangles=staticTriangles)
189
190
191#put here, such that it is transparent in background
192oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
193 visualization=VObjectGround(graphicsData=[gSilo]+gDataList)))
194
195
196mbs.Assemble()
197print("finish gContact")
198# print(gContact.GetPythonObject())
199
200items=gContact.GetItemsInBox(pMin=[-4,-4,0], pMax=[4,4,20])
201print('n spheres=',len(items['MarkerBasedSpheres']), ', ss=',ss)
202
203
204tEnd = 10
205simulationSettings = exu.SimulationSettings()
206simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
207#simulationSettings.solution.file.write = True
208simulationSettings.solution.file.write = True
209simulationSettings.solution.file.writePeriod = 0.02
210simulationSettings.solution.precision = 5 #make files smaller
211simulationSettings.solution.file.export.accelerations = False
212simulationSettings.solution.file.export.velocities = False
213simulationSettings.solution.file.name = 'solution/test.txt'
214simulationSettings.show.computationTime = True
215#simulationSettings.show.statistics = True
216simulationSettings.timeIntegration.verboseMode = 1
217simulationSettings.timeIntegration.stepInformation += 32 #show time to go
218simulationSettings.parallel.numberOfThreads = 8 #this should not be higher than the number of real cores (not threads)
219
220simulationSettings.timeIntegration.explicit.computeEndOfStepAccelerations = False
221simulationSettings.timeIntegration.explicit.computeMassMatrixInversePerBody = True
222
223SC.visualizationSettings.general.graphicsUpdateInterval=2
224SC.visualizationSettings.general.circleTiling=20
225SC.visualizationSettings.view0.scene.drawCoordinateSystem=True
226SC.visualizationSettings.loads.show=False
227SC.visualizationSettings.bodies.show=True
228SC.visualizationSettings.markers.show=False
229
230SC.visualizationSettings.nodes.show=True
231SC.visualizationSettings.nodes.drawNodesAsPoint = False
232SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
233SC.visualizationSettings.nodes.tiling = 8
234
235SC.visualizationSettings.view0.window.renderWindowSize=[1200,1200]
236#SC.visualizationSettings.view0.window.renderWindowSize=[1024,1400]
237SC.visualizationSettings.openGL.multiSampling = 4
238#improved OpenGL rendering
239
240SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
241SC.visualizationSettings.exportImages.saveImageTimeOut=10000 #5000 is too shot sometimes!
242
243if False:
244 simulationSettings.solution.recordImagesInterval = 0.005
245
246
247simulate=True
248if simulate:
249 if useGraphics:
250 SC.visualizationSettings.general.autoFitScene = False
251 SC.renderer.Start()
252 SC.renderer.RestoreSavedState()
253 #SC.renderer.DoIdleTasks()
254
255 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
256 simulationSettings.timeIntegration.endTime = tEnd
257 mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.VelocityVerlet)
258 #print(gContact)
259
260 if useGraphics:
261 SC.renderer.DoIdleTasks()
262 SC.renderer.Stop() #safely close rendering window!
263
264if not simulate:
265 SC.visualizationSettings.general.autoFitScene = False
266 SC.visualizationSettings.general.graphicsUpdateInterval=0.5
267
268 print('load solution file')
269 #sol = LoadSolutionFile(OutputFilePath('solution/test2.txt'), safeMode=False)
270 sol = LoadSolutionFile(OutputFilePath('solution/test.txt'), safeMode=True, verbose = True)#, maxRows=100)
271 print('start SolutionViewer')
272 mbs.SolutionViewer(sol)
273
274
275#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
276#timings (best of 3), a=0.25, Core i9-7940X@3.10GHz:
277#1.9.66, ss=32, staticTriangles = False
278# n spheres= 8000 , ss= 16
279# +++++++++++++++++++++++++++++++
280# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
281# Start multi-threading with 8 threads
282# STEP526, t = 0.0526s, timeToGo = 5.61318s, Nit/step = 0
283# STEP1031, t = 0.1031s, timeToGo = 3.75953s, Nit/step = 0
284# STEP1530, t = 0.153s, timeToGo = 1.84418s, Nit/step = 0
285# STEP2000, t = 0.2s, timeToGo = 2.26949e-13s, Nit/step = 0
286# Solver terminated successfully after 7.96149 seconds.
287# ====================
288# CPU-time statistics:
289# total time = 7.96 seconds
290# measured time= 7.65 seconds (=96.1%)
291# non-zero timer [__ sub-timer]:
292# newtonIncrement = 2.97%
293# integrationFormula= 6.03%
294# ODE2RHS = 84.2%
295# writeSolution = 5.01%
296# overhead = 1.74%
297# visualization/user= 0.00322%
298# special timers:
299# Contact:BoundingBoxes = 0.92391 (12.1%)s
300# Contact:SearchTree = 0.59905 (7.83%)s
301# Contact:Overall = 5.1361 (67.2%)s
302
303# n spheres= 8000 , ss= 32
304# +++++++++++++++++++++++++++++++
305# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
306# Start multi-threading with 8 threads
307# STEP562, t = 0.0562s, timeToGo = 5.12578s, Nit/step = 0
308# STEP1126, t = 0.1126s, timeToGo = 3.10595s, Nit/step = 0
309# STEP1688, t = 0.1688s, timeToGo = 1.10932s, Nit/step = 0
310# STEP2000, t = 0.2s, timeToGo = 2.03936e-13s, Nit/step = 0
311# Solver terminated successfully after 7.15063 seconds.
312# ====================
313# CPU-time statistics:
314# total time = 7.15 seconds
315# measured time= 6.86 seconds (=95.9%)
316# non-zero timer [__ sub-timer]:
317# newtonIncrement = 3.16%
318# integrationFormula= 6.44%
319# ODE2RHS = 84.1%
320# writeSolution = 4.47%
321# overhead = 1.81%
322# visualization/user= 0.00365%
323# special timers:
324# Contact:BoundingBoxes = 0.89462 (13%)s
325# Contact:SearchTree = 1.142 (16.7%)s
326# Contact:Overall = 4.5078 (65.8%)s
327
328# n spheres= 8000 , ss= 48
329# +++++++++++++++++++++++++++++++
330# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
331# Start multi-threading with 8 threads
332# STEP440, t = 0.044s, timeToGo = 7.10356s, Nit/step = 0
333# STEP861, t = 0.0861s, timeToGo = 5.29512s, Nit/step = 0
334# STEP1294, t = 0.1294s, timeToGo = 3.27367s, Nit/step = 0
335# STEP1707, t = 0.1707s, timeToGo = 1.37385s, Nit/step = 0
336# STEP2000, t = 0.2s, timeToGo = 2.71236e-13s, Nit/step = 0
337# Solver terminated successfully after 9.50154 seconds.
338# ====================
339# CPU-time statistics:
340# total time = 9.5 seconds
341# measured time= 9.17 seconds (=96.5%)
342# non-zero timer [__ sub-timer]:
343# newtonIncrement = 2.45%
344# integrationFormula= 4.94%
345# ODE2RHS = 87.5%
346# writeSolution = 3.6%
347# overhead = 1.54%
348# visualization/user= 0.00405%
349# special timers:
350# Contact:BoundingBoxes = 0.92908 (10.1%)s
351# Contact:SearchTree = 2.3468 (25.6%)s
352# Contact:Overall = 6.676 (72.8%)s
353
354
355
356#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
357#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
358#1.9.66, staticTriangles = True
359
360# n spheres= 8000 , ss= 16
361# +++++++++++++++++++++++++++++++
362# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
363# Start multi-threading with 8 threads
364# WARNING: VelocityVerlet: still under development
365# STEP526, t = 0.0526s, timeToGo = 5.60837s, Nit/step = 0
366# STEP1091, t = 0.1091s, timeToGo = 3.33273s, Nit/step = 0
367# STEP1655, t = 0.1655s, timeToGo = 1.25102s, Nit/step = 0
368# STEP2000, t = 0.2s, timeToGo = 2.0821e-13s, Nit/step = 0
369# Solver terminated successfully after 7.30395 seconds.
370# ====================
371# CPU-time statistics:
372# total time = 7.3 seconds
373# measured time= 7.01 seconds (=96%)
374# non-zero timer [__ sub-timer]:
375# newtonIncrement = 3.19%
376# integrationFormula= 6.58%
377# ODE2RHS = 83.6%
378# writeSolution = 4.71%
379# overhead = 1.93%
380# visualization/user= 0.003%
381# special timers:
382# Contact:BoundingBoxes = 0.98863 (14.1%)s
383# Contact:SearchTree = 0.4626 (6.6%)s
384# Contact:Overall = 4.3393 (61.9%)s
385
386#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
387#update to 1.9.68 (exact bin computation for triangles)
388# start create: number of masses = 8000
389# create mass 0
390# finish gContact
391# n spheres= 8000 , ss= 32
392# +++++++++++++++++++++++++++++++
393# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
394# Start multi-threading with 8 threads
395# STEP639, t = 0.0639s, timeToGo = 4.26062s, Nit/step = 0
396# STEP1309, t = 0.1309s, timeToGo = 2.11214s, Nit/step = 0
397# STEP2000, t = 0.2s, timeToGo = 1.70989e-13s, Nit/step = 0
398# Solver terminated successfully after 5.99597 seconds.
399# ====================
400# CPU-time statistics:
401# total time = 6 seconds
402# measured time= 5.72 seconds (=95.4%)
403# non-zero timer [__ sub-timer]:
404# newtonIncrement = 3.65%
405# integrationFormula= 7.59%
406# ODE2RHS = 81.1%
407# writeSolution = 5.57%
408# overhead = 2.07%
409# visualization/user= 0.00279%
410# special timers:
411# Contact:BoundingBoxes = 0.86635 (15.1%)s
412# Contact:SearchTree = 0.59121 (10.3%)s
413# Contact:Overall = 3.3974 (59.4%)s
414
415
416# runfile('C:/DATA/cpp/EXUDYN_git/main/pythonDev/Examples/particlesSilo.py', wdir='C:/DATA/cpp/EXUDYN_git/main/pythonDev/Examples')
417# start create: number of masses = 8000
418# create mass 0
419# finish gContact
420# n spheres= 8000 , ss= 48
421# +++++++++++++++++++++++++++++++
422# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
423# Start multi-threading with 8 threads
424# STEP517, t = 0.0517s, timeToGo = 5.74121s, Nit/step = 0
425# STEP1051, t = 0.1051s, timeToGo = 3.614s, Nit/step = 0
426# STEP1601, t = 0.1601s, timeToGo = 1.49891s, Nit/step = 0
427# STEP2000, t = 0.2s, timeToGo = 2.14357e-13s, Nit/step = 0
428# Solver terminated successfully after 7.51458 seconds.
429# ====================
430# CPU-time statistics:
431# total time = 7.51 seconds
432# measured time= 7.21 seconds (=95.9%)
433# non-zero timer [__ sub-timer]:
434# newtonIncrement = 3.12%
435# integrationFormula= 6.21%
436# ODE2RHS = 84.7%
437# writeSolution = 4.11%
438# overhead = 1.84%
439# visualization/user= 0.00239%
440# special timers:
441# Contact:BoundingBoxes = 0.88801 (12.3%)s
442# Contact:SearchTree = 0.93017 (12.9%)s
443# Contact:Overall = 4.8152 (66.8%)s
444
445#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
446#large test:
447# 500000 particles, stepSize=2.5e-5
448# STEP40046, t = 1.00115s, timeToGo = 5.81135 days, tCPU=2.84971h, Nit/step = 0
449# STEP121643, t = 3.04107s, timeToGo = 6.78003 days, tCPU=10.5378h, Nit/step = 0
450# STEP121644 (stopped), t = 3.04107s, tCPU=37936.2s, Nit/step = 0
451# solver stopped by user after 37935.8 seconds.
452# ====================
453# CPU-time statistics:
454# total time = 3.79e+04 seconds
455# measured time= 3.63e+04 seconds (=95.8%)
456# non-zero timer [__ sub-timer]:
457# newtonIncrement = 3.02%
458# integrationFormula= 5.06%
459# ODE2RHS = 86.1%
460# writeSolution = 0.974%
461# overhead = 1.18%
462# visualization/user= 3.61%
463# special timers:
464# Contact:BoundingBoxes = 3156.8 (8.68%)s
465# Contact:SearchTree = 3350.5 (9.22%)s
466# Contact:Overall = 25633 (70.5%)s