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.graphics as graphics #only import if it does not conflict
 17
 18import numpy as np
 19
 20useGraphics = True #without test
 21#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 22#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 23try: #only if called from test suite
 24    from modelUnitTests import exudynTestGlobals #for globally storing test results
 25    useGraphics = exudynTestGlobals.useGraphics
 26except:
 27    class ExudynTestGlobals:
 28        pass
 29    exudynTestGlobals = ExudynTestGlobals()
 30#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 31
 32SC = exu.SystemContainer()         #container of systems
 33mbs = SC.AddSystem()               #add a new system to work with
 34
 35mass = mbs.CreateMassPoint(referencePosition=[1,0,0],
 36                           initialVelocity=[1,0,0],
 37                           #gravity=[10,0,0],
 38                           physicsMass=10,
 39                           graphicsDataList=[graphics.Sphere(radius=0.2,color=graphics.color.red)])
 40ground = mbs.CreateGround()
 41mbs.CreateSpringDamper(bodyNumbers=[ground, mass], referenceLength=1,
 42                       stiffness = 1000, damping = 0,
 43                       drawSize = 0.1,
 44                       )
 45sMass = mbs.AddSensor(SensorBody(bodyNumber = mass, storeInternal=True,
 46                                 outputVariableType=exu.OutputVariableType.Position))
 47
 48mbs.Assemble()                     #assemble system and solve
 49simulationSettings = exu.SimulationSettings()
 50
 51if useGraphics: #only start graphics once, but after background is set
 52    exu.StartRenderer()
 53    mbs.WaitForUserToContinue()
 54
 55sumSol = 0
 56for i in range(4):
 57    h = 0.001 / 2**i
 58    tEnd = 2
 59    simulationSettings.timeIntegration.endTime = tEnd
 60    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
 61    simulationSettings.timeIntegration.verboseMode=1 #provide some output
 62    # simulationSettings.timeIntegration.simulateInRealtime = True
 63
 64    #mbs.WaitForUserToContinue()
 65    mbs.SolveDynamic(simulationSettings, solverType = exudyn.DynamicSolverType.VelocityVerlet)
 66    #mbs.SolveDynamic(simulationSettings, solverType = exudyn.DynamicSolverType.ExplicitEuler)
 67    # mbs.SolveDynamic(simulationSettings, solverType = exudyn.DynamicSolverType.ExplicitMidpoint)
 68    #mbs.SolveDynamic(simulationSettings, solverType = exudyn.DynamicSolverType.RK67)
 69    p = mbs.GetSensorValues(sMass)
 70    #ref solution with RK67 !
 71    #err = p - [1.0332409398209812, 0.0, 0.0] #init vel
 72    err = p - [1.091294525072764, 0.0, 0.0] #init vel, no damping
 73    #err = p - [1.082490077681814, 0.0, 0.0]  #gravity
 74    #err = p - [1.059191793818665, 0.0, 0.0]  #gravity, no damping
 75    exu.Print('p  =',list(p))
 76    exu.Print('h=', round(h,8), ', err=', err[0])
 77
 78    sumSol += p[0]
 79
 80
 81    #Euler:
 82    # p  = [1.0368839392961136, 0.0, 0.0]
 83    # p  = [1.035017478371429, 0.0, 0.0]
 84    # p  = [1.0341182324516134, 0.0, 0.0]
 85    # p  = [1.033676874325782, 0.0, 0.0]
 86    # h= 0.001 , err= 0.003642999475132358
 87    # h= 0.0005 , err= 0.0017765385504477926
 88    # h= 0.00025 , err= 0.0008772926306321871
 89    # h= 0.000125 , err= 0.00043593450480083895
 90
 91    #VelocityVerlet:
 92    # p  = [1.0333130215415933, 0.0, 0.0]
 93    # p  = [1.0332766497073687, 0.0, 0.0]
 94    # p  = [1.033258711414605, 0.0, 0.0]
 95    # p  = [1.0332498047052472, 0.0, 0.0]
 96    # h= 0.001 , err= 7.208172061212714e-05
 97    # h= 0.0005 , err= 3.570988638745831e-05
 98    # h= 0.00025 , err= 1.777159362381653e-05
 99    # h= 0.000125 , err= 8.86488426599108e-06
100
101    #ExplicitMidpoint:
102    # h= 0.001 , err= 3.6605276949597254e-06
103    # h= 0.0005 , err= 9.043581823409141e-07
104    # h= 0.00025 , err= 2.247182926407021e-07
105    # h= 0.000125 , err= 5.6006646875772503e-08
106
107    #RK67:    p= [1.0332409398209816, 0.0, 0.0]
108    #RK67/2:  p= [1.0332409398209812, 0.0, 0.0]
109if useGraphics: #only start graphics once, but after background is set
110    exu.StopRenderer()
111    mbs.PlotSensor(sMass,components=[0])
112
113exu.Print("velocityVerletTest result=",sumSol)
114
115exudynTestGlobals.testError = sumSol - (4.365184132226787) #2024-10-06
116exudynTestGlobals.testResult = sumSol
117exu.Print("velocityVerletTest error=", exudynTestGlobals.testError)