Segmentation fault(core dumped)

We need a MWE to reproduce this. Not part of it. https://yade-dev.gitlab.io/trunk/user.html#questions-and-answers

Hi, Janek,

Sorry, my MWE is as follow:

# -*- coding: utf-8 -*-
from yade import pack
import math

mn,mx=Vector3(0,0,0),Vector3(1,2,1)
young = 0.208e9
poisson = 0.3
compFricDegree = 15.0
finalFricDegree = 30.0
porosity = 0.95
num_spheres = 10000
confiningS = 50e3
isoconfining = 5e3
damp = 0.2
rate = -0.02
AxialStrainLimit = 0.3
key="triax_"
cd = 1 # cd=1 - consolidate drained cd=0- undrained
stabilityThreshold = 0.001
thick = 0.01

O.materials.append(FrictMat(
    young=young,
    poisson=poisson,
    density=2650,
    frictionAngle=radians(compFricDegree),
    label='spheres')
    )

O.materials.append(FrictMat(
    young=15e6,
    poisson=0.4,
    frictionAngle=0,
    density=0,
    label='walls')
    )

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

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,1.0/7.0,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

O.dt=.1*utils.PWaveTimeStep()
O.usesTimeStepper=True

triax=TriaxialStressController(
    maxMultiplier=1.+2e4/young,
    finalMaxMultiplier=1.+2e3/young,
    thickness = thick,
    stressMask = 7,
    internalCompaction=True,
)

newton = NewtonIntegrator(damping=damp)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys()],
        [Law2_ScGeom_MindlinPhys_Mindlin()]
    ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
    triax,
    TriaxialStateRecorder(iterPeriod=1000,file='TriaxialRecorder'+key,truncate=1),
    newton,
]

triax.goal1=triax.goal2=triax.goal3=-isoconfining #(5kpa)

while 1:
    O.run(1000, True)
    unb=unbalancedForce()
    print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
    if unb < stabilityThreshold and abs(-isoconfining-triax.meanStress)/isoconfining<0.001:
        break

O.save('confinedState'+key+'.yade.gz')
print "### Isotropic state saved ###"

####### mini particle radius#####
minr=0
for b0 in O.bodies:
    if not isinstance(b0.shape,Sphere):continue
    if (b0.shape.radius)<=minr or (minr==0.0):minr=(b0.shape.radius)

#################################
##### Particle breakage #####
#################################
import math
def selectBodiesToCrush():
    breakparticle = []
    for b in O.bodies:
        if not isinstance(b.shape,Sphere):continue # exclude walls
        sigma_lim0 = 400e4
        diamter0 = 2.0
        diamter = 2.0 * b.shape.radius
        m = 10.0
        sigma_lim = sigma_lim0 * math.pow((diamter/diamter0), (-3.0/m))

        Fn_max = 0
        for i in b.intrs():
            Fn_abs = math.sqrt(i.phys.normalForce[0]*i.phys.normalForce[0]+i.phys.normalForce[1]*i.phys.normalForce[1]+i.phys.normalForce[2]*i.phys.normalForce[2])
            if Fn_abs > Fn_max:Fn_max = Fn_abs
            if Fn_abs == Fn_max:
                b1,b2=[O.bodies[ii] for ii in (i.id1,i.id2)]
                if not isinstance(b1.shape,Sphere):continue
                if not isinstance(b2.shape,Sphere):continue
                r_apostrophe = math.pow(1/b1.shape.radius + 1/b2.shape.radius,-1.0)
                E_apostrophe = math.pow(((1-b1.material.poisson*b1.material.poisson)/b1.material.young)+((1-b2.material.poisson*b2.material.poisson)/b2.material.young),-1.0)
                rH = math.pow((3*r_apostrophe)/(4*E_apostrophe),1.0/3.0)
                F_lim = math.pow((sigma_lim * math.pi * rH * rH),3.0)
                if Fn_max >= F_lim:
                    breakparticle.append(b)
    return breakparticle

def createChildren(bb):
    p = bb.state.pos
    r = bb.shape.radius
    s1 = sphere(p+(0,0.5359*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s2 = sphere(p+(-0.4641*r,-0.2679*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s3 = sphere(p+(0.4641*r,-0.2679*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s4 = sphere(p+(-0.6720*r,0.3880*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s5 = sphere(p+(0.6720*r,0.3880*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s6 = sphere(p+(0,-0.7760*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s7 = sphere(p+(0,0,0.6339*r),0.3659*r,material='spheres',color=(0,0,1))
    s8 = sphere(p+(0,0,-0.6339*r),0.3659*r,material='spheres',color=(0,0,1))
    s9 = sphere(p+(0,-0.5942*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s10 = sphere(p+(0.5146*r,0.2671*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s11 = sphere(p+(-0.5146*r,0.2971*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s12 = sphere(p+(0,-0.5942*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s13 = sphere(p+(0.5146*r,0.2971*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s14 = sphere(p+(-0.5146*r,0.2971*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    return s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14

def replace():
    bodiesToDelete = selectBodiesToCrush()
    for bb in bodiesToDelete:
        if isinstance(bb.shape,Sphere) and bb.shape.radius > minr/4:
            children = createChildren(bb)
            O.bodies.erase(bb.id)
            O.bodies.append(children)

triax.internalCompaction=False

setContactFriction(radians(finalFricDegree))

if cd:
    triax.stressMask = 5
    triax.goal2 = rate
    triax.goal1 = -confiningS
    triax.goal3 = -confiningS
else:
    triax.stressMask = 0
    triax.goal2 = rate
    triax.goal1 = triax.goal3 = -rate/2

newton.damping=0.1

O.saveTmp()

from yade import plot

def history():
    plot.addData(
            e22=-triax.strain[1],
            e11=-triax.strain[0],
            e33=-triax.strain[2],
       ev=triax.volumetricStrain,
      s11=-triax.stress(triax.wall_right_id)[0],
      s22=-triax.stress(triax.wall_top_id)[1],
      s33=-triax.stress(triax.wall_front_id)[2],
            args=(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_top_id)[1]-triax.stress(triax.wall_front_id)[2])/3,
            eta = 3*(-triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_right_id)[0])/(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_top_id)[1]-triax.stress(triax.wall_front_id)[2]),
            devi = -triax.stress(triax.wall_top_id# -*- coding: utf-8 -*-
from yade import pack
import math

nRead=utils.readParamsFromTable(
        young=0.208e9,
        poisson = 0.3,
        compFricDegree=15.0,
        alphaKr = 0.2,
        Eta = 1.0,
        finalFricDegree=30.0,
        porosity = 0.95,
        num_spheres=10000,
        confiningS=50e3,
        isoconfining = 5e3,
        AxialStrainLimit=0.3,
        rate=-0.02,
        damp=0.2,
        unknownOk = True
)
from yade.params import table

mn,mx=Vector3(0,0,0),Vector3(1,2,1)
young = table.young
poisson = table.poisson
compFricDegree = table.compFricDegree
alphaKr = table.alphaKr
Eta = table.Eta
finalFricDegree = table.finalFricDegree
porosity = table.porosity
num_spheres = table.num_spheres
confiningS = table.confiningS
isoconfining = table.isoconfining
damp = table.damp
rate = table.rate
AxialStrainLimit=table.AxialStrainLimit
key="triax_"
cd = 1 # cd=1 - consolidate drained cd=0- undrained
stabilityThreshold = 0.001
thick = 0.01

O.materials.append(FrictMat(
    young=young,
    poisson=poisson,
    density=2650,
    frictionAngle=radians(compFricDegree),
    label='spheres')
    )

O.materials.append(FrictMat(
    young=15e6,
    poisson=0.4,
    frictionAngle=0,
    density=0,
    label='walls')
    )

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

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,1.0/7.0,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

O.dt=.1*utils.PWaveTimeStep()
O.usesTimeStepper=True

triax=TriaxialStressController(
    maxMultiplier=1.+2e4/young,
    finalMaxMultiplier=1.+2e3/young,
    thickness = thick,
    stressMask = 7,
    internalCompaction=True,
)

newton = NewtonIntegrator(damping=damp)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys()],
        [Law2_ScGeom_MindlinPhys_Mindlin()]
    ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
    triax,
    TriaxialStateRecorder(iterPeriod=1000,file='TriaxialRecorder'+key,truncate=1),
    newton,
]

triax.goal1=triax.goal2=triax.goal3=-isoconfining #(5kpa)

while 1:
    O.run(1000, True)
    unb=unbalancedForce()
    print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
    if unb < stabilityThreshold and abs(-isoconfining-triax.meanStress)/isoconfining<0.001:
        break

O.save('confinedState'+key+'.yade.gz')
print "### Isotropic state saved ###"

####### mini particle radius#####
minr=0
for b0 in O.bodies:
    if not isinstance(b0.shape,Sphere):continue
    if (b0.shape.radius)<=minr or (minr==0.0):minr=(b0.shape.radius)

#################################
##### Particle breakage #####
#################################
import math
def selectBodiesToCrush():
    breakparticle = []
    for b in O.bodies:
        if not isinstance(b.shape,Sphere):continue # exclude walls
        sigma_lim0 = 400e4
        diamter0 = 2.0
        diamter = 2.0 * b.shape.radius
        m = 10.0
        sigma_lim = sigma_lim0 * math.pow((diamter/diamter0), (-3.0/m))

        Fn_max = 0
        for i in b.intrs():
            Fn_abs = math.sqrt(i.phys.normalForce[0]*i.phys.normalForce[0]+i.phys.normalForce[1]*i.phys.normalForce[1]+i.phys.normalForce[2]*i.phys.normalForce[2])
            if Fn_abs > Fn_max:Fn_max = Fn_abs
            if Fn_abs == Fn_max:
                b1,b2=[O.bodies[ii] for ii in (i.id1,i.id2)]
                if not isinstance(b1.shape,Sphere):continue
                if not isinstance(b2.shape,Sphere):continue
                r_apostrophe = math.pow(1/b1.shape.radius + 1/b2.shape.radius,-1.0)
                E_apostrophe = math.pow(((1-b1.material.poisson*b1.material.poisson)/b1.material.young)+((1-b2.material.poisson*b2.material.poisson)/b2.material.young),-1.0)
                rH = math.pow((3*r_apostrophe)/(4*E_apostrophe),1.0/3.0)
                F_lim = math.pow((sigma_lim * math.pi * rH * rH),3.0)
                if Fn_max >= F_lim:
                    breakparticle.append(b)
    return breakparticle

def createChildren(bb):
    p = bb.state.pos
    r = bb.shape.radius
    s1 = sphere(p+(0,0.5359*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s2 = sphere(p+(-0.4641*r,-0.2679*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s3 = sphere(p+(0.4641*r,-0.2679*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s4 = sphere(p+(-0.6720*r,0.3880*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s5 = sphere(p+(0.6720*r,0.3880*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s6 = sphere(p+(0,-0.7760*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s7 = sphere(p+(0,0,0.6339*r),0.3659*r,material='spheres',color=(0,0,1))
    s8 = sphere(p+(0,0,-0.6339*r),0.3659*r,material='spheres',color=(0,0,1))
    s9 = sphere(p+(0,-0.5942*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s10 = sphere(p+(0.5146*r,0.2671*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s11 = sphere(p+(-0.5146*r,0.2971*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s12 = sphere(p+(0,-0.5942*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s13 = sphere(p+(0.5146*r,0.2971*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s14 = sphere(p+(-0.5146*r,0.2971*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    return s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14

def replace():
    bodiesToDelete = selectBodiesToCrush()
    for bb in bodiesToDelete:
        if isinstance(bb.shape,Sphere) and bb.shape.radius > minr/4:
            children = createChildren(bb)
            O.bodies.erase(bb.id)
            O.bodies.append(children)

triax.internalCompaction=False

setContactFriction(radians(finalFricDegree))

if cd:
    triax.stressMask = 5
    triax.goal2 = rate
    triax.goal1 = -confiningS
    triax.goal3 = -confiningS
else:
    triax.stressMask = 0
    triax.goal2 = rate
    triax.goal1 = triax.goal3 = -rate/2

newton.damping=0.1

O.saveTmp()

from yade import plot

def history():
    plot.addData(
            e22=-triax.strain[1],
            e11=-triax.strain[0],
            e33=-triax.strain[2],
       ev=triax.volumetricStrain,
      s11=-triax.stress(triax.wall_right_id)[0],
      s22=-triax.stress(triax.wall_top_id)[1],
      s33=-triax.stress(triax.wall_front_id)[2],
            args=(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_top_id)[1]-triax.stress(triax.wall_front_id)[2])/3,
            eta = 3*(-triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_right_id)[0])/(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_top_id)[1]-triax.stress(triax.wall_front_id)[2]),
            devi = -triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_right_id)[0],
            e=triax.porosity/(1-triax.porosity),
            i=O.iter,
            )

if 1:
    O.engines=O.engines[0:5]+[PyRunner(iterPeriod=500,command='history()',label='recorder')]+O.engines[5:7]
    O.engines=O.engines+[PyRunner(command='replace()',iterPeriod=100,label ='replace')]

############################
####### Start loading ##########
############################
while 1:
 O.run(1000,True)
 if O.iter%100 == 0:
    print "step = ", O.iter,"number of particles:",len(O.bodies)
         print "axialstrain = ",abs(triax.strain[1])
        if abs(triax.strain[1]) > AxialStrainLimit:
  break

plot.plots={'e22':('ev',),' e22 ':('eta')}
plot.plot()

plot.saveDataTxt('results'+key)
utils.waitIfBatch()
)[1]+triax.stress(triax.wall_right_id)[0],
            e=triax.porosity/(1-triax.porosity),
            i=O.iter,
            )

if 1:
    O.engines=O.engines[0:5]+[PyRunner(iterPeriod=500,command='history()',label='recorder')]+O.engines[5:7]
    O.engines=O.engines+[PyRunner(command='replace()',iterPeriod=100,label ='replace')]

###################################
####### Start loading ##########
###################################
while 1:
 O.run(1000,True)
 if O.iter%100 == 0:
    print "step = ", O.iter,"number of particles:",len(O.bodies)
         print "axialstrain = ",abs(triax.strain[1])
        if abs(triax.strain[1]) > AxialStrainLimit:
  break

plot.plots={'e22':('ev',),' e22 ':('eta')}
plot.plot()

plot.saveDataTxt('results'+key)

please make sure that yuor MWE reproduces the issue.

```
  File "./test.py", line 180
    from yade import pack
```
SyntaxError: invalid syntax

also make sure that it has no mixed indentation (use either tabs or spaces) and use brackets ( ) for your print function.

It may be a time-step problem. After making some modification to time-step, it runs normally.