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                                         #numberOfThreads=2,
162                                         showProgress=True,
163                                         )
164
165    print("--- %s seconds ---" % (time.time() - start_time))
166    print('values[-1]=', values[-1]) # values[-1] = 3.8418270115351496
167
168    from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 unused import
169    import matplotlib.pyplot as plt
170    #from matplotlib import cm
171    #from matplotlib.ticker import LinearLocator, FormatStrFormatter
172    import numpy as np
173    colorMap = plt.cm.get_cmap('jet') #finite element colors
174
175    plt.close('all')
176    fig = plt.figure()
177    ax = fig.add_subplot(projection='3d')
178
179    #reshape output of parametervariation to fit plot_surface
180    X = np.array(pDict['mass']).reshape((n,n))
181    Y = np.array(pDict['spring']).reshape((n,n))
182    Z = np.array(values).reshape((n,n))
183
184    surf = ax.plot_surface(X, Y, Z,
185                           cmap=colorMap, linewidth=2,
186                           antialiased=True,
187                           shade = True)
188    plt.colorbar(surf, shrink=0.5, aspect=5)
189    plt.tight_layout()
190
191    #++++++++++++++++++++++++++++++++++++++++++++++++++
192    #now add a refined parameter variation
193    #visualize results with scatter plot
194    [pDict2, values2] = ParameterVariation(parameterFunction=ParameterFunction,
195                                         parameters={'mass':(1.5,1.7,n), 'spring':(3000,5000,n)},
196                                         debugMode=False,
197                                         addComputationIndex=True,
198                                         useMultiProcessing=True,
199                                         showProgress=True,
200                                         )
201
202    print('values2[-1]=', values2[-1]) # values2[-1]=1.8943208246113492
203    fig = plt.figure()
204    ax = fig.add_subplot(projection='3d')
205
206    X = np.concatenate((pDict['mass'],pDict2['mass']))
207    Y = np.concatenate((pDict['spring'],pDict2['spring']))
208    Z = np.concatenate((values, values2))
209
210    #plt.scatter(pDict['mass'], pDict['spring'], values, c='b', marker='o')
211    ps = ax.scatter(X, Y, Z, c=Z, marker='o', cmap = colorMap)
212    plt.colorbar(ps)
213    plt.tight_layout()
214
215    plt.show()
216
217
218
219#for mpi parallelization use the following example: