mininum timesteps after change to calculate the convergene accuracy

Hello Luc,

A full minimal working example (MWE) is usually the best way to get help on any programming question forum. If you are unfamiliar, an MWE is a small barebones script that we can copy and paste into our machines, run, and experience exactly what you are experiencing. MWEs give us a much better understanding of what you are trying to achieve, what you are doing wrong, and how you can change your script to obtain your desired results.

I am having trouble understanding your question, please help me clarify:

>>
In my simulation it is not the case as I guess the number of timestep before doing the calculation of the new unbalancedforce minium is too small :
>>

Why do you say the "number of timestep before doing the calculation of the new unbalanced force" is too small? Do is your simulation behaving in an unexpected manner? Is Yade crashing? Throwing an error? Please elaborate.

And maybe even more difficult for me to understand:

>>
my question is how to define or calculate or adjust this nimimal value depending of the parameter of my expriments: number of grains, friciton coefficient, and so on...
>>

Are you asking which parameter you should change to reduce the required number of timesteps to force equilibrium? Tough to answer such a question without knowing what material, physics, and law you are using.

Cheers,

rc

Hello Robert
thanks for trying to help me

in more detail, I am filling abox by classical gravity depositon than slowly titlt the box and look at the sphere displacements.

as the first 'step' made a correct converging force balance, the change of the gravity direction is not playing immediately enough change in the balancedforcecalculation to reach the 'new equilibrium', so I have used the same idea as in your oedometrixc sample case
i hope that I am enough clear?

 So here one lineof the table for my batch run then the full code py.
table:
pas_box theta_max converg_min nb_layer init_seed
.050 35. .01 20 10

py code:
# gravity deposition, continuing with oedometric test after stabilization
# shows also how to run parametric studies with yade-batch

# The components of the batch are:
# 1. table with parameters, one set of parameters per line (ccc.table)
# 2. readParamsFromTable which reads respective line from the parameter file
# 3. the simulation muse be run using yade-batch, not yade
#
# $ yade-batch --job-threads=1 03-oedometric-test.table 03-oedometric-test.py
#

# load parameters from file if run in batch
# default values are used if not run from batch
# gravity deposition in box, showing how to plot and save history of data,
# and how to control the simulation while it is running by calling
# python functions from within the simulation loop
#pas_box theta_max converg_min
#.05 35. .01
readParamsFromTable(pas_box=0.05, theta_max=35.0, converg_min=0.01,nb_layer=20.0,init_seed=10)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *

# import yade modules that we will use below
from yade import pack, plot, export

box_size = 5.0
# create rectangular box from facets
O.bodies.append(geom.facetBox((box_size/2.0,box_size/2.0,box_size*0.62),(box_size/2.0,box_size/2.0,box_size*0.62),wallMask=31))

# create empty sphere packing
# sphere packing is not equivalent to particles in simulation, it contains only the pure geometry
packing_fraction = 60./100.
global i_pas, pas_box,step0,gravity_y,gravity_z,theta_max,converg_min,nb_layer,init_seed
i_pas = 0
pas_box=.0250
step0 = 0
nb_layer=20.0
gravity_y = 9.81*sin( pas_box*i_pas*3.14/180.0)
gravity_z = 9.81*cos( pas_box*i_pas*3.14/180.0)

rayon = 0.025
height = 5. * rayon*nb_layer
nombre= int(packing_fraction*(box_size*box_size*nb_layer)/(rayon*rayon*3.14159) )
print(nombre)
print(i_pas*pas_box)
sp=pack.SpherePack()
# generate randomly spheres with uniform radius distribution
sp.makeCloud((0,0,0),(box_size,box_size,height),rMean=rayon,rRelFuzz=.005,num= nombre,seed=init_seed)
# add the sphere pack to the simulation
sp.toSimulation()

# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
 global i_pas, pas_box, step0,gravity_y,gravity_z,theta_max,converg_min
   # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
 if O.iter<(step0+5000): return
 if utils.unbalancedForce()<converg_min:
        i_pas = i_pas+ 1
        gravity_y = 9.81*sin( pas_box*i_pas*3.14/180.0)
        gravity_z = 9.81*cos( pas_box*i_pas*3.14/180.0)
        Newton_integrator.gravity=Vector3(0,-gravity_y,-gravity_z)
        step0 = O.iter
# print('%8.3f %8.3f %8.3f %8.3f' step0,pas_box*i_pas,sin( pas_box*i_pas*3.14/180.0),cos( pas_box*i_pas*3.14/180.0))
        print (step0, (pas_box*i_pas), gravity_y, gravity_z)
        if (pas_box*i_pas)>theta_max :
            O.pause()
# plot.saveDataTxt('bbb.txt.bz2')
        nom_file=str(i_pas)+"_"+str(pas_box)+"_"+str(nb_layer)+"_"+str(converg_min)
        filename1='./output_'+nom_file+'.txt'
        # sauvegarde fichier txt
        export.textExt(filename1,format='x_y_z_r',comment='x y z r')
        filename2='./output_'+nom_file+'.vtk'
        export.text2vtk(filename1,filename2)

# simulation loop
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
  # handle sphere+sphere and facet+sphere collisions
     [Ig2_Sphere_Sphere_ScGeom(),
      Ig2_Facet_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys(frictAngle=radians(15.0))],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),

 NewtonIntegrator(damping=0.4,gravity=(0,-gravity_y,-gravity_z),label='Newton_integrator'),
 # call the checkUnbalanced function (defined below) every 2 seconds
 PyRunner(command='checkUnbalanced()',realPeriod=5),
 # call the addPlotData function every 200 steps
 PyRunner(command='addPlotData()',iterPeriod=200)
]

#without friction
O.materials[0].frictionAngle=radians(15.0)
# timesteps
O.dt=.5*PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy=True

# collect history of data which will be plotted
def addPlotData():
 global i_pas, pas_box
# print( pas_box*i_pas,sin( pas_box*i_pas*3.14/180.0),cos( pas_box*i_pas*3.14/180.0)),
 # each item is given a names, by which it can be the unsed in plot.plots
 # the **O.energy converts dictionary-like O.energy to plot.addData arguments
 plot.addData(i=O.iter,unbalanced=utils.unbalancedForce(),**O.energy)

# define how to plot data: 'i' (step number) on the x-axis, unbalanced force
# on the left y-axis, all energies on the right y-axis
# (O.energy.keys is function which will be called to get all defined energies)
# None separates left and right y-axis
plot.plots={'i':('unbalanced',None,O.energy.keys)}

# show the plot on the screen, and update while the simulation runs
plot.plot()

# save tmp pour rerun convergence
O.saveTmp()

O.run()
# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()

I am traveling, but I will look at this when I arrive home!

Hi,

I did not completely get you :-)

This being said, assuming your concern is about the line https://github.com/yade/trunk/blob/master/doc/sphinx/tutorial/03-oedometric-test.py#L51 (and the comment L50) that is reflected in https://yade-dem.org/doc/tutorial-examples.html#oedometric-test:

This line and the choice and of a 5000 step number is only motivated by the fact unbalancedForce is not a correct equilibrium (in the sense resultant force on particles ~ 0) indicator at the start of a simulation where no / low contact forces exist.

Then some iterations (5000 here) are performed regardless of unbalancedForce() to allow some non-zero contact force / non-zero stress state appear in the numerical sample, making unbalancedForce() useable from this point onwards.

I'd advice you try to skip this "if O.iter<5000: return " kind-of exception in your general workflow (once your simulation exhibits some stress state).

Note running a (unknown beforehand) number of DEM steps until equilibrium can also be done with
while unbalancedForce() > yourThresholdValue:
    O.run(500,True) # you can use plenty of other values instead of 500

Dear Jerome,

When I tilted the gravity field, the full simulation needs 'some timesteps" before giving a new correct unbalancedforce calculation (otherwise the threshold test is immediatly still valid); checking every 500 steps will produce the same problem of initial recalculation of the new equilibrium of all the forces.
so this is the reason of the waiting of 5000 timesteps BEFOREchecking the unbalance, but it looks like that it is too small when I check the quality of the results!

so your comment about "you can use plenty of other values instead of 500" is what exactly I was looking for : how to approach the correct order of magnitude without doing a lot of test runs depending of the number of spheres, the friction coefificent, the stiffness and so on....
doing that at random is right now a waste of time.

Hi,
The following assertion is wrong, and consequently I don't understand the question:
"When I tilted the gravity field, the full simulation needs 'some timesteps" before giving a new correct unbalancedforce calculation"

If gravity is changed the unbalanced force is different and correct instantaneously, even after 0 iteration.
Proof below.
Cheers
Bruno

Yade [1]: tt=TriaxialTest()
Yade [2]: tt.generate('test.yade')
Yade [3]: O.load("test.yade")
Yade [4]: O.run(10000,1)
Yade [5]: unbalancedForce()
 -> [5]: 0.036382094190800705
Yade [6]: typedEngine('NewtonIntegrator').gravity=(100,0,0)
Yade [7]: unbalancedForce()
 -> [7]: 0.18847344871180757

this is not exactly my case:
for me, the spheres are only subjected to a gravity field : no external compression at all and I am awaiting for a "quasistatic equilibrium" before changing and the imposed change of the gravity is of the order of 10-4
which meant that the calcul of the 'new unbalancedForce is only changed immediatly by 10-6 order of magnitude before changing more after some time....
i have tried gravity (0.01,0 ,0) then (0.02,0,0) and looked at the immediat difference and after few timesteps

Hi,
In that case the question is very theoretical and I don't think a
theoretical answer is known.
I would think the amplitude of the gravity change in itself is a key
parameter.
But in fact the question is not so /important/interesting since it should
be enough to increase by very very small increments over time and check the
evolution of [unbalancedForce|kineticEnergy|momentum|whatever] as a
function of gravity.
In other words: instead of prescribing a change of magnitude, prescribe a
rate of change of magnitude.
Bruno