Errors with PeriIsoCompressor's stress computation and unbalancedForce update

Following this question in the launchpad, I open up this issue concerning PeriIsoCompressor.

The problem

It looks like there is two errors in PeriIsoCompressor:

  1. PeriIsoCompressor.sigma doesn't hold the same main stresses' values as the ones computed with getStress.
  2. PeriIsoCompressor.currUnbalanced isn't recomputed every PeriIsoCompressor.globalUpdateInt iterations, as a consequence its value stays -1. The engine then almost always think that the unbalanced force is lower than PeriIsoCompressor.maxUnbalanced, while it isn't.

To illustrate these problems, here is a MWE:

O.periodic = True
O.cell.setBox((1, 1, 1))

sp=pack.SpherePack()
sp.makeCloud(minCorner=(0,0,0), maxCorner=(1, 1, 1), rMean=5e-3, rRelFuzz=5e-4, num=5000, periodic=True, seed=42)
sp.toSimulation()

iso_comp = PeriIsoCompressor(charLen=5e-3, stresses=[-1e6], doneHook="O.pause()")
O.engines = O.engines + [iso_comp]

O.run(wait=True)

print("\nPIC mean stress: ", iso_comp.sigma.mean())
print("getStress mean stress: ", getStress().trace()/3)

print("\nPIC max unbalancedForce: ", iso_comp.maxUnbalanced)
print("PIC unbalancedForce: ", iso_comp.currUnbalanced)
print("Actual unbalancedForce: ", unbalancedForce())

And its output when executed with Yade 2021-06-08.git-4852d560:

Welcome to Yade 2021-06-08.git-4852d56 
Using python version: 3.6.9 (default, Jan 26 2021, 15:33:00) 
[GCC 8.4.0]
TCP python prompt on localhost:9001, auth cookie `csudye'
XMLRPC info provider on http://localhost:21001
Running script mwe.py

PIC mean stress:  -995101.1166703544
getStress mean stress:  -37873.57983158174

PIC max unbalancedForce:  0.0001
PIC unbalancedForce:  -1.0
Actual unbalancedForce:  0.01334929567748287

An incomplete fix

I tried to make a few changes in the code to fix this (see correct_PIC branch), but it is incomplete:

  • the variable avgStiffness is no longer computed since totalForceInVolume is no longer called. I am not sure why it is used, but maybe calling totalForceInVolume without using its output would be enough to solve the problem. Or another function (computeAvgStiffness ?) could be created just for that, maybe that would be more clean.
  • the velocity of the boundaries is decreased using an hard-coded coefficient on velGrad (I did that so my MWE would work, otherwise the sample doesn't stabilize). I think that adding a parameter to tune the compaction's velocity would be nice.

Also, I have put the if (allStressesOK) block inside a if ((step % globalUpdateInt) == 0) (the stress is recomputed only every globalUpdateInt steps anyway).

Here is the output of the MWE after these changes:

Welcome to Yade 2021-06-22.git-c37c18e 
Using python version: 3.6.9 (default, Jan 26 2021, 15:33:00) 
[GCC 8.4.0]
TCP python prompt on localhost:9000, auth cookie `yesuka'
XMLRPC info provider on http://localhost:21000
Running script mwe.py

PIC mean stress:  -998941.7264993662
getStress mean stress:  -998941.7264993662

PIC max unbalancedForce:  0.0001
PIC unbalancedForce:  8.793531162646358e-05
Actual unbalancedForce:  8.793531162646358e-05