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