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; also those of items whose signature an earlier item has (#2664)
  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
 18
 19testIsActive = exu.sys.get('testIsActive', True)
 20
 21esym = exu.symbolic
 22import numpy as np
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27useSymbolicUF = True
 28
 29#variable can be used to switch behavior
 30#esym.variables.Set('flag',1)
 31
 32if useSymbolicUF:
 33    def springForceUserFunction(mbs, t, itemNumber, deltaL, deltaL_t, stiffness, damping, force):
 34        #f0 = damping*deltaL_t + stiffness*deltaL + force #linear
 35        #fact = esym.variables.Get('flag') #potential way to couple behaviour to external triggers
 36        f0 = 10*damping*deltaL_t + stiffness*esym.sign(deltaL) * (esym.abs(deltaL))**1.2 + force
 37        return f0
 38    
 39    def UFload(mbs, t, load):
 40        return load*esym.sin(10*(2*pi)*t)
 41else:
 42    def springForceUserFunction(mbs, t, itemNumber, deltaL, deltaL_t, stiffness, damping, force):
 43        f0 = 10*damping*deltaL_t + stiffness*np.sign(deltaL) * (np.abs(deltaL))**1.2 + force
 44        return f0
 45    
 46    def UFload(mbs, t, load):
 47        return load*np.sin(10*(2*pi)*t)
 48
 49# no user function:    
 50# UFload=0
 51# springForceUserFunction=0
 52
 53oGround = mbs.CreateGround()
 54
 55oMassPoint = mbs.CreateMassPoint(referencePosition=[1.+0.05,0,0], mass=1)
 56co = mbs.CreateSpringDamper(itemNumbers=[oGround, oMassPoint],
 57                            referenceLength = 0.1, stiffness = 100, 
 58                            damping = 1,
 59                            springForceUserFunction = springForceUserFunction,
 60                            )
 61
 62mc = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=0, coordinate=0))
 63load = mbs.AddLoad(LoadCoordinate(markerNumber=mc, load=10,
 64                                  loadUserFunction=UFload))
 65
 66if useSymbolicUF:
 67    symbolicFunc = CreateSymbolicUserFunction(mbs, springForceUserFunction, 'springForceUserFunction', co)
 68    mbs.SetObjectParameter(co, 'springForceUserFunction', symbolicFunc)
 69    
 70    symFuncLoad = CreateSymbolicUserFunction(mbs, UFload, 'loadUserFunction', load)
 71    mbs.SetLoadParameter(load, 'loadUserFunction', symFuncLoad)
 72    
 73    #check function:
 74    exu.Print('spring user function:',symbolicFunc)
 75    exu.Print('load user function:  ',symFuncLoad)
 76
 77
 78#assemble and solve system for default parameters
 79mbs.Assemble()
 80
 81endTime = 50
 82stepSize = 0.005
 83
 84simulationSettings = exu.SimulationSettings()
 85#simulationSettings.solution.file.writePeriod = 0.01
 86simulationSettings.solution.file.write = False
 87simulationSettings.timeIntegration.verboseMode = 1
 88
 89simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
 90simulationSettings.timeIntegration.endTime = endTime
 91simulationSettings.timeIntegration.newton.useModifiedNewton = True
 92
 93if not testIsActive:
 94    SC.renderer.Start()
 95    # SC.renderer.DoIdleTasks()
 96
 97import time
 98ts = time.time()
 99exu.Print('start simulation')
100mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.RK44)
101exu.Print('finished: ', time.time()-ts, 'seconds')
102
103if not testIsActive:
104    SC.renderer.Stop() #safely close rendering window!
105
106n = mbs.GetObject(oMassPoint)['nodeNumber']
107p = mbs.GetNodeOutput(n, exu.OutputVariableType.Position)
108u = np.linalg.norm(p)
109
110#every user function can be made symbolic, also one whose signature an earlier item has (#2664): the coordinate
111#spring-damper shares it with the spring-damper, the torque and the mass-proportional load with the force vector
112def UFcoordinate(mbs, t, itemNumber, displacement, velocity, stiffness, damping, offset):
113    return stiffness*(displacement-offset) + damping*velocity
114def UFvector(mbs, t, loadVector):
115    return [loadVector[0]*t, loadVector[1], loadVector[2]]
116symbolicCoordinate = CreateSymbolicUserFunction(mbs, UFcoordinate, 'springForceUserFunction',
117                                                itemTypeName='ObjectConnectorCoordinateSpringDamper')
118uErrors = abs(symbolicCoordinate.Evaluate(mbs, 0., 0, 0.1, 0.2, 100., 2., 0.05) - UFcoordinate(mbs, 0., 0, 0.1, 0.2, 100., 2., 0.05))
119for itemTypeName in ['LoadTorqueVector', 'LoadMassProportional']:
120    symbolicVector = CreateSymbolicUserFunction(mbs, UFvector, 'loadVectorUserFunction', itemTypeName=itemTypeName)
121    uErrors += np.linalg.norm(np.array(symbolicVector.Evaluate(mbs, 2., [1., 2., 3.])) - UFvector(mbs, 2., [1., 2., 3.]))
122exu.Print('symbolic user functions with a shared signature, errors:', uErrors)
123u += uErrors
124
125exu.Print('u=',u)
126exu.Print('solution of symbolicUserFunctionTest=',u)
127
128# result for 10000 steps; identical for both UF cases
129exu.sys['testResult'] = u
130
131#++++++++++++++++++++++++++++++
132#i7-1390, boost
133#results for ExplicitMidpoint, 1e6 steps, best of three, exudynFast=True:
134#C++:                    1.71s  #1710ns/step
135#Python user function:  13.56s  #
136#Symbolic user function: 2.28s  #570ns overhead for 2 x user function
137# => speedup of user function part 20.8
138
139#++++++++++++++++++++++++++++++
140#OLD/original timings, one user function spring-damper
141#timings for 2e6 steps
142#i7-1390, power saving
143#results for ExplicitMidpoint, 2e6 steps, best of three, exudynFast=False:
144#C++:                   2.99s
145#Python user function:  16.01s
146#Symbolic user function:4.37s
147
148#i7-1390, boost
149#results for ExplicitMidpoint, 2e6 steps, best of three, exudynFast=True:
150#C++:                   1.14s  #570ns / step
151#Python user function:  5.62s
152#Symbolic user function:1.34s  #100ns overhead for user function
153# => speedup 22.4
154
155#i9
156#results for ExplicitMidpoint, 2e6 steps, best of three, exudynFast=True:
157#C++:                   1.80s  #570ns / step
158#Python user function:  9.52s
159#Symbolic user function:2.09s  #100ns overhead for user function