# GlobalStiffnessTimestepper gives dt of 0.0 when the saved file is restored

Asked by JINSUN LEE

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?

====

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

====

The MWE contains local path

==== 1st script, CBT_T1.py

from yade import pack, plot
import sys

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

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

==== 2nd script, CBT_T2.py

from yade import pack, plot
import sys

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

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='checkUnbalanced()', iterPeriod=1000),
# DomainLimiter(lo=(-width,-width,-height), hi=(width,width,5*height), iterPeriod = 10000, label = 'Domain') # destroy balls outside domain in every 100 steps
]

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")

# plotting

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

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='checkUnbalanced()', iterPeriod=1000),
# DomainLimiter(lo=(-width,-width,-height), hi=(width,width,5*height), iterPeriod = 10000, label = 'Domain') # destroy balls outside domain in every 100 steps
]

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
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")
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()

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
JINSUN LEE
Solved:
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2023-03-07: #1

Hello,
Do you need full O.save/O.load of full simulation?
Wouldn't be just spheres position export/import be sufficient for continuation?
Cheers
Jan

 Revision history for this message JINSUN LEE (ppk21) said on 2023-03-07: #2

Actually, the final goal of this code is to perform cyclic simple shear test.

Thus, the K0 stress condition overburdened by using top-cap plate should be transferred to the following stages (CBT_T4.py)
The minimum informations for the final test are position of bodies(spheres, facets), material properties, interactions (internal stresses and stresses acting between spheres and facets)

Regards.

 Revision history for this message Bruno Chareyre (bruno-chareyre) said on 2023-03-07: #3

Hello,
Could you please try an concatenate the files as much as possible (ideally just one file)? or find a way to reproduce the zero timestep which wouldn't require us to go through 3 scripts?
Bruno

 Revision history for this message Bruno Chareyre (bruno-chareyre) said on 2023-03-07: #4

Something rings a bell:

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

Does body "-1" have a mass, at least?...
I'm unsure save/load needs to be in the picture. That's why we ask to not send multiple scripts, it raises false issues in many cases.
Just load+run, nothing else, and see if the 0 timestep happens?

Do you really need to keep "z" unblocked? Blocking everything would probably fix the problem.

Bruno

 Revision history for this message JINSUN LEE (ppk21) said on 2023-03-08: #5

Thank you for your interest,

The first reason why I use O.save/O.load is to plan parametric study which is designed to find the effect of confining pressure, stress path and loading cycles(frequency, amplitude etc.) for the identical sphere sample.
Also, It is time consuming process to generate sphere sample by using pluviation process in DEM code when the number of particle exceed 10k.

The second reason is to see "what is going on?" during the simulation. If I use O.wait() then the graphical interface becomes to black

The last reason is trivial, because I got used to handle PFC code :)

===

DOF of the body "-1" is only allowed to move in "Z" direction. Because, the script is intended to apply overburden pressure using the mass of the top plate (facets). If DOF of the other directions "X" and "Y" is unblocked, the top plate will rotate or translate when it touches the spheres.

I was tried to use "ServoPIDController" but, it cannot keep pressure acting on "Z" direction when the top plate is subjected to translational sinusoidal motions in "Y" direction which is required to simulate cyclic simple shear test.

Regards.

 Revision history for this message Bruno Chareyre (bruno-chareyre) said on 2023-03-09: #6

I was not questioning why you use multiple scripts to produce results, this is your business (although none of the reasons you list are applicable see [1,2], there are better solutions). The point is that it is better to avoid multiple scripts at least when asking questions, namely because nobody wants to go through 3 scripts to try and understand your problem.
And, in fact, by trying to provide a unique script for this question you would have find that reloading wasn't the cause of your problem, at all.
The only problem is to unblock the DOFs of bodies with null mass (facets).
Bruno

[1] According to https://www.yade-dem.org/wiki/Howtoask: "If the input of one script is the output of another script, then both scripts should be merged into one. That is: avoid sending multiple scripts just like you avoid sending input files to a script. If needed, the method to save and reload in one single script is illustrated here [2]."