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