explicitSolversPostNewtonTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Objects whose state is updated in the PostNewton step, solved with every explicit integrator
  5#           against the implicit generalized-alpha solver (#2754): a ball dropping on the ground through
  6#           ObjectContactCoordinate and ObjectContactSphereSphere (contact switches on and off), a body
  7#           turning several times against an ObjectConnectorTorsionalSpringDamper (the rotation counter)
  8#           and a mass under a periodic force with the stick-slip friction of
  9#           ObjectConnectorCoordinateSpringDamperExt. The explicit solvers perform the PostNewton step
 10#           after every step, through the discontinuous iteration all solvers share; measured: every
 11#           explicit integrator follows the implicit solution to the accuracy of its step size, e.g. the
 12#           deepest point of the drop 0.09639 for all, the counted rotation angle 0.2616 at the end.
 13#           Last, the drop with DOPRI5 and a large maximum step size (300 steps for 0.3 s), with and without
 14#           the step size the contacts recommend where they begin and end (#2109): the height after the
 15#           rebound 0.11615 with the recommendation for both contacts, against 0.11625 of the implicit solver
 16#           with 20 000 steps; without it the step size control alone gives 0.1147 to 0.1196.
 17#
 18# Author:   Johannes Gerstmayr
 19# Date:     2026-09-30
 20#
 21# 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.
 22#
 23#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 24
 25import exudyn as exu
 26from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 27
 28import numpy as np
 29
 30testIsActive = exu.sys.get('testIsActive', False)
 31
 32solverTypes = ['GeneralizedAlpha', 'ExplicitEuler', 'ExplicitMidpoint', 'RK33', 'RK44', 'RK67', 'ODE23',
 33               'DOPRI5', 'VelocityVerlet']
 34
 35
 36def Solve(mbs, endTime, numberOfSteps, solverType, useRecommendedStepSize=True):
 37    mbs.Assemble()
 38    simulationSettings = exu.SimulationSettings()
 39    simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
 40    simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = useRecommendedStepSize
 41    simulationSettings.timeIntegration.numberOfSteps = numberOfSteps
 42    simulationSettings.timeIntegration.endTime = endTime
 43    simulationSettings.solution.file.write = False
 44    simulationSettings.solution.sensors.writePeriod = endTime/numberOfSteps
 45    mbs.SolveDynamic(simulationSettings, solverType=getattr(exu.DynamicSolverType, solverType))
 46
 47
 48def Drop(kind, solverType, numberOfSteps=6000, useRecommendedStepSize=True, lastValue=False):
 49    """the deepest point of the center of a ball of radius 0.1 dropped from 0.2 onto the ground, or the
 50    height at the end"""
 51    r = 0.1; m = 1; g = 9.81; z0 = 0.2
 52    SC = exu.SystemContainer(); mbs = SC.AddSystem()
 53    oGround = mbs.AddObject(ObjectGround())
 54    if kind == 'coordinate':
 55        node = mbs.AddNode(Node1D(referenceCoordinates=[0], initialCoordinates=[z0]))
 56        mbs.AddObject(ObjectMass1D(nodeNumber=node, mass=m))
 57        mBall = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=node, coordinate=0))
 58        mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
 59        nData = mbs.AddNode(NodeGenericData(numberOfDataCoordinates=1, initialCoordinates=[z0 - r]))
 60        mbs.AddObject(ObjectContactCoordinate(markerNumbers=[mGround, mBall], nodeNumber=nData, offset=r,
 61                                              contactStiffness=1e5, contactDamping=100))
 62        mbs.AddLoad(LoadCoordinate(markerNumber=mBall, load=-m*g))
 63        sensor = mbs.AddSensor(SensorNode(nodeNumber=node, outputVariableType=exu.OutputVariableType.Coordinates,
 64                                          storeInternal=True, writeToFile=False))
 65    else:
 66        #Lie group node: the explicit integrators do not take the constraint of Euler parameters
 67        node = mbs.AddNode(NodeRigidBodyRotVecLG(referenceCoordinates=[0.1, 0.1, z0, 0, 0, 0]))
 68        ball = mbs.AddObject(ObjectRigidBody(nodeNumber=node, mass=m,
 69                                             inertia=InertiaSphere(mass=m, radius=r).GetInertia6D()))
 70        mBall = mbs.AddMarker(MarkerBodyRigid(bodyNumber=ball))
 71        mbs.AddLoad(LoadForceVector(markerNumber=mBall, loadVector=[0, 0, -m*g]))
 72        nData = mbs.AddNode(NodeGenericData(numberOfDataCoordinates=4, initialCoordinates=[0, 0, 0, 0]))
 73        mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0.1, 0.1, -1]))
 74        mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mGround, mBall], nodeNumber=nData,
 75                                                spheresRadii=[1, r], contactStiffness=1e5, contactDamping=100))
 76        sensor = mbs.AddSensor(SensorNode(nodeNumber=node, outputVariableType=exu.OutputVariableType.Position,
 77                                          storeInternal=True, writeToFile=False))
 78    Solve(mbs, 0.3, numberOfSteps, solverType, useRecommendedStepSize)
 79    if lastValue:
 80        return mbs.GetSensorStoredData(sensor)[-1, -1]
 81    return mbs.GetSensorStoredData(sensor)[:, -1].min()
 82
 83
 84def Torsional(solverType):
 85    """the rotation angle counted by the data node, after several turns back and forth"""
 86    SC = exu.SystemContainer(); mbs = SC.AddSystem()
 87    node = mbs.AddNode(NodeRigidBodyRotVecLG(referenceCoordinates=[0]*6, initialVelocities=[0, 0, 0, 0, 0, 20]))
 88    body = mbs.AddObject(ObjectRigidBody(nodeNumber=node, mass=1, inertia=[1, 1, 1, 0, 0, 0]))
 89    ground = mbs.AddObject(ObjectGround())
 90    nData = mbs.AddNode(NodeGenericData(numberOfDataCoordinates=1, initialCoordinates=[0]))
 91    mbs.AddObject(ObjectConnectorTorsionalSpringDamper(markerNumbers=[mbs.AddMarker(MarkerBodyRigid(bodyNumber=ground)),
 92                  mbs.AddMarker(MarkerBodyRigid(bodyNumber=body))], nodeNumber=nData, stiffness=10, damping=0))
 93    Solve(mbs, 2, 2000, solverType)
 94    return float(np.atleast_1d(mbs.GetNodeOutput(nData, exu.OutputVariableType.Coordinates))[0])
 95
 96
 97def StickSlip(solverType):
 98    """the position of a mass with stick-slip friction under a periodic force"""
 99    SC = exu.SystemContainer(); mbs = SC.AddSystem()
100    node = mbs.AddNode(Node1D(referenceCoordinates=[0]))
101    mbs.AddObject(ObjectMass1D(nodeNumber=node, mass=1))
102    mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
103    mMass = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=node, coordinate=0))
104    nData = mbs.AddNode(NodeGenericData(numberOfDataCoordinates=3, initialCoordinates=[0, 0, 0]))
105    mbs.AddObject(ObjectConnectorCoordinateSpringDamperExt(markerNumbers=[mGround, mMass], nodeNumber=nData, stiffness=0,
106                  damping=0, dynamicFrictionForce=2, staticFrictionOffsetForce=0.5, stickingStiffness=1e4,
107                  stickingDamping=100, frictionProportionalZone=0))
108    mbs.AddLoad(LoadCoordinate(markerNumber=mMass, load=0, loadUserFunction=lambda mbs, t, load: 3*np.sin(2*np.pi*t)))
109    Solve(mbs, 2, 4000, solverType)
110    return float(np.atleast_1d(mbs.GetNodeOutput(node, exu.OutputVariableType.Coordinates))[0])
111
112
113testResult = 0
114cases = [('drop, coordinate', lambda s: Drop('coordinate', s)), ('drop, sphere', lambda s: Drop('sphere', s)),
115         ('torsional', Torsional), ('stick-slip', StickSlip)]
116for (name, Case) in cases:
117    values = dict((solverType, Case(solverType)) for solverType in solverTypes)
118    reference = values['GeneralizedAlpha']
119    deviation = max(abs(v - reference) for v in values.values())
120    exu.Print(name + ': implicit', round(reference, 6), ', largest deviation of an explicit integrator', deviation)
121    testResult += sum(values.values())
122
123#DOPRI5 with a large maximum step size; without the recommended step size it did not finish (#2109)
124for kind in ['coordinate', 'sphere']:
125    for useRecommendedStepSize in [True, False]:
126        height = Drop(kind, 'DOPRI5', 300, useRecommendedStepSize, lastValue=True)
127        exu.Print('drop with DOPRI5, 300 steps,', kind, 'contact, recommended step size', useRecommendedStepSize,
128                  ': height at the end', round(height, 5))
129        testResult += height
130
131exu.Print('solution of explicitSolversPostNewtonTest=', testResult)
132exu.sys['testResult'] = testResult