dispyParameterVariationExample.py

You can view and download this file on Github: dispyParameterVariationExample.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  This example performs a parameter variation of a simple
  5#           mass-spring-damper system; varying mass, spring, using the python module dispy
  6#           In this example, we return the computed sensor values to the ParameterVariation funtion
  7#           and visualize a set of results
  8#           NOTE: the real speedup using multiprocessing or cluster is about 3 on a 4-core machine, but not the value given by dispy!
  9#
 10# Author:   Stefan Holzinger, Johannes Gerstmayr
 11# Date:     2022-04-11
 12#
 13# 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.
 14#
 15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 16
 17import numpy as np
 18from exudyn.basicUtilities import OutputFilePath
 19
 20#this is the function that is called via dispy; parameterSet contains the dictionary of current parameters and 
 21#optional 'functionData' with additional data passed to this function (e.g. numpy arrays)
 22#parameters are written into the 'P' structure, such that they can be overwritten easily
 23def ParameterFunction(parameterSet): 
 24    # modules such as exudyn etc. need to be imported inside parameter function
 25    import socket
 26    import numpy as np
 27    import exudyn as exu
 28    import exudyn.itemInterface as eii
 29    
 30    SC = exu.SystemContainer()
 31    mbs = SC.AddSystem()
 32
 33    #++++++++++++++++++++++++++++++++++++++++++++++
 34    #++++++++++++++++++++++++++++++++++++++++++++++
 35    #store default parameters in structure (all these parameters can be varied!)
 36    class P: pass #create emtpy structure for parameters; simplifies way to update parameters
 37
 38    #default values
 39    P.mass = 1.6          #mass in kg
 40    P.spring = 4000       #stiffness of spring-damper in N/m
 41    P.damper = 8          #damping constant in N/(m/s)
 42    P.u0=-0.08            #initial displacement
 43    P.v0=1                #initial velocity
 44    P.f =80               #force applied to mass
 45    P.L=0.5               #spring length (for drawing)
 46    P.computationIndex = 'Ref' #if computationIndex not provided
 47    
 48    # #now update parameters with parameterSet (will work with any parameters in structure P)
 49    for key,value in parameterSet.items():
 50        setattr(P,key,value)
 51
 52     #++++++++++++++++++++++++++++++++++++++++++++++
 53     #++++++++++++++++++++++++++++++++++++++++++++++
 54     #START HERE: create parameterized model, using structure P, which is updated in every computation
 55        
 56    # +++++++++++++++++++++++++++++++++++++++
 57    # create mechanical system
 58        
 59    #mass-spring-damper system
 60    P.L=0.5               #spring length (for drawing)
 61    
 62    x0=P.f/P.spring         #static displacement
 63    
 64    #node for 3D mass point:
 65    n1=mbs.AddNode(eii.Point(referenceCoordinates = [P.L,0,0], 
 66                                initialCoordinates = [P.u0,0,0], 
 67                                initialVelocities= [P.v0,0,0]))
 68    
 69    
 70    #ground node
 71    nGround=mbs.AddNode(eii.NodePointGround(referenceCoordinates = [0,0,0]))
 72    
 73    #add mass point (this is a 3D object with 3 coordinates):
 74    massPoint = mbs.AddObject(eii.MassPoint(mass = P.mass, nodeNumber = n1))
 75    
 76    #marker for ground (=fixed):
 77    groundMarker=mbs.AddMarker(eii.MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
 78    #marker for springDamper for first (x-)coordinate:
 79    nodeMarker  =mbs.AddMarker(eii.MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
 80    
 81    #spring-damper between two marker coordinates
 82    nC = mbs.AddObject(eii.CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker], 
 83                                                  stiffness = P.spring, damping = P.damper)) 
 84    
 85    #add load:
 86    mbs.AddLoad(eii.LoadCoordinate(markerNumber = nodeMarker, load = P.f))
 87    
 88    #add sensor:
 89    fileName = 'solution/paramVarDisplacement'+str(P.computationIndex)+'.txt'
 90    mbs.AddSensor(eii.SensorObject(objectNumber=nC, fileName=fileName, 
 91                                outputVariableType=exu.OutputVariableType.Force))    
 92    
 93   
 94    #print(mbs)
 95    mbs.Assemble()
 96    
 97    steps = 1000000  #number of steps to show solution; use many steps to see speedup!
 98    tEnd = 1     #end time of simulation
 99    
100    simulationSettings = exu.SimulationSettings()
101    simulationSettings.solution.file.write = False
102    simulationSettings.solution.sensors.writePeriod = 5e-3  #output interval of sensors
103    simulationSettings.timeIntegration.numberOfSteps = steps
104    simulationSettings.timeIntegration.endTime = tEnd
105    
106    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #no damping
107    
108    #start solver:
109    mbs.SolveDynamic(simulationSettings)
110
111    # +++++++++++++++++++++++++++++++++++++++++++++++++++++
112    # evaluate difference between reference and optimized solution
113    
114    # get reference solution
115    if 'functionData' in parameterSet: # load reference solution from function data
116        functionData = parameterSet['functionData']
117        refSol = functionData['refSol'] 
118    else: # load from file
119        refSol = np.loadtxt(fileName, comments='#', delimiter=',')
120        
121    # computation solution from current simulation
122    sol = np.loadtxt(fileName, comments='#', delimiter=',')
123
124    # compute error
125    diff = refSol[:,1] - sol[:,1] 
126    errorNorm = np.sqrt(np.dot(diff,diff))/steps*tEnd
127        
128        
129    return [errorNorm, sol]  #we can also return any pickle-able data in return value, which is then stored in values of ParameterVariation
130    #return errorNorm #for debug: (socket.gethostname(), errorNorm)
131
132
133 
134#%%
135if __name__ == '__main__':
136    from exudyn.processing import ParameterVariation
137    import time
138
139    
140    #++++++++++++++++++++++++++++++++++++++++++++++++
141    # generate reference solution
142    refval = ParameterFunction({}) # compute reference solution
143    #print("refval =", refval)
144
145    # add reference solution to function data
146    referenceSolution = np.loadtxt(OutputFilePath('solution/paramVarDisplacementRef.txt'), comments='#', delimiter=',')
147    functionData = {'refSol': referenceSolution}
148
149
150    
151    #++++++++++++++++++++++++++++++++++++++++++++++++
152    # process inputs to form array with each parameter values etc. see ParameterVariation()
153    useCluster = True # if false, multiprocessing is used
154    if useCluster: # specify nodes by providing host names or IPv4-addresses
155        if False: #put your dispy hosts here
156            activeHosts = ['111.112.113.114'] # enter list of IPv4 addresses or host names here
157        else: #just use local host (for tests only!); 
158            import socket
159            hostName = socket.gethostname()
160            print('your hostname is: ',hostName)
161            activeHosts = [hostName] 
162        clusterHostNames = activeHosts
163        useDispyWebMonitor = 'useDispyWebMonitor' # if given to ParameterVariation(), a web browser is started which can be used to manage the cluster computation
164        
165    else:
166        clusterHostNames = []
167        useDispyWebMonitor = ''
168    
169    
170
171    #++++++++++++++++++++++++++++++++++++++++++++++++
172    # perform parameter variation
173    n = 2 # n*n = number of variations
174    start_time = time.time()
175    [pDict, values] = ParameterVariation(parameterFunction = ParameterFunction, 
176                                         parameters = {'mass':(1,2,n), 
177                                                       'spring':(2000,8000,n),
178                                                       },
179                                         parameterFunctionData = functionData,
180                                         debugMode = False,
181                                         addComputationIndex = True,
182                                         useMultiProcessing = True,
183                                         showProgress = True, 
184                                         clusterHostNames = clusterHostNames, #not that there is a significant overhead, thus cluster only makes sense for computations that take longer than 1 second
185                                         #useDispyWebMonitor = useDispyWebMonitor, #slows down a little and waits 5 seconds in the end to finish
186                                         )
187    print("--- %s seconds ---" % (time.time() - start_time))
188    
189   
190    
191    
192#%% post processing of all computed variations
193    if True:
194        import exudyn as exu
195        SC = exu.SystemContainer()
196        mbs = SC.AddSystem() #dummy mbs for PlotSensor (and in order to skip PlotSensor in test suite)
197        
198        #extract solution from return values
199        valuesPlot = []
200        for item in values: 
201            valuesPlot +=[item[0]]
202        
203        #plot first 8 results:
204        for i in range(len(values)):
205            mbs.PlotSensor(sensorNumbers=[values[i][1]], labels='result'+str(i), newFigure=(i==0),
206                       colorCodeOffset=i, fontSize=12, closeAll=(i==0))