solverExplicitODE1ODE2test.py

You can view and download this file on Github: solverExplicitODE1ODE2test.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  test model for ODE1 and ODE2 coordinates
  5#           test explicit integrators
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2020-06-19
  9#
 10# 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.
 11#
 12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 13
 14import exudyn as exu
 15from exudyn.itemInterface import *
 16
 17import numpy as np
 18
 19useGraphics = True #without test
 20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 21#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 22try: #only if called from test suite
 23    from modelUnitTests import exudynTestGlobals #for globally storing test results
 24    useGraphics = exudynTestGlobals.useGraphics
 25except:
 26    class ExudynTestGlobals:
 27        pass
 28    exudynTestGlobals = ExudynTestGlobals()
 29#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 30
 31SC = exu.SystemContainer()
 32mbs = SC.AddSystem()
 33exu.Print('EXUDYN version='+exu.GetVersionString())
 34
 35tEnd = 1     #end time of simulation
 36steps = 100
 37
 38if True:
 39    m=2 #mass
 40    nODE1 = mbs.AddNode(NodeGenericODE1(referenceCoordinates=[0,0,0,0],
 41                                        initialCoordinates=[1,0,0,0],
 42                                        numberOfODE1Coordinates=4))
 43
 44    #build system matrix and force vector
 45    A = np.array([[0,0,1,0],
 46                  [0,0,0,1],
 47                  [-200/m,0,0,0],
 48                  [0,-100/m,0,0]])
 49    b = np.array([0,0,0,1/m]) #includes m!
 50
 51    oGenericODE1 = mbs.AddObject(ObjectGenericODE1(nodeNumbers=[nODE1],
 52                                                   systemMatrix=A,
 53                                                   rhsVector=b))
 54
 55    sCoords1 = mbs.AddSensor(SensorNode(nodeNumber = nODE1,
 56                                       fileName='solution/testODE1.txt',
 57                                       writeToFile = False,
 58                                       outputVariableType=exu.OutputVariableType.Coordinates))
 59
 60nODE2 = mbs.AddNode(NodeGenericODE2(referenceCoordinates=[0,0],
 61                                    initialCoordinates=[1,0],
 62                                    initialCoordinates_t=[0,0],
 63                                    numberOfODE2Coordinates=2))
 64
 65#build system matrices and force vector
 66M = np.diag([m,m])
 67K = np.diag([200,100])
 68f = np.array([0,1])
 69
 70oGenericODE2 = mbs.AddObject(ObjectGenericODE2(nodeNumbers=[nODE2],
 71                                               massMatrix=M,stiffnessMatrix=K,forceVector=f,
 72                                               visualization=VObjectGenericODE2(show=False)))
 73
 74sCoords2 = mbs.AddSensor(SensorNode(nodeNumber = nODE2,
 75                                    storeInternal=True,#fileName='solution/testODE2.txt',
 76                                    writeToFile = False,
 77                                    outputVariableType=exu.OutputVariableType.Coordinates))
 78sCoords2_t = mbs.AddSensor(SensorNode(nodeNumber = nODE2,
 79                                    storeInternal=True,#fileName='solution/testODE2_t.txt',
 80                                    writeToFile = False,
 81                                    outputVariableType=exu.OutputVariableType.Coordinates_t))
 82
 83
 84#assemble and solve system for default parameters
 85mbs.Assemble()
 86
 87# exu.Print(mbs.systemData.GetObjectLTGODE1(0))
 88# exu.Print(mbs.systemData.GetObjectLTGODE2(1))
 89
 90sims=exu.SimulationSettings()
 91tEnd = 2 #2000000 steps in 1.28s on Python3.7 64bits
 92sims.timeIntegration.endTime = tEnd
 93sims.solutionSettings.writeSolutionToFile = False
 94sims.solutionSettings.sensorsWritePeriod = 10
 95sims.timeIntegration.verboseMode = 0
 96
 97
 98err = 0
 99ref = 0.40808206181339224 #reference value computed with RK67, h=5e-4
100
101fullConvergence = False #compute larger range of step size and tolerances
102offset = 1
103if fullConvergence:
104    offset = 0
105printResults = False
106
107if True: #check automatic step size control
108    #sims.timeIntegration.verboseMode = 1
109    for i in range(6-2*offset):
110        #sims.solutionSettings.writeSolutionToFile = True
111        #sims.solutionSettings.solutionWritePeriod = 0
112        tEnd = 2
113        h=1
114        sims.timeIntegration.numberOfSteps = int(tEnd/h)
115        sims.timeIntegration.endTime = tEnd
116        #sims.timeIntegration.initialStepSize = 1e-5
117
118        sims.timeIntegration.absoluteTolerance = 10**(-2*i)
119        sims.timeIntegration.relativeTolerance = 1e-4
120
121        solverType = exu.DynamicSolverType.DOPRI5
122        mbs.SolveDynamic(solverType=solverType, simulationSettings=sims)
123        if printResults:
124            exu.Print(str(solverType)+", h="+str(h)+", tol="+str(sims.timeIntegration.absoluteTolerance),
125                  ", val=", mbs.GetSensorValues(sCoords1)[0]-ref)
126        #exu.Print(mbs.GetSensorValues(sCoords1),mbs.GetSensorValues(sCoords2))
127        err += abs(mbs.GetSensorValues(sCoords1)[0]-ref) + abs(mbs.GetSensorValues(sCoords2)[0]-ref)
128
129if True: #check orders:
130    if printResults:
131        exu.Print("RK67:")
132    for i in range(3-offset):
133        h=1e-1*10**-i
134        sims.timeIntegration.numberOfSteps = int(tEnd/h)
135        mbs.SolveDynamic(solverType=exu.DynamicSolverType.RK67, simulationSettings=sims)
136        err += abs(mbs.GetSensorValues(sCoords1)[0]-ref) + abs(mbs.GetSensorValues(sCoords2)[0]-ref)
137        if printResults:
138            exu.Print("h=10**-"+str(i+1), ", val=", mbs.GetSensorValues(sCoords1)[0]-ref)
139    h=5e-4
140    sims.timeIntegration.numberOfSteps = int(tEnd/h)
141    mbs.SolveDynamic(solverType=exu.DynamicSolverType.RK67, simulationSettings=sims)
142    if printResults:
143        exu.Print("h=5e-4  ",", val=", mbs.GetSensorValues(sCoords1)[0]-ref)
144
145    if printResults:
146        exu.Print("RK44:")
147    for i in range(4-offset):
148        h=1e-1*10**-i
149        sims.timeIntegration.numberOfSteps = int(tEnd/h)
150        mbs.SolveDynamic(solverType=exu.DynamicSolverType.RK44, simulationSettings=sims)
151        err += abs(mbs.GetSensorValues(sCoords1)[0]-ref) + abs(mbs.GetSensorValues(sCoords2)[0]-ref)
152        if printResults:
153            exu.Print("h=10**-"+str(i+1), ", val=", mbs.GetSensorValues(sCoords1)[0]-ref)
154
155    if printResults:
156        exu.Print("RK33:")
157    for i in range(5-offset*2):
158        h=1e-1*10**-i
159        sims.timeIntegration.numberOfSteps = int(tEnd/h)
160        mbs.SolveDynamic(solverType=exu.DynamicSolverType.RK33, simulationSettings=sims)
161        err += abs(mbs.GetSensorValues(sCoords1)[0]-ref) + abs(mbs.GetSensorValues(sCoords2)[0]-ref)
162        if printResults:
163            exu.Print("h=10**-"+str(i+1), ", val=", mbs.GetSensorValues(sCoords1)[0]-ref)
164
165    if printResults:
166        exu.Print("RK22:")
167    for i in range(4-offset*2):
168        h=1e-2*10**-i
169        sims.timeIntegration.numberOfSteps = int(tEnd/h)
170        mbs.SolveDynamic(solverType=exu.DynamicSolverType.ExplicitMidpoint, simulationSettings=sims)
171        err += abs(mbs.GetSensorValues(sCoords1)[0]-ref) + abs(mbs.GetSensorValues(sCoords2)[0]-ref)
172        if printResults:
173            exu.Print("h=10**-"+str(i+1), ", val=", mbs.GetSensorValues(sCoords1)[0]-ref)
174
175    #OUTPUT for convergence tests (2021-01-24)
176    #RK67:
177    #h=10**-1 , val= -0.0024758038186170617
178    #h=10**-2 , val= -1.1607280747671922e-08
179    #h=10**-3 , val= -1.2934098236883074e-14
180    #h=5e-4   , val= 0.0
181    #RK44:
182    #h=10**-1 , val= 0.040740959793283626
183    #h=10**-2 , val= 1.4887435646093738e-05
184    #h=10**-3 , val= 1.515853220723784e-09
185    #h=10**-4 , val= 1.5693002453076588e-13
186    #RK33:
187    #h=10**-1 , val= -0.5131558622402683
188    #h=10**-2 , val= -0.0004058849181840518
189    #h=10**-3 , val= -3.461431334894627e-07
190    #h=10**-4 , val= -3.406751547530007e-10
191    #h=10**-5 , val= 5.202355213285159e-11
192    #RK22:
193    #h=10**-1 , val= -9.614173942611721
194    #h=10**-2 , val= -0.029910241095991663
195    #h=10**-3 , val= -0.0003033091724236603
196    #h=10**-4 , val= -3.0421319984763606e-06
197    #h=10**-5 , val= -3.037840295982974e-08
198
199
200#mbs.GetSensorValues(sCoords2),
201#h=1e-4, tEnd = 2:
202#Expl. Euler: [0.4121895  0.01004991]
203#Trapez. rule:[0.40808358 0.01004968]
204#h=1e-5, tEnd = 2:
205#Expl. Euler: [0.40849041 0.01004971]
206#Trapez. rule:[0.40808208 0.01004969]
207#h=1e-6, tEnd = 2:
208#Expl. Euler: [0.40812287 0.01004969]
209#RK4:         [0.40808206 0.01004969]
210#exu.Print("ODE1 end values=",mbs.GetSensorValues(sCoords1))
211
212
213exu.Print("solverExplicitODE1ODE2 err=",err)
214
215exudynTestGlobals.testError = err - (3.3767933275918964) #2021-01-25: 3.3767933275918964
216exudynTestGlobals.testResult = err
217
218#+++++++++++++++++++++++++++++++++++++++++++++++++++++
219if False:
220
221    mbs.PlotSensor([sCoords1],[1])
222    mbs.PlotSensor([sCoords2],[1])