Problems in 3D simulation with PSC

Dear all :

    I have tested this simulation again and find that in Yadedaily 1.12, the problem is remained. So I give all the information about it as follow:

(1). I have checked that there is nothing different between my script and psd.py in yadedaily.doc in particle generation with a given psd. The key words in both are:

1. mincorner and maxcorner to define the generation box.
2. particles number.
3. psdSizes
4. psdCumm
5. distributeMass = 0 or 1

all above are just in sp0.makeCloud() classes.

(2). There are two psd.txt for the test. As mentioned above, one has a small value of the ratio of min size to the max and the other has the value been approximately equal to 1.0.

(3). With no change with the codes, the value error warning ocurrs very soon after the simulation starts when the first data is used as the psd generation but there is nothing appearing when the second psd data is employed.

I am lost in it and need your help. The simulation codes and condition, as well as the error warning are as follow.

codes:
#########################################################################################
# unicode: UTF-8
Filename='3D-triaxial compression test with a given PSD'
from yade import pack,os,utils
import math

'''
The purpose of this simulation is to run sucessfully a 2D-biaxial compression
with a pre-define PSD in Yade DEM platform and then to compare its results to
lab test achievements.

The whole processing of this simulation can be divided into 3 parts, which can
be described as follow:

[1].material defination. Of course, there are just two types of materials needed
to be defined, the first one is the particles-material, cohefrict-cohefrict material
(incohesive but can consider rolling mechanism)
and the other is the walls-material, frict-frict material(incohesive).

[2].Bodies generation. In this simulation, just two types of bodies should be considered,
particles which here are simplified as idealistic spheres and walls which are introduced
for representing the sample's boundary. These two here are treated as rigid body. And the
particles are generated with a given PSD.

[3].Running the simulation. The whole process of this simulation composes of just two parts,
which can be called, of course, the isotropic consolidation process and the triaxial
compression process. So this step is to define the calculation engines for consolidation.
Despite that the algorithms of contact detection in Yade (it can be said that in all of
the DEM software here) and of random packing in certain space make it difficult to satisfy
that the initial porosity of a given simulation sample at the start of the consolidtaion
is as the same as it of one in lab, some mathematics methods can be employed as a approximate
solution, such as radius expansion menthod, to obtain a appropriate porosity with a certain
confine pressure.Some more details are in the below section.
'''

################################## materials defination ###############################
'''
O.materials.append(FrictMat(young=1e9,poisson=0.5,frictionAngle=radians(20),
                    density=2600,label='particles_mat'))
O.materials.append(FrictMat(young=1e9,poisson=0.5,frictionAngle=0,
                    density=0,label='walls_mat'))
'''
# Read the simulation condition.
condition_data = []
f = file('condition.txt', 'r')
for line in f.readlines():
    condition_data.append(line.strip('\n').split(':'))
for i in xrange(len(condition_data)):
    if i < 2 :
        condition_data[i][1] = condition_data[i][1].split(',')
        for j in xrange(len(condition_data[i][1])):
            condition_data[i][1][j] = float(condition_data[i][1][j])
        condition_data[i][1] = tuple(condition_data[i][1])
    elif i == 2 :
        condition_data[i][1] = int(condition_data[i][1])
    elif i == 3 :
        condition_data[i][1] = bool(float(condition_data[i][1]))
    else:
        condition_data[i][1] = float(condition_data[i][1])
f.close()
# Particles-material defination
material_list_1 = []
material_list_1.append(8e8) # list[0]='young module'
material_list_1.append(0.5) # list[1]='poisson ratio'
material_list_1.append(0.59) # list[2]='friction Angle'
material_list_1.append(2650) # list[3]='density'
material_list_1.append(0.3) # list[4]='rolling coefficient'
material_list_1.append(1.0) # list[5]='considering plastic moment of rolling'
compFricDegree = material_list_1[2] # for radius expansion
particles_mat = O.materials.append(CohFrictMat( young = material_list_1[0],
                                                poisson = material_list_1[1],
                                                frictionAngle = material_list_1[2],
                                                density =material_list_1[3],
                                                isCohesive = False,
                                                momentRotationLaw = True,
                                                alphaKr = material_list_1[4],
                                                etaRoll = material_list_1[5]))

# walls-material defination
material_list_2 = []
material_list_2.append(1e9) # list[0]='young module'
material_list_2.append(0.5) # list[1]='poisson ratio'
material_list_2.append(radians(0)) # list[2]='friction Angle'
material_list_2.append(0) # list[3]='density'
walls_mat = O.materials.append(FrictMat(young=material_list_2[0],
                                        poisson=material_list_2[1],
                                        frictionAngle=material_list_2[2],
                                        density=material_list_2[3]))

################################## bodies defination ###############################

# walls generation.
mincorner = condition_data[0][1]
maxcorner = condition_data[1][1]
walls_width = 1e-9
# In 2D simulation, spheres generation corners are:
# genmincorner = [mincorner[0],mincorner[1],(maxcorner[2]-mincorner[2])/2.0]
# genmaxcorner = [maxcorner[0],maxcorner[1],(maxcorner[2]-mincorner[2])/2.0]
# In 3D simulation, spheres generation corners are:
genmincorner = [mincorner[0]+walls_width,mincorner[1]+walls_width,mincorner[2]+walls_width]
genmaxcorner = [maxcorner[0]-walls_width,maxcorner[1]-walls_width,maxcorner[2]-walls_width]
walls = aabbWalls([mincorner,maxcorner],thickness=walls_width,material=walls_mat)
wallids=O.bodies.append(walls)
left_wall_posi = O.bodies[0].state.pos[0] # pos[0]--x direction
right_wall_posi = O.bodies[1].state.pos[0]
bottom_wall_posi = O.bodies[2].state.pos[1] # pos[1]--y direction
top_wall_posi = O.bodies[3].state.pos[1]
back_wall_posi = O.bodies[4].state.pos[2] # pos[2]--z direction
front_wall_posi = O.bodies[5].state.pos[2]
x_length = abs(right_wall_posi - left_wall_posi)
y_length = abs(top_wall_posi - bottom_wall_posi)
z_length = abs(front_wall_posi - back_wall_posi)
# In 3D simulation:
sample_volume = x_length * y_length * z_length
# In 2D simualtion:
# sample_volume = x_length * y_length
print 'The x_length of sample is: ',x_length
print 'The y_length of sample is: ',y_length
print 'The z_length of sample is: ',z_length
print 'The volume of sample is: ',sample_volume,'\n'

# particles generation with a given PSD curve. The objective is to ensure
# that all the particles generated in packbox can have different unique
# sizes and all these sizes can fix approximately the PSD data of realistic
# granular materials. The programs for this purpose are composed of 7 parts,
# which are as follows:

# The first step: read PSD data for pre-generation.
f = file('psd.txt','r')
read_PSD_data_str = []; read_PSD_data_flo = []
for line in f.readlines():
    datalist = line.split()
    read_PSD_data_str.append(datalist)
for i in xrange(1,len(read_PSD_data_str),1):
    read_PSD_data_flo.append([])
    for j in read_PSD_data_str[i]:
        j = float(j)
        read_PSD_data_flo[i-1].append(j)

# The second step: particles generating in makecloud with this given PSD.
PSD_diameters_list = []; PSD_frequencies_list = [];PSD = []
for i in xrange(len(read_PSD_data_flo)):
    PSD_diameters_list.append(read_PSD_data_flo[i][0])
    PSD_frequencies_list.append(read_PSD_data_flo[i][1])
PSD.append(genmincorner)
PSD.append(genmaxcorner)
PSD.append(condition_data[2][1])
PSD.append(PSD_diameters_list)
PSD.append(PSD_frequencies_list)
sp0 = pack.SpherePack();
sp0.makeCloud(PSD[0],PSD[1],num = PSD[2],psdSizes = PSD[3],
              psdCumm = PSD[4],distributeMass=True)

# The third step: particles generating from makecloud to yade.
sp0.toSimulation(material=particles_mat)
radius_list = []
for b in O.bodies:
    if isinstance(b.shape,Sphere):
        init_particles_rad = b.shape.radius
        radius_list.append(init_particles_rad)
'''
# The fourth step: blockedDOFs for 2D simulation tasks.
for b in O.bodies:
    if isinstance(b.shape,Sphere):
        b.state.blockedDOFs='zXY'
'''
# Making an output to test the particles number and the min or max radii.
min_diameter = 2*min(radius_list);
max_diameter = 2*max(radius_list)
print 'The total number of particles in simulation is: ',len(radius_list)
print 'The initial max diameter is: ',max_diameter
print 'The initial min diameter is: ',min_diameter,'\n'

# Making post-processing for PSD, which may include:
# [1]. particle diameter VS relative frequency curve by munber accumulation.
# [2]. particle diameter VS cumulative frequency curve by munber accumulation.
# [3]. particle diameter VS cumulative frequency curve by volume(mass) accumulation.
particle_dia_VS_num_rel = []
particle_dia_VS_num_cum = []
particle_dia_VS_mass_cum = []
maxd = max_diameter
mind = min_diameter
mund = len(PSD_diameters_list)-1
for i in xrange(mund):
    particle_dia_VS_num_rel.append([])
    particle_dia_VS_num_cum.append([])
    particle_dia_VS_mass_cum.append([])
    lower_level_limit = mind + (maxd-mind)*i/mund
    upper_level_limit = mind + (maxd-mind)*(i+1)/mund
    for j in radius_list:
        if lower_level_limit <= 2*j and 2*j <= upper_level_limit:
            particle_dia_VS_num_rel[i].append(2*j)
        if 2*j <= upper_level_limit:
            particle_dia_VS_num_cum[i].append(2*j)
            particle_dia_VS_mass_cum[i].append((4/3.0)*(math.pi)*j*j*j)

# Making output of PSD data, which may include:
# [1]. the data of diameter VS relative frequency curve in number based.
# [2]. the data of diameter VS commulative frequency curve in number based.
# [3]. mean radii of particles according to cummulative number.
# [4]. the data of diameter VS commulative frequency curve in mass based.
# [5]. mean radii of particles according to cummulative mass.
# [6]. the porosity of this sample.
print 'the data of diameter VS relative frequency curve in number based are: '
for i in xrange(mund):
    print float(len(particle_dia_VS_num_rel[i]))/len(radius_list)
print '\n'

print 'the data of diameter VS commulative frequency curve in number based are: '
for i in xrange(mund):
    print float(len(particle_dia_VS_num_cum[i]))/len(radius_list)
print '\n'

print 'mean radii of particles according to cummulative number is: '
mean_radii_size = sum(radius_list)/len(radius_list)
print mean_radii_size,'\n'

ori_particles_sum_volume = []
for i in xrange(len(radius_list)):
    orirad = radius_list[i]
    ori_particle_volume = ((4/3.0)*(math.pi)*(orirad**3))
    ori_particles_sum_volume.append(ori_particle_volume)
initial_particles_vol = float(sum(ori_particles_sum_volume))
print 'the data of diameter VS commulative frequency curve in mass based are: '
for i in xrange(mund):
    print float(sum(particle_dia_VS_mass_cum[i]))/initial_particles_vol
print '\n'

print 'mean radii of particles according to cummulative mass is: '
mean_radius_vol = ((initial_particles_vol)/((4/3.0)*(math.pi)*len(radius_list)))**(1.0/3.0)
print mean_radius_vol,'\n'

print 'the porosity of this sample is: '
initial_porosity = 1-(initial_particles_vol)/sample_volume
print initial_porosity,'\n'

# Making an output for getting all the desired data,which include:
# [1]. PSD-data.dat
# [2]. radius-data.dat
f=file("PSD-information.csv",'a')
f.write('%-12s,%-12s,%-12s,%-12s,\n'%('Diameter','Frequency',
                                    'Cum-number','Cum-mass'))
for i in xrange(mund):
    lower_level_limit = mind + (maxd-mind)*i/mund
    upper_level_limit = mind + (maxd-mind)*(i+1)/mund
    a = ( lower_level_limit + upper_level_limit )/2
    b = float(len(particle_dia_VS_num_rel[i]))/len(radius_list)
    c = float(len(particle_dia_VS_num_cum[i]))/len(radius_list)
    d = float(sum(particle_dia_VS_mass_cum[i]))/initial_particles_vol
    f.write('%12.5E,%12.5E,%12.5E,%12.5E,\n'%(a,b,c,d))
f.close()

############################## consolidation engines defination #####################
sample_density = condition_data[3][1]
confine_pressure = condition_data[4][1]
strain_rate = condition_data[5][1]
final_axial_strain = condition_data[6][1]
print '#####################################################'

consolidation = TriaxialStressController(
# defining the ids of walls, of which the normal directions are respectively parallel to x,y,z-coordinate.
    wall_left_id=wallids[0],wall_right_id=wallids[1],
    wall_bottom_id=wallids[2],wall_top_id=wallids[3],
    wall_back_id=wallids[4],wall_front_id=wallids[5],
    wall_front_activated = False,wall_back_activated = False,
# In consolidation process, choosing the 'radius expansion' method which can provide the chance that
# all the radius of particles may increase with the certain Multiplier until the target confine-pressure
# or porosity are attained. At this time, we choose the confine-pressure as the target.
    maxMultiplier = 1.01, # spheres growing factor (fast growth)
    finalMaxMultiplier = 1.001, # spheres growing factor (slow growth)
    thickness=walls_width,
    internalCompaction=True, # 'False' here may close the task

# switch stress/strain control using a bitmask. What is a bitmask ?
# Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
# Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
# to put it differently, the mask is the integer whose binary representation is xyz, i.e.
# "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
    stressMask=7,
        goal1 = - confine_pressure, # confine-pressure applied in x-derection.
        goal2 = - confine_pressure, # confine-pressure applied in y-derection.
        goal3 = - confine_pressure, # confine-pressure applied in z-derection.
# Defining the max velocity of walls to make sure the isotropic conpression in the whole process
# max_vel = 10,
)

# Define consolidation results output
f=file("consolidation.csv",'a')
f.write('%-12s,%-12s,%-12s,\n'%('Iter','confine-prs','porosity'))
f.close()

def consolidation_results():
# Define mean sigma as the confinment pressure.
    sigma_1 = -consolidation.stress(consolidation.wall_top_id)[1]
    sigma_2 = -consolidation.stress(consolidation.wall_right_id)[0]
    sigma_3 = -consolidation.stress(consolidation.wall_front_id)[2]
    meanS = (sigma_1 + sigma_2 + sigma_3)/3.0
# Define variable porosity in the process of consolidation.
    sum_particles_volume = []; max_particles_rad = []
    for b in O.bodies:
        if isinstance(b.shape,Sphere):
            rad = b.shape.radius
            single_particles_volume = (4.0/3.0)*(math.pi)*(rad**3)
            sum_particles_volume.append(single_particles_volume)
            max_particles_rad.append(rad)
    total_particles_volume = sum(sum_particles_volume)
    max_particles_radii = max(max_particles_rad)
    print 'the max radii of particles is: ',max_particles_radii
    print 'the total particle volume is: ',total_particles_volume
    porosity = 1 - (total_particles_volume/sample_volume)
# Make all the things above outputing
    f=file("consolidation.csv",'a')
    f.write('%12.5E,%12.5E,%12.5E,\n'%(O.iter,meanS,porosity)),
    f.close()

# Define simulation engine
O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
    InteractionLoop(
# Box-sphere interactions will be the simple normal-shear law, we use ScGeom for them
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
# Boxes will be frictional (FrictMat), so the sphere-box physics is FrictMat vs. CohFrictMat.
        [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
# Finally, two different contact laws for sphere-box and sphere-sphere
        [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
            useIncrementalForm=True,
            always_use_moment_law=True,
            label='cohesiveLaw')]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    consolidation,
    NewtonIntegrator(damping=.2),
# PyRunner(command='consolidation_results()',iterPeriod=200)
]

# The class 'PWaveTimeStep' is used as 'dt'.
O.dt = utils.PWaveTimeStep()
O.usesTimeStepper = True

for i in xrange(10):
    O.run(100,True)
    consolidation_results()
O.wait()

simulation condition: (Please copy all the condition in the 'condition.txt')
########################################################################################
The min corner of the box is: 0,0,0
The max corner of the box is: 0.2,0.6,0.2
The total number of particles is: 3000
Shall we use the dense sample: 1
Confine-pressure is: 500000
Strain rate in y-direction is: 0.25
The final axial strain is: 0.2

psd-1 and psd-2 (Please make sure that they are renamed as 'psd.txt' when they are used)

########################################################################################
psd-1:

diameters percentage
0.250 0
0.297 0.0004
0.344 0.0027
0.391 0.0114
0.438 0.0384
0.484 0.1032
0.531 0.2242
0.578 0.4002
0.625 0.5997
0.672 0.7757
0.719 0.8967
0.766 0.9615
0.813 0.9885
0.860 0.9972
0.906 0.9995
0.953 0.9999
1.000 1

###########################################################################

psd-2:

diameters percentage
0.984 0
0.985 0.0004
0.986 0.0027
0.987 0.0114
0.988 0.0384
0.989 0.1032
0.990 0.2242
0.991 0.4002
0.992 0.5997
0.993 0.7757
0.994 0.8967
0.995 0.9615
0.996 0.9885
0.997 0.9972
0.998 0.9995
0.999 0.9999
1.000 1

#########################################################

error information:
###########################################################
ValueError Traceback (most recent call last)
/usr/lib/x86_64-linux-gnu/yadedaily/py/yade/__init__.py in refreshEvent(self)
    196 def zxySlot(self): self.setViewAxes((0,-1,0),(1,0,0))
    197 def refreshEvent(self):
--> 198 self.refreshValues()
    199 self.activateControls()
    200 def deactivateControls(self):

/usr/lib/x86_64-linux-gnu/yadedaily/py/yade/__init__.py in refreshValues(self)
    269 self.iterLabel.setText('#%ld / %ld, %.1f/s %s'%(O.iter,stopAtIter,self.iterPerSec,subStepInfo))
    270 if t!=float('inf'):
--> 271 s=int(t); ms=int(t*1000)%1000; us=int(t*1000000)%1000; ns=int(t*1000000000)%1000
    272 self.virtTimeLabel.setText(u'%03ds%03dm%03dμ%03dn'%(s,ms,us,ns))
    273 else: self.virtTimeLabel.setText(u'[ ∞ ] ?!')

ValueError: cannot convert float NaN to integer

Dear Prof. Chareyre and all users:

    As far as I concerned, there maybe something error in the combination of PSD generation algorithm and cohesive contact Law.

    In order to check this problem, I added a group of 'if' sentences in the whole codes above for the definition of particle material and O. engine, which are as follow:

##############################################################################
materials:

coFrictMat = 1

if coFrictMat == 1:
    print 'rolling is considered'
    particles_mat = O.materials.append(CohFrictMat( young = material_list_1[0],
                                                    poisson = material_list_1[1],
                                                    frictionAngle = material_list_1[2],
                                                    density =material_list_1[3],
                                                    isCohesive = False,
                                                    momentRotationLaw = True,
                                                    alphaKr = material_list_1[4],
                                                    etaRoll = material_list_1[5]))
if coFrictMat == 0:
    print 'no rolling'
    particles_mat = O.materials.append(FrictMat( young = material_list_1[0],
                                                 poisson = material_list_1[1],
                                                 frictionAngle = material_list_1[2],
                                                 density =material_list_1[3]))
#############################################################################

contact law:

if coFrictMat == 1:
    O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
        InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
            [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
                useIncrementalForm=True,
                always_use_moment_law=True,
                label='cohesiveLaw')]
        ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        consolidation,
        NewtonIntegrator(damping=.2),
    ]

if coFrictMat == 0:
    O.engines=[
     ForceResetter(),
     InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
     InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
     ),
     GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
     consolidation,
     NewtonIntegrator(damping=.2)
    ]
###############################################################################

Then use the psd-1 data as generation. When ' coFrictMat == 0' , there is nothing error but when ' coFrictMat == 1' , error warning may appear. I am lost in it.

Seeking your help!

# unicode: UTF-8
Filename='3D-triaxial compression test with a given PSD'
from yade import pack,os,utils
import math

'''
The purpose of this simulation is to run sucessfully a 2D-biaxial compression
with a pre-define PSD in Yade DEM platform and then to compare its results to
lab test achievements.

The whole processing of this simulation can be divided into 3 parts, which can
be described as follow:

[1].material defination. Of course, there are just two types of materials needed
to be defined, the first one is the particles-material, cohefrict-cohefrict material
(incohesive but can consider rolling mechanism)
and the other is the walls-material, frict-frict material(incohesive).

[2].Bodies generation. In this simulation, just two types of bodies should be considered,
particles which here are simplified as idealistic spheres and walls which are introduced
for representing the sample's boundary. These two here are treated as rigid body. And the
particles are generated with a given PSD.

[3].Running the simulation. The whole process of this simulation composes of just two parts,
which can be called, of course, the isotropic consolidation process and the triaxial
compression process. So this step is to define the calculation engines for consolidation.
Despite that the algorithms of contact detection in Yade (it can be said that in all of
the DEM software here) and of random packing in certain space make it difficult to satisfy
that the initial porosity of a given simulation sample at the start of the consolidtaion
is as the same as it of one in lab, some mathematics methods can be employed as a approximate
solution, such as radius expansion menthod, to obtain a appropriate porosity with a certain
confine pressure.Some more details are in the below section.
'''

################################## materials defination ###############################
'''
O.materials.append(FrictMat(young=1e9,poisson=0.5,frictionAngle=radians(20),
                    density=2600,label='particles_mat'))
O.materials.append(FrictMat(young=1e9,poisson=0.5,frictionAngle=0,
                    density=0,label='walls_mat'))
'''
# Read the simulation condition.
condition_data = []
f = file('condition.txt', 'r')
for line in f.readlines():
    condition_data.append(line.strip('\n').split(':'))
for i in xrange(len(condition_data)):
    if i < 2 :
        condition_data[i][1] = condition_data[i][1].split(',')
        for j in xrange(len(condition_data[i][1])):
            condition_data[i][1][j] = float(condition_data[i][1][j])
        condition_data[i][1] = tuple(condition_data[i][1])
    elif i == 2 :
        condition_data[i][1] = int(condition_data[i][1])
    elif i == 3 :
        condition_data[i][1] = bool(float(condition_data[i][1]))
    else:
        condition_data[i][1] = float(condition_data[i][1])
f.close()
# Particles-material defination
material_list_1 = []
material_list_1.append(8e8) # list[0]='young module'
material_list_1.append(0.5) # list[1]='poisson ratio'
material_list_1.append(0.59) # list[2]='friction Angle'
material_list_1.append(2650) # list[3]='density'
material_list_1.append(0.3) # list[4]='rolling coefficient'
material_list_1.append(1.0) # list[5]='considering plastic moment of rolling'
compFricDegree = material_list_1[2] # for radius expansion

coFrictMat = 0

if coFrictMat == 1:
    print 'rolling is considered'
    particles_mat = O.materials.append(CohFrictMat( young = material_list_1[0],
                                                    poisson = material_list_1[1],
                                                    frictionAngle = material_list_1[2],
                                                    density =material_list_1[3],
                                                    isCohesive = False,
                                                    momentRotationLaw = True,
                                                    alphaKr = material_list_1[4],
                                                    etaRoll = material_list_1[5]))
if coFrictMat == 0:
    print 'no rolling'
    particles_mat = O.materials.append(FrictMat( young = material_list_1[0],
                                                 poisson = material_list_1[1],
                                                 frictionAngle = material_list_1[2],
                                                 density =material_list_1[3]))
# walls-material defination

material_list_2 = []
material_list_2.append(1e9) # list[0]='young module'
material_list_2.append(0.5) # list[1]='poisson ratio'
material_list_2.append(radians(0)) # list[2]='friction Angle'
material_list_2.append(0) # list[3]='density'
walls_mat = O.materials.append(FrictMat(young=material_list_2[0],
                                        poisson=material_list_2[1],
                                        frictionAngle=material_list_2[2],
                                        density=material_list_2[3]))

################################## bodies defination ###############################

# walls generation.
mincorner = condition_data[0][1]
maxcorner = condition_data[1][1]
walls_width = 1e-9
# In 2D simulation, spheres generation corners are:
# genmincorner = [mincorner[0],mincorner[1],(maxcorner[2]-mincorner[2])/2.0]
# genmaxcorner = [maxcorner[0],maxcorner[1],(maxcorner[2]-mincorner[2])/2.0]
# In 3D simulation, spheres generation corners are:
genmincorner = [mincorner[0]+walls_width,mincorner[1]+walls_width,mincorner[2]+walls_width]
genmaxcorner = [maxcorner[0]-walls_width,maxcorner[1]-walls_width,maxcorner[2]-walls_width]
walls = aabbWalls([mincorner,maxcorner],thickness=walls_width,material=walls_mat)
wallids=O.bodies.append(walls)
left_wall_posi = O.bodies[0].state.pos[0] # pos[0]--x direction
right_wall_posi = O.bodies[1].state.pos[0]
bottom_wall_posi = O.bodies[2].state.pos[1] # pos[1]--y direction
top_wall_posi = O.bodies[3].state.pos[1]
back_wall_posi = O.bodies[4].state.pos[2] # pos[2]--z direction
front_wall_posi = O.bodies[5].state.pos[2]
x_length = abs(right_wall_posi - left_wall_posi)
y_length = abs(top_wall_posi - bottom_wall_posi)
z_length = abs(front_wall_posi - back_wall_posi)
# In 3D simulation:
sample_volume = x_length * y_length * z_length
# In 2D simualtion:
# sample_volume = x_length * y_length
print 'The x_length of sample is: ',x_length
print 'The y_length of sample is: ',y_length
print 'The z_length of sample is: ',z_length
print 'The volume of sample is: ',sample_volume,'\n'

# particles generation with a given PSD curve. The objective is to ensure
# that all the particles generated in packbox can have different unique
# sizes and all these sizes can fix approximately the PSD data of realistic
# granular materials. The programs for this purpose are composed of 7 parts,
# which are as follows:

# The first step: read PSD data for pre-generation.
f = file('psd.txt','r')
read_PSD_data_str = []; read_PSD_data_flo = []
for line in f.readlines():
    datalist = line.split()
    read_PSD_data_str.append(datalist)
for i in xrange(1,len(read_PSD_data_str),1):
    read_PSD_data_flo.append([])
    for j in read_PSD_data_str[i]:
        j = float(j)
        read_PSD_data_flo[i-1].append(j)

# The second step: particles generating in makecloud with this given PSD.
PSD_diameters_list = []; PSD_frequencies_list = [];PSD = []
for i in xrange(len(read_PSD_data_flo)):
    PSD_diameters_list.append(read_PSD_data_flo[i][0])
    PSD_frequencies_list.append(read_PSD_data_flo[i][1])
PSD.append(genmincorner)
PSD.append(genmaxcorner)
PSD.append(condition_data[2][1])
PSD.append(PSD_diameters_list)
PSD.append(PSD_frequencies_list)
sp0 = pack.SpherePack();
sp0.makeCloud(PSD[0],PSD[1],num = PSD[2],psdSizes = PSD[3],
              psdCumm = PSD[4],distributeMass=True)

# The third step: particles generating from makecloud to yade.
sp0.toSimulation(material=particles_mat)
radius_list = []
for b in O.bodies:
    if isinstance(b.shape,Sphere):
        init_particles_rad = b.shape.radius
        radius_list.append(init_particles_rad)

# Making an output to test the particles number and the min or max radii.
min_diameter = 2*min(radius_list);
max_diameter = 2*max(radius_list)
print 'The total number of particles in simulation is: ',len(radius_list)
print 'The initial max diameter is: ',max_diameter
print 'The initial min diameter is: ',min_diameter,'\n'

# Making post-processing for PSD, which may include:
# [1]. particle diameter VS relative frequency curve by munber accumulation.
# [2]. particle diameter VS cumulative frequency curve by munber accumulation.
# [3]. particle diameter VS cumulative frequency curve by volume(mass) accumulation.
particle_dia_VS_num_rel = []
particle_dia_VS_num_cum = []
particle_dia_VS_mass_cum = []
maxd = max_diameter
mind = min_diameter
mund = len(PSD_diameters_list)-1
for i in xrange(mund):
    particle_dia_VS_num_rel.append([])
    particle_dia_VS_num_cum.append([])
    particle_dia_VS_mass_cum.append([])
    lower_level_limit = mind + (maxd-mind)*i/mund
    upper_level_limit = mind + (maxd-mind)*(i+1)/mund
    for j in radius_list:
        if lower_level_limit <= 2*j and 2*j <= upper_level_limit:
            particle_dia_VS_num_rel[i].append(2*j)
        if 2*j <= upper_level_limit:
            particle_dia_VS_num_cum[i].append(2*j)
            particle_dia_VS_mass_cum[i].append((4/3.0)*(math.pi)*j*j*j)

# Making output of PSD data, which may include:
# [1]. the data of diameter VS relative frequency curve in number based.
# [2]. the data of diameter VS commulative frequency curve in number based.
# [3]. mean radii of particles according to cummulative number.
# [4]. the data of diameter VS commulative frequency curve in mass based.
# [5]. mean radii of particles according to cummulative mass.
# [6]. the porosity of this sample.
print 'the data of diameter VS relative frequency curve in number based are: '
for i in xrange(mund):
    print float(len(particle_dia_VS_num_rel[i]))/len(radius_list)
print '\n'

print 'the data of diameter VS commulative frequency curve in number based are: '
for i in xrange(mund):
    print float(len(particle_dia_VS_num_cum[i]))/len(radius_list)
print '\n'

print 'mean radii of particles according to cummulative number is: '
mean_radii_size = sum(radius_list)/len(radius_list)
print mean_radii_size,'\n'

ori_particles_sum_volume = []
for i in xrange(len(radius_list)):
    orirad = radius_list[i]
    ori_particle_volume = ((4/3.0)*(math.pi)*(orirad**3))
    ori_particles_sum_volume.append(ori_particle_volume)
initial_particles_vol = float(sum(ori_particles_sum_volume))
print 'the data of diameter VS commulative frequency curve in mass based are: '
for i in xrange(mund):
    print float(sum(particle_dia_VS_mass_cum[i]))/initial_particles_vol
print '\n'

print 'mean radii of particles according to cummulative mass is: '
mean_radius_vol = ((initial_particles_vol)/((4/3.0)*(math.pi)*len(radius_list)))**(1.0/3.0)
print mean_radius_vol,'\n'

print 'the porosity of this sample is: '
initial_porosity = 1-(initial_particles_vol)/sample_volume
print initial_porosity,'\n'

# Making an output for getting all the desired data,which include:
# [1]. PSD-data.dat
# [2]. radius-data.dat
f=file("PSD-information.csv",'a')
f.write('%-12s,%-12s,%-12s,%-12s,\n'%('Diameter','Frequency',
                                    'Cum-number','Cum-mass'))
for i in xrange(mund):
    lower_level_limit = mind + (maxd-mind)*i/mund
    upper_level_limit = mind + (maxd-mind)*(i+1)/mund
    a = ( lower_level_limit + upper_level_limit )/2
    b = float(len(particle_dia_VS_num_rel[i]))/len(radius_list)
    c = float(len(particle_dia_VS_num_cum[i]))/len(radius_list)
    d = float(sum(particle_dia_VS_mass_cum[i]))/initial_particles_vol
    f.write('%12.5E,%12.5E,%12.5E,%12.5E,\n'%(a,b,c,d))
f.close()

############################## consolidation engines defination #####################
sample_density = condition_data[3][1]
confine_pressure = condition_data[4][1]
strain_rate = condition_data[5][1]
final_axial_strain = condition_data[6][1]
print '#####################################################'

consolidation = TriaxialStressController(
# defining the ids of walls, of which the normal directions are respectively parallel to x,y,z-coordinate.
    wall_left_id=wallids[0],wall_right_id=wallids[1],
    wall_bottom_id=wallids[2],wall_top_id=wallids[3],
    wall_back_id=wallids[4],wall_front_id=wallids[5],
    wall_front_activated = False,wall_back_activated = False,
# In consolidation process, choosing the 'radius expansion' method which can provide the chance that
# all the radius of particles may increase with the certain Multiplier until the target confine-pressure
# or porosity are attained. At this time, we choose the confine-pressure as the target.
    maxMultiplier = 1.01, # spheres growing factor (fast growth)
    finalMaxMultiplier = 1.001, # spheres growing factor (slow growth)
    thickness=walls_width,
    internalCompaction=True, # 'False' here may close the task

# switch stress/strain control using a bitmask. What is a bitmask ?
# Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
# Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
# to put it differently, the mask is the integer whose binary representation is xyz, i.e.
# "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
    stressMask=7,
        goal1 = - confine_pressure, # confine-pressure applied in x-derection.
        goal2 = - confine_pressure, # confine-pressure applied in y-derection.
        goal3 = - confine_pressure, # confine-pressure applied in z-derection.
# Defining the max velocity of walls to make sure the isotropic conpression in the whole process
# max_vel = 10,
)

# Define consolidation results output
f=file("consolidation.csv",'a')
f.write('%-12s,%-12s,%-12s,\n'%('Iter','confine-prs','porosity'))
f.close()

def consolidation_results():
# Define mean sigma as the confinment pressure.
    sigma_1 = -consolidation.stress(consolidation.wall_top_id)[1]
    sigma_2 = -consolidation.stress(consolidation.wall_right_id)[0]
    sigma_3 = -consolidation.stress(consolidation.wall_front_id)[2]
    meanS = (sigma_1 + sigma_2 + sigma_3)/3.0
# Define variable porosity in the process of consolidation.
    sum_particles_volume = []; max_particles_rad = []
    for b in O.bodies:
        if isinstance(b.shape,Sphere):
            rad = b.shape.radius
            single_particles_volume = (4.0/3.0)*(math.pi)*(rad**3)
            sum_particles_volume.append(single_particles_volume)
            max_particles_rad.append(rad)
    total_particles_volume = sum(sum_particles_volume)
    max_particles_radii = max(max_particles_rad)
    print 'the max radii of particles is: ',max_particles_radii
    print 'the total particle volume is: ',total_particles_volume
    porosity = 1 - (total_particles_volume/sample_volume)
# Make all the things above outputing
    f=file("consolidation.csv",'a')
    f.write('%12.5E,%12.5E,%12.5E,\n'%(O.iter,meanS,porosity)),
    f.close()

# Define simulation engine
if coFrictMat == 1:
    O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
        InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
            [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
                useIncrementalForm=True,
                always_use_moment_law=True,
                label='cohesiveLaw')]
        ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        consolidation,
        NewtonIntegrator(damping=.2),
    ]
if coFrictMat == 0:
    O.engines=[
     ForceResetter(),
     InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
     InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
     ),
     GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
     consolidation,
     NewtonIntegrator(damping=.2)
    ]

# The class 'PWaveTimeStep' is used as 'dt'.
O.dt = utils.PWaveTimeStep()
O.usesTimeStepper = True

for i in xrange(5):
    O.run(100,True)
    consolidation_results()

I don't understand the problem. I'm sorry.

Hi,

    All the codes above are about to model the shearing bands in PSC with a given PSD . In this case, particle rolling should be considerd as one of the most important mechanics in each contact for those macro bands formation and evolution. So:

(1). I employed the 'coFrictMat' as the particle material and 'cohesiveLaw' as the contact law to consider relative rolling. The detials of engine definition in my aodes above are as follow:
########################################################################
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
        InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
            [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
                useIncrementalForm=True,
                always_use_moment_law=True,
                label='cohesiveLaw')]
        ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        consolidation,
        NewtonIntegrator(damping=.2),
    ]
######################################################################

(2). particle generation with a given PSD. In Yade, the psd.py in examples suggested me to use sp0.makeCloud() class
for this purpose. In this class, There are 5 keys:

1. mincorner and maxcorner to define the generation box.
2. particles number.
3. psdSizes
4. psdCumm
5. distributeMass = 0 or 1

Now the formation of psdSizes and psdCumm in final codes are inputed from a 'psd.txt' data, which is #4.

(3). I have tested that if I need to consider the two things above simultaneously, there are something wrong when the running starts.

(4). In my own testing process, I choose the two different contact laws, one is labeled 'cohesiveLaw' as mentioned above and the other is the most common law 'Ip2_FrictMat_FrictMat_FrictPhys()' and two groups of psd data(psd-1 and psd-2 in #4). Other testing conditions are all the same. The summary of results are that

1. When using 'cohesiveLaw' with psd-1, there is an error mentioned in #5 but with psd-2 there is nothing wrong.
2. When using 'Ip2_FrictMat_FrictMat_FrictPhys()' with the two psd data. All of them can run completely.

So I doubt that there may be something difficult in the combination of the two considering tasks.

Send my whole running codes with the condition data is not easy in this website, so I copied all the testing codes in this page.
I am sorry for bringing you some troubles in your reading and need your help to check the reason.

Cheers

Dear Fu,
It seems to me than one hour would not be enough to only understand all the above messages and the corresponding scripts.
Add to this, the need to download input data files and rename them as you suggest (why input txt files for python scripts?!), put them in the right place, find which version of this long script is the one we are supposed to try (two complete scripts plus many code fragments in different messages), and you will anticipate that nobody will ever try to reply.

To give people a chance to help you I suggest the following:
- Close this question. Open a new one with a simpler problem statement. "There is something wong" is helpless.
- Don't mix physical considerations (e.g. how shear bands are developping in lab experiments) with technical problems (e.g. "ValueError: cannot convert float NaN to integer", why is this error not mentionned before #5?!).
- Prepare a script as simple as possible in order to show the problem.

Bruno

This question is too confused to read.