Different stress values from different Yade versions when running PeriodicFlowEngine

I’m working on a Fluid-DEM coupling simulation script and trying to export two things : (1) the total stress acting at the boundary and (2) the effective stress with changing pore water pressure.

What the script will do :
State I : Execute an iso-compression using PeriTriaxController engine
State II : The PeriTriaxController engine is switched off, the PeriodicFlowEngine is then applied with fluid pressure increasing by stepping value Delta_P = 100 kPa after every 1000 iterations. Once the fluid pressure reaches 500 kPa , the simulation stops and a stress plot will be given.

The script runs fine without any technical error ; however, when running with YADE 2018.02b (installed from Ubuntu repo), the output stress is calculated as it’s supposed to be : the boundary stress is constant while the effective stress decreases according to the imposed fluid pressure; but when running with Yadedaily (version Yade 20200701-4022~5cf8771~bionic1 at spoken time) , I receive a different behaviour: the boundary stress increases gradually.

Can you please tell me what changes have been made in YADE between the two versions that produces such different results on the same script? What modification should be done in latest Yade version to obtain the result given in YADE 2018.2b?

More info:
1) I’m running Ubuntu 18.04.04
2) useSolver was set to zero because Yade 2018.2b gives me error about not-compiled CHOLMOD when I used 3, but I think I should ask in another question about this issue.

The MWE script is below:

from __future__ import print_function
from yade import pack,plot,export,utils
import math

sp=pack.SpherePack()
O.periodic=True

# dimensions of sample (fixed by particle size such as L/D~X)
DIAMETER =1.e-1
RADIUS =0.5*DIAMETER
length =10*(DIAMETER)
height =length
width =length
thickness =length/100.

# microproperties
DENS =2600
E =1.e9
P =0.25
compFRIC =1.
FRIC =30.
TENS =0.
COH =0.

# boundary conditions
PI =1.e5 # for sample preparation: isotropic compaction up to PI pressure
PC =PI # confining pressure for assembly preparation
SN =5.e6 # normal stress

# simulation control
DAMPSHEAR =0.
ITER =1e6
OUT ='test'

#### create sample and loading boxes

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

O.materials.append(CohFrictMat(isCohesive=True,density=DENS,young=E,poisson=P,frictionAngle=radians(0.),normalCohesion=1e100,shearCohesion=1e100,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=DENS,young=E,poisson=P,frictionAngle=radians(compFRIC),normalCohesion=TENS,shearCohesion=COH,label='sphereMat'))

upBox = utils.box(center=(length/2.0,2*height+thickness/2.0,width/2.0),orientation=Quaternion(1,0,0,0),extents=(length,thickness/2.,width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2.0,height-thickness/2.0,width/2.0),orientation=Quaternion(1,0,0,0),extents=(length,thickness/2.,width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
O.bodies.append([upBox,lowBox])

sp.makeCloud((0,height+1.5*RADIUS,0),(length,2*height-1.5*RADIUS,width),rMean=RADIUS,rRelFuzz=0.2,periodic=True)
O.bodies.append([utils.sphere(s[0],s[1],color=(0,0.5,1),material='sphereMat') for s in sp])

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

#### engines

flow=PeriodicFlowEngine(
    isActivated=0 # no flow calculation if this is False
    ,useSolver=0 # value 3 not work with YADE 2018.2b package version
    ,defTolerance=+1 # with this value, triangulation is done according to meshUpdateInterval
    ,meshUpdateInterval=1000
    ,duplicateThreshold=0.5
    ,boundaryUseMaxMin=[0,0,1,1,0,0]
    ,wallIds=[-1,-1,1,0,-1,-1]
    ,wallThickness=thickness
    ,bndCondIsPressure=[0,0,0,0,0,0]
        ,bndCondValue=[0,0,0,0,0,0]
    ,permeabilityFactor=1
    ,viscosity=1e-3
    ,fluidBulkModulus=2.2e9
    ,label='flowEng'
    )

O.engines=[
    ForceResetter()
    ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
    ,InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
    )
    ,flow
    ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
    ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep(),label='timeStepper')
    ,NewtonIntegrator(damping=0.3,label='newton')
    ,PyRunner(command='dataRecorder()',iterPeriod=1,label='recData',dead=True)
]

inputP=0
def dataRecorder():
    global inputP
    h=vol=vol_s=nb_s=0.
    h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
    vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
    contactStress=getStress(vol)
    for o in O.bodies:
        if isinstance(o.shape,Sphere):
            nb_s += 1
            vol_s += 4.*pi/3.*(o.shape.radius)**3
        n = 1-vol_s/vol
        nbFrictCont=0.
    for i in O.interactions:
        if i.isReal and i.phys.cohesionBroken:
            nbFrictCont+=1
    plot.addData(
        iter=O.iter
        ,stress_upWall0=abs(O.forces.f(0)[0]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall1=abs(O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall2=abs(O.forces.f(0)[2]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
        ,contactStress00=(contactStress[0,0]),contactStress01=(contactStress[0,1]),contactStress02=(contactStress[0,2]),contactStress10=(contactStress[1,0]),contactStress11=abs(contactStress[1,1]),contactStress12=(contactStress[1,2]),contactStress20=(contactStress[2,0]),contactStress21=(contactStress[2,1]),contactStress22=(contactStress[2,2])
        ,x=O.bodies[0].state.pos[0]
        ,pf=inputP
        ,height=h
        ,volume=vol
        ,porosity=n
        ,k=2.0*nbFrictCont/nb_s
        ,Ek=kineticEnergy()
        ,unbF=unbalancedForce()
    )

phase=0
def triaxDone():
    global phase
    volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])
    if phase==0:
        h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
        vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
        contactStress=getStress(vol)
        vol_s=Rmean=Rmax=nbSph=0
        Rmin=1e6
        for o in O.bodies:
            if isinstance(o.shape,Sphere):
                nbSph+=1
                Rmean+=o.shape.radius
                if o.shape.radius>Rmax: Rmax=o.shape.radius
                if o.shape.radius<Rmin: Rmin=o.shape.radius
                vol_s += 4.*pi/3.*(o.shape.radius)**3
        Rmean=Rmean/nbSph
        n = 1-vol_s/vol
        print('DONE! iter =',O.iter,'| sample generated: nb spheres =',nbSph,', Rmean =',Rmean,', Rratio =',Rmax/Rmin,', porosity =',n)
        print('Changing contact properties now')
        utils.setContactFriction(radians(FRIC))
        print('APPLYING CONFINING PRESSURE: sx,sy and sz will go to PC =',PC)
        triax.goal=(-PC/volRatio,-PC/volRatio,-PC/volRatio) # - - -
        phase+=1
    elif phase==1:
        print('DONE! iter =',O.iter,'| isotropic confinement done: stresses =',volRatio*triax.stress)
        triax.dead=True
        O.pause()

#### Initialization
print('SAMPLE PREPARATION!')

O.run(1000000,1)
O.step()

#### Applying normal stress
print('NORMAL LOADING! iter =',O.iter)

stage=0
stiff=fnPlaten=currentSN=0.
def servo():
    global stage,stiff,fnPlaten,currentSN
    if stage==0:
        currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
        unbF=unbalancedForce()
        boundaryVel=copysign(min(0.1,abs(0.5*(currentSN-SN))),currentSN-SN)
        O.bodies[0].state.vel[1]=boundaryVel
        if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
            stage+=1
            fnPlaten=O.forces.f(0)[1]
            print('Normal stress =',currentSN,' | unbF =',unbF)
            for i in O.interactions.withBody(O.bodies[0].id):
                stiff+=i.phys.kn
            print('Normal stiffness =',stiff)
            print('DONE! iter =',O.iter)
            O.pause()
    if stage==1:
        fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
        boundaryVel=copysign(min(0.1,abs(0.35*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired) # why 0.35?
        O.bodies[0].state.vel[1]=boundaryVel

O.engines = O.engines[:4]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[4:]
O.run(1000000,1)

recData.dead=False
O.run(1,1)

#### injecting fluid
print('FLUID! iter =',O.iter)
flowEng.isActivated=1
flowEng.boundaryUseMaxMin=[1,1,0,0,1,1]
flowEng.bndCondIsPressure=[0,0,1,0,0,0]
newton.damping=0.

O.run(1,1)

# automatic
iter0=O.iter
iterMax=5e3
deltaP=1e5
while O.iter<(iter0+iterMax):
  O.run(1000,True)
  inputP+=deltaP
  flowEng.bndCondValue=[0.,0.,+inputP,0.,0.,0.]
  flowEng.updateBCs()
  print('updateBCs! inputP=',inputP)
  O.run(1,1)

plot.plots={'iter':('pf','stress_upWall1','contactStress11')}

plot.plot()