> 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

-- 
You received this bug notification because you are a member of Yade
developers, which is subscribed to Yade.
https://bugs.launchpad.net/bugs/1790167

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

Status in Yade:
  New

Bug description:
  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()

To manage notifications about this bug go to:
https://bugs.launchpad.net/yade/+bug/1790167/+subscriptions

_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to     : yade-dev@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-dev
More help   : https://help.launchpad.net/ListHelp

Reply via email to