JCFPM: "neverErase" modifies the simulated behavior while it should not

[Luc Scholtès]

I noticed that running the exact same simulation (with same initial packing) gives different behaviors (stress-strain response in, e.g., a compression test) when neverErase is True or False. Given the purpose of neverErase (keep record of broken contacts, primarily for DFNFlow), it should not. The difference can be more or less important depending on the situation but it always exists. I could not figure out the cause of this yet but it seems that it comes from the treatment of broken contacts (obviously).

Here is a simulation (uniaxial compression) that illustrates the problem. Running the same script using the exact same sample (to make sure the error does not come from a difference in the packings used) with either neverErase=True or neverErase=False produces 2 stress-strain curves which deviate at some point during the simulation. I made sure that the error is only due to neverErase by running the same simulations several times. The curves obtained with neverErase=True are always identical, as the curves obtained with neverErase=False. For those who would be interested, I also attach a packing and the python script to plot the curves (below the simulation script).

### yade script ###

from yade import ymport, pack, plot

#### material definition
def sphereMat(): return JCFpmMat(type=1,density=3000,young=1e9,poisson=0.2,tensileStrength=1e6,cohesion=10e6,frictionAngle=radians(30))

##### create the specimen
#L=0.10
#D=0.05
#pred=pack.inCylinder((0,0,0),(0,0,L),D/2.)
#O.bodies.append(pack.regularHexa(pred,radius=D/20.,gap=0.,material=sphereMat))
#O.bodies.append(pack.randomDensePack(pred,radius=D/20.,rRelFuzz=0.4,spheresInCell=1000,memoizeDb='/tmp/gts-triax-packings.sqlite',returnSpherePack=False,color=(0.9,0.8,0.6),material=sphereMat))

#### import the specimen
O.bodies.append(ymport.text('121_3k.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

#### help define boundary conditions (see utils.uniaxialTestFeatures)
bb=utils.uniaxialTestFeatures()
negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

################# DEM loop + ENGINES DEFINED HERE

O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.2,label='Saabb')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.2,label='SSgeom')],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(neverErase=False,label='interactionLaw')]
 ),
 UniaxialStrainer(strainRate=-0.01,axis=longerAxis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,dead=1,label='strainer'),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.4,label='newton'),
 PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
]

################# RECORDER DEFINED HERE

def recorder():
    yade.plot.addData({'i':O.iter,
                       'eps':strainer.strain,
                       'sigma':strainer.avgStress,
                       'tc':interactionLaw.nbTensCracks,
                       'sc':interactionLaw.nbShearCracks,
                       'te':interactionLaw.totalTensCracksE,
                       'se':interactionLaw.totalShearCracksE,
                       'unbF':utils.unbalancedForce()})
    plot.saveDataTxt('compressionTest_1')

################# PREPROCESSING

#### manage interaction detection factor during the first timestep and then set default interaction range
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

################# SIMULATION REALLY STARTS HERE
strainer.dead=0
O.run(50000)

### python script ###

# -*- coding: utf-8 -*-
from pylab import *

### processing function
def store(var,textFile):
    data=loadtxt(textFile,skiprows=1)
    it=[]
    e=[]
    s=[]
    tc=[]
    sc=[]
    uf=[]
    for i in range(0,len(data)):
      it.append(float(data[i,1]))
      e.append(-float(data[i,0]))
      s.append(-float(data[i,4]))
      tc.append(float(data[i,5]))
      sc.append(float(data[i,2]))
      uf.append(float(data[i,7]))
    var.append(it)
    var.append(e)
    var.append(s)
    var.append(tc)
    var.append(sc)
    var.append(uf)

### data input
dataFile1='compressionTest'
a1=[]
store(a1,dataFile1)

dataFile2='compressionTest_neverErase'
a2=[]
store(a2,dataFile2)

rcParams.update({'legend.numpoints':1,'font.size':20,'axes.labelsize':28,'xtick.major.pad':10,'ytick.major.pad':10,'legend.fontsize':18})

figure(1,figsize=(10,10))
xlabel(r'$\epsilon_1$ [millistrain]')
#axis(xmax=0.1)
plot([x*1e3 for x in a1[1]],[x/1e6 for x in a1[2]],'-k',linewidth=2)
plot([x*1e3 for x in a2[1]],[x/1e6 for x in a2[2]],'-r',linewidth=2)
ylabel(r'$\sigma_1$ [MPa]')
axis(ymin=0)
#savefig(dataFile1+'_qVSeps.eps',dpi=1000,format='eps',transparent=False)

### show
show()

[Chareyre (bruno-chareyre-9)]

«neverErase» is only ok if another functor is taking care of erasing. It
doesn't seem to be the case here, so it should just be «false».
Bruno

Le ven. 31 août. 2018 17:01, Luc Scholtès <email address hidden> a écrit :

> Public bug reported:
>
> I noticed that running the exact same simulation (with same initial
> packing) gives different behaviors (stress-strain response in, e.g., a
> compression test) when neverErase is True or False. Given the purpose of
> neverErase (keep record of broken contacts, primarily for DFNFlow), it
> should not. The difference can be more or less important depending on
> the situation but it always exists. I could not figure out the cause of
> this yet but it seems that it comes from the treatment of broken
> contacts (obviously).
>
> Here is a simulation (uniaxial compression) that illustrates the
> problem. Running the same script using the exact same sample (to make
> sure the error does not come from a difference in the packings used)
> with either neverErase=True or neverErase=False produces 2 stress-strain
> curves which deviate at some point during the simulation. I made sure
> that the error is only due to neverErase by running the same simulations
> several times. The curves obtained with neverErase=True are always
> identical, as the curves obtained with neverErase=False. For those who
> would be interested, I also attach a packing and the python script to
> plot the curves (below the simulation script).
>
>
> ### yade script ###
>
>
> from yade import ymport, pack, plot
>
> #### material definition
> def sphereMat(): return
> JCFpmMat(type=1,density=3000,young=1e9,poisson=0.2,tensileStrength=1e6,cohesion=10e6,frictionAngle=radians(30))
>
> ##### create the specimen
> #L=0.10
> #D=0.05
> #pred=pack.inCylinder((0,0,0),(0,0,L),D/2.)
> #O.bodies.append(pack.regularHexa(pred,radius=D/20.,gap=0.,material=sphereMat))
>
>
> #O.bodies.append(pack.randomDensePack(pred,radius=D/20.,rRelFuzz=0.4,spheresInCell=1000,memoizeDb='/tmp/gts-triax-packings.sqlite',returnSpherePack=False,color=(0.9,0.8,0.6),material=sphereMat))
>
> #### import the specimen
>
> O.bodies.append(ymport.text('121_3k.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))
>
> #### help define boundary conditions (see utils.uniaxialTestFeatures)
> bb=utils.uniaxialTestFeatures()
>
> negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
>
> ################# DEM loop + ENGINES DEFINED HERE
>
> O.engines=[
> ForceResetter(),
>
> InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.2,label='Saabb')]),
> InteractionLoop(
>
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.2,label='SSgeom')],
>
> [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
>
> [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(neverErase=False,label='interactionLaw')]
> ),
>
> UniaxialStrainer(strainRate=-0.01,axis=longerAxis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,dead=1,label='strainer'),
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0...

[Luc Scholtès (luc)]

OK... We need neverErase=True for DFNFlow so what should we do? Create a functor to erase the interactions?

[Robert Caulk (rcaulk)]

> I noticed that running the exact same simulation (with same initial packing) gives different behaviors

I ran your script multiple times but I am having trouble identifying the different stress-strain behaviors you mention. What exactly is different? Is it the post failure part of the curve (I guess it has to be if neverErase is responsible)? Are my post-failure curves more consistent than yours [see attached plots]?

> The curves obtained with neverErase=True are always identical

This is not the case when I run your script (or am I missing something?) [see attached figures].

[Robert Caulk (rcaulk)]

I am also curious as to why we see changes in these curves between replicates (even for neverErase=False), despite our use of a constant sphere packing...

[Bruno Chareyre (bruno-chareyre)]

Random results (we are speaking of "yade -j1" right?) usually come from using null pointers or uninitialized variables.
Running yade with valgrind (and debug symbols) would probably tell more.
In the case neverErase=true, I would double-check what happens if the interaction becomes active again.
Bruno

[Luc Scholtès (luc)]

Hi guys,

As illustrated in attached figures, I always get the same curve when running the simulation with the same packing (with either neverErase=True or neverErase=False). I don't face the problem you mentioned Robert...

I do have different curves for neverErase=True and neverErase=False though (initial point of the bug report). The difference is "minimal" in this case (post peak) but I have other results where it is pretty important (even pre peak).

Just to explain a bit further:

I try to compute the change of permeability due to stress induced damage. To do so, I run a compression and compute the permeability of the sample at different stages of the loading (the simulation is not coupled, I just call for DFNFlow to compute the permeability). I thus need to keep track of the induced cracks (hence the neverErase=True) so as to take into account the evolution of their aperture during the loading.

Luc

[Robert Caulk (rcaulk)]

>>
I always get the same curve when running the simulation with the same packing (with either neverErase=True or neverErase=False). I don't face the problem you mentioned Robert...
>>

Ok, I re-ran your scripts with -j1 (instead of -j6), and sure enough, I am able to replicate results (neverErase=True and False). Looking back at the comparison of two replicate runs with -j6, the difference is almost as staggering as comparing neverErase=True vs neverErase=False...how could floating point differences add up to almost 3% difference for peak strength?! Also, my curves do not match yours exactly (they are close though), why? I guess this is a not the bug we are focused on here...

W.r.t the bug in this thread: is it possible that the interaction between two particles gets "neverErased" (i.e. shearForce = normalForce = FnMax=FsMax = 0), but then the particles later regain contact and "want" to interact frictionally? In this case they but cannot interact frictionally since we have assigned all the interaction values to 0, right? In contrast, the neverErase=False case would allow two particles to lose cohesive contact and later regain frictional contact since the collider handles interaction disposal and creation, right?

>>
The difference is "minimal" in this case (post peak) but I have other results where it is pretty important (even pre peak).
>>

I would say any difference is significant considering we are able to get perfect replicates. Also, you mention pre peak differences: have you ever noticed any differences BEFORE a single bond has failed?

Robert

[Bruno Chareyre (bruno-chareyre)]

> OK... We need neverErase=True for DFNFlow so what should we do? Create a functor to erase the interactions?

No there is no need for an extra functor. The place where erasing is triggered is Law2::action() (based on the bool returned).

I don't think we really need neverErase=True for DFNFlow. What we really need is to keep interactions alive until a condition is met. What about erasing (i.e. returning false) beyond a certain distance?
This way it would be possible to keep interactions for a sufficient period without interfering with neverErase (which really means to defer deletion to another functor).
This being said, using neverErase=True is probably a solution, in the short term.

Again, I guess some contacts forces could go wrong when two particles detach-touch-detach-..., if not-erased virtual interactions don't receive correct updates. A two sphere example could help check this.

Bruno

[Luc Scholtès (luc)]

"Also, my curves do not match yours exactly (they are close though), why?"

-> Probably because I ran these with the option to delete the contact as soon as it fails in shear (JCFPM.cpp line 249: option 2: delete contact if in tension). That's not the issue here.

"is it possible that the interaction between two particles gets "neverErased" (i.e. shearForce = normalForce = FnMax=FsMax = 0), but then the particles later regain contact and "want" to interact frictionally?"

-> As suggested Bruno, I tend to think so... My first guess was that it is due to the interaction range: when the same particles come into contact after failure, interaction range is different for erased contact than for "neverErased" contact (the neverErased contacts keeps record of the initial equilibrium distance which corresponds to the centre to centre distance at t=0 while the erased one initializes it to the default value which is the sum of the 2 radii). I tried to run the simulation with the default interaction range (bonds created for strict geometric contact) and the discrepancy tends to reduce (but it still exists). I am working on the 2 sphere examples but I am also wondering:

Can we have 2 interactions (in the interaction container") which concerns the same particles pair? For instance could we have, at the same time, a "non erased / broken /frictional" and a "new / frictional" interaction between the 2 same particles?

"What we really need is to keep interactions alive until a condition is met. What about erasing (i.e. returning false) beyond a certain distance?"

-> we need to keep record of the crack aperture for the computation of the local permeability. Ideally, we need this during the whole duration of the simulation.

[Bruno Chareyre (bruno-chareyre)]

>Can we have 2 interactions (in the interaction container") which concerns the same particles pair?

Impossible.

> keep record of the crack aperture [...] during the whole duration of the simulation.

Saying that we erase when distance > 10*radius would surely enough keep it alive for the whole (typical) DFNFlow simulation.
On the other hand if one day some particles are moving apart by more than 10*radius (possibly without DFNFlow) because there is sliding along a fracture plane it makes no sense to keep the corresponding interaction. The "10" factor should be a user defined variable, of course.

Bruno

[Janek Kozicki (cosurgi)]

moved from https://bugs.launchpad.net/yade/+bug/1790167 to https://gitlab.com/yade-dev/trunk/issues/43