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

Asked by gaoxuesong on 2019-05-06

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.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:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2019-05-08
Last query:
2019-05-08
2019-05-07
 Jan Stránský (honzik) said on 2019-05-06: #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 
- please use a 'seed' parameter to sp.makeCloud

thanks
Jan

 gaoxuesong (260582472-9) said on 2019-05-07: #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
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()

 gaoxuesong (260582472-9) said on 2019-05-07: #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. Jan Stránský (honzik) said on 2019-05-07: #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 
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
harmEngine.A = ...
harmEngine.f = ...
elif ...
###

could you modify your code and test if it works?

cheers
Jan

 gaoxuesong (260582472-9) said on 2019-05-08: #5

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

 gaoxuesong (260582472-9) said on 2019-05-08: #6

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

To post a message you must log in.