Why particle velocity varies with time step

Asked by Aoxi Zhang

Hi,
I was trying to simulate a bender element test using DEM, and I'd like to ask that why the evolution of particle velocity varies with time step, even dt is small enough compared to the dt obtained from PWaveTimeStep().

The attached MWE can illustrate my question, the MWE consists of 3 scripts, which can be found from [1].

(1) phase1.py is for sample generation and stress relaxation, it takes around 7 minutes if you want to reproduce the results.

The idea is: firstly having a sample according to Bruno's example[2], then closing triax engine and increasing damping to quickly reduce unb. phase1 stops when unb reaches a small value (e.g.,1e-7).

(2) phase2.py is the test under different dt, if you want to reproduce the resutls, you may use yade-batch phase2_batchParam.txt phase2.py, it takes around 8 minutes in general.

The idea of this test is: firstly load the sample generated from phase1, then turn damping to 0, close time stepper, select some particles whose velocity will be recorded. phase2 will terminate when the time reaches 3e-4.

(3) phase3.py is for plotting the results to get a figure which compared the results under different dt. You can directly find it in the attached files.

The MWE above leads to two questions:

(a) When I run phase2 in batch mode, not all simulations terminate when the time reaches 3e-4. Some of them terminated in advance. This can be seen from the figure.

(b) As you can see from the figure, different dt leads to different velocity of the same group of particles. If I understood correctly, when dt is small enough, different dt should not affect the results. Do you have any ideas about it?

Thank you very much!
Aoxi

[1]https://we.tl/t-B5s6rc6dRP
[2]https://gitlab.com/yade-dev/trunk/blob/master/examples/triax-tutorial/script-session1.py

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#1

Hi,

Could you just just paste your code? Please find enclosed some good reasons provided for such a practice [1].

Best wishes,
KB

[1] https://www.yade-dem.org/wiki/Howtoask

PS I know that it is a little bit ironic that I am providing link saying why you shouldn't incorporate links in your posts. But mostly refers to external file repositories :)

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

Hello,

please read [3], section "Please, no external links!".
MWE can be put directly in the question (assuming satisfying M=minimal condition).

> even dt is small enough compared to the dt obtained from PWaveTimeStep()
> If I understood correctly, when dt is small enough, different dt should not affect the results. Do you have any ideas about it?

Straightforward answer could be that it is not small enough to neglect the effect of its value

But, more importantly:
You are examining a few "random" particle's velocities in a simulation with "many" particles.
Such DEM simulation is "chaotic", i.e. a small change in input may lead to big changes for such quantities like velocity or position of one individual particle.
You can try running exactly the same simulation (same O.dt, same initial configuration, everything the same) in parallel, let's say with "-j 4" option and you should get different result for each run [4].
So what you are experiencing with different velocities is the feature of DEM and its computer solution.

So you should be interested in some "overall" quantities, like average velocity of some band or region of particles

Concerning premature exit, I have no idea.
The log files of those simulations does not contain the line "Running finished". But it has status: 0 (OK) result. So it must have been stopped somehow differently than by stopSimulation() function..

Cheers
Jan

[3] https://www.yade-dem.org/wiki/Howtoask
[4] https://yade-dem.org/doc/formulation.html#result-indeterminism

Revision history for this message
Aoxi Zhang (aoxizhang) said :
#3

Hi Karol and Jan,

Thanks for your feedback!

Below are the scripts of the MWE:

#### phase1.py####
from __future__ import division
from yade import pack, plot
import math
import numpy as np
import random
from random import gauss
import timeit
import pickle

num_spheres=5000
targetPorosity = 0.405
compFricDegree = 30
finalFricDegree = 19
rate=-0.01
damp=0.6
stabilityThreshold=0.001

confinement=100e3

mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)

MatWall=O.materials.append(FrictMat(young=2e9,poisson=0.3,frictionAngle=0,density=0,label='walls'))

MatSand = O.materials.append(CohFrictMat(isCohesive=True,young=2e8,alphaKr=0.15,alphaKtw=0.15,\
           poisson=0.3,frictionAngle=radians(30),etaRoll=0.15,etaTwist=0.15,\
           density=2650.0,normalCohesion=0, shearCohesion=0,\
           momentRotationLaw=True,label='sand'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()

sp.makeCloud(mn,mx,-1,0.3,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='sand') for center,rad in sp])

Gl1_Sphere.quality=3

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="Ip2Coh"),
   Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,label="TimeStepper"),
 TriaxialStressController(
  maxMultiplier=1.4,
  finalMaxMultiplier=1.004,
  thickness = 0,
  stressMask = 7,
  internalCompaction=True,
  label='triax'),
 newton,
]

Gl1_Sphere.stripes=0
triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1]
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001:
   break

print "### state 1 Isotropic completed ###"
triax.internalCompaction=False
import sys
while triax.porosity - targetPorosity>0.0001:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)

print "### state 2 Reach target porosity completed ###"
print "### state 3 - click run to start stress relaxation ###"
setContactFriction(radians(finalFricDegree))

NewtonIntegrator.damping=0.9
O.engines=O.engines+[PyRunner(command='stopStressRelaxation()', iterPeriod=1,label="StopRelaxation")]

for i in range(6):
 O.bodies[i].state.vel = Vector3(0,0,0)
 O.bodies[i].state.blockedDOFs='xyzXYZ'

triax.dead=True

def stopStressRelaxation():
 print('unb',unbalancedForce())
 if unbalancedForce()<1e-7:
  O.pause()
  NewtonIntegrator.damping=0.0#newton.damping=0
  print "stress relaxation finished, damping has been set as 0"
  print "Please go to phase2"
  O.save('sample_generatedFromPhase1_afterRelax.yade.gz')

#####
##### phase2.py#####
from __future__ import division
from yade import pack, plot
import math
import numpy as np
import random
import pickle
from yade import ymport,export

utils.readParamsFromTable(dt=1e-7)
from yade.params import table

finalFricDegree=19

O.load('sample_generatedFromPhase1_afterRelax.yade.gz')

RunTimeLength=3e-4 # how long for the time of data recording
Odt=table.dt
O.dt=Odt

O.engines[3].dead=True## turn of GlobalStiffnessTimeStepper
# TimeStepper.dead=True## Note, use label cannot successfully close time stepper
StopRelaxation.dead=True## turn of StopRelaxation in previous phase
NewtonIntegrator.damping=0.0

monitoredSand=[1000,2001,3002,4003,5004]

def avgVel(idList):
 vel=0.0
 avg=0.0
 for i in idList:
   vel+=O.bodies[i].state.vel[1]
 avg=vel/len(idList)
 return avg

print("O.dt=",O.dt)
O.engines=O.engines+[PyRunner(command='outputData()', iterPeriod=1,label="outputData")]
O.engines=O.engines+[PyRunner(command='stopSimulation()', iterPeriod=1,label="stopSimu")]

def stopSimulation():
 print ('Running, O.time,',O.time)
 if O.time > RunTimeLength:
  O.pause()
  print ('Running finished, O.realtime,',O.realtime)

def outputData():
 f = open("ResultsOfPhase2_dt{}.txt".format(Odt), "a+")
 f.write('{} {} {}\n'.format(O.time,O.dt,avgVel(monitoredSand)))
 f.close()

O.resetTime()
O.run()
utils.waitIfBatch()
######
###### phase3.py#####
from yade import pack, plot
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator
import matplotlib.ticker as ticker
import numpy as np
import pandas as pd
font1 = {'family': 'Arial', 'fontweight': 'normal', 'size': 14, 'style':'normal'}

linewidth=2
fig , ax = plt.subplots(figsize=(6.25,3))
bwith = 0.5# here control the width of axis
ax.spines['bottom'].set_linewidth(bwith)
# ax.spines['bottom'].set_color('black')
ax.spines['left'].set_linewidth(bwith)
# ax.spines['left'].set_color('black')
ax.spines['top'].set_linewidth(bwith)
# ax.spines['top'].set_color('black')
ax.spines['right'].set_linewidth(bwith)
# ax.spines['right'].set_color('black')

def plotResults(dt):
 f=pd.read_csv("ResultsOfPhase2_dt{}.txt".format(dt),sep=" ",header=None)
 f.columns=["Otime","Odt","avgVel"]
 Otime,Odt,avgVel= np.array(f["Otime"]),np.array(f["Odt"]),np.array(f["avgVel"])
 ax.plot(Otime,avgVel,label="dt={}".format(dt))

plotResults(1e-7)
plotResults(1e-8)
plotResults(1e-9)
plotResults(5e-9)
# =============================================================================
# Figure setting
#1 set legend
ax.legend(loc=0,frameon=False,prop={'family': 'Arial','size': 11})

ax.set_xlabel('Time (s)',font1)
ax.set_ylabel('Velocity (m/s)',font1)

ax.tick_params(which='both',labelsize=10,direction='in')
plt.xticks(fontproperties = 'Arial', size = 10)
plt.yticks(fontproperties = 'Arial', size = 10)
# plt.ylim(-3e-5,-2e-5)
plt.xlim(0,3e-4)

plt.tight_layout(pad=0.4)
plt.savefig('Fig_compare_dt.jpg', dpi=600)
plt.show()
#####

#### For run phase2.py in batch mode, another txt file could be as below:
dt
1e-7
5e-8
1e-8
5e-9
1e-9

Thanks
Aoxi

Revision history for this message
Aoxi Zhang (aoxizhang) said :
#4

Hi Jan,

Thanks for your comments.

>> it is not small enough to neglect the effect of its value.

I further tested on other smaller dt values, e.g. 1e-10, but it has the same problem.

>> You are examining a few "random" particle's velocities in a simulation with "many" particles. Such DEM simulation is "chaotic".

That is a good point. I updated phase2.py (attached below) by monitoring a layer of sands, thereby more sands are monitored. But I still got this problem.

#### updated phase2.py in which more sands are monitored###
from __future__ import division
from yade import pack, plot
import math
import numpy as np
import random
import pickle
from yade import ymport,export

utils.readParamsFromTable(dt=1e-7)
from yade.params import table

finalFricDegree=19

O.load('sample_generatedFromPhase1_afterRelax.yade.gz')

RunTimeLength=3e-4 # how long for the time of data recording
Odt=table.dt
O.dt=Odt

O.engines[3].dead=True## turn of GlobalStiffnessTimeStepper
# TimeStepper.dead=True## Note, use label cannot successfully close time stepper
StopRelaxation.dead=True## turn of StopRelaxation in previous phase
NewtonIntegrator.damping=0.0

# monitoredSand=[1000,2001,3002,4003,5004]
monitoredSand=[]
def getMonitoredSand():
 halfHeigh=aabbDim()[1]/2
 for i in O.bodies:
  if isinstance(i.shape, Sphere):
   if halfHeigh<i.state.pos[1]<halfHeigh+0.0015:
    i.shape.color=[1,0,0]
    monitoredSand.append(i.id)
   else:
    i.shape.color=[1,1,1]
 print ("number of monitoredSand is",len(monitoredSand))

getMonitoredSand()

def avgVel(idList):
 vel=0.0
 avg=0.0
 for i in idList:
   vel+=O.bodies[i].state.vel[1]
 avg=vel/len(idList)
 return avg

print("O.dt=",O.dt)
O.engines=O.engines+[PyRunner(command='outputData()', iterPeriod=1,label="outputData")]
O.engines=O.engines+[PyRunner(command='stopSimulation()', iterPeriod=1,label="stopSimu")]
# print("O.dt=",O.dt)
# NewtonIntegrator.damping=0

def stopSimulation():
 print ('Running, O.time,',O.time)
 if O.time > RunTimeLength:
  O.pause()
  print ('Running finished, O.realtime,',O.realtime)
  # plot.saveDataTxt('historyData_calmdownVel{}.txt'.format(calmDownVel)) # should revise the monitored data at history

def outputData():
 f = open("ResultsOfPhase2_dt{}.txt".format(Odt), "a+")
 f.write('{} {} {}\n'.format(O.time,O.dt,avgVel(monitoredSand)))
 f.close()

O.resetTime()
O.run()
utils.waitIfBatch()
#####
Do you have any ideas about this?

Thanks!
Aoxi

Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#5

Hi,

I didn't run your code yet, but I am looking at the example. Are you sure that you modify the damping by

> NewtonIntegrator.damping=0.9

earlier you created an instance called 'newton', so maybe it should be rather

>newton.damping=0.9

I prefer using labels for engines.

Cheers,
Karol

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

@Karol: I was surprised, too, but it does works also this way :-D
It seems that somehow using class directly passes the argument to some instances.
Maybe something about how boost.python and YADE macros deals with it?
Anyway, I would not call it the best practice..

### MWE
newton1 = NewtonIntegrator()
print("1",newton1.damping) # 0.2
newton1.damping = 0.3
print("1",newton1.damping) # 0.3
NewtonIntegrator.damping = 0.4
print("1",newton1.damping) # 0.4 (!)
newton2 = NewtonIntegrator
print("1",newton1.damping) # 0.4
print("2",newton2.damping) # 0.4
NewtonIntegrator.damping = 0.5
print("1",newton1.damping) # 0.5 (!)
print("2",newton2.damping) # 0.5 (!)
newton1.damping = 0.6
newton2.damping = 0.7
print("1",newton1.damping) # 0.6
print("2",newton2.damping) # 0.7
NewtonIntegrator.damping = 0.8
print("1",newton1.damping) # 0.6 (!)
print("2",newton2.damping) # 0.8 (!)
###

Cheers
Jan

Revision history for this message
Aoxi Zhang (aoxizhang) said :
#7

Hi,

I'd like also use labels for engines.
But sometimes I found that using label cannot successfully control the engine. For instance, in phase2.py:

O.engines[3].dead=True## turn of GlobalStiffnessTimeStepper
# TimeStepper.dead=True## Note, use label cannot successfully close time stepper

In above codes, I want to close GlobalStiffnessTimeStepper, if I use TimeStepper.dead=True, it doesn't work. Thereby, I use O.engines[3].dead=True to turn of GlobalStiffnessTimeStepper. I don't know why.

Regarding the original question, I use newton.damping=0.9 and/or NewtonIntegrator.damping=0.9 to set the damping. The problem is still there :(

Any ideas or comments are welcome.

Best,
Aoxi

Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#8

Just a comment regarding changing the damping. I loaded the save from phase1, and inspected damping in the following way

> O.engines[5].damping

Output is 0.6, which confirms that damping didn't change during the simulation at all. Also, please note that you set damping in function stopStressRelaxation. After setting damping there as
>newton.damping = 0
it changes as expected (please inspect the save file after this). Unfortunately, it doesn't answer your question.

Also, I think that saving data every iteration may cause some efficiency issues (in phase 2). Instead, I propose the same time interval for all the cases (virtPeriod = RunTimeLength/1000).

Cheers,
Karol

Revision history for this message
Aoxi Zhang (aoxizhang) said :
#9

Hi Karol,

I have (1) updated the way of setting daping in phase2.py, now it works well, (2) changed iterPeriod=1 to virtPeriod = RunTimeLength/1000 as you suggested.

Thanks,
Aoxi

Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#10

Hi Zhang,

Do you mean 'work well', so the damping is changing properly, but you still need help or work well and the problem is solved?

Cheers,
Karol

Revision history for this message
Aoxi Zhang (aoxizhang) said :
#11

Hi Karol,

Sorry for not making it clear.
I mean the damping is changing properly, but the problem is not solved and I still need help on it.

Thanks,
Aoxi

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

Hi,
I guess you want to simulate elastic waves or something, but as I was reading phase 2 I did not see anything special exerted on the boundaries.
If you analyze velocities in a situation of static equilibrium, what you get is just some random noise. It is not deeply surprizing that random noise depends on timestep (and, basically, on anything else).
Did I miss something?
Bruno

Can you help with this problem?

Provide an answer of your own, or ask Aoxi Zhang for more information if necessary.

To post a message you must log in.