taskmanagerTest.py
You can view and download this file on Github: taskmanagerTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Systematic tests with several types of models, object types, and constraints;
5# test for implicit, explicit and static solvers for number of threads in task manager
6#
7# Author: Johannes Gerstmayr
8# Date: 2024-10-08
9#
10# 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.
11#
12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
13
14import exudyn as exu
15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
16import exudyn.graphics as graphics #only import if it does not conflict
17
18import numpy as np
19
20useGraphics = True #without test
21#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
23try: #only if called from test suite
24 from modelUnitTests import exudynTestGlobals #for globally storing test results
25 useGraphics = exudynTestGlobals.useGraphics
26except:
27 class ExudynTestGlobals:
28 pass
29 exudynTestGlobals = ExudynTestGlobals()
30#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31
32useGraphics = False
33
34SC = exu.SystemContainer()
35mbs = SC.AddSystem()
36
37resultTotal = 0
38#general dimentions
39a = 1
40b = 0.5
41g = [0,-9.81,0]
42mass = 0.2
43stiffness = 1000
44damping = 2
45load = 2
46totalCases = 0
47
48#user function for spring force
49def springForce(mbs, t, itemIndex, u, v, k, d, offset): #changed 2023-01-21:, mu, muPropZone):
50 return k*u+k*u**3+v*d
51
52#user function for load
53def userLoadVector(mbs, t, loadVector):
54 return np.sin(4*pi*t)*np.array(loadVector)
55
56nMasses = 1
57caseConstraints = 0
58caseODE1 = 0
59caseUser = 0
60caseStatic = 0
61nThreads = 1
62
63for nMasses in [0, 1, 2, 4]:
64 for caseConstraints in range(2):
65 for caseStatic in range(2):
66 for caseUser in range(2):
67 for caseODE1 in range(2):
68 for nThreads in [1,2,5]:
69 mbs.Reset()
70
71 if caseStatic and caseODE1: continue
72 if nMasses == 0 and not caseODE1: continue
73
74 useExplicitSolver = (nMasses == 2 and caseConstraints == 0 and caseStatic == 0)
75
76 if not caseUser:
77 springForce = 0
78
79 case = [nThreads, caseODE1, nMasses, caseConstraints, caseUser,
80 caseStatic, int(useExplicitSolver)]
81
82 ground0 = mbs.CreateGround()
83 ground1 = mbs.CreateGround(referencePosition=[0,b,0])
84 bLast0 = ground0
85 bLast1 = ground1
86
87 for i in range(nMasses):
88 xOff = i*a
89 b0=mbs.CreateMassPoint(referencePosition=[xOff+a, 0, 0],
90 physicsMass=mass,
91 drawSize = 0.1*a, color=graphics.color.red,
92 create2D=True,
93 gravity = g)
94 b1=mbs.CreateMassPoint(referencePosition=[xOff+a, b, 0],
95 physicsMass=mass,
96 drawSize = 0.1*a, color=graphics.color.red,
97 create2D=True,
98 gravity = g)
99
100 mbs.CreateSpringDamper(bodyNumbers=[bLast1, b1], stiffness=stiffness, damping=damping,
101 springForceUserFunction=springForce)
102 if caseConstraints:
103 mbs.CreateDistanceConstraint(bodyNumbers=[bLast0, b0])
104 mbs.CreateDistanceConstraint(bodyNumbers=[bLast1, b0])
105 mbs.CreateDistanceConstraint(bodyNumbers=[b0, b1])
106 else:
107 mbs.CreateSpringDamper(bodyNumbers=[bLast0, b0], stiffness=stiffness, damping=damping)
108 mbs.CreateSpringDamper(bodyNumbers=[bLast1, b0], stiffness=stiffness, damping=damping)
109 mbs.CreateSpringDamper(bodyNumbers=[b0, b1], stiffness=stiffness, damping=damping)
110
111 bLast0 = b0
112 bLast1 = b1
113
114 if caseUser:
115 mbs.CreateForce(bodyNumber=bLast1,
116 loadVector=[10,20,0],
117 loadVectorUserFunction=userLoadVector)
118
119 sLast = mbs.AddSensor(SensorBody(bodyNumber=bLast0, storeInternal=True,
120 outputVariableType=exu.OutputVariableType.Displacement))
121
122 sODE1 = None
123 if caseODE1:
124 #set up a 2-DOF system
125 nODE1 = mbs.AddNode(NodeGenericODE1(referenceCoordinates=[0,0,1],
126 initialCoordinates=[1,0,0],
127 numberOfODE1Coordinates=3))
128
129 #build system matrix and force vector
130 #undamped mechanical system with m=1, K=100, f=1
131 #additionally 1 ODE1 coordinate
132 Amat = np.array([[0,1,0],
133 [-100,0,0],
134 [0,0,-10]])
135 bVec = np.array([0,1,50])
136
137 oGenericODE1 = mbs.AddObject(ObjectGenericODE1(nodeNumbers=[nODE1],
138 systemMatrix=Amat,
139 rhsVector=bVec))
140 sODE1 = mbs.AddSensor(SensorNode(nodeNumber=nODE1, storeInternal=True,
141 outputVariableType=exu.OutputVariableType.Coordinates))
142
143
144 mbs.Assemble()
145
146 tEnd = 0.1*20
147 h = 0.01
148 simulationSettings = exu.SimulationSettings()
149 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
150 simulationSettings.timeIntegration.endTime = tEnd
151 # simulationSettings.displayStatistics = True
152 # simulationSettings.displayComputationTime = True
153 simulationSettings.timeIntegration.verboseMode = 0
154 simulationSettings.staticSolver.verboseMode = 0
155
156 simulationSettings.parallel.numberOfThreads = nThreads
157
158 SC.visualizationSettings.nodes.drawNodesAsPoint = False
159 SC.visualizationSettings.nodes.tiling = 32
160 SC.visualizationSettings.window.alwaysOnTop = True
161
162 #start solver:
163 if useGraphics:
164 exu.StartRenderer()
165 mbs.WaitForUserToContinue()
166
167 if caseStatic:
168 mbs.SolveStatic(simulationSettings)
169 else:
170 if not useExplicitSolver:
171 mbs.SolveDynamic(simulationSettings)
172 else:
173 solverType = exu.DynamicSolverType.RK44
174 if caseUser and not caseODE1:
175 solverType = exu.DynamicSolverType.VelocityVerlet
176
177 mbs.SolveDynamic(simulationSettings,
178 solverType=solverType)
179
180 if useGraphics:
181 exu.StopRenderer()
182
183 #+++++++++++++++++++++++++++++++++++++
184 #evaluate results
185 x = mbs.GetSensorValues(sLast)
186 resultTotal += sum(x)
187 output = 'tip-disp='+str(list(x))
188
189 if caseODE1:
190 x = mbs.GetSensorValues(sODE1)
191 output += ', ODE1='+str(list(x))
192 resultTotal += sum(x)
193
194 if True:
195 exu.Print('** CASE '+str(totalCases)+' **:',case,
196 '; RESULTS:',output)
197
198 totalCases += 1
199
200#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
201#evaluate final (=current) output values
202resultTotal *= 0.001 #to avoid too large non-deterministic round-off errors due to multithreading
203exu.Print('\ntotal cases:', totalCases)
204exu.Print('result taskmanagerTest=', resultTotal)
205
206exudynTestGlobals.testResult = resultTotal #-0.17210771618100057
207
208# mbs.SolutionViewer()
209
210#%%++++++++++++++++
211if False:
212 mbs.PlotSensor(sLast, components=[1])
213 if sODE1 != None:
214 mbs.PlotSensor(sODE1, components=[0,1,2])