restartFileTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  A simulation continued from its restart file computes what the uninterrupted one computes
  5#           (#2850). A double pendulum of rigid bodies with revolute joints (algebraic coordinates,
  6#           Lie group nodes) and a mass point in sphere-sphere contact (data coordinates)
  7#           runs to t = 1 in one go, and again stopped by a post-step function at t = 0.6 and started
  8#           a second time with solution.restart.continueIfAvailable, after the restart file of t = 0.6
  9#           was replaced by the one of t = 0.4 - as if the job had been killed after 0.4. The rows of
 10#           the solution and sensor files and the final state are identical to the uninterrupted run;
 11#           the same for the explicit RK44 solver, with springs instead of the joints. A restart file of another system raises; without a
 12#           file the simulation starts at its start time.
 13#
 14# Author:   Johannes Gerstmayr
 15# Date:     2026-10-05
 16#
 17# 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.
 18#
 19#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 20
 21import exudyn as exu
 22from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 23import numpy as np
 24import os
 25import shutil
 26
 27testIsActive = exu.sys.get('testIsActive', False)
 28
 29directory = 'solution/restartFileTest/'
 30errors = 0
 31def Check(condition, text):
 32    global errors
 33    if not condition:
 34        errors += 1
 35        exu.Print('restartFileTest failed:', text)
 36
 37def Model(explicit, addNode=False):
 38    SC = exu.SystemContainer()
 39    mbs = SC.AddSystem()
 40    oGround = mbs.CreateGround()
 41    inertia = InertiaCuboid(density=1000, sideLengths=[0.4,0.05,0.05])
 42    nodeType = exu.NodeType.RotationRotationVector if explicit else exu.NodeType.RotationEulerParameters
 43    b0 = mbs.CreateRigidBody(inertia=inertia, referencePosition=[0.2,0,0], gravity=[0,-9.81,0], nodeType=nodeType)
 44    b1 = mbs.CreateRigidBody(inertia=inertia, referencePosition=[0.6,0,0], gravity=[0,-9.81,0], initialAngularVelocity=[0,0,2],
 45                             nodeType=nodeType)
 46    if explicit: #an explicit solver takes no joints: stiff springs
 47        for (items, position) in [([oGround, b0], [-0.2,0,0]), ([b0, b1], [0.2,0,0])]:
 48            mbs.CreateCartesianSpringDamper(itemNumbers=items, localPosition1=position, stiffness=[1e4]*3, damping=[10]*3)
 49    else:
 50        mbs.CreateRevoluteJoint(itemNumbers=[oGround, b0], position=[0,0,0], axis=[0,0,1])
 51        mbs.CreateRevoluteJoint(itemNumbers=[b0, b1], position=[0.4,0,0], axis=[0,0,1])
 52    oBall = mbs.CreateMassPoint(referencePosition=[1,0.3,0], initialVelocity=[0,-1,0], mass=0.5, gravity=[0,-9.81,0])
 53    mbs.CreateSphereSphereContact(itemNumbers=[oGround, oBall], localPosition0=[1,-0.6,0], spheresRadii=[0.5,0.1],
 54                                  contactStiffness=1e4, contactDamping=20)
 55    if addNode:
 56        mbs.AddNode(NodePointGround())
 57    nBall = mbs.GetObject(oBall)['nodeNumber']
 58    mbs.AddSensor(SensorBody(bodyNumber=b1, localPosition=[0.2,0,0], fileName=directory+'sensorTip.txt',
 59                             outputVariableType=exu.OutputVariableType.Position))
 60    mbs.AddSensor(SensorNode(nodeNumber=nBall, fileName=directory+'sensorBall.txt',
 61                             outputVariableType=exu.OutputVariableType.Velocity))
 62    mbs.Assemble()
 63    return (SC, mbs)
 64
 65def Settings(solverType, solutionFile, stopTime=None, continueIfAvailable=False):
 66    simulationSettings = exu.SimulationSettings()
 67    simulationSettings.timeIntegration.numberOfSteps = 500
 68    simulationSettings.timeIntegration.endTime = 1
 69    simulationSettings.timeIntegration.verboseMode = 0
 70    simulationSettings.solution.file.name = directory + solutionFile
 71    simulationSettings.solution.file.writePeriod = 0.01
 72    simulationSettings.solution.sensors.writePeriod = 0.02
 73    simulationSettings.solution.precision = 16
 74    simulationSettings.solution.restart.write = True
 75    simulationSettings.solution.restart.name = directory + 'restart.txt'
 76    simulationSettings.solution.restart.writePeriod = 0.2
 77    simulationSettings.solution.restart.continueIfAvailable = continueIfAvailable
 78    simulationSettings.timeIntegration.newton.useModifiedNewton = False #the Jacobian of an earlier step is not in the restart file
 79    return simulationSettings
 80
 81def Run(solverType, solutionFile, stopTime=None, continueIfAvailable=False, addNode=False):
 82    (SC, mbs) = Model(solverType != exu.DynamicSolverType.GeneralizedAlpha, addNode)
 83    if stopTime is not None:
 84        mbs.SetPostStepUserFunction(lambda mbs, t: t < stopTime - 1e-10)
 85    mbs.SolveDynamic(Settings(solverType, solutionFile, stopTime, continueIfAvailable), solverType=solverType)
 86    return (mbs, mbs.systemData.GetODE2Coordinates())
 87
 88def Path(fileName):
 89    outputDirectory = exu.config.outputDirectory if hasattr(exu.config, 'outputDirectory') else ''
 90    return os.path.join(outputDirectory, fileName) if outputDirectory else fileName
 91
 92def Rows(fileName):
 93    with open(Path(fileName)) as f:
 94        return [line for line in f if line[0] != '#']
 95
 96def RemoveRestartFiles():
 97    for ending in ['', '.bck', '.tmp']:
 98        if os.path.isfile(Path(directory + 'restart.txt' + ending)):
 99            os.remove(Path(directory + 'restart.txt' + ending))
100
101u = 0
102for solverType in [exu.DynamicSolverType.GeneralizedAlpha, exu.DynamicSolverType.RK44]:
103    RemoveRestartFiles()
104    (mbsFull, qFull) = Run(solverType, 'full.txt')
105    rowsFull = [Rows(directory + name) for name in ['full.txt', 'sensorTip.txt', 'sensorBall.txt']]
106
107    RemoveRestartFiles()
108    (mbsStopped, qStopped) = Run(solverType, 'continued.txt', stopTime=0.6)
109    Check(0.6 - 1e-10 <= mbsStopped.systemData.GetTime() < 0.62, 'stopped at ' + str(mbsStopped.systemData.GetTime())) #the step size adapts at contact
110    #as if killed after the restart state of t = 0.4: the files hold rows up to 0.6 that are written again
111    shutil.copyfile(Path(directory + 'restart.txt.bck'), Path(directory + 'restart.txt'))
112    (mbsContinued, qContinued) = Run(solverType, 'continued.txt', continueIfAvailable=True)
113    restartTime = mbsContinued.sys['dynamicSolver'].output.restartTime
114    Check(0.4 - 1e-10 <= restartTime < 0.42, 'restart time ' + str(restartTime))
115    rowsContinued = [Rows(directory + name) for name in ['continued.txt', 'sensorTip.txt', 'sensorBall.txt']]
116
117    Check(np.linalg.norm(qContinued - qFull) == 0, str(solverType) + ': final state differs by ' + str(np.linalg.norm(qContinued - qFull)))
118    for (full, continued, name) in zip(rowsFull, rowsContinued, ['solution', 'sensor tip', 'sensor ball']):
119        Check(full == continued, str(solverType) + ': ' + name + ' file: ' + str(len(full)) + ' rows against ' + str(len(continued)))
120    u += np.sum(qFull)
121    exu.Print(solverType, 'final state', np.sum(qFull), ', rows', [len(r) for r in rowsFull])
122
123#a restart file of another system raises; continuing at the end time computes nothing
124try:
125    Run(exu.DynamicSolverType.RK44, 'other.txt', continueIfAvailable=True, addNode=True)
126    Check(False, 'a restart file of another system is taken')
127except Exception as error:
128    Check('does not fit' in str(error), str(error))
129(mbs, q) = Run(exu.DynamicSolverType.RK44, 'atEnd.txt', continueIfAvailable=True)
130Check(np.linalg.norm(q - qFull) == 0 and abs(mbs.sys['dynamicSolver'].output.restartTime - 1) < 1e-12, 'continued at the end time')
131
132#without a restart file, the simulation starts at its start time
133RemoveRestartFiles()
134(mbs, q) = Run(exu.DynamicSolverType.RK44, 'fresh.txt', continueIfAvailable=True)
135Check(mbs.sys['dynamicSolver'].output.restartTime == -1 and np.linalg.norm(q - qFull) == 0, 'without a restart file')
136
137exu.Print('restartFileTest: errors', errors)
138#what the test checks is the equality within one platform (errors); the final state itself differs between platforms
139#in the 10th digit through the contact (Linux, #2874), so it enters rounded
140u = round(u, 6) + errors
141exu.Print('solution of restartFileTest=', u)
142
143exu.sys['testResult'] = u