Bug in Calculating getSaturation in TwoPhaseFlow Engine

hi,
I had a problem using "TwoPhaseFlow Engine" in 'Yadedaily 2018.02b-290bf6a54e~xenial' and 'Yade-2019.01a'.
The value of Saturation Degree 'getSaturation(False)' was shown properly while doing drainage,but after the dry completed, I close the pressure related to water in bottom of box by 'bndCondIsPressure=[0,0,0,1,0,0]' . I provide Oedometer conditions for the loading.

According to the formula:
Porosity=VolumeOfVoids/BoxVolume
VoidVolumes=Porosity*BoxVolume

The volume of Voids containing the air and Water.
We have blocked the way out water by 'bndCondIsPressure=[0,0,0,1,0,0]' While Loading, So Volume of Water while loading is Constant.
And only the volume of air in the ٰVoid decreases by increasing loading.

According to the formula:
DegreeOfSaturation=VolumeOfWater/VolumeOfVoids

while Oedometer loading:
VolumeOfWater > Constant
VolumeOfVoids > Decresed

So DegreeOfSaturation must be Increased that The same is true with the laboratory results and the literature
but in TwoPhaseFlow Engine Degree of Saturation No changed while loading!!!
-----------------------------------
CODE:

import matplotlib; matplotlib.rc('axes',grid=True)
from yade import pack
import pylab
from numpy import *

compFricDegree = 30.0
confiningS=-1e5
psdSizes,psdCumm=[.599,0.6,0.849,0.85,1.0,1.40],[0.,0.35,0.35,0.70,.70,1.]
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(10,10,10)
sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=1000,seed=1)

O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=0,density=0,label='frictionless'))

walls=aabbWalls((mn,mx),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(
 internalCompaction=True,
 goal1=confiningS,
 goal2=confiningS,
 goal3=confiningS,
 max_vel=10,
 label="triax"
)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 NewtonIntegrator(damping=0.4,label="newton"),
 TwoPhaseFlowEngine(label="unsat")
]

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.01:
    break

##################################
## Drainage Test under oedometer conditions #
##################################
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=0

meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]
unsat.isPhaseTrapped=True
unsat.initialization()
unsat.surfaceTension = 10

for pg in arange(1.0e-5,5.0,0.1):
  unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]
  unsat.invasion()
  unsat.computeCapillaryForce()
  for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))
  while 1:
    O.run(1000,True)
    unb=unbalancedForce()
    if unb<0.01:
      break
  print 'Drainage:::','PC:',-unsat.bndCondValue[2],' Sw:',unsat.getSaturation(False),' VoidVolume:',triax.porosity*triax.boxVolume

print '----------LOADING--------------'
## Oedometer conditions
unsat.bndCondIsPressure=[0,0,0,1,0,0]
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=False
loadingMatrix=[-500e3,-1000e3,-1500e3,-2000e3,-2500e3,-3000e3]
for i in arange(0,len(loadingMatrix),1):
 triax.goal2=loadingMatrix[i]
 O.run(2000,True)
 print 'Load:',loadingMatrix[i],' Sw:',unsat.getSaturation(False),' VoidVolume:',triax.porosity*triax.boxVolume
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-------------------------------------------------------------------------
end of code
-------------------------------------------------------------------------

The results I see is wrong:

Load: -500000.0 Sw: 0.268726543553 VoidVolume: 436.590507582
Load: -1000000.0 Sw: 0.268726543553 VoidVolume: 433.836864791
Load: -1500000.0 Sw: 0.268726543553 VoidVolume: 431.014178921
Load: -2000000.0 Sw: 0.268726543553 VoidVolume: 428.158278656
Load: -2500000.0 Sw: 0.268726543553 VoidVolume: 425.283191863
Load: -3000000.0 Sw: 0.268726543553 VoidVolume: 422.394919453
---------------------------------------------------------

It is observed that due to the closure of the water outlet, the volume of Voids decreased but the amount of saturation remained constant during loading, which is incorrect.
It seems that when loading and changing the sample size, the effect of the saturation change is not calculated.