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
 51if n >= 4000*64:
 52    #a*=0.75
 53    SHsilo = 8*L
 54    row=int(1.3*row)
 55    stepSize *= 0.25
 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(physicsMass=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        dictMass = mbs.CreateRigidBody(
151                      inertia=RBinertia,
152                      nodeType=exu.NodeType.RotationRotationVector,
153                      referencePosition=pRef,
154                      initialVelocity=v0,
155                      returnDict=True)
156        [nMass, oMass] = [dictMass['nodeNumber'], dictMass['bodyNumber']]
157
158        mbs.SetNodeParameter(nMass, 'VdrawSize', 2*gRad)
159        mbs.SetNodeParameter(nMass, 'Vcolor', color4node)
160        mThis = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
161
162    mbs.AddLoad(Force(markerNumber=mThis, loadVector= [0,0,-m*9.81]))
163
164    gContact.AddSphereWithMarker(mThis, radius=gRad, contactStiffness=k, contactDamping=d,
165                                     frictionMaterialIndex=0)
166
167
168
169if True: #add Silo
170    SR = 3.1*L
171    SH2 = 1*L #hole
172    SR2 = holeRad   #hole
173    ST = 0.25*L
174    #contour=8*np.array([[0,0.2],[0.3,0.2],[0.5,0.3],[0.7,0.4],[1,0.4],[1,0.]])
175    contour=np.array([[0,SR2],[0,SR2+ST],[SH2-ST,SR2+ST],[2*SH2-ST,SR+ST],[2*SH2+SHsilo,SR+ST],
176                      [2*SH2+SHsilo,SR],[2*SH2,SR],[SH2,SR2],[0,SR2]])
177    contour = list(contour)
178    contour.reverse()
179    gSilo = graphics.SolidOfRevolution(pAxis=[0,0,3*L], vAxis=[0,0,1],
180            contour=contour, color=[0.8,0.1,0.1,0.5], nTiles = 64)
181
182    [meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gSilo)
183    gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
184        pointList=meshPoints,  triangleList=meshTrigs, staticTriangles=staticTriangles)
185
186
187#put here, such that it is transparent in background
188oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
189                                   visualization=VObjectGround(graphicsData=[gSilo]+gDataList)))
190
191
192mbs.Assemble()
193print("finish gContact")
194# print(gContact.GetPythonObject())
195
196items=gContact.GetItemsInBox(pMin=[-4,-4,0], pMax=[4,4,20])
197print('n spheres=',len(items['MarkerBasedSpheres']), ', ss=',ss)
198
199
200tEnd = 10
201simulationSettings = exu.SimulationSettings()
202simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
203#simulationSettings.solutionSettings.writeSolutionToFile = True
204simulationSettings.solutionSettings.writeSolutionToFile = True
205simulationSettings.solutionSettings.solutionWritePeriod = 0.02
206simulationSettings.solutionSettings.outputPrecision = 5 #make files smaller
207simulationSettings.solutionSettings.exportAccelerations = False
208simulationSettings.solutionSettings.exportVelocities = False
209simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/test.txt'
210simulationSettings.displayComputationTime = True
211#simulationSettings.displayStatistics = True
212simulationSettings.timeIntegration.verboseMode = 1
213simulationSettings.timeIntegration.stepInformation += 32 #show time to go
214simulationSettings.parallel.numberOfThreads = 8 #this should not be higher than the number of real cores (not threads)
215
216simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False
217simulationSettings.timeIntegration.explicitIntegration.computeMassMatrixInversePerBody = True
218
219SC.visualizationSettings.general.graphicsUpdateInterval=2
220SC.visualizationSettings.general.circleTiling=200
221SC.visualizationSettings.general.drawCoordinateSystem=True
222SC.visualizationSettings.loads.show=False
223SC.visualizationSettings.bodies.show=True
224SC.visualizationSettings.markers.show=False
225
226SC.visualizationSettings.nodes.show=True
227SC.visualizationSettings.nodes.drawNodesAsPoint = False
228SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
229SC.visualizationSettings.nodes.tiling = 4
230
231SC.visualizationSettings.window.renderWindowSize=[1200,1200]
232#SC.visualizationSettings.window.renderWindowSize=[1024,1400]
233SC.visualizationSettings.openGL.multiSampling = 4
234#improved OpenGL rendering
235
236SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
237SC.visualizationSettings.exportImages.saveImageTimeOut=10000 #5000 is too shot sometimes!
238
239if False:
240    simulationSettings.solutionSettings.recordImagesInterval = 0.005
241
242
243simulate=True
244if simulate:
245    if useGraphics:
246        SC.visualizationSettings.general.autoFitScene = False
247        exu.StartRenderer()
248        if 'renderState' in exu.sys:
249            SC.SetRenderState(exu.sys['renderState'])
250        #mbs.WaitForUserToContinue()
251
252    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
253    simulationSettings.timeIntegration.endTime = tEnd
254    mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.VelocityVerlet)
255    #print(gContact)
256
257    if useGraphics:
258        SC.WaitForRenderEngineStopFlag()
259        exu.StopRenderer() #safely close rendering window!
260
261if not simulate:
262    SC.visualizationSettings.general.autoFitScene = False
263    SC.visualizationSettings.general.graphicsUpdateInterval=0.5
264
265    print('load solution file')
266    #sol = LoadSolutionFile('solution/test2.txt', safeMode=False)
267    sol = LoadSolutionFile('solution/test.txt', safeMode=True, verbose = True)#, maxRows=100)
268    print('start SolutionViewer')
269    mbs.SolutionViewer(sol)
270
271
272#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
273#timings (best of 3), a=0.25, Core i9-7940X@3.10GHz:
274#1.9.66, ss=32, staticTriangles = False
275# n spheres= 8000 , ss= 16
276# +++++++++++++++++++++++++++++++
277# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
278# Start multi-threading with 8 threads
279# STEP526, t = 0.0526s, timeToGo = 5.61318s, Nit/step = 0
280# STEP1031, t = 0.1031s, timeToGo = 3.75953s, Nit/step = 0
281# STEP1530, t = 0.153s, timeToGo = 1.84418s, Nit/step = 0
282# STEP2000, t = 0.2s, timeToGo = 2.26949e-13s, Nit/step = 0
283# Solver terminated successfully after 7.96149 seconds.
284# ====================
285# CPU-time statistics:
286#   total time   = 7.96 seconds
287#   measured time= 7.65 seconds (=96.1%)
288#   non-zero timer [__ sub-timer]:
289#   newtonIncrement   = 2.97%
290#   integrationFormula= 6.03%
291#   ODE2RHS           = 84.2%
292#   writeSolution     = 5.01%
293#   overhead          = 1.74%
294#   visualization/user= 0.00322%
295# special timers:
296#   Contact:BoundingBoxes = 0.92391 (12.1%)s
297#   Contact:SearchTree = 0.59905 (7.83%)s
298#   Contact:Overall = 5.1361 (67.2%)s
299
300# n spheres= 8000 , ss= 32
301# +++++++++++++++++++++++++++++++
302# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
303# Start multi-threading with 8 threads
304# STEP562, t = 0.0562s, timeToGo = 5.12578s, Nit/step = 0
305# STEP1126, t = 0.1126s, timeToGo = 3.10595s, Nit/step = 0
306# STEP1688, t = 0.1688s, timeToGo = 1.10932s, Nit/step = 0
307# STEP2000, t = 0.2s, timeToGo = 2.03936e-13s, Nit/step = 0
308# Solver terminated successfully after 7.15063 seconds.
309# ====================
310# CPU-time statistics:
311#   total time   = 7.15 seconds
312#   measured time= 6.86 seconds (=95.9%)
313#   non-zero timer [__ sub-timer]:
314#   newtonIncrement   = 3.16%
315#   integrationFormula= 6.44%
316#   ODE2RHS           = 84.1%
317#   writeSolution     = 4.47%
318#   overhead          = 1.81%
319#   visualization/user= 0.00365%
320# special timers:
321#   Contact:BoundingBoxes = 0.89462 (13%)s
322#   Contact:SearchTree = 1.142 (16.7%)s
323#   Contact:Overall = 4.5078 (65.8%)s
324
325# n spheres= 8000 , ss= 48
326# +++++++++++++++++++++++++++++++
327# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
328# Start multi-threading with 8 threads
329# STEP440, t = 0.044s, timeToGo = 7.10356s, Nit/step = 0
330# STEP861, t = 0.0861s, timeToGo = 5.29512s, Nit/step = 0
331# STEP1294, t = 0.1294s, timeToGo = 3.27367s, Nit/step = 0
332# STEP1707, t = 0.1707s, timeToGo = 1.37385s, Nit/step = 0
333# STEP2000, t = 0.2s, timeToGo = 2.71236e-13s, Nit/step = 0
334# Solver terminated successfully after 9.50154 seconds.
335# ====================
336# CPU-time statistics:
337#   total time   = 9.5 seconds
338#   measured time= 9.17 seconds (=96.5%)
339#   non-zero timer [__ sub-timer]:
340#   newtonIncrement   = 2.45%
341#   integrationFormula= 4.94%
342#   ODE2RHS           = 87.5%
343#   writeSolution     = 3.6%
344#   overhead          = 1.54%
345#   visualization/user= 0.00405%
346# special timers:
347#   Contact:BoundingBoxes = 0.92908 (10.1%)s
348#   Contact:SearchTree = 2.3468 (25.6%)s
349#   Contact:Overall = 6.676 (72.8%)s
350
351
352
353#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
354#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
355#1.9.66, staticTriangles = True
356
357# n spheres= 8000 , ss= 16
358# +++++++++++++++++++++++++++++++
359# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
360# Start multi-threading with 8 threads
361# WARNING: VelocityVerlet: still under development
362# STEP526, t = 0.0526s, timeToGo = 5.60837s, Nit/step = 0
363# STEP1091, t = 0.1091s, timeToGo = 3.33273s, Nit/step = 0
364# STEP1655, t = 0.1655s, timeToGo = 1.25102s, Nit/step = 0
365# STEP2000, t = 0.2s, timeToGo = 2.0821e-13s, Nit/step = 0
366# Solver terminated successfully after 7.30395 seconds.
367# ====================
368# CPU-time statistics:
369#   total time   = 7.3 seconds
370#   measured time= 7.01 seconds (=96%)
371#   non-zero timer [__ sub-timer]:
372#   newtonIncrement   = 3.19%
373#   integrationFormula= 6.58%
374#   ODE2RHS           = 83.6%
375#   writeSolution     = 4.71%
376#   overhead          = 1.93%
377#   visualization/user= 0.003%
378# special timers:
379#   Contact:BoundingBoxes = 0.98863 (14.1%)s
380#   Contact:SearchTree = 0.4626 (6.6%)s
381#   Contact:Overall = 4.3393 (61.9%)s
382
383#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
384#update to 1.9.68 (exact bin computation for triangles)
385# start create: number of masses = 8000
386# create mass 0
387# finish gContact
388# n spheres= 8000 , ss= 32
389# +++++++++++++++++++++++++++++++
390# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
391# Start multi-threading with 8 threads
392# STEP639, t = 0.0639s, timeToGo = 4.26062s, Nit/step = 0
393# STEP1309, t = 0.1309s, timeToGo = 2.11214s, Nit/step = 0
394# STEP2000, t = 0.2s, timeToGo = 1.70989e-13s, Nit/step = 0
395# Solver terminated successfully after 5.99597 seconds.
396# ====================
397# CPU-time statistics:
398#   total time   = 6 seconds
399#   measured time= 5.72 seconds (=95.4%)
400#   non-zero timer [__ sub-timer]:
401#   newtonIncrement   = 3.65%
402#   integrationFormula= 7.59%
403#   ODE2RHS           = 81.1%
404#   writeSolution     = 5.57%
405#   overhead          = 2.07%
406#   visualization/user= 0.00279%
407# special timers:
408#   Contact:BoundingBoxes = 0.86635 (15.1%)s
409#   Contact:SearchTree = 0.59121 (10.3%)s
410#   Contact:Overall = 3.3974 (59.4%)s
411
412
413# runfile('C:/DATA/cpp/EXUDYN_git/main/pythonDev/Examples/particlesSilo.py', wdir='C:/DATA/cpp/EXUDYN_git/main/pythonDev/Examples')
414# start create: number of masses = 8000
415# create mass 0
416# finish gContact
417# n spheres= 8000 , ss= 48
418# +++++++++++++++++++++++++++++++
419# EXUDYN V1.9.68.dev1 solver: explicit time integration (VelocityVerlet)
420# Start multi-threading with 8 threads
421# STEP517, t = 0.0517s, timeToGo = 5.74121s, Nit/step = 0
422# STEP1051, t = 0.1051s, timeToGo = 3.614s, Nit/step = 0
423# STEP1601, t = 0.1601s, timeToGo = 1.49891s, Nit/step = 0
424# STEP2000, t = 0.2s, timeToGo = 2.14357e-13s, Nit/step = 0
425# Solver terminated successfully after 7.51458 seconds.
426# ====================
427# CPU-time statistics:
428#   total time   = 7.51 seconds
429#   measured time= 7.21 seconds (=95.9%)
430#   non-zero timer [__ sub-timer]:
431#   newtonIncrement   = 3.12%
432#   integrationFormula= 6.21%
433#   ODE2RHS           = 84.7%
434#   writeSolution     = 4.11%
435#   overhead          = 1.84%
436#   visualization/user= 0.00239%
437# special timers:
438#   Contact:BoundingBoxes = 0.88801 (12.3%)s
439#   Contact:SearchTree = 0.93017 (12.9%)s
440#   Contact:Overall = 4.8152 (66.8%)s
441
442#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
443#large test:
444# 500000 particles, stepSize=2.5e-5
445# STEP40046, t = 1.00115s, timeToGo = 5.81135 days, tCPU=2.84971h, Nit/step = 0
446# STEP121643, t = 3.04107s, timeToGo = 6.78003 days, tCPU=10.5378h, Nit/step = 0
447# STEP121644 (stopped), t = 3.04107s, tCPU=37936.2s, Nit/step = 0
448# solver stopped by user after 37935.8 seconds.
449# ====================
450# CPU-time statistics:
451#   total time   = 3.79e+04 seconds
452#   measured time= 3.63e+04 seconds (=95.8%)
453#   non-zero timer [__ sub-timer]:
454#   newtonIncrement   = 3.02%
455#   integrationFormula= 5.06%
456#   ODE2RHS           = 86.1%
457#   writeSolution     = 0.974%
458#   overhead          = 1.18%
459#   visualization/user= 3.61%
460# special timers:
461#   Contact:BoundingBoxes = 3156.8 (8.68%)s
462#   Contact:SearchTree = 3350.5 (9.22%)s
463#   Contact:Overall = 25633 (70.5%)s