parameterVariationExample.py
You can view and download this file on Github: parameterVariationExample.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, ...
6# The value computed in every parameter variation is the error compared to
7# a reference solution using reference/nominal values
8#
9# Author: Johannes Gerstmayr
10# Date: 2020-11-18
11#
12# 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.
13#
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16import exudyn as exu
17from exudyn.itemInterface import *
18from exudyn.processing import ParameterVariation
19
20import numpy as np #for postprocessing
21
22SC = exu.SystemContainer()
23mbs = SC.AddSystem()
24#this is the function which is repeatedly called from ParameterVariation
25#parameterSet contains dictinary with varied parameters
26def ParameterFunction(parameterSet):
27 global mbs
28 mbs.Reset()
29
30 #++++++++++++++++++++++++++++++++++++++++++++++
31 #++++++++++++++++++++++++++++++++++++++++++++++
32 #store default parameters in structure (all these parameters can be varied!)
33 class P: pass #create emtpy structure for parameters; simplifies way to update parameters
34
35 #default values
36 P.mass = 1.6 #mass in kg
37 P.spring = 4000 #stiffness of spring-damper in N/m
38 P.damper = 8 #damping constant in N/(m/s)
39 P.u0=-0.08 #initial displacement
40 P.v0=1 #initial velocity
41 P.f =80 #force applied to mass
42 P.L=0.5 #spring length (for drawing)
43 P.computationIndex = 'Ref'
44
45 # #now update parameters with parameterSet (will work with any parameters in structure P)
46 for key,value in parameterSet.items():
47 setattr(P,key,value)
48
49 #++++++++++++++++++++++++++++++++++++++++++++++
50 #++++++++++++++++++++++++++++++++++++++++++++++
51 #START HERE: create parameterized model, using structure P, which is updated in every computation
52
53 x0=P.f/P.spring #static displacement
54
55 # print('resonance frequency = '+str(np.sqrt(spring/mass)))
56 # print('static displacement = '+str(x0))
57
58 #node for 3D mass point:
59 n1=mbs.AddNode(Point(referenceCoordinates = [P.L,0,0],
60 initialCoordinates = [P.u0,0,0],
61 initialVelocities= [P.v0,0,0]))
62
63 #ground node
64 nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
65
66 #add mass point (this is a 3D object with 3 coordinates):
67 massPoint = mbs.AddObject(MassPoint(physicsMass = P.mass, nodeNumber = n1))
68
69 #marker for ground (=fixed):
70 groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
71 #marker for springDamper for first (x-)coordinate:
72 nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
73
74 #spring-damper between two marker coordinates
75 nC = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker],
76 stiffness = P.spring, damping = P.damper))
77
78 #add load:
79 mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
80 load = P.f))
81 #add sensor:
82 #not needed, if file not written:
83 fileName = ''
84 if P.computationIndex == 'Ref':
85 fileName = 'solution/paramVarDisplacementRef.txt'
86 sForce = mbs.AddSensor(SensorObject(objectNumber=nC, fileName=fileName,
87 storeInternal = True,
88 outputVariableType=exu.OutputVariableType.Force))
89
90 #print(mbs)
91 mbs.Assemble()
92
93 steps = 1000 #number of steps to show solution
94 tEnd = 1 #end time of simulation
95
96 simulationSettings = exu.SimulationSettings()
97 #simulationSettings.solutionSettings.solutionWritePeriod = 5e-3 #output interval general
98 simulationSettings.solutionSettings.writeSolutionToFile = False
99 simulationSettings.solutionSettings.sensorsWritePeriod = 5e-3 #output interval of sensors
100 simulationSettings.timeIntegration.numberOfSteps = steps
101 simulationSettings.timeIntegration.endTime = tEnd
102
103 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #no damping
104
105 #exu.StartRenderer() #start graphics visualization
106 #mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
107
108 #start solver:
109 mbs.SolveDynamic(simulationSettings)
110
111 #SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
112 #exu.StopRenderer() #safely close rendering window!
113
114 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
115 #evaluate difference between reference and optimized solution
116 #reference solution:
117 dataRef = np.loadtxt('solution/paramVarDisplacementRef.txt', comments='#', delimiter=',')
118 #data = np.loadtxt(fileName, comments='#', delimiter=',')
119 data = mbs.GetSensorStoredData(sForce)
120 diff = data[:,1]-dataRef[:,1]
121
122 errorNorm = np.sqrt(np.dot(diff,diff))/steps*tEnd
123
124 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
125 #compute exact solution:
126 if False:
127 from matplotlib import plt
128
129 plt.close('all')
130 plt.plot(data[:,0], data[:,1], 'b-', label='displacement (m)')
131
132 ax=plt.gca() # get current axes
133 ax.grid(True, 'major', 'both')
134 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
135 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
136 plt.legend() #show labels as legend
137 plt.tight_layout()
138 plt.show()
139
140 return errorNorm
141
142
143#for mpi parallelization see below
144#now perform parameter variation
145if __name__ == '__main__': #include this to enable parallel processing
146 import time
147
148 refval = ParameterFunction({}) # compute reference solution
149 #print("refval =", refval)
150
151 n = 16
152 start_time = time.time()
153 [pDict, values] = ParameterVariation(parameterFunction = ParameterFunction,
154 parameters = {'mass':(1,2,n),
155 'spring':(2000,8000,n),
156 #'test':(1,3,4)
157 },
158 debugMode = False,
159 addComputationIndex = True,
160 useMultiProcessing = True,
161 showProgress = True,
162 )
163
164 print("--- %s seconds ---" % (time.time() - start_time))
165 print('values[-1]=', values[-1]) # values[-1] = 3.8418270115351496
166
167 from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import
168 import matplotlib.pyplot as plt
169 from matplotlib import colormaps
170 import numpy as np
171 colorMap = colormaps.get_cmap('jet') #finite element colors
172
173 plt.close('all')
174 fig = plt.figure()
175 ax = fig.add_subplot(projection='3d')
176
177 #reshape output of parametervariation to fit plot_surface
178 X = np.array(pDict['mass']).reshape((n,n))
179 Y = np.array(pDict['spring']).reshape((n,n))
180 Z = np.array(values).reshape((n,n))
181
182 surf = ax.plot_surface(X, Y, Z,
183 cmap=colorMap, linewidth=2,
184 antialiased=True,
185 shade = True)
186 plt.colorbar(surf, shrink=0.5, aspect=5)
187 plt.tight_layout()
188
189 #++++++++++++++++++++++++++++++++++++++++++++++++++
190 #now add a refined parameter variation
191 #visualize results with scatter plot
192 [pDict2, values2] = ParameterVariation(parameterFunction = ParameterFunction,
193 parameters = {'mass':(1.5,1.7,n), 'spring':(3000,5000,n)},
194 debugMode = False,
195 addComputationIndex = True,
196 useMultiProcessing = True,
197 showProgress = True,
198 )
199
200 print('values2[-1]=', values2[-1]) # values2[-1]=1.8943208246113492
201 fig = plt.figure()
202 ax = fig.add_subplot(projection='3d')
203
204 X = np.concatenate((pDict['mass'],pDict2['mass']))
205 Y = np.concatenate((pDict['spring'],pDict2['spring']))
206 Z = np.concatenate((values, values2))
207
208 #plt.scatter(pDict['mass'], pDict['spring'], values, c='b', marker='o')
209 ps = ax.scatter(X, Y, Z, c=Z, marker='o', cmap = colorMap)
210 plt.colorbar(ps)
211 plt.tight_layout()
212
213 plt.show()
214
215
216
217#for mpi parallelization use the following example: