velocityVerletTest.py
You can view and download this file on Github: velocityVerletTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for Velocity Verlet integrator;
5# Note: only second order accurate for conservative systems
6#
7# Author: Johannes Gerstmayr
8# Date: 2024-10-06
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.utilities import * #includes itemInterface and rigidBodyUtilities
16import exudyn
17import exudyn.graphics as graphics
18
19import numpy as np
20
21testIsActive = exu.sys.get('testIsActive', False)
22
23SC = exu.SystemContainer() #container of systems
24mbs = SC.AddSystem() #add a new system to work with
25
26mass = mbs.CreateMassPoint(referencePosition=[1,0,0],
27 initialVelocity=[1,0,0],
28 #gravity=[10,0,0],
29 mass=10,
30 graphicsDataList=[graphics.Sphere(radius=0.2,color=graphics.color.red)])
31ground = mbs.CreateGround()
32mbs.CreateSpringDamper(itemNumbers=[ground, mass], referenceLength=1,
33 stiffness = 1000, damping = 0,
34 drawSize = 0.1,
35 )
36sMass = mbs.AddSensor(SensorBody(bodyNumber = mass, storeInternal=True,
37 outputVariableType=exu.OutputVariableType.Position))
38
39mbs.Assemble() #assemble system and solve
40simulationSettings = exu.SimulationSettings()
41
42if not testIsActive: #only start graphics once, but after background is set
43 SC.renderer.Start()
44 SC.renderer.DoIdleTasks()
45
46sumSol = 0
47for i in range(4):
48 h = 0.001 / 2**i
49 tEnd = 2
50 simulationSettings.timeIntegration.endTime = tEnd
51 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
52 simulationSettings.timeIntegration.verboseMode=1 #provide some output
53 # simulationSettings.timeIntegration.realtime.active = True
54
55 #SC.renderer.DoIdleTasks()
56 mbs.SolveDynamic(simulationSettings, solverType = exudyn.DynamicSolverType.VelocityVerlet)
57 #mbs.SolveDynamic(simulationSettings, solverType = exudyn.DynamicSolverType.ExplicitEuler)
58 # mbs.SolveDynamic(simulationSettings, solverType = exudyn.DynamicSolverType.ExplicitMidpoint)
59 #mbs.SolveDynamic(simulationSettings, solverType = exudyn.DynamicSolverType.RK67)
60 p = mbs.GetSensorValues(sMass)
61 #ref solution with RK67 !
62 #err = p - [1.0332409398209812, 0.0, 0.0] #init vel
63 err = p - [1.091294525072764, 0.0, 0.0] #init vel, no damping
64 #err = p - [1.082490077681814, 0.0, 0.0] #gravity
65 #err = p - [1.059191793818665, 0.0, 0.0] #gravity, no damping
66 exu.Print('p =',list(p))
67 exu.Print('h=', round(h,8), ', err=', err[0])
68
69 sumSol += p[0]
70
71
72 #Euler:
73 # p = [1.0368839392961136, 0.0, 0.0]
74 # p = [1.035017478371429, 0.0, 0.0]
75 # p = [1.0341182324516134, 0.0, 0.0]
76 # p = [1.033676874325782, 0.0, 0.0]
77 # h= 0.001 , err= 0.003642999475132358
78 # h= 0.0005 , err= 0.0017765385504477926
79 # h= 0.00025 , err= 0.0008772926306321871
80 # h= 0.000125 , err= 0.00043593450480083895
81
82 #VelocityVerlet:
83 # p = [1.0333130215415933, 0.0, 0.0]
84 # p = [1.0332766497073687, 0.0, 0.0]
85 # p = [1.033258711414605, 0.0, 0.0]
86 # p = [1.0332498047052472, 0.0, 0.0]
87 # h= 0.001 , err= 7.208172061212714e-05
88 # h= 0.0005 , err= 3.570988638745831e-05
89 # h= 0.00025 , err= 1.777159362381653e-05
90 # h= 0.000125 , err= 8.86488426599108e-06
91
92 #ExplicitMidpoint:
93 # h= 0.001 , err= 3.6605276949597254e-06
94 # h= 0.0005 , err= 9.043581823409141e-07
95 # h= 0.00025 , err= 2.247182926407021e-07
96 # h= 0.000125 , err= 5.6006646875772503e-08
97
98 #RK67: p= [1.0332409398209816, 0.0, 0.0]
99 #RK67/2: p= [1.0332409398209812, 0.0, 0.0]
100if not testIsActive: #only start graphics once, but after background is set
101 SC.renderer.Stop()
102 mbs.PlotSensor(sMass,components=[0])
103
104exu.Print("velocityVerletTest result=",sumSol)
105
106exu.sys['testResult'] = sumSol