yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #19639
Re: [Question #680642]: Segmentation fault(core dumped)
Question #680642 on Yade changed:
https://answers.launchpad.net/yade/+question/680642
Status: Needs information => Open
Cloud gave more information on the question:
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)
--
You received this question notification because your team yade-users is
an answer contact for Yade.
