negative volume for an ordinary pore, stack smashing detected

Hello,

I'm trying to duplicate the hydraulic fracturing process in Papachristos et al., 2017[1]. As a very newbie of Yade, I just start from the existing example. The fractured rock is modified from example/jointedCohesiveFrictionalPM: remove the gravity engine. Then I put a point Flux on the fracture plane (actually I'm not quite sure it is the right way to model the hydraulic fracture). However, when I use the FlowEngine, it gives me the following error:

negative volume for an ordinary pore (temp warning, should still be safe)
negative volume for an ordinary pore (temp warning, should still be safe)
negative volume for an ordinary pore (temp warning, should still be safe)
negative volume for an ordinary pore (temp warning, should still be safe)
negative volume for an ordinary pore (temp warning, should still be safe)
negative volume for an ordinary pore (temp warning, should still be safe)
negative volume for an ordinary pore (temp warning, should still be safe)
*** stack smashing detected ***: <unknown> terminated
Aborted (core dumped)

I searched around the forum, even though there were similar questions posted before, but I have no luck finding a effective solution for me. Here is the script (the parallellepiped_10_persistentPlane30Deg and jointDip30_jointFrict20 files are same as example/jointedCohesiveFrictionalPM):

==================================
from __future__ import print_function
O=Omega()
from yade import plot, pack,utils,ymport

#### controling parameters
packing='parallellepiped_10_persistentPlane30Deg'
smoothContact=True
jointFrict=radians(20)
jointDil=radians(0)
output='jointDip30_jointFrict20'
maxIter=10000

#### Import of the sphere assembly
def sphereMat(): return JCFpmMat(type=1,young=1e8,frictionAngle=radians(30),density=3000,poisson=0.3,tensileStrength=1e6,cohesion=1e6,jointNormalStiffness=1e7,jointShearStiffness=1e7,jointCohesion=1e6,jointFrictionAngle=jointFrict,jointDilationAngle=jointDil)
print("\n In case of errors please look at README about generating parallellepiped_10_persistentPlane30Deg.spheres file\n")
O.bodies.append(ymport.text(packing+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

## preprocessing to get dimensions of the packing
dim=utils.aabbExtrema()
dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

## preprocessing to get spheres dimensions
R=0
Rmax=0
numSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   numSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
Rmean=R/numSpheres

inFile=open(packing+'_jointedPM.spheres','r')
for line in inFile:
 if '#' in line : continue
 id = int(line.split()[0])
 onJ = int(line.split()[1])
 nj = int(line.split()[2])
 j11 = float(line.split()[3])
 j12 = float(line.split()[4])
 j13 = float(line.split()[5])
 j21 = float(line.split()[6])
 j22 = float(line.split()[7])
 j23 = float(line.split()[8])
 j31 = float(line.split()[9])
 j32 = float(line.split()[10])
 j33 = float(line.split()[11])
 O.bodies[id].state.onJoint=onJ
 O.bodies[id].state.joint=nj
 O.bodies[id].state.jointNormal1=(j11,j12,j13)
 O.bodies[id].state.jointNormal2=(j21,j22,j23)
 O.bodies[id].state.jointNormal3=(j31,j32,j33)
inFile.close

#### joint strength degradation
for i in O.interactions:
 if i.phys.isOnJoint :
  if i.phys.isCohesive:
   i.phys.isCohesive=False
   i.phys.FnMax=0.
   i.phys.FsMax=0.

newton=NewtonIntegrator(damping=0.2)

interactionRadius=1.
O.engines=[

 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='is2aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=smoothContact,label='interactionLaw'),Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
    FlowEngine(dead=1,label="flow"),
 GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.8),
 newton,
]

from yade import qt
v=qt.Controller()
v=qt.View()

flow.dead=0
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=1
flow.boundaryUseMaxMin=[1,1,1,1,1,1]
flow.fluidBulkModulus=2.2e9
flow.fluidRho=1000
O.dt=0.1e-4
O.dynDt=False
flow.updateTriangulation=True
flow.bndCondIsPressure=[1,1,1,1,1,1]
flow.bndCondValue=[0,0,0,0,0,0]
flow.imposeFlux=[(0,1,0),8.3e-5]
O.run(1,1)
=================================
Can anyone help me with this? Or any tips or tricks.

Thank you