Errors with PeriIsoCompressor's stress computation and unbalancedForce update
Following this question in the launchpad, I open up this issue concerning PeriIsoCompressor.
It looks like there is two errors in PeriIsoCompressor:
PeriIsoCompressor.sigmadoesn't hold the same main stresses' values as the ones computed with
PeriIsoCompressor.currUnbalancedisn't recomputed every
PeriIsoCompressor.globalUpdateIntiterations, 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
avgStiffnessis no longer computed since
totalForceInVolumeis no longer called. I am not sure why it is used, but maybe calling
totalForceInVolumewithout 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