segmentation fault after erasing some particles

[azim]

Dear all,
I built a simple script according to script-session1.py (B.Chareyre). I tried to erase some particles, it was done without any problem. but after finishing the simulation, as soon as I tried to run more iterations, it showed segmentation fault and quit yade.
The length of O.bodies is similar to The one before erasing.(Does not change after erasing)
you can see my question in:
https://answers.launchpad.net/yade/+question/668274
I'm using Ubuntu 14.04 LTS and yadedaily version: 2018.02b-85-f861843~trusty

And here is my script:

# -*- coding: utf-8 -*-
from yade import pack, plot
young=1e8
compFricDegree = 5
finalFricDegree = 35
mn,mx=Vector3(0,0,0),Vector3(.005,.005,.005)

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
psdSizes=[.00001,.00006,.00008,.0002,.0004,.0005,.0008,.001]
psdCumm=[0,.0175,.025,.4,.5,.7,.85,1]
sp=pack.SpherePack()
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=500)
sp.toSimulation(material='spheres')
triax=TriaxialStressController(
 maxMultiplier=1.001,
 finalMaxMultiplier=1.0001,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)

newton=NewtonIntegrator(damping=0.2)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 triax,
 newton,
]
O.dt=PWaveTimeStep()
triax.goal1=triax.goal2=triax.goal3=-10000

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  print 'porosity', triax.porosity
  if unb<0.01 and abs(-10000-triax.meanStress)/10000<.01:
    break

bodiesToBeDeleted=[]
for b in O.bodies:
 if b.id in range(6):
  continue
 else:
  if b.state.pos[0]<.002:
   bodiesToBeDeleted.append(b)

for b in bodiesToBeDeleted:
 O.bodies.erase(b.id)

[Jérôme Duriez (jduriez)]

Hi,

I launched this script with Yade 2018-04-19.git-79b2edb.

I got the repeated output of "porosity X" with the porosity X-value soon getting negativ, and then decreasing againt until ~ -1e100, then I stopped manually the simulation.

I did not get any segmentation fault. Could you clarify the situation, please ?

[azim (mirzavand)]

Hi jerome,
>>I did not get any segmentation fault. Could you clarify the situation, please ?
that confused me too, I run the same simulation and watch different results. Some times segmentation fault and the other times decreasing the porosity over and over to infinity.

you mentioned porosity, I always see an irrational difference between triaxialStressController.porosity and the porosity(),(in many cases not only in one here) I thought it may be resulted from breaking of particles according to low Modules of elasticity (micro!!) I increased it, for the same simulation, in some cases it works well and other cases leads to segmentaion fault. i referred to the system that these two parameters(porosities)are calculated to figure it out but got no results.

here is the same script writen by Brunu Chareyre,scipt-session.py
I only did some cleaning, please run it. it leads to segmentation fault.
Thanks Jerome.

# -*- coding: utf-8 -*-
from yade import pack
num_spheres=1000
compFricDegree = 30
targetPorosity = 0.43
finalFricDegree = 30
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=5e6
mn,mx=Vector3(0,0,0),Vector3(1,1,1)
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()

clumps=False
if clumps:
 volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
 mean_rad = pow(0.09*volume/num_spheres,0.3333)
 c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
 sp.makeClumpCloud(mn,mx,[c1],periodic=False)
 sp.toSimulation(material='spheres')
 O.bodies.updateClumpProperties()
else:
 sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1)
 O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]
triax.goal1=triax.goal2=triax.goal3=-10000

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
    break
internalCompaction=False
bodiesToBeDeleted=[]
for b in O.bodies:
 if b.id in range(6):
  continue
 else:
  if b.state.pos[0]<.1:
   bodiesToBeDeleted.append(b)

for b in bodiesToBeDeleted:
 O.bodies.erase(b.id)

[Jérôme Duriez (jduriez)]

Thanks, I could reproduce the crash with this script.

Solution
********

The change [1] should fix your problems. Obviously you need to update your yade version to enjoy this change.

In terms of package versions, I do not know when this change will be reflected in yadedaily, considering in particular we currently have compiling problems. [2]

If you use a local source code version, the best is for you to apply locally the same change as [1] (and to recompile.. :-) )

Reason
******

As told by Jan in the Launchpad question thread, the bug came from growParticles() where interactions operations were performed before checking whether the interaction actually is real. Because these soon-to-be-erased interactions (involving deleted bodies) were accessed, deleted bodies, were tried to be accessed as well, leading to the crash.

Advice
******

The MWE below was enough to trigger the exact same bug:
-------------
O.bodies.append(sphere(Vector3(0,0,0),1))
O.bodies.append(sphere(Vector3(0,0,2),1))
O.step()
O.bodies.erase(1)
growParticles(1.5)
-------------
and would have maybe (??) led to a faster bug resolution, enticing more people to look faster at the issue.

I understand showing real MWE may require some YADE practice.. Nevertheless, it's worth trying it, and here you had some elements from Jan narrowing down the problem to Shop::growParticles()

This being said, I just hope this will motivate you to dig into the source code, happy YADE-ing ! ;-)

[1] https://github.com/yade/trunk/commit/bd1089e77eae83c41c4ce0c60873e1e234ba9c89
[2] https://yade-dem.org/buildbot/one_line_per_build

[azim (mirzavand)]

Thank you Jerome,
>>This being said, I just hope this will motivate you to dig into the source code, happy YADE-ing ! ;-)
I really want to, but I don't know what to do.I'm searching and studying in order to learn this things.
I will appreciate you if you show me the way or instruct me in the procedure.
thanks again for the help.

[Jérôme Duriez (jduriez)]

For playing with source code, you may refer to https://yade-dem.org/doc/installation.html#source-code for source code download and compilation (i.e. YADE installation), then the Programmer's manual https://yade-dem.org/doc/prog.html.

[azim (mirzavand)]

Hi Jerome
I can delete particles without any problem in my newest yadedaily version.
it is fixed, am I right?
sorry for asking such a beginner question.!
thanks

[Jérôme Duriez (jduriez)]

Hi,

If it now works while it used not to, I think we can say it's fixed ;-)

A recently updated yadedaily should indeed reflect the above-mentioned commit https://github.com/yade/trunk/commit/bd1089e77eae83c41c4ce0c60873e1e234ba9c89
(the compilation problems I mentioned in #3 are solved now --- and, actually, they never affected yadedaily packages anyway)

Enjoy,

Jérôme

[azim (mirzavand)]

Hi Jérôme
Thanks.