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
17#a local run must not open a plot window when windows are suppressed (#2477)
18if exu.special.userInterface.suppressPlots: import matplotlib; matplotlib.use('Agg')
19from exudyn.basicUtilities import OutputFilePath
20from exudyn.itemInterface import *
21from exudyn.processing import ParameterVariation
22
23import numpy as np #for postprocessing
24
25SC = exu.SystemContainer()
26mbs = SC.AddSystem()
27#this is the function which is repeatedly called from ParameterVariation
28#parameterSet contains dictinary with varied parameters
29def ParameterFunction(parameterSet):
30 global mbs
31 mbs.Reset()
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'
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 x0=P.f/P.spring #static displacement
57
58 # print('resonance frequency = '+str(np.sqrt(spring/mass)))
59 # print('static displacement = '+str(x0))
60
61 #node for 3D mass point:
62 n1=mbs.AddNode(Point(referenceCoordinates = [P.L,0,0],
63 initialCoordinates = [P.u0,0,0],
64 initialVelocities= [P.v0,0,0]))
65
66 #ground node
67 nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
68
69 #add mass point (this is a 3D object with 3 coordinates):
70 massPoint = mbs.AddObject(MassPoint(mass = P.mass, nodeNumber = n1))
71
72 #marker for ground (=fixed):
73 groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
74 #marker for springDamper for first (x-)coordinate:
75 nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
76
77 #spring-damper between two marker coordinates
78 nC = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker],
79 stiffness = P.spring, damping = P.damper))
80
81 #add load:
82 mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
83 load = P.f))
84 #add sensor:
85 #not needed, if file not written:
86 fileName = ''
87 if P.computationIndex == 'Ref':
88 fileName = 'solution/paramVarDisplacementRef.txt'
89 sForce = mbs.AddSensor(SensorObject(objectNumber=nC, fileName=fileName,
90 storeInternal = True,
91 outputVariableType=exu.OutputVariableType.Force))
92
93 #print(mbs)
94 mbs.Assemble()
95
96 steps = 1000 #number of steps to show solution
97 tEnd = 1 #end time of simulation
98
99 simulationSettings = exu.SimulationSettings()
100 #simulationSettings.solution.file.writePeriod = 5e-3 #output interval general
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 #SC.renderer.Start() #start graphics visualization
109 #SC.renderer.DoIdleTasks() #wait for pressing SPACE bar to continue
110
111 #start solver:
112 mbs.SolveDynamic(simulationSettings)
113
114 #SC.renderer.DoIdleTasks()#wait for pressing 'Q' to quit
115 #SC.renderer.Stop() #safely close rendering window!
116
117 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
118 #evaluate difference between reference and optimized solution
119 #reference solution:
120 dataRef = np.loadtxt(OutputFilePath('solution/paramVarDisplacementRef.txt'), comments='#', delimiter=',')
121 #data = np.loadtxt(fileName, comments='#', delimiter=',')
122 data = mbs.GetSensorStoredData(sForce)
123 diff = data[:,1]-dataRef[:,1]
124
125 errorNorm = np.sqrt(np.dot(diff,diff))/steps*tEnd
126
127 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
128 #compute exact solution:
129 if False:
130 from matplotlib import plt
131
132 plt.close('all')
133 plt.plot(data[:,0], data[:,1], 'b-', label='displacement (m)')
134
135 ax=plt.gca() # get current axes
136 ax.grid(True, 'major', 'both')
137 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
138 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
139 plt.legend() #show labels as legend
140 plt.tight_layout()
141 plt.show()
142
143 return errorNorm
144
145
146#for mpi parallelization see below
147#now perform parameter variation
148if __name__ == '__main__': #include this to enable parallel processing
149 import time
150
151 refval = ParameterFunction({}) # compute reference solution
152 #print("refval =", refval)
153
154 n = 16
155 start_time = time.time()
156 [pDict, values] = ParameterVariation(parameterFunction = ParameterFunction,
157 parameters = {'mass':(1,2,n),
158 'spring':(2000,8000,n),
159 #'test':(1,3,4)
160 },
161 debugMode = False,
162 addComputationIndex = True,
163 useMultiProcessing = True,
164 showProgress = True,
165 )
166
167 print("--- %s seconds ---" % (time.time() - start_time))
168 print('values[-1]=', values[-1]) # values[-1] = 3.8418270115351496
169
170 from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import
171 import matplotlib.pyplot as plt
172 from matplotlib import colormaps
173 import numpy as np
174 colorMap = colormaps.get_cmap('jet') #finite element colors
175
176 plt.close('all')
177 fig = plt.figure()
178 ax = fig.add_subplot(projection='3d')
179
180 #reshape output of parametervariation to fit plot_surface
181 X = np.array(pDict['mass']).reshape((n,n))
182 Y = np.array(pDict['spring']).reshape((n,n))
183 Z = np.array(values).reshape((n,n))
184
185 surf = ax.plot_surface(X, Y, Z,
186 cmap=colorMap, linewidth=2,
187 antialiased=True,
188 shade = True)
189 plt.colorbar(surf, shrink=0.5, aspect=5)
190 plt.tight_layout()
191
192 #++++++++++++++++++++++++++++++++++++++++++++++++++
193 #now add a refined parameter variation
194 #visualize results with scatter plot
195 [pDict2, values2] = ParameterVariation(parameterFunction = ParameterFunction,
196 parameters = {'mass':(1.5,1.7,n), 'spring':(3000,5000,n)},
197 debugMode = False,
198 addComputationIndex = True,
199 useMultiProcessing = True,
200 showProgress = True,
201 )
202
203 print('values2[-1]=', values2[-1]) # values2[-1]=1.8943208246113492
204 fig = plt.figure()
205 ax = fig.add_subplot(projection='3d')
206
207 X = np.concatenate((pDict['mass'],pDict2['mass']))
208 Y = np.concatenate((pDict['spring'],pDict2['spring']))
209 Z = np.concatenate((values, values2))
210
211 #plt.scatter(pDict['mass'], pDict['spring'], values, c='b', marker='o')
212 ps = ax.scatter(X, Y, Z, c=Z, marker='o', cmap = colorMap)
213 plt.colorbar(ps)
214 plt.tight_layout()
215
216 plt.show()
217
218
219
220#for mpi parallelization use the following example: