Yade

whats the difference between “GlobalStiffnessTimeStepper”and "O.dt"

Asked by mikexue

Hi，

My simulation is about a rotating rod inserted into a chamber that contains about 470,000 particles.
When I use "GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8)", I get a dynamic dt around 3.64e-6.
The program proceeded with a speed of 3 iter/s. It takes a long time to run the whole simulation.
For some reasons I want to set a fixed dt to calculate the position.
I wonder if I can set the dt equals to maybe 4e-6 or 1e-5 or even a larger dt to save some time?
Would it affect the stability or the data such as the contact force between the particles?

I appreciate your help,
Mikexue

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said : 		#1

Hello,

> whats the difference between “GlobalStiffnessTimeStepper”and "O.dt"

GlobalStiffnessTimeStepper sets the time step (O.dt actually) "automatically" based on simulation conditions.
It can be very handy e.g. if you insert particles with different material etc., you just add them and do not take care of O.dt, which is set to a safe value.

Apart from letting GlobalStiffnessTimeStepper to set O.dt automatically, you can set it manually.

You can combine both approaches (e.g. part of simulation this, part that).

> It takes a long time to run the whole simulation.

yes, large-scale simulations might take long time..

> I wonder if I can set the dt equals to maybe 4e-6 or 1e-5 or even a larger dt to save some time?

Yes, you can.
But probably it is not a good idea. You would save time in exchange for (possibly, likely) non-sense results.
Higher value than determined by the GlobalStiffnessTimeStepper could (most likely would) make the simulation unstable.

> Would it affect the stability or the data such as the contact force between the particles?

Most likely yes.

To make the simulation faster, you can try to run it in parallel (multi-core).

Critical time step is determined by mass/density, stiffness and size of particles [1].
Based on your problem*, you can do some mass/stiffness/size tricks, but each of them changes the physics (!!!). Some more, some less..

* please provide more information, e.g. what particle size distribution you use, what materials and laws, ...

Cheers
Jan

[1] https://yade-dem.org/doc/formulation.html#stability-considerations

Revision history for this message
mikexue (mikexue) said (last edit ): 		#2

 Hi, Jan.
Thank you for your answer. I find that "O.dt=SafetyCoefficient*utils.PWaveTimeStep() " is also commonly used in the Yade documentation.

Here is another question related to dt:
Besides the fact that GlobalStiffnessTimeStepper is dynamic and utils.PWaveTimeStep() is fixed, should they be of the same value at the very beginning of the simulation?
I mean dt="GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8) "equals to dt="0.8*utils.PWaveTimeStep()" at the beginning of the simulation?
I am confused coz I tried these two ways and they are different in my simulation.

best,
Mikexue

Revision history for this message
Best Jan Stránský (honzik) said : 		#3

Instant replay:

please provide more information, e.g. what particle size distribution you use, what materials and laws, ...
Basically a MWE [2].

utils.PWaveTimeStep() is simply min(radius/sqrt(young/density))
GlobalStiffnessTimeStepper [3] computes actual stiffness, includes viscosity if needed etc. (basically is much more sophisticated) so a different result is expected.
You can e.g. try to put two spheres and then four spheres at almost same position, the GlobalStiffnessTimeStepper is different, whereas the PWaveTimeStep is constant:
###
# GlobalStiffnessTimeStepper is included in the default engines
xs = [0, 1]
#xs = [0, 1, 0.1, 1.1] # un-comment this line for the second run
O.bodies.append([sphere((x,0,0),1) for x in xs])
O.step()
print(O.dt,PWaveTimeStep())
###

Cheers
Jan

[2] https://www.yade-dem.org/wiki/Howtoask
[3] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/GlobalStiffnessTimeStepper.cpp

Revision history for this message
mikexue (mikexue) said : 		#4

Thanks Jan Stránský, that solved my question.

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#5

Hi,
> I am confused coz I tried these two ways and they are different in my simulation

They are derived with slightly different logics and so they have no reason to be exactly the same. Usually they should be of the same order.

More generally I would give a response even stronger that Jan's: if GSTS gives a certain O.dt with timestepSafetyCoefficient=0.8, then the maximum O.dt you can theoretically use is 1.25 times greater (1/0.8). If you exceed that it will not just "affect the stability". In most cases it will go completely wrong (division by zero, NaN and such things) and at some point it will segfault.

Better experiment with 4000 particles first, it will be easier to understand what happens.

Also note that there is an example rotating drum using MPI parallelization:
https://gitlab.com/yade-dev/trunk/-/blob/master/examples/DEM2020Benchmark/Case2_rotating_drum_mpi.py

Cheers
Bruno