PeriTriaxController not initiating next stage

Hello,

I'm trying to do a uniaxial compression of a periodic box with a 3x3x3 grid of particles.

I'm implementing periodic boundary conditions for level-set particles. I had some unexpected results (level set spheres gave a much lower density). Now I want to do a unit test by comparing a very simple simulation using analytical spheres and level-set spheres.

Whenever I run the simulation (see code below), the pressure increases to the target pressure properly. However, hereafter the simulation keeps iterating but does not initiate the decompression stage defined by compactionFinished(). I have attempted to change maxUnbalanced and relStressTol to see if the issue is there but had no success.

Would any of you happen to know why the next stage is not initiated by PeriTriaxController?

With kind regards,
Danny

###########################################################
# encoding: utf-8

from __future__ import print_function
sigmaIso = -1e6
ydim = 1.0
xdim = 1.0
zdim = 1.0
r = 1.0/3.0/2.0

#import matplotlib
#matplotlib.use('Agg')

# generate loose packing
from yade import pack, qt, plot

O.periodic = True

O.bodies.append(
        [
                # In xy-plane
                sphere(center=(r, r, r), radius=r, fixed=True),
                sphere(center=(3*r, r, r), radius=r, fixed=True),
                sphere(center=(5*r, r, r), radius=r, fixed=True),
                sphere(center=(r, 3*r, r), radius=r, fixed=True),
                sphere(center=(3*r, 3*r, r), radius=r, fixed=True),
                sphere(center=(5*r, 3*r, r), radius=r, fixed=True),
                sphere(center=(r, 5*r, r), radius=r, fixed=True),
                sphere(center=(3*r, 5*r, r), radius=r, fixed=True),
                sphere(center=(5*r, 5*r, r), radius=r, fixed=True),
                # Repeat in z-direction
                sphere(center=(r, r, 3*r), radius=r, fixed=True),
                sphere(center=(3*r, r, 3*r), radius=r, fixed=True),
                sphere(center=(5*r, r, 3*r), radius=r, fixed=True),
                sphere(center=(r, 3*r, 3*r), radius=r, fixed=True),
                sphere(center=(3*r, 3*r, 3*r), radius=r, fixed=True),
                sphere(center=(5*r, 3*r, 3*r), radius=r, fixed=True),
                sphere(center=(r, 5*r, 3*r), radius=r, fixed=True),
                sphere(center=(3*r, 5*r, 3*r), radius=r, fixed=True),
                sphere(center=(5*r, 5*r, 3*r), radius=r, fixed=True),
                # Repeat in z-direction
                sphere(center=(r, r, 5*r), radius=r, fixed=True),
                sphere(center=(3*r, r, 5*r), radius=r, fixed=True),
                sphere(center=(5*r, r, 5*r), radius=r, fixed=True),
                sphere(center=(r, 3*r, 5*r), radius=r, fixed=True),
                sphere(center=(3*r, 3*r, 5*r), radius=r, fixed=True),
                sphere(center=(5*r, 3*r, 5*r), radius=r, fixed=True),
                sphere(center=(r, 5*r, 5*r), radius=r, fixed=True),
                sphere(center=(3*r, 5*r, 5*r), radius=r, fixed=True),
                sphere(center=(5*r, 5*r, 5*r), radius=r, fixed=True)
        ]
)

O.cell.hSize = Matrix3(xdim, 0.0, 0.0, 0.0, ydim, 0.0, 0.0, 0.0, zdim)

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
        PeriTriaxController(
                label='triax',
                # specify target values and whether they are strains or stresses
                goal=(0, 0, sigmaIso),
                stressMask=4,
                # type of servo-control
                dynCell=True,
                maxStrainRate=(0, 0, 0.005),
                # wait until the unbalanced force goes below this value
                maxUnbalanced=.05,
                relStressTol=1e-4,
                # call this function when goal is reached and the packing is stable
                doneHook='compactionFinished()'
        ),
        NewtonIntegrator(damping=.2),
        PyRunner(command='addPlotData()', iterPeriod=100),
]
O.dt = .25 * PWaveTimeStep()

def compactionFinished():
 # set the current cell configuration to be the reference one
 O.cell.trsf = Matrix3.Identity
 # change control type: keep constant confinement in x,y, 20% compression in z
 triax.goal = (0, 0, .2)
 triax.stressMask = 3
 # allow faster deformation along x,y to better maintain stresses
 triax.maxStrainRate = (0, 0, 0.005)
 # next time, call triaxFinished instead of compactionFinished
 triax.doneHook = 'triaxFinished()'
 # do not wait for stabilization before calling triaxFinished
 triax.maxUnbalanced = 10

def triaxFinished():
 print('Finished')
        plot.saveDataTxt('results_ptt3.txt',vars=('compDens','sz'))
 O.pause()

totalMass = 0
for b in O.bodies:
        totalMass+=b.state.mass

def addPlotData():
 plot.addData(
         unbalanced=unbalancedForce(),
         i=O.iter,
         sz=abs(triax.stress[2]),
                compDens=totalMass/O.cell.volume
 )

# ez=abs(triax.strain[2]),

print("The total mass is", totalMass)
print("The number of particles is",len(O.bodies))
#print("The volume is", O.cell.volume)

# define what to plot
plot.plots = {
        'i': ('unbalanced'),
        'compDens': ('sz'),
}
# show the plot
plot.plot()

O.run()