Can I control the seed in randomdensepack?

Asked by kalogeropoulos on 2018-11-23

Hello to everyone,

I am trying to make a simulation of cutting of a synthetic sample (i will post the code above).
I am creating the sample with randomdensepack, but every time i am making the simulation , i take different diagramm of the force. This is due to the seed, that is in randomdensepack (which i can control).

My question is: If i want to put a constant value of the seed number, how i do it?

I haven't build yade from trunk, i install it via < sudo ....>
My ubuntu distribution is Bionic 18.04.

Thanks a lot.

The code is:

!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division

from yade import plot,pack,timing, utils, geom
import time, sys, os, copy

# default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
    density = 2674.3,
    frictionAngle= 44, #frictionAngle
    tensileStrength = 65e6, #Tensile Strength
 poisson=.20,
 cohesion = 65e6, #Shear Strength
 young = 60e9,

 intRadius=1.00, # to set initial contacts to larger neighbours
 dtSafety=.8,
 damping=0.6,
 # Characteristic of sample:
 x = 75e-3, #75e-3
 y = 27e-3,
 z= 6.35e-3, # edw to exoume balei na einai idio me to platos tou koptikou
 # Characteristic of shperes
 sphereRadius=1e-3, #0.35e-3
 rRelFuzz=.34,
 #characteristic of cutter
 x_cutter = 5e-3, # mhkos kata ton x-axona
 y_cutter = 1e-3, # ipsos koptikou (y-axona)
 z_cutter = 6.35e-3, # platos koptikou
 DOC = 5e-3, #Bathos kopis
)

from yade.params.table import *

if 'description' in O.tags.keys(): O.tags['id']=O.tags['id']+O.tags['description']

#material properties
sample = O.materials.append(JCFpmMat(young=young, poisson=poisson, frictionAngle=radians(frictionAngle), cohesion=cohesion, tensileStrength=tensileStrength,density = density, label='spheres'))

sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((-x,-y,-z),(x,y,z)),radius=sphereRadius,rRelFuzz=rRelFuzz,spheresInCell=1000,returnSpherePack=True) #spheresInCell=2000
sp.toSimulation(color=(1.09,1.0,0.0),wire=False,material=sample)

#kataskeuh koptikou
bx = geom.facetBox(center=(-80e-3,27e-3+y_cutter,0.0),extents=(x_cutter,-(5e-3+DOC),z_cutter),orientation=Quaternion((0,0,1),-pi/36), wallMask=(2+8))
O.bodies.append(bx)

for facet in bx:
 facet.state.blockedDOFs='xyzXYZ' #block e.w na metakineitai mono kata ton x-axona
 facet.state.vel=(5.0,0,0) # taxuthta kata ton x-axona

O.dt=dtSafety*PWaveTimeStep()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),Bo1_Facet_Aabb()],verletDist=.05*sphereRadius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc'), Ig2_Facet_Sphere_ScGeom()],
        [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)], # Ip2_FrictMat_FrictMat_FrictPhys(),
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(Key='cracks',recordCracks=True,cracksFileExist=True,label='interactionLaw')],
    ),
    GlobalStiffnessTimeStepper(),
    VTKRecorder(fileName='3d-vtk-',recorders=['spheres','intr','all','cracks','jcfpm','facets'],Key='cracks',iterPeriod=150),
    NewtonIntegrator(damping=damping,label='damper'),
    PyRunner(command='addPlotData()',realPeriod=2),
    PyRunner(command='StopSimulation()',iterPeriod=50),
]

def addPlotData():
   global Fx
   global dx
   Fx = 0.0
   Fx = abs(sum(O.forces.f(facet.id)[0] for facet in bx))
   plot.addData(i = O.iter , Fx=Fx, dx = bx[0].state.pos[0], tc=interactionLaw.nbTensCracks, sc=interactionLaw.nbShearCracks)
   plot.saveDataTxt(O.tags['d.id']+'.txt')

def StopSimulation():
    if bx[0].state.pos[0] >= 0.0:
        O.pause()

plot.plots={'dx':('Fx'), ' dx':('tc', 'sc')}
plot.plot()

print " Number of particles: %s" %(len(O.bodies))
print " Number of Interactions per Particle: %s" %(utils.avgNumInteractions())
print " Porosity of sample is: %s" %(utils.porosity())
print "\n"
O.step()
is2aabb.aabbEnlargeFactor = 1.00
ss2sc.interactionDetectionFactor = 1.00
print " Number of particles: %s" %(len(O.bodies))
print " Number of Interactions per Particle: %s" %(utils.avgNumInteractions())
print " Porosity of sample is: %s" %(utils.porosity())

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2018-11-23
Last query:
2018-11-23
Last reply:
2018-11-23
Best Jan Stránský (honzik) said : #1

Hello,
randomDensePack currently does not support seed.
options:
1) rewrite (copy-paste) the function with seed parameter
2) run randomDensePack once, save the results* and the next time use the saved results

* using memoizeDB parameter of randomDensePack [1] or using save and load functions of SpherePack [2,3] or ... (plenty of ways)

cheers
Jan

[1] https://yade-dem.org/doc/yade.pack.html#yade.pack.randomDensePack
[2] https://yade-dem.org/doc/yade.pack.html#yade._packSpheres.SpherePack.save
[3] https://yade-dem.org/doc/yade.pack.html#yade._packSpheres.SpherePack.load

kalogeropoulos (antoniskal) said : #2

Thank you a lot Jan Stránský (honzik).

Well i need to make a parametric analysis of how random generation of spherical elements influence the force during the rock cutting simulatio, so i must rewrite the function.

i make that but it doesnt run.

The code with the seed:

def randomDensePack(predicate,radius,material=-1,dim=None,cropLayers=0,rRelFuzz=0.,spheresInCell=0,memoizeDb=None,useOBB=False,memoDbg=False,color=None,seed=0,returnSpherePack=None):
   """Generator of random dense packing with given geometry properties, using TriaxialTest (aperiodic)
   or PeriIsoCompressor (periodic). The periodicity depens on whether the spheresInCell parameter is given.

   *O.switchScene()* magic is used to have clean simulation for TriaxialTest without deleting the original simulation.
   This function therefore should never run in parallel with some code accessing your simulation.

   :param predicate: solid-defining predicate for which we generate packing
   :param spheresInCell: if given, the packing will be periodic, with given number of spheres in the periodic cell.
   :param radius: mean radius of spheres
   :param rRelFuzz: relative fuzz of the radius -- e.g. radius=10, rRelFuzz=.2, then spheres will have radii 10 ± (10*.2)), with an uniform distribution.
      0 by default, meaning all spheres will have exactly the same radius.
   :param cropLayers: (aperiodic only) how many layers of spheres will be added to the computed dimension of the box so that there no
      (or not so much, at least) boundary effects at the boundaries of the predicate.
   :param dim: dimension of the packing, to override dimensions of the predicate (if it is infinite, for instance)
   :param memoizeDb: name of sqlite database (existent or nonexistent) to find an already generated packing or to store
      the packing that will be generated, if not found (the technique of caching results of expensive computations
      is known as memoization). Fuzzy matching is used to select suitable candidate -- packing will be scaled, rRelFuzz
      and dimensions compared. Packing that are too small are dictarded. From the remaining candidate, the one with the
      least number spheres will be loaded and returned.
   :param useOBB: effective only if a inGtsSurface predicate is given. If true (not default), oriented bounding box will be
      computed first; it can reduce substantially number of spheres for the triaxial compression (like 10× depending on
      how much asymmetric the body is), see examples/gts-horse/gts-random-pack-obb.py
   :param memoDbg: show packings that are considered and reasons why they are rejected/accepted
   :param returnSpherePack: see the corresponding argument in :yref:`yade.pack.filterSpherePack`

   :return: SpherePack object with spheres, filtered by the predicate.
   """
   import sqlite3, os.path, cPickle, time, sys, _packPredicates, numpy
   from math import pi
   wantPeri=(spheresInCell>0)
   if 'inGtsSurface' in dir(_packPredicates) and type(predicate)==inGtsSurface and useOBB:
      center,dim,orientation=gtsSurfaceBestFitOBB(predicate.surf)
      print "Best-fit oriented-bounding-box computed for GTS surface, orientation is",orientation
      dim*=2 # gtsSurfaceBestFitOBB returns halfSize
   else:
      if not dim: dim=predicate.dim()
      if max(dim)==float('inf'): raise RuntimeError("Infinite predicate and no dimension of packing requested.")
      center=predicate.center()
      orientation=None
   if not wantPeri: fullDim=tuple([dim[i]+4*cropLayers*radius for i in 0,1,2])
   else:
      # compute cell dimensions now, as they will be compared to ones stored in the db
      # they have to be adjusted to not make the cell to small WRT particle radius
      fullDim=dim
      cloudPorosity=0.25 # assume this number for the initial cloud (can be underestimated)
      beta,gamma=fullDim[1]/fullDim[0],fullDim[2]/fullDim[0] # ratios β=y₀/x₀, γ=z₀/x₀
      N100=spheresInCell/cloudPorosity # number of spheres for cell being filled by spheres without porosity
      x1=radius*(1/(beta*gamma)*N100*(4/3.)*pi)**(1/3.)
      y1,z1=beta*x1,gamma*x1; vol0=x1*y1*z1
      maxR=radius*(1+rRelFuzz)
      x1=max(x1,8*maxR); y1=max(y1,8*maxR); z1=max(z1,8*maxR); vol1=x1*y1*z1
      N100*=vol1/vol0 # volume might have been increased, increase number of spheres to keep porosity the same
      sp=_getMemoizedPacking(memoizeDb,radius,rRelFuzz,x1,y1,z1,fullDim,wantPeri,fillPeriodic=True,spheresInCell=spheresInCell,memoDbg=False)
      if sp:
         if orientation:
            sp.cellSize=(0,0,0) # resetting cellSize avoids warning when rotating
            sp.rotate(*orientation.toAxisAngle())
         return filterSpherePack(predicate,sp,material=material,returnSpherePack=returnSpherePack)
      else: print "No suitable packing in database found, running",'PERIODIC compression' if wantPeri else 'triaxial'
      sys.stdout.flush()
   O.switchScene(); O.resetThisScene() ### !!
   if wantPeri:
      # x1,y1,z1 already computed above
      sp=SpherePack()
      O.periodic=True
      #O.cell.refSize=(x1,y1,z1)
      O.cell.setBox((x1,y1,z1))
      #print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.cell.refSize
      #print x1,y1,z1,radius,rRelFuzz
      O.materials.append(FrictMat(young=3e10,density=2400))
      num=sp.makeCloud(Vector3().Zero,O.cell.refSize,radius,rRelFuzz,spheresInCell,True,seed=seed)
      O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=.05*radius),InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5,keepProportions=True),NewtonIntegrator(damping=.6)]
      O.materials.append(FrictMat(young=30e9,frictionAngle=.5,poisson=.3,density=1e3))
      for s in sp: O.bodies.append(utils.sphere(s[0],s[1]))
      O.dt=utils.PWaveTimeStep()
      O.run(); O.wait()
      sp=SpherePack(); sp.fromSimulation()
      #print 'Resulting cellSize',sp.cellSize,'proportions',sp.cellSize[1]/sp.cellSize[0],sp.cellSize[2]/sp.cellSize[0]
      # repetition to the required cell size will be done below, after memoizing the result
   else:
      assumedFinalDensity=0.6
      V=(4.0/3.0)*pi*radius**3.0; N=assumedFinalDensity*fullDim[0]*fullDim[1]*fullDim[2]/V;
      TriaxialTest(
         numberOfGrains=int(N),radiusMean=radius,radiusStdDev=rRelFuzz,
         # upperCorner is just size ratio, if radiusMean is specified
         upperCorner=fullDim,
         ## no need to touch any the following
         noFiles=True,lowerCorner=[0,0,0],sigmaIsoCompaction=-1e7,sigmaLateralConfinement=-1e5,compactionFrictionDeg=1,StabilityCriterion=.05,strainRate=.2,thickness=-1,maxWallVelocity=.1,wallOversizeFactor=1.5,autoUnload=True,autoCompressionActivation=False).load()
      while ( numpy.isnan(utils.unbalancedForce()) or utils.unbalancedForce()>0.005 ) :
         O.run(100,True)
      sp=SpherePack(); sp.fromSimulation()
   O.switchScene() ### !!
   _memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri,fullDim)
   if wantPeri: sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2]))
   if orientation:
      sp.cellSize=(0,0,0); # reset periodicity to avoid warning when rotating periodic packing
      sp.rotate(*orientation.toAxisAngle())
   return filterSpherePack(predicate,sp,material=material,color=color,returnSpherePack=returnSpherePack)

kalogeropoulos (antoniskal) said : #3

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

Robert Caulk (rcaulk) said : #4

In case you need non-periodic seed functionality, I've added the seed functionality to randomDensePack and TriaxialTest functions. Simply pass a seed argument:

sp=pack.randomDensePack(pack.inAlignedBox((-x,-y,-z),(x,y,z)),radius=sphereRadius,rRelFuzz=rRelFuzz,spheresInCell=1000,returnSpherePack=True,seed=1)

or

TriaxialTest(seed=1)

Cheers,

Robert