.. _examples-parametervariationexample: **************************** parameterVariationExample.py **************************** You can view and download this file on Github: `parameterVariationExample.py `_ .. code-block:: python :linenos: #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # This is an EXUDYN example # # Details: This example performs a parameter variation of a simple # mass-spring-damper system; varying mass, spring, ... # The value computed in every parameter variation is the error compared to # a reference solution using reference/nominal values # # Author: Johannes Gerstmayr # Date: 2020-11-18 # # 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. # #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ import exudyn as exu from exudyn.itemInterface import * from exudyn.processing import ParameterVariation import numpy as np #for postprocessing SC = exu.SystemContainer() mbs = SC.AddSystem() #this is the function which is repeatedly called from ParameterVariation #parameterSet contains dictinary with varied parameters def ParameterFunction(parameterSet): global mbs mbs.Reset() #++++++++++++++++++++++++++++++++++++++++++++++ #++++++++++++++++++++++++++++++++++++++++++++++ #store default parameters in structure (all these parameters can be varied!) class P: pass #create emtpy structure for parameters; simplifies way to update parameters #default values P.mass = 1.6 #mass in kg P.spring = 4000 #stiffness of spring-damper in N/m P.damper = 8 #damping constant in N/(m/s) P.u0=-0.08 #initial displacement P.v0=1 #initial velocity P.f =80 #force applied to mass P.L=0.5 #spring length (for drawing) P.computationIndex = 'Ref' # #now update parameters with parameterSet (will work with any parameters in structure P) for key,value in parameterSet.items(): setattr(P,key,value) #++++++++++++++++++++++++++++++++++++++++++++++ #++++++++++++++++++++++++++++++++++++++++++++++ #START HERE: create parameterized model, using structure P, which is updated in every computation x0=P.f/P.spring #static displacement # print('resonance frequency = '+str(np.sqrt(spring/mass))) # print('static displacement = '+str(x0)) #node for 3D mass point: n1=mbs.AddNode(Point(referenceCoordinates = [P.L,0,0], initialCoordinates = [P.u0,0,0], initialVelocities= [P.v0,0,0])) #ground node nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0])) #add mass point (this is a 3D object with 3 coordinates): massPoint = mbs.AddObject(MassPoint(physicsMass = P.mass, nodeNumber = n1)) #marker for ground (=fixed): groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0)) #marker for springDamper for first (x-)coordinate: nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0)) #spring-damper between two marker coordinates nC = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker], stiffness = P.spring, damping = P.damper)) #add load: mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker, load = P.f)) #add sensor: #not needed, if file not written: fileName = '' if P.computationIndex == 'Ref': fileName = 'solution/paramVarDisplacementRef.txt' sForce = mbs.AddSensor(SensorObject(objectNumber=nC, fileName=fileName, storeInternal = True, outputVariableType=exu.OutputVariableType.Force)) #print(mbs) mbs.Assemble() steps = 1000 #number of steps to show solution tEnd = 1 #end time of simulation simulationSettings = exu.SimulationSettings() #simulationSettings.solutionSettings.solutionWritePeriod = 5e-3 #output interval general simulationSettings.solutionSettings.writeSolutionToFile = False simulationSettings.solutionSettings.sensorsWritePeriod = 5e-3 #output interval of sensors simulationSettings.timeIntegration.numberOfSteps = steps simulationSettings.timeIntegration.endTime = tEnd simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #no damping #exu.StartRenderer() #start graphics visualization #mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue #start solver: mbs.SolveDynamic(simulationSettings) #SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit #exu.StopRenderer() #safely close rendering window! #+++++++++++++++++++++++++++++++++++++++++++++++++++++ #evaluate difference between reference and optimized solution #reference solution: dataRef = np.loadtxt('solution/paramVarDisplacementRef.txt', comments='#', delimiter=',') #data = np.loadtxt(fileName, comments='#', delimiter=',') data = mbs.GetSensorStoredData(sForce) diff = data[:,1]-dataRef[:,1] errorNorm = np.sqrt(np.dot(diff,diff))/steps*tEnd #+++++++++++++++++++++++++++++++++++++++++++++++++++++ #compute exact solution: if False: from matplotlib import plt plt.close('all') plt.plot(data[:,0], data[:,1], 'b-', label='displacement (m)') ax=plt.gca() # get current axes ax.grid(True, 'major', 'both') ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) plt.legend() #show labels as legend plt.tight_layout() plt.show() return errorNorm #for mpi parallelization see below #now perform parameter variation if __name__ == '__main__': #include this to enable parallel processing import time refval = ParameterFunction({}) # compute reference solution #print("refval =", refval) n = 16 start_time = time.time() [pDict, values] = ParameterVariation(parameterFunction = ParameterFunction, parameters = {'mass':(1,2,n), 'spring':(2000,8000,n), #'test':(1,3,4) }, debugMode = False, addComputationIndex = True, useMultiProcessing = True, showProgress = True, ) print("--- %s seconds ---" % (time.time() - start_time)) print('values[-1]=', values[-1]) # values[-1] = 3.8418270115351496 from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import import matplotlib.pyplot as plt from matplotlib import colormaps import numpy as np colorMap = colormaps.get_cmap('jet') #finite element colors plt.close('all') fig = plt.figure() ax = fig.add_subplot(projection='3d') #reshape output of parametervariation to fit plot_surface X = np.array(pDict['mass']).reshape((n,n)) Y = np.array(pDict['spring']).reshape((n,n)) Z = np.array(values).reshape((n,n)) surf = ax.plot_surface(X, Y, Z, cmap=colorMap, linewidth=2, antialiased=True, shade = True) plt.colorbar(surf, shrink=0.5, aspect=5) plt.tight_layout() #++++++++++++++++++++++++++++++++++++++++++++++++++ #now add a refined parameter variation #visualize results with scatter plot [pDict2, values2] = ParameterVariation(parameterFunction = ParameterFunction, parameters = {'mass':(1.5,1.7,n), 'spring':(3000,5000,n)}, debugMode = False, addComputationIndex = True, useMultiProcessing = True, showProgress = True, ) print('values2[-1]=', values2[-1]) # values2[-1]=1.8943208246113492 fig = plt.figure() ax = fig.add_subplot(projection='3d') X = np.concatenate((pDict['mass'],pDict2['mass'])) Y = np.concatenate((pDict['spring'],pDict2['spring'])) Z = np.concatenate((values, values2)) #plt.scatter(pDict['mass'], pDict['spring'], values, c='b', marker='o') ps = ax.scatter(X, Y, Z, c=Z, marker='o', cmap = colorMap) plt.colorbar(ps) plt.tight_layout() plt.show() #for mpi parallelization use the following example: