>> 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 -- 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