use loadGeometry and ConnectionFinder together - it sort of works but doesn't?

Hi,

At the end there is a question, I promise, but first I thought I describe how I got there.

I have a block of unbonded particles of random size that fall down on a surface. I then wait a bit until all the particles have come to rest so that I have an unmoving heap of particles, and then I do other stuff with that heap, like slowly moving it into a barrier to form a wedge. One thing I wanted to do was create bonds between the particles once they have formed the heap, just to see what might happen.

You can't create new bonds in the middle of the simulation (right?), so I take a checkpointer file from when the heap has come to rest and turn it into a .geo file. The .geo-file then gets loaded in with sim.loadGeometry(...) and I add the bonds with sim.createConnections(ConnectionFinder(...))... wait, that doesn't work. Dion said at https://answers.launchpad.net/esys-particle/+question/222478

"It is not possible to combine readGeometry and ConnectionFinder. ConnectionFinder is only designed to work with geometries created in-simulation, via RandomBoxPacker for example."

But there are ways to work around that. We can make (just make, NOT create!) SimpleSpheres that have the exact same attributes as the particles we just loaded in with loadGeometry, put the Simple Spheres in a list and give that list to ConnectionFinder:

sim.readGeometry("geometry.geo")

particleList = sim.getParticleList()
SphereList = []

for particle in particleList:
 SSphere = SimpleSphere(id=particle.getId(), posn=particle.getPosn(), radius=particle.getRadius(), mass=particle.getMass())
 SSphere.setTag(particle.getTag())
 SphereList.append(SSphere)

Does this work? Yes! Is it a good idea? I have no idea.
Instead of loading a .geo file (which automatically creates particles) we can load the checkpointer .txt file and create SimpleSpheres based on what is written there:

particleList = []
with open("test1_chkptr_t=0_1.txt") as inputfile:
 for line in inputfile:
  numbers = line.split()

  if len(numbers) > 5:
   particle = SimpleSphere(id=int(numbers[4]), posn=Vec3(float(numbers[0]), float(numbers[1]), float(numbers[2])), radius=float(numbers[3]), mass=float(numbers[6]))
   particle.setTag(int(numbers[5]))
   particleList.append(particle)

sim.createParticles(particleList)

This also works, and we don't have to worry about mixing SimpleSpheres and NRotParticles. Both methods deliver the exact same results. But do they deliver the exact same results as if I created the spheres with RandomBoxPacker instead of loading them in?

Sort of. To investigate I made a small simulation: 36 particles created by RandomBoxPacker fall from a height of ~ 8m onto a surface over a time of 5s and 5000 time steps. To control particle mass I use setParticleDensity() after the particles are created or loaded. And for some reason the mass is different between creating particles with RandomBoxPacker or loading them. The difference is tiny, but it is there. For example RBP m = 0.198232504 and .geo m = 0.198232503. Where does this difference come from?

Anyway, there is an easy way to get around this. When you use loadGeometry you have to use setParticleDensity afterwards, because .geo files don't include mass information. If you load the particle information from the checkpinter .txt file you can just use that mass and don't have to recalculate it with setParticleDensity.

So now I have two simulations that have identical starting conditions, so they should do the exact same thing. The only difference is that one has particles created by RandomBoxPacker, and the other loads them in.

The start looks good: the checkpointer files for both simulations at t=0 are identical. However, at the end of the simulation there is a significant difference. For example the position of one particle: RBP (0.725675323 0.479624327 0), load (0.699173192 0.611043039 0)

If the particles are all the same and the forces are all the same the simulations should be identical, right? There should be absolutely no randomness in a DEM simulation, so where does the difference come from?

Here is the code for the test simulation. Just comment/uncomment the proper block and change checkpoint name to switch between simulations

from esys.lsm import *
from esys.lsm.util import *
from esys.lsm.geometry import *

# create simulation container object
sim = LsmMpi(numWorkerProcesses = 1, mpiDimList=[1,1,1])
sim.initNeighbourSearch(
 particleType = "NRotSphere",
 gridSpacing = 1.25,
 verletDist = 0.05)

# force 2D computation
sim.force2dComputations(True)

# specify number of timesteps an increment
sim.setNumTimeSteps(5000)
sim.setTimeStepSize(0.001)

# specify spatial domain
sim.setSpatialDomain(
 bBox = BoundingBox(Vec3(0,0,0), Vec3(10,10,0)),
 circDimList = [False,False,False])

# construct rectangle of unbonded particles with non-uniform size
packer = RandomBoxPacker(
 minRadius = 0.25,
 maxRadius = 0.5,
 cubicPackRadius = 1.1,
 maxInsertFails = 1000,
 bBox = BoundingBox(Vec3(0,8,0),Vec3(10,10,0)),
 circDimList=[False,False,False],
 tolerance = 1.0e-5)
packer.generate()
particleList = packer.getSimpleSphereCollection()

for particle in particleList:
 particle.setTag(1)
 sim.createParticle(particle)

"""
# load particles
particleList = []
with open("bondtest1_chkptr_t=0_1.txt", "r") as inputfile:
 for line in inputfile:
  numbers = line.split()

  if len(numbers) > 5:
   particle = SimpleSphere(id=int(numbers[4]), posn=Vec3(float(numbers[0]), float(numbers[1]), float(numbers[2])), radius=float(numbers[3]), mass=float(numbers[6]))
   particle.setTag(int(numbers[5]))
   particleList.append(particle)

sim.createParticles(particleList)
"""

# bind particles
sim.createConnections(
 ConnectionFinder(
  maxDist = 0.05,
  bondTag = 1,
  pList = particleList))

# specify the bond between particles
BondGrp = sim.createInteractionGroup(
 NRotBondPrms(
  name = "particleBonds",
  normalK = 100.0,
  breakDistance = 3,
  tag = 1,
  scaling = True))

# add wall to the bottom of the simulation
sim.createWall(
 name = "bottomWall",
 posn = Vec3(0,0,0),
 normal = Vec3(0,1,0))

# add repulsion to floor
sim.createInteractionGroup(
 NRotElasticWallPrms(
  name = "bottomWallRepulsion",
  wallName = "bottomWall",
  normalK = 10000))

# friction between particles
sim.createInteractionGroup(
 NRotFrictionPrms(
  name = "pp-friction",
  normalK = 1000,
  dynamicMu = 0.6,
  shearK = 100,
  scaling = True))

# set particle density
# sim.setParticleDensity(tag = 1, mask = 0, Density = 1.0)

# add gravity
sim.createInteractionGroup(
 GravityPrms(name="gravity", acceleration = Vec3(0,-9.81,0)))

# add damping
sim.createInteractionGroup(
 LinDampingPrms(
  name = "damping",
  viscosity = 1,
  maxIterations = 100))

# create check pointer as data output
sim.createCheckPointer(
 CheckPointPrms(
  fileNamePrefix = "bondtest1_chkptr",
  beginTimeStep = 0,
  endTimeStep = 5000,
  timeStepIncr = 1000))

# execute simulation
sim.run()