symbolicUserFunctionTest.py
You can view and download this file on Github: symbolicUserFunctionTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Tests for symbolic user function
5#
6# Author: Johannes Gerstmayr
7# Date: 2023-11-28
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#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12import sys
13sys.exudynFast = True
14
15import exudyn as exu
16from exudyn.utilities import *
17
18useGraphics = False #without test
19#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
20#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
21try: #only if called from test suite
22 from modelUnitTests import exudynTestGlobals #for globally storing test results
23 useGraphics = exudynTestGlobals.useGraphics
24except:
25 class ExudynTestGlobals:
26 pass
27 exudynTestGlobals = ExudynTestGlobals()
28#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
29
30esym = exu.symbolic
31import numpy as np
32
33SC = exu.SystemContainer()
34mbs = SC.AddSystem()
35
36useSymbolicUF = True
37
38#variable can be used to switch behavior
39#esym.variables.Set('flag',1)
40
41if useSymbolicUF:
42 def springForceUserFunction(mbs, t, itemNumber, deltaL, deltaL_t, stiffness, damping, force):
43 #f0 = damping*deltaL_t + stiffness*deltaL + force #linear
44 #fact = esym.variables.Get('flag') #potential way to couple behaviour to external triggers
45 f0 = 10*damping*deltaL_t + stiffness*esym.sign(deltaL) * (esym.abs(deltaL))**1.2 + force
46 return f0
47
48 def UFload(mbs, t, load):
49 return load*esym.sin(10*(2*pi)*t)
50else:
51 def springForceUserFunction(mbs, t, itemNumber, deltaL, deltaL_t, stiffness, damping, force):
52 f0 = 10*damping*deltaL_t + stiffness*np.sign(deltaL) * (np.abs(deltaL))**1.2 + force
53 return f0
54
55 def UFload(mbs, t, load):
56 return load*np.sin(10*(2*pi)*t)
57
58# no user function:
59# UFload=0
60# springForceUserFunction=0
61
62oGround = mbs.CreateGround()
63
64oMassPoint = mbs.CreateMassPoint(referencePosition=[1.+0.05,0,0], physicsMass=1)
65co = mbs.CreateSpringDamper(bodyList=[oGround, oMassPoint],
66 referenceLength = 0.1, stiffness = 100,
67 damping = 1,
68 springForceUserFunction = springForceUserFunction,
69 )
70
71mc = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=0, coordinate=0))
72load = mbs.AddLoad(LoadCoordinate(markerNumber=mc, load=10,
73 loadUserFunction=UFload))
74
75if useSymbolicUF:
76 symbolicFunc = CreateSymbolicUserFunction(mbs, springForceUserFunction, 'springForceUserFunction', co)
77 mbs.SetObjectParameter(co, 'springForceUserFunction', symbolicFunc)
78
79 symFuncLoad = CreateSymbolicUserFunction(mbs, UFload, 'loadUserFunction', load)
80 mbs.SetLoadParameter(load, 'loadUserFunction', symFuncLoad)
81
82 #check function:
83 exu.Print('spring user function:',symbolicFunc)
84 exu.Print('load user function: ',symFuncLoad)
85
86
87#assemble and solve system for default parameters
88mbs.Assemble()
89
90endTime = 50
91stepSize = 0.005
92
93simulationSettings = exu.SimulationSettings()
94#simulationSettings.solutionSettings.solutionWritePeriod = 0.01
95simulationSettings.solutionSettings.writeSolutionToFile = False
96simulationSettings.timeIntegration.verboseMode = 1
97
98simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
99simulationSettings.timeIntegration.endTime = endTime
100simulationSettings.timeIntegration.newton.useModifiedNewton = True
101
102if useGraphics:
103 exu.StartRenderer()
104 # mbs.WaitForUserToContinue()
105
106import time
107ts = time.time()
108exu.Print('start simulation')
109mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.RK44)
110exu.Print('finished: ', time.time()-ts, 'seconds')
111
112if useGraphics:
113 exu.StopRenderer() #safely close rendering window!
114
115n = mbs.GetObject(oMassPoint)['nodeNumber']
116p = mbs.GetNodeOutput(n, exu.OutputVariableType.Position)
117u = NormL2(p)
118
119exu.Print('u=',u)
120exu.Print('solution of symbolicUserFunctionTest=',u)
121
122# result for 10000 steps; identical for both UF cases
123exudynTestGlobals.testError = u - (0.10039884426884882)
124exudynTestGlobals.testResult = u
125
126#++++++++++++++++++++++++++++++
127#i7-1390, boost
128#results for ExplicitMidpoint, 1e6 steps, best of three, exudynFast=True:
129#C++: 1.71s #1710ns/step
130#Python user function: 13.56s #
131#Symbolic user function: 2.28s #570ns overhead for 2 x user function
132# => speedup of user function part 20.8
133
134#++++++++++++++++++++++++++++++
135#OLD/original timings, one user function spring-damper
136#timings for 2e6 steps
137#i7-1390, power saving
138#results for ExplicitMidpoint, 2e6 steps, best of three, exudynFast=False:
139#C++: 2.99s
140#Python user function: 16.01s
141#Symbolic user function:4.37s
142
143#i7-1390, boost
144#results for ExplicitMidpoint, 2e6 steps, best of three, exudynFast=True:
145#C++: 1.14s #570ns / step
146#Python user function: 5.62s
147#Symbolic user function:1.34s #100ns overhead for user function
148# => speedup 22.4
149
150#i9
151#results for ExplicitMidpoint, 2e6 steps, best of three, exudynFast=True:
152#C++: 1.80s #570ns / step
153#Python user function: 9.52s
154#Symbolic user function:2.09s #100ns overhead for user function