GlobalStiffnessTimestepper gives dt of 0.0 when the saved file is restored

As noted in the previous question (https://answers.launchpad.net/yade/+question/261239),
I have encountered similar situations.

The program code consists of three Yade sccript.
The first one (CBT_T1.py) makes boxes and walls using facets and generate save file "CDSS_box_1.yade"
The second one (CBT_T2.py) generate spheres using free falling procedures
at this stage, the critical time step is calculated using "GlobalStiffnessTimestepper" and generate save file "CDSS_ball_2.yade"

The problem occurs when running the third stage script (CBT_T3.py)
At this stage, the saved file (CDSS_ball_2.yade) is restored and the critical time step is calculated using "GlobalStiffnessTimestepper"

The simulation goes well until it reaches the iteration of 1000 which re-calculate critical time step using c

The "GlobalStiffnessTimestepper" gives critical time step as 0.0.

What is the problem?

Thanks in advance.

====

Working environment :
Ubuntu 22.10 kinetic
Yade 2022.01a
Intel Xeon Gold 6252 2EA

====

The MWE contains local path

==== 1st script, CBT_T1.py

from yade import pack, plot
import sys
sys.path.append('/home/jinsun/Dropbox/Yadee/CDSS_BOX_TEST2')

def setGeomVars (): # initialize variables
     width = 0.05
     height = 0.03
     margin = 30
     # Calculate extra length
     dx = width/100/2*margin
     dy = width/100/2*margin
     dz = height/100/2*margin
     saveVars('geoms',loadnow=True,**locals())

setGeomVars ()
from yade.params.geoms import * # load initilized variables

# side pannel
p1s = (-width/2,-width/2-dy,-height/2-dz)
p5s = (-width/2,-width/2-dy,height/2+dz)
p6s = (-width/2,width/2+dy,height/2+dz)
p2s = (-width/2,width/2+dy,-height/2-dz)

side1_1 = utils.facet(vertices=[p1s,p5s,p2s], wire=True, highlight=False) #p1 p5 p2
side1_2 = utils.facet(vertices=[p2s,p5s,p6s], wire=True, highlight=False) #p1 p5 p2
O.bodies.append(side1_1)
O.bodies.append(side1_2)

p4s = (width/2,-width/2-dy,-height/2-dz)
p8s = (width/2,-width/2-dy,height/2+dz)
p7s = (width/2,width/2+dy,height/2+dz)
p3s = (width/2,width/2+dy,-height/2-dz)

side2_1 = utils.facet(vertices=[p4s,p8s,p3s], wire=True, highlight=False) #p1 p5 p2
side2_2 = utils.facet(vertices=[p3s,p8s,p7s], wire=True, highlight=False) #p1 p5 p2
O.bodies.append(side2_1)
O.bodies.append(side2_2)

# front pannel
p1f = (-width/2-dx,-width/2,-height/2-dz)
p5f = (-width/2-dx,-width/2,height/2+dz)
p4f = (width/2+dx,-width/2,-height/2-dz)
p8f = (width/2+dx,-width/2,height/2+dz)

front1 = utils.facet(vertices=[p1f,p5f,p8f], wire=True, highlight=False)
front2 = utils.facet(vertices=[p1f,p8f,p4f], wire=True, highlight=False)
O.bodies.append(front1)
O.bodies.append(front2)

# back pannel
p2b = (-width/2-dx,width/2,-height/2-dz)
p6b = (-width/2-dx,width/2,height/2+dz)
p3b = (width/2+dx,width/2,-height/2-dz)
p7b = (width/2+dx,width/2,height/2+dz)

back1 = utils.facet(vertices=[p2b,p6b,p7b], wire=True, highlight=False)
back2 = utils.facet(vertices=[p2b,p7b,p3b], wire=True, highlight=False)
O.bodies.append(back1)
O.bodies.append(back2)

#bottom pannel
p1bt = (-width/2-dx,-width/2-dy,-height/2)
p2bt = (-width/2-dx,width/2+dy,-height/2)
p3bt = (width/2+dx,width/2+dy,-height/2)
p4bt = (width/2+dx,-width/2-dy,-height/2)

bot1 = utils.facet(vertices=[p1bt,p2bt,p3bt], wire=True, highlight=False)
bot2 = utils.facet(vertices=[p1bt,p4bt,p3bt], wire=True, highlight=False)
O.bodies.append(bot1)
O.bodies.append(bot2)

collar = geom.facetBox(center=(0,0,(2*height)), extents=(width/2,width/2,3*height/2),wallMask=15)
O.bodies.append(collar)

#save
O.save("CDSS_box_1.yade")

==== 2nd script, CBT_T2.py

from yade import pack, plot
import sys
sys.path.append('/home/jinsun/Dropbox/Yadee/CDSS_BOX_TEST2')
O.load("CDSS_box_1.yade")
loadVars('geoms')

from yade.params.geoms import * # load initilized variables

O.materials.append(FrictMat(young=20e6, poisson=0.17, density=2700, frictionAngle=0.523, label='frictmat'))
#O.materials.append(FrictMat(young=70e9, poisson=0.17, density=2300, frictionAngle=0.523, label='frictmat')) # silica

cp1 = (-width/2,-width/2,3*height) #corner point #1
cp2 = (width/2,width/2,2*height+3*height/2) #corner point #2
radius_mean = 0.002

sp = pack.SpherePack()
sp.makeCloud(cp1, cp2, rMean=radius_mean, rRelFuzz=0.0, num = 10000)
sp.toSimulation(color=(0,0,1)) # pure blue

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
    InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()], # collision geometry
            [Ip2_FrictMat_FrictMat_FrictPhys()], # collision "physics"
            [Law2_ScGeom_FrictPhys_CundallStrack()] # contact law -- apply forces
    ),
    # Apply gravity force to particles. damping: numerical dissipation of energy.
    GlobalStiffnessTimeStepper(active=True, timestepSafetyCoefficient=0.8, timeStepUpdateInterval=1000, label = 'timeStepper'),
    NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.5),
    PyRunner(command='addData()', iterPeriod=1000),
    PyRunner(command='checkUnbalanced()', iterPeriod=1000),
    PyRunner(command='addPlotData()', iterPeriod=1000),
# DomainLimiter(lo=(-width,-width,-height), hi=(width,width,5*height), iterPeriod = 10000, label = 'Domain') # destroy balls outside domain in every 100 steps
]

def addData():
     for b in O.bodies:
         b.shape.color=scalarOnColorScale(b.state.vel.norm(),0.0,1e-1)

elapsed_time=0.0
def checkUnbalanced():
# print(unbalancedForce())
     global elapsed_time
     if (O.time-elapsed_time) > 1.0:
        ball_highest_z = numpy.max([b.state.pos[2] for b in O.bodies if isinstance(b.shape,Sphere)]) # make list of ball z position if body is sphere and find max value
        print("The highest ball position is = ", ball_highest_z)
        if ball_highest_z < height/2:
             sp.makeCloud(cp1, cp2, rMean=radius_mean, rRelFuzz=0.0, num = 10000)
             sp.toSimulation(color=(0,0,1)) # pure blue
             print("Total number of balls =", len(O.bodies))
             elapsed_time=O.time
        else:
             O.pause(); print("simulation paused")
             O.save("CDSS_ball_2.yade")

# plotting
def addPlotData():
    plot.addData(i=O.iter, unbalanced=unbalancedForce())

plot.plots={'i':('unbalanced')}
plot.plot()

O.dt=0.5*PWaveTimeStep()
O.run() #; O.wait()

==== 3rd script, CBT_T3.py

from yade import pack, plot
import sys
sys.path.append('/home/jinsun/Dropbox/Yadee/CDSS_BOX_TEST2')
O.load("CDSS_ball_2.yade")
loadVars('geoms')

from yade.params.geoms import * # load initilized variables

for i in range(10, 18):
    O.bodies.erase(i)

# create topcap
p5t = (-width/2-dx,-width/2-dy,height/2+2*dz)
p6t = (-width/2-dx,width/2+dy,height/2+2*dz)
p7t = (width/2+dx,width/2+dy,height/2+2*dz)
p8t = (width/2+dx,-width/2-dy,height/2+2*dz)

topcap1 = utils.facet(vertices=[p5t,p6t,p7t], wire=True, highlight=False)
topcap2 = utils.facet(vertices=[p5t,p8t,p7t], wire=True, highlight=False)
O.bodies.append(topcap1)
O.bodies.append(topcap2)

topcap_area = width**2
target_press = 100000 # unit in Pa
mass = target_press*topcap_area/2/9.81
target_force = target_press * topcap_area
O.bodies[-1].state.mass=mass
O.bodies[-2].state.mass=mass
print('Each facet mass is', mass)
print('Each facet target force is', mass)

#O.engines=O.engines[:6]+O.engines[7:]
O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
    InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()], # collision geometry
            [Ip2_FrictMat_FrictMat_FrictPhys()], # collision "physics"
            [Law2_ScGeom_FrictPhys_CundallStrack()] # contact law -- apply forces
    ),
    # Apply gravity force to particles. damping: numerical dissipation of energy.
    GlobalStiffnessTimeStepper(active=True, timestepSafetyCoefficient=0.8, timeStepUpdateInterval=1000, label = 'timeStepper'),
    NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.5),
    PyRunner(command='addData()', iterPeriod=1000),
# PyRunner(command='checkUnbalanced()', iterPeriod=1000),
    PyRunner(command='addPlotData()', iterPeriod=1000),
# DomainLimiter(lo=(-width,-width,-height), hi=(width,width,5*height), iterPeriod = 10000, label = 'Domain') # destroy balls outside domain in every 100 steps
]

def addData():
     for b in O.bodies:
         b.shape.color=scalarOnColorScale(b.state.vel.norm(),0.0,1e-1)

def cal_stress(listup):
    for i in listup:
        i+=1
    return i

# plotting
i_times=0
def addPlotData():
    topcap_zforce=(O.forces.f(O.bodies[-1].id)[2] + O.forces.f(O.bodies[-2].id)[2])
    topcap_zdisp=O.bodies[-1].state.displ()[2] # direction in z
    topcap_zvel=O.bodies[-1].state.vel[2] # direction in z
    #Calculate shear stress
    global Box_Vol
    Box_Vol = (height + dz + O.bodies[-1].state.displ()[2])*width**2 # correct topcap z position
    Force_Matrix=Matrix3.Zero
    Cauchy_stress=Matrix3.Zero

    for i in O.interactions:
        if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
            f = i.phys.normalForce + i.phys.shearForce
            l = O.bodies[i.id2].state.pos - O.bodies[i.id1].state.pos
            Force_Matrix+=f.outer(l)
    Cauchy_stress=Force_Matrix/Box_Vol
    #Plot data
    plot.addData(i=O.iter, topcap_zdisp=topcap_zdisp, topcap_zvel= topcap_zvel, topcap_zstress= topcap_zforce/topcap_area,\
         Sxx=Cauchy_stress[0][0], Sxy=Cauchy_stress[0][1], Sxz=Cauchy_stress[0][2],\
            Syx=Cauchy_stress[1][0], Syy=Cauchy_stress[1][1], Syz=Cauchy_stress[1][2],\
                Szx=Cauchy_stress[2][0], Szy=Cauchy_stress[2][1], Szz=Cauchy_stress[2][2], box_vol=Box_Vol)
    global i_times
    servo_error=abs(topcap_zforce - target_force)/target_force*100.0
    if servo_error < 0.1:
        i_times += 1
        print("Number of servo condition satisfied =", i_times, servo_error)
        if i_times > 100:
            O.pause(); print("Target servo force is reached")
            O.save("CDSS_ovp_3.yade")
    else:
        i_times = 0 # zero setting

O.bodies[-1].shape.highlight = True
O.bodies[-2].shape.highlight = True
O.bodies[-1].shape.wire = False
O.bodies[-2].shape.wire = False

O.bodies[-1].state.blockedDOFs='xyXYZ'
O.bodies[-2].state.blockedDOFs='xyXYZ'

plot.plots={'i':('topcap_zdisp',), 'i ':('box_vol'), 'i ':('topcap_zstress'), 'i ':('Szz')}
plot.plot()

O.dt=0.5*PWaveTimeStep()
O.run()