fix particledistribution/discrete not working correctly for ntemp > 2

Submitted by mattkesseler on Fri, 01/06/2017 - 14:23

Hi all, I am currently trying to run a LIGGGHTS script of granular particles sliding down a chute. I have tested the code with a monodisperse size distribution and now want to match the size distribution given for a specific sand that I'm using in my experiments. My understanding is that I should be able to use the following code to produce a distribution with 3 different particle sizes:

fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2650 radius constant 0.012
fix pts2 all particletemplate/sphere 15485867 atom_type 2 density constant 2650 radius constant 0.014
fix pts3 all particletemplate/sphere 15485867 atom_type 3 density constant 2650 radius constant 0.016
fix pdd1 all particledistribution/discrete 32452843 3 pts1 0.3 pts2 0.4 pts3 0.3

However, when I run my simulation and look at the VTK files, it appears to have only generated particles from pts2 and pts3. I am wondering whether this is due to an error in the function or perhaps due to the particle generation method. Here is the rest of the script in case it's relevant.

atom_style granular
atom_modify map array
boundary f f f
newton off

communicate single vel yes

units si

region reg block -0.01 1.01 -2.11 2.51 -0.01 2.51 units box
create_box 2 reg

neighbor 0.003 bin
neigh_modify delay 0

fix m1 all property/global youngsModulus peratomtype 65.e9 65.e9 65.e9 210.e9
fix m2 all property/global poissonsRatio peratomtype 0.25 0.25 0.25 0.25
fix m3 all property/global coefficientRestitution peratomtypepair 4 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893
fix m4 all property/global coefficientFriction peratomtypepair 4 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
fix m5 all property/global k_finnie peratomtypepair 4 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

pair_style gran model hertz tangential history
pair_coeff * *

timestep 5e-7

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix ramp all mesh/surface file meshes/biglabtest40.stl element_exclusion_list read biglabtest40fixed.stl type 4 move 0.0 -2.08981 -0.003 curvature_tolerant yes
fix shut all mesh/surface file meshes/shutterbiglabtest40.stl type 4
fix xwalls1 all wall/gran model hertz tangential history primitive type 4 xplane 0.0
fix xwalls2 all wall/gran model hertz tangential history primitive type 4 xplane 1.0
fix ywalls1 all wall/gran model hertz tangential history primitive type 4 yplane 1.74652
fix granwalls all wall/gran model hertz tangential history mesh n_meshes 2 meshes ramp shut

fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2650 radius constant 0.012
fix pts2 all particletemplate/sphere 15485867 atom_type 2 density constant 2650 radius constant 0.014
fix pts3 all particletemplate/sphere 15485867 atom_type 3 density constant 2650 radius constant 0.016
fix pdd1 all particledistribution/discrete 32452843 3 pts1 0.3 pts2 0.4 pts3 0.3

group nve_group region reg

region bc block 0.0 1.0 1.74652 2.08981 1.60035 1.79186 units box
#region del1 prism 0.0 1.0 0.0 1.77866 1.75356 1.75356 0.0 0.0 0.0 side out
region del2 prism 0.0 1.0 1.74652 1.90722 1.60035 1.79186 0.0 0.0 -0.16070 side out
region del3 prism 0.0 1.0 1.90722 2.08981 1.60035 1.79186 0.0 0.0 0.18259 side out
region wedge intersect 3 bc del2 del3

fix ins nve_group insert/pack seed 32452867 distributiontemplate pdd1 insert_every 20000 overlapcheck yes all_in yes vel constant 0.0 0.0 0.0 volumefraction_region 0.5 region wedge ntry_mc 10000

fix integr nve_group nve/sphere

thermo_style custom spcpu step time cpu cpuremain atoms ke
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes

label loop
variable a loop 20000
variable b equal ke
run 20000
if "${b} < 0.01" then "jump SELF break"
next a
jump SELF loop
label break
print "Settlement complete, granular slide triggered."

unfix ins
unfix granwalls
unfix ywalls1
fix trueramp all wall/gran model hertz tangential history mesh n_meshes 1 meshes ramp

dump dmp all custom 2000 test2/dump*.test2 id type x y z radius vx vy vz
dump dmp2 all mesh/gran/VTK 2000 test2/dump*.vtk id ramp

label loop2
variable c loop 20000
variable d equal ke
run 20000
if "${d} < 0.01" then "jump SELF break2"
next c
jump SELF loop2
label break2
print "Granular slide finished."

There are probably other inefficiencies in the code since this is my first go at this.

Kind regards,

Matt.

Daniel Queteschiner | Mon, 01/09/2017 - 09:34

First of all, in the script you posted you are creating a simulation with 2 types (create_box 2 reg) but want to use 4 types. However, I guess that's just a typo, otherwise you wouldn't be able to run the simulation at all. Second, I suppose that at some point you run out of space for inserting particles. A volume fraction of 0.5 is already rather high for the random particle insertion algorithm used in LIGGGHTS (simple sequential inhibition). In that the case you will see a warning like this:
WARNING: Particle insertion: Less insertions than requested
This means that at some point of the insertion process LIGGGHTS was not able to find an empty space for the particle to be inserted, at which point LIGGGHTS just stops the insertion process (even though smaller particles in the list might still fit into the volume). Since LIGGGHTS sorts the particles from largest to smallest for the insertion procedure, it can easily happen that it just omits all particles from the smallest template. In other words: as soon as you get the warning above, your size distribution will be messed up. In the worst case entire size classes may be missing from your distribution.