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