i.geom.penetrationDepth = nan after load yade.gz file

Hi Jan,
Many thanks for your reply.

With your suggestion, I got the i.geom.penetrationDepth=someValue. It is really good compared to Nan before.

However, I found another two problems that I'd like to discuss:
(1) In the Script2(i.e., loading sample). No matter I set O.dt=0 or not, the i.geom.penetrationDepth (for a certain interaction) is the same. You may see it by skipping O.dt=0 in Script2 below.

(2) At the end of Script1 (i.e., making sample), I record a penetrationDepth of a certain interaction, for example O.interactions[2154,3930].geom.penetrationDepth=1.54203739497e-06.
And when I go to Script2, I check the penetrationDepth of that interaction , the value I got in Script2 is slightly different with Script1 (i.e., O.interactions[2154,3930].geom.penetrationDepth=1.5419131404164482e-06). which means that the state has changed although it is slight.

You may reproduce the problem using the following MWE (Script1 and Script2) below if you are interested. (if O.interactions[2154,3930] isn't existed, You may need to select a new one from script1):

##################Script 1#######################
from __future__ import division
from yade import pack, plot
num_spheres=7000#
targetPorosity = 0.45
compFricDegree = 30
finalFricDegree = 18
rate=-0.01
damp=0.6
stabilityThreshold=0.001
young=2e8
confinement=100e3

mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)

MatWall=O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
MatSand = O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3,num_spheres,False, 0.95,seed=1)

O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
triax=TriaxialStressController(
    maxMultiplier=1.+3e7/young,
    finalMaxMultiplier=1.+16e4/young,
    thickness = 0,
    stressMask = 7,
    internalCompaction=True,
)
newton=NewtonIntegrator(damping=damp)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    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='WallStresses'),
    newton,
]

Gl1_Sphere.stripes=0
triax.goal1=triax.goal2=triax.goal3=-confinement

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]
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001:
      break

print "### Isotropic state saved ###"

import sys
while triax.porosity>targetPorosity:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)

O.save('sample.yade.gz')
print "### Compacted state saved ###"

for i in O.interactions:
    print i.id1,i.id2,i.geom.penetrationDepth

##################Script 2#######################
from __future__ import division
from yade import pack, plot

O.load('sample.yade.gz')
confinement=100e3
triax=TriaxialStressController(
    thickness = 0,
    stressMask = 5,
    internalCompaction=False,
)
newton=NewtonIntegrator(damping=0.4)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),
         Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    TriaxialStateRecorder(iterPeriod=1000,file='WallStresses'),
    newton,
]
triax.goal2=-1
triax.goal1=-confinement
triax.goal3=-confinement

print "penetrationDepth:", O.interactions[2154,3930].geom.penetrationDepth
print "O.dt:", O.dt
# O.dt=0 #(or skip this step for my question 1)
O.step()
print "penetrationDepth:", O.interactions[2154,3930].geom.penetrationDepth
#############################
Do you have any ideas about this?

Thanks
Leonard

These are limitations of the (explicit? not sure) discrete model, where the values in one time step may differ in "actual" time

1) have a look at the O.interactionLoop. Ig2 (where penetrationDepth is computed) is before NewtonIntegrator (where the position are updated).
After one simulation step (one execution of O.interactionLoop engines) you have "old" penetrationDepth even the particle position are "new".

2) is very similar.
After script1, you already have "new" positons, from which penetrationDepth in script2 is computed, but if you just print penetrationDepth in script1, you are getting "old" values.

You would have to set O.dt=0 also for the last step in script1 to get the same values in both script1 and script2

cheers
Jan

Hi Jan,

> After one simulation step (one execution of O.interactionLoop engines)
you have "old" penetrationDepth even the particle position are "new".

From this insightful comment is it ok to conclude that re-running after
reload will give same result as with no save/load (I guess so), since
interaction loop will update wrt new positions before computing forces?
If so the issue is nearly the same as "if I put penetrationDepth recorder
before/after interactionLoop I get different values at iteration N".
Cheers
B

@Bruno:

more like:
re-running after reload is equivalent to continuing with simulation (further step(s)) in the original simulation.
I.e. the values are intrinsically different that the saved state

(I did not really get your comment :-)

cheers
Jan

Thanks Jan and Bruno for your comments,

I tried a few more times in different ways (e.g., set O.dt=0 or not in both scripts) and still can't find the "exact state (i.e., the same penetrationDepth)'' between script1 and script2. What Jan explained make sense.

It is good that I can use O.step() to re-gain the penetrationDepth after loading the sample. I will check further and hope this one additional step will not affect the final results compared to the case which doesn't involve this O.step.

Thanks
Leonard

Thanks Jan Stránský, that solved my question.

> (I did not really get your comment :-)

I was pointing the fact that the notion of "value of X at iteration N" is
vague in the first place, since X can change at any point within one
time-loop.
Virtually you could insert multiple O.save() or any form of snapshot within
an iteration, interleaved with the usual engines, and you would get
different numbers in each of them depending on their relative positions in
the sequence of engines.
The above save/load issue is nearly the same thing I would say.
B

@Leonard
The rationale for not saving/loading Un is that it is always unambiguously
deduced from positions and sizes.
We try to store only independent variables, not derived variables.
What you report is not actually a problem, or it is only a problem of the
way you output penetrationDepth.
Insert a PyRunner to print Un, after NewtonIntegrator. The apparent issue
will disappear.
B

On Thu, 20 Aug 2020 at 17:48, Bruno Chareyre <bruno.chareyre@3sr-grenoble.fr>
wrote:

> > (I did not really get your comment :-)
>
> I was pointing the fact that the notion of "value of X at iteration N" is
> vague in the first place, since X can change at any point within one
> time-loop.
> Virtually you could insert multiple O.save() or any form of snapshot
> within an iteration, interleaved with the usual engines, and you would get
> different numbers in each of them depending on their relative positions in
> the sequence of engines.
> The above save/load issue is nearly the same thing I would say.
> B
>

Great, thanks Bruno!

Cheers
Leonard