GSTS gives an unstable dt?

Hello everyone,

I'm Alessandro and I'm a new Yade user. I have to simulate the behaviour of a given set of spheres (radii, material, accelerations are given) inside a cylinder. I use yadedaily, xenial version.

The given set of accelerations translates into amplitudes of the HarmonicMotionEngine up to 7 mm. And here I find my problem: up to 1-2 mm all is ok, no unstability problems. But when I try, for example, with an amplitude of 3.5 mm, the spheres begin to overpass the facets of the cylinder and, as far as I know, it could be an issue linked to instabilities.
Firstly I tried to use a smaller dt, then I read about GSTS and decided to move to that. It wolud be perfect for my purpose: in the first part of the script I need a sort of "gravity deposition" and so a bigger dt would be nice, then the GSTS should assure the right dt in order to avoid instabilities (if I understood correctly).

But using GSTS, even with a very small safety coefficient, the spheres keep overpass the facets. I have to admit that with a lower coefficient the simulation seems to be improved, but the problem is not solved.

So: do you think that's an instability issue? Is the decreasing of the safety coefficient the only way to improve it? And shouldn't GSTS give a pretty good dt even without this very small coefficient?

I've used both CoundallStack (now commented) and ViscElPhys, and my problem is present in both cases.

Here I attach my script:

### MOVIMENTO ARMONICO ALTO BASSO

from yade import pack,utils,plot,ymport,export,qt,timing

## Physical Parameters
density=7950
en = 0.5
es = 0.5
cyl_r = 12.7
cyl_h = 24.6
sp_p_r = 1.5
fill_size = 11.5
young = 200e9
poisson = 0.305

## Creation Cylinder
facetMat=O.materials.append(ViscElMat(density=density,young=young, poisson=poisson,en=en,et=es))
Cylinder=O.bodies.append(geom.facetCylinder((0,0,cyl_h/2),cyl_r,height=cyl_h,segmentsNumber=15,wallMask=7,material=facetMat))

## Creation Spheres
sphereMat=O.materials.append(ViscElMat(density=density,young=young,poisson=poisson,en=en,et=es))
sp=pack.SpherePack()

sp.makeCloud((-(cyl_r-2),-(cyl_r-2),1.5),(cyl_r-2,cyl_r-2,fill_size-1.5),rMean=sp_p_r,num=50,distributeMass=False)
sp.toSimulation(material=sphereMat)

sp2=pack.SpherePack()
sp2.makeCloud((-(cyl_r-2),-(cyl_r-2),fill_size),(cyl_r-2,cyl_r-2,cyl_h),rMean=3,num=3,distributeMass=False)
sp2.toSimulation(material=sphereMat)

## Engine
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
      #[Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_ViscElPhys_Basic()]
      #[Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   GlobalStiffnessTimeStepper(timeStepUpdateInterval=1,label='tss',timestepSafetyCoefficient=0.001,viscEl=True),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.3),
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker')
   ]

def checkUnbalanced():
 if O.iter<15000: return
        if kineticEnergy()>18000000: return
        O.engines=O.engines+[HarmonicMotionEngine(A=[0,0,3.5],f=[0,0,60.0],fi=[0,0,pi],ids=Cylinder), PyRunner(command='addPlotData()',iterPeriod=10)]
 checker.command='nothing()'

def nothing():
        return

def addPlotData():
 Ek,maxId=kineticEnergy(findMaxId=True)
 plot.addData(i=O.iter,t=O.time,Ninter=O.interactions.countReal(),total=O.energy.total(),maxId=maxId,**O.energy)

O.trackEnergy=True
O.saveTmp()

plot.plots={'i':('Ninter',),'t':(O.energy.keys,None,'total')}

plot.plot()

qt.Controller()
qt.View()