TriaxialStressController: I cannot guarantee isotropic pressure in all directions.

I use the “TriaxialStressController” engine for compaction. Expand the particle diameter and try to control the pressure at the boundaries.

The settings are as follows.
##################################################
triax = TriaxialStressController(
    maxMultiplier = 1.004,
    finalMaxMultiplier = 1.0001,
    thickness = 0.05,
    stressMask = 7,
    internalCompaction = True,
)
triax.goal1 = triax.goal2 = triax.goal3 = -100e3

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()], verletDist=.05 * .00075),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    TriaxialStateRecorder(iterPeriod=100,file=outPath+"/WallStresses.sand.expand"),
    NewtonIntegrator(damping = damp),
]

while 1:
    O.run(1000, True)
    #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
    unb = unbalancedForce()
    print("unbF:",unb," meanStress:",triax.meanStress," top:",triax.stress(triax.wall_top_id)[1]," porosity:",porosity())
    if unb < stabilityThreshold and abs(1-triax.meanStress/sigmaIso) < stabilityThreshold:
        break
##################################################

However, I found that the pressure in the y-direction is consistently greater than in the x- and z-directions, as follows

unbF: 2.06184e-05 meanStress: -100000.19435 top: -127238.75836 porosity: 0.406501