PeriodicFlowEngine - periodicity is broken

Other problems I found:
- I was assuming that the walls were ids 0 and 1 in your problem but it is not the case. You need to give their ids in flow.wallIds=[-1,-1,1,0,-1,-1]; or you append them before the spheres so that they really get ids 0 and 1 (that's what I did here)
- the duplicateThreshold value is a length, and I see that your period size is around 12, so that flow.duplicateThreshold=0.5 is still a small thing, when I set it equal to 2.5 the problem disappears. This parameter (the thickness of a layer in which periodic images of the particles are inserted all around the main periodic cell) should be calculated automatically in a robust way but I still don't know how exactly...

What I end up with is below.
Did you get an initial script by someone else, by the way?

Bruno

# -*- coding: utf-8 -*-

from yade import pack, plot

O.periodic=True

radius=0.3
rRelFuzz=0.1
numParticles=2000
numCohParticles=100
porosity=0.3
height=40*radius
length=25*radius
width=25*radius
youngSpheres=1e6
poissonSpheres=0.5
frictDegreeSpheres=30
densitySpheres=2000
youngWalls=1e6
poissonWalls=0.5
frictDegreeWalls=30
densityWalls=7000
thickness=0.01
maxLoad=1000.*length*width
coh=1000000.
gap=0.01*length

O.cell.hSize=Matrix3(length,0.,0.,
       0.,3*height,0.,
       0.,0.,width)

O.materials.append(CohFrictMat(young=youngSpheres,poisson=poissonSpheres,frictionAngle=radians(frictDegreeSpheres),density=densitySpheres,isCohesive=False,label='spheres'))
O.materials.append(CohFrictMat(young=youngSpheres,poisson=poissonSpheres,frictionAngle=radians(frictDegreeSpheres),density=densitySpheres,fragile=False,normalCohesion=coh,shearCohesion=poissonSpheres*coh,label='cohSpheres'))
O.materials.append(CohFrictMat(young=youngWalls,poisson=poissonWalls,frictionAngle=radians(frictDegreeWalls),density=densityWalls,fragile=False,normalCohesion=coh,shearCohesion=poissonWalls*coh,label='walls'))

bottomWall = utils.box(center=(length/2.0,height-thickness/2.0,width/2.0), extents=(length*1000.0,thickness/2.0,width*1000.0) ,fixed=True,wire=False,material='walls')
topWall = utils.box(center=(length/2.0,2.1*height-thickness/2.0,width/2.0), extents=(length*1000.0,thickness/2.0,width*1000.0) ,fixed=True,wire=False,material='walls')
O.bodies.append([bottomWall,topWall])

sp1=pack.SpherePack()
sp1.makeCloud((0.,height+1.1*radius,0.),(length,2.*height,width),radius,rRelFuzz,numParticles,periodic=True,porosity=porosity)
O.bodies.append([utils.sphere(s[0],s[1],color=(1,1,1),material='spheres') for s in sp1])

#pred=pack.inAlignedBox(Vector3(0.,2.0*height,0.),Vector3(length,2.1*height,width))
#sp2=pack.regularHexa(pred,radius=radius,gap=gap,material='cohSpheres')
#O.bodies.append(sp2)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],allowBiggerThanPeriod=True),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=False,setCohesionOnNewContacts=True)],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 PeriodicFlowEngine(dead=1,label="flow"),
 NewtonIntegrator(gravity=(0,0,0),damping=0.2),
 PyRunner(command='loading()',iterPeriod=100,label='stop')
]

def loading():
   topWall.state.vel=(1,0,0)
   if abs(O.forces.f(topWall.id)[1])<maxLoad:
    topWall.state.vel=(1,-0.1,0)

O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=100)]

def addPlotData():
   force=O.forces.f(topWall.id)[1]
   plot.addData(force=force,i=O.iter)

plot.plots={'i':('force',)}
plot.plot()

flow.dead=0
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.eps=0.04*radius
flow.viscosity=1
#flow.debug=1
flow.bndCondIsPressure=[0,0,0,1,0,0]
flow.bndCondValue=[0,0,0,0,0,0]
flow.duplicateThreshold=1.2
flow.boundaryUseMaxMin=[0,0,1,1,0,0]
flow.wallIds=[-1,-1,1,0,-1,-1]
O.dt=0.1e-3
O.dynDt=True