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