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