Small density sphere in micron size makes segment fault (core dumped)

Asked by gaoxuesong

In my case, the sphere used is in micron size. When i use the real material property, especially the density, the calculation crashes after some time(Segment fault (core dumped)). Then i try to magnify the density and find the run goes well when the density amplified by one million times.
Following is the testing code where the parameter masscof is the density amplification factor,

from yade import geom,pack
import math
from yade import qt
import os

os.system('rm -rf plt;mkdir ./plt/')
os.system('rm -rf out;mkdir ./out/')

# PhysicalParameters
muS = 0.57735 # Friction coefficient (30 degree)
FricAngleS = math.atan(muS)

### mass amplification factor ###
masscof = 1e6

# PhysicalParameters
matSph = CohFrictMat(
 density = 7990*masscof,
 young = 193e9,
 poisson = 0.3,
 frictionAngle = FricAngleS,
 momentRotationLaw = True)
SMat = O.materials.append(matSph)

# create empty sphere packing
sp = pack.SpherePack()
# generate spheres1 with a certain diameter distribution
sp.makeCloud((0,0,0.35e-3), (1.0e-3,1.0e-3,1.4e-3), psdSizes=[0.012e-3,0.0186e-3,0.0312e-3,0.0484e-3,0.053e-3], psdCumm=[0,0.1,0.5,0.9,1.0])
# add the sphere pack to the simulation
sp.toSimulation(material=SMat)

# create rectangular box1(Left) from facets
O.bodies.append(geom.facetBox((0.5e-3,0.5e-3,0.5e-3),(0.5e-3,0.5e-3,0.5e-3), wallMask=31, material=SMat))
O.bodies.append(geom.facetBox((0.5e-3,0.5e-3,0.15e-3),(0.5e-3,0.5e-3,0.2e-3), wallMask=63, color=(0,1,0), wire=False, material=SMat))

O.dt = 0.85*utils.PWaveTimeStep()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
       # handle sphere+sphere and facet+sphere collisions
         [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
         [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
          [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
    ),
 NewtonIntegrator(damping=0.75, exactAsphericalRot=True, gravity=(0,0,-9.81)),
    qt.SnapshotEngine(fileBase="./plt/",iterPeriod=200,label='snapshot'),
    VTKRecorder(iterPeriod=200, recorders=['spheres','colors'], fileName='./out/'),
]

#show geometry:
qtr = qt.Renderer()
qtr.bgColor = [1,1,1]
qt.Controller()
v = qt.View()
#v.center()
v.axes=False
v.viewDir=Vector3(0,1,0)
v.eyePosition=Vector3(1.8e-3,-5.2e-3,1.0e-3)

#O.run()

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,

thanks for the code. I have tried the script for a while without any problem. Please provide more info:
- At which iteration the problem occur?
- is the SnapshotEngine and VTKRecorder needed to reproduce the problem? If no, please remove them to make the script a true MWE [1]
- please use a 'seed' parameter to sp.makeCloud

thanks
Jan

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

Revision history for this message
gaoxuesong (260582472-9) said :
#2

Hi, Jan. Yes, you are right. I have tried the above code and no segment fault occurred. Then I add some lines which i thought were unimportant in the past and the problems happened at about 341300 iteration. As you suggested, i leave the snap and VTK code out. Also, to make the crash easy to happen, i amplify the density by a small factor as 1e3.
The following is the code,

from yade import geom,pack
import math
from yade import qt
import os

#os.system('rm -rf plt;mkdir ./plt/')
#os.system('rm -rf out;mkdir ./out/')

# PhysicalParameters
muS = 0.57735 # Friction coefficient (30 degree)
muLB = 0.0001 # Friction coefficient
muF = 0.17632 # Friction coefficient (10 degree)
FricAngleS = math.atan(muS)
FricAngleLB = math.atan(muLB)
FricAngleF = math.atan(muF)

### mass amplification factor ###
masscof = 1e3

# PhysicalParameters
matSph = CohFrictMat(
 density = 7990*masscof,
 young = 193e9,
 poisson = 0.3,
 frictionAngle = FricAngleS,
 momentRotationLaw = True)
SMat = O.materials.append(matSph)

# create empty sphere packing
sp = pack.SpherePack()
# generate spheres1 with a certain diameter distribution
sp.makeCloud((0,0,0.35e-3), (1.0e-3,1.0e-3,1.4e-3), psdSizes=[0.012e-3,0.0186e-3,0.0312e-3,0.0484e-3,0.053e-3], psdCumm=[0,0.1,0.5,0.9,1.0])
# add the sphere pack to the simulation
sp.toSimulation(material=SMat)

# create rectangular box1(Left) from facets
O.bodies.append(geom.facetBox((0.5e-3,0.5e-3,0.5e-3),(0.5e-3,0.5e-3,0.5e-3), wallMask=31, material=SMat))
Cylinder1IDs=O.bodies.append(geom.facetBox((0.5e-3,0.5e-3,0.15e-3),(0.5e-3,0.5e-3,0.2e-3), wallMask=63, color=(0,1,0), wire=False, material=SMat))

################## add these lines, problem happens #####################
################## add these lines, problem happens #####################

def changeValues():

# Standstill for stabilization
   if O.time < 0.2:
       O.engines = O.engines + [
               TranslationEngine(dead=False, translationAxis= [0,0,1], velocity=0.0, ids=Cylinder1IDs)
   ]
       # Vibrate Piston1
   elif 0.2 <= O.time < 0.4:
       O.engines = O.engines + [
               HarmonicMotionEngine(A=(0.0,0.0,1.5e-5), f=(0.0,0.0,3600/60.0), ids=Cylinder1IDs)
   ]
# Vibrate Stop
   elif 0.4 <= O.time < 0.5129:
       O.engines = O.engines + [
               HarmonicMotionEngine(A=(0.0,0.0,0.0), f=(0.0,0.0,0.0), ids=Cylinder1IDs)
   ]

   O.dt = 0.85*utils.PWaveTimeStep()

################## add these lines, problem happens #####################
################## add these lines, problem happens #####################

def print_ite():
    print("ite is %d" % O.iter)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
       # handle sphere+sphere and facet+sphere collisions
         [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
         [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
          [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
    ),
 NewtonIntegrator(damping=0.75, exactAsphericalRot=True, gravity=(0,0,-9.81)),
     PyRunner(iterPeriod=100,command='changeValues()'),
    PyRunner(iterPeriod=100,command='print_ite()'),
# qt.SnapshotEngine(fileBase="./plt/",iterPeriod=200,label='snapshot'),
# VTKRecorder(iterPeriod=200, recorders=['spheres','colors'], fileName='./out/'),
]

#show geometry:
qtr = qt.Renderer()
qtr.bgColor = [1,1,1]
qt.Controller()
v = qt.View()
#v.center()
v.axes=False
v.viewDir=Vector3(0,1,0)
v.eyePosition=Vector3(1.8e-3,-5.2e-3,1.0e-3)

O.run()

Revision history for this message
gaoxuesong (260582472-9) said :
#3

Sorry to add the ubuntu and yade version used.
The ubuntu is Ubuntu 18.04.1 LTS, and the yade is 2018.02b.
Much thanks.

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

Hello, just running your code and waiting for the result. However, the changeValues() function should not be used in this form (!)

You are adding new TranslationEngine every 100 iterations. Now I have O.iter=288300 and if I print len(O.engines), I got 2893.. Furthermore, what is the meaning of the TranslationEngine with velocity=0.0?

A proper way could be something like:
###
O.engines = [
   ...
   TranslationEngine(...,dead=False,label='transEngine'), # just one engine
   HarmonicMotionEngine(...,dead=True,label='harmEngine'), # just one engine
   NewtonIntegrator(...), # consider pacing it after kinematic engines [1]
   PyRunner(iterPeriod=100,command='changeValues()'),
]

def changeValues():
   if O.time < 0.2:
      pass # probably initialized above
   elif O.time < 0.4: # no need of 0.2 <= O.time, because if so, it would have finished in above branch
      transEngine.dead = True
      harmEngine.dead = False
      harmEngine.A = ...
      harmEngine.f = ...
   elif ...
###

could you modify your code and test if it works?

cheers
Jan

[1] https://yade-dev.gitlab.io/trunk/user.html#base-engines

Revision history for this message
gaoxuesong (260582472-9) said :
#5

Thanks. Jan. The adapted code runs well. I think the problem is the increasing length of O.engines.

Revision history for this message
gaoxuesong (260582472-9) said :
#6

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