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:
-
PeriIsoCompressor.sigma
doesn't hold the same main stresses' values as the ones computed withgetStress
. -
PeriIsoCompressor.currUnbalanced
isn't recomputed everyPeriIsoCompressor.globalUpdateInt
iterations, as a consequence its value stays -1. The engine then almost always think that the unbalanced force is lower thanPeriIsoCompressor.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 sincetotalForceInVolume
is no longer called. I am not sure why it is used, but maybe callingtotalForceInVolume
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