Segmentation Fault - PFV2
Dear all,
I am facing an error of segmentation fault with this script.
Not sure if it is hardware limitation or if it is related to tessellation issue, since this packing is very peculiar.
I have checked similar questions but couldn't find a way out.
Do you have an idea of size ration between particle diameter and aggregate diameter.
Cheers,
# encoding: utf-8
# This script demonstrates a simple case of drainage simulation using the "2PFV" two-phase model implemented in UnsaturatedEngine.
# The script was used to generate the result and supplementary material (video) of [1]. The only difference is the problem size (40k particles in the paper vs. 1k (default) in this version)
# [1] Yuan, C., & Chareyre, B. (2017). A pore-scale method for hydromechanical coupling in deformable granular media. Computer Methods in Applied Mechanics and Engineering, 318, 1066-1079. (http://
from yade import plot
from yade import ymport
from yade import bodiesHandling
from yade import export
from yade import utils
import matplotlib; matplotlib.
from yade import pack
import pylab
from numpy import *
utils.readParam
from yade.params import table
seed=table.seed
#num_spheres=
compFricDegree = table.compFricD
confiningS=-1e6
## creat a packing with a specific particle side distribution (PSD)
#psdSizes,
sp=pack.
sp1=pack.
mn,mx=Vector3(
mnn,mxx=
mnc,mxc=
#Coating
sp.makeCloud(
#Matrix
sp1.makeCloud(
O.materials.
O.materials.
O.materials.
## create walls around the packing
walls=aabbWalls
wallIds=
#O.bodies.
#for center,rad in sp:
# print (center,rad)
#O.bodies.
#O.bodies.
for p in sp1:
a=p[0]
# print (a)
f1=O.bodies.
triax=TriaxialS
# wall_back_
# wall_bottom_
# wall_front_
# wall_left_
# wall_right_
# wall_top_
internalCompac
goal1=confiningS,
goal2=confiningS,
goal3=confiningS,
max_vel=5,
label="triax"
)
newton=
O.engines=[
ForceResetter(),
InsertionSortC
InteractionLoop(
[Ig2_
[Ip2_
[Law2_
# [Ig2_Sphere_
# [Ip2_FrictMat_
# [Law2_ScGeom_
),
# GlobalStiffness
triax,
# VTKRecorder(
newton
]
O.dt=0.
O.dynDt=False
while 1:
O.run(1000,True)
unb=unbalance
# print(triax.
if unb<0.01 and abs(triax.
break
#######
## REACH NEW EQU. STATE ###
#######
finalFricDegree = 30 # contact friction during the deviatoric loading
#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalC
# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFrict
while 1:
O.run(1000,True)
unb=unbalance
if unb<0.01 and abs(triax.
break
triax.depth0=
triax.height0=
triax.width0=
O.save(
O.run(1000,True)
ei0=-triax.
si0=-triax.
from yade import plot
O.engines=
def history():
plot.addData(
s11=-
s22=-
s33=-
pc=-unsat.
sw=unsat.
i=O.iter
)
plot.plots=
plot.plot()
#######
## Drainage Test under oedometer conditions ###
#######
##oedometer conditions
triax.stressMask=2
triax.goal1=
goalTop=
#print(goalTop)
triax.goal2=goalTop
triax.wall_
recorder.dead=0
##Instantiate a two-phase engine
unsat=TwoPhaseF
meanDiameter=
##set boundary conditions, the drainage is controlled by decreasing W-phase pressure and keeping NW-phase pressure constant
unsat.bndCondIs
unsat.bndCondVa
unsat.isPhaseTr
unsat.initializ
unsat.surfaceTe
#start invasion, the data of normalized pc-sw-strain will be written into pcSwStrain.txt
f1=open(
f2=open(
f3=open(
f4=open(
f5=open(
for pg in arange(
#increase gaz pressure at the top boundary
unsat.
#compute the evolution of interfaces
unsat.invasion()
#save the phases distribution in vtk format, to be displayed by paraview
# unsat.savePhase
#compute and apply the capillary forces on each particle
unsat.
for b in O.bodies:
O.forces.
if pg==0.60001:
cels=
celsW1 = [0.0]*cels
celsV1 = [0.0]*cels
celsBar1 = [0.0]*cels
for ii in range(cels):
celsW1=
celsV1=
# celsBar1=
celsBar1=
f1.
f1.close()
if pg==2.30001:
cels2=
celsW2 = [0.0]*cels2
celsV2 = [0.0]*cels2
celsBar2 = [0.0]*cels2
for jj in range(cels2):
celsW2=
celsV2=
celsBar2=
f2.
f2.close()
if pg==3.10001:
cels3=
celsW3 = [0.0]*cels3
celsV3 = [0.0]*cels3
celsBar3 = [0.0]*cels3
for gg in range(cels3):
celsW3=
celsV3=
celsBar3=
f3.
f3.close()
if pg==9.90001:
cels4=
celsW4 = [0.0]*cels4
celsV4 = [0.0]*cels4
celsBar4 = [0.0]*cels4
for hh in range(cels4):
celsW4=
celsV4=
celsBar4=
f4.
f4.close()
#reac
while 1:
O.run(
unb=
if unb<0.001:
break
f5.write(
f5.close()
Question information
- Language:
- English Edit question
- Status:
- Expired
- For:
- Yade Edit question
- Assignee:
- No assignee Edit question
- Last query:
- Last reply: