yade-users team mailing list archive

[Luc OGER <question692681@xxxxxxxxxxxxxxxxxxxxx>]

Question #692681 on Yade changed:
https://answers.launchpad.net/yade/+question/692681

    Status: Needs information => Open

Luc OGER gave more information on the question:
here the cmake command

and the py code
#Result indeterminism
#It is naturally expected that running the same simulation several times will give exactly the same results: although the computation is done 
#with finite precision, round-off errors would be deterministically the same at every run. While this is true for single-threaded computation 
#where exact order of all operations is given by the simulation itself, it is not true anymore in multi-threaded computation which is described in detail in later sections.
#The straight-forward manner of parallel processing in explicit DEM is given by the possibility of treating interactions in arbitrary order. 
#Strain and stress is evaluated for each interaction independently, but forces from interactions have to be summed up. If summation order is also arbitrary 
#(in Yade, forces are accumulated for each thread in the order interactions are processed, then summed together), then the results can be slightly different. For instance

#(1/10.)+(1/13.)+(1/17.)=0.23574660633484162
#(1/17.)+(1/13.)+(1/10.)=0.23574660633484165

#As forces generated by interactions are assigned to bodies in quasi-random order, summary force Fi
#on the body can be different between single-threaded and multi-threaded computations, but also between different runs of multi-threaded computation with exactly the same parameters. 
#Exact thread scheduling by the kernel is not predictable since it depends on asynchronous events (hardware interrupts) and other unrelated tasks running on the system; 
#and it is thread scheduling that ultimately determines summation order of force contributions from interactions.


# The components of the batch are:
# 1. table with parameters, one set of parameters per line (ccc.table)
# 2. readParamsFromTable which reads respective line from the parameter file
# 3. the simulation muse be run using yade-batch, not yade
#
# $ yade-batch --job-threads=1 03-oedometric-test.table 03-oedometric-test.py
#

# load parameters from file if run in batch
# default values are used if not run from batch
# gravity deposition in box, showing how to plot and save history of data,
# and how to control the simulation while it is running by calling
# python functions from within the simulation loop


#pas_box theta_max converg_min cover_pack_fraction init_seed friction ratio
#model_type bottom_cover pas_box theta_max converg_min cover_pack_fraction init_seed friction ratio
#.1 55. .002 0.1 210 0.4 4
readParamsFromTable(model_type=1, packing_fraction=0.5,pas_box=0.5, theta_max=5.0, nb_cycles=3,converg_min=0.00005,nb_layers=10,init_seed=10, friction=0.15,ratio=1)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *

# import yade modules that we will use below
from yade import pack, plot, export,math
global ratio,nombre,position_prec,position_init # size ratio between the glued spheres and the moving ones
global converg_min, init_seed,nombre_moving,z_max # coverage percent for moving spheres
global  i_pas, box_size,pas_box,rayon,step0,step1,step2,step_precedent,gravity_y,gravity_z,theta_max,theta_max,nom_file,filename_yade,traitement_file,str_Angle, nb_cycles,numero_cycle,num_Angle


#some parameters passed by batch_table:

####batch : friction		= 0.5
####batch : theta_max       = 30.0
####batch : pas_box = 0.1
####batch : ratio = 3	# size ratio between the glued spheres and the moving ones
####batch : cover_pack_fraction = 0.2 # coverage percent for moving spheres
####batch : init_seed=10
####batch : packing_fraction = 70./100.
####batch : model_type = 1   

#some parameters:
shear_modulus 	= 1e5
poisson_ratio 	= 0.3
young_modulus	= 2*shear_modulus*(1+poisson_ratio)
local_damping 	= 0.01
viscous_normal	= 0.021
viscous_shear	= 0.8*viscous_normal
angle			= math.atan(friction)
# initialisation coordonnees initiales
#newTable("position_init",600,4) # Create a new table with 5 rows and 3 column
#creating a material (FrictMat):
id_SphereMat=O.materials.append(FrictMat(young=young_modulus,poisson=poisson_ratio,density=2500,frictionAngle=angle,label="glass_beads"))
SphereMat=O.materials[id_SphereMat]
 
box_size = 1.0
i_pas = 0 
step0=0
mask1=0b01
num_Angle= pas_box*i_pas
gravity_y = 9.81*sin( num_Angle*3.14/180.0)
gravity_z = 9.81*cos( num_Angle*3.14/180.0)

O.reset

nom_file=str(model_type)+"_"+str(packing_fraction)+"_"+str(theta_max)+"_"+str(nb_cycles)+"_"+str(nb_layers)+"_"+str(ratio)+'_'+str(friction)+'_'+str(pas_box)+'_'+str(converg_min)+'_'+str(init_seed)
filename_yade=nom_file+'.yade'
rayon = 0.025
height = 5.0*rayon*nb_layers
# create rectangular box from facets
O.bodies.append(geom.facetBox((box_size/2.0,box_size/2.0,height),(box_size/2.0,box_size/2.0,height),wallMask=31))

# create empty sphere packing
# sphere packing is not equivalent to particles in simulation, it contains only the pure geometry

nombre= int(packing_fraction*(box_size*box_size*height)/(4.0*rayon*rayon*rayon*3.14159)/3.0 )
print("nombre = ",nombre)
sp=pack.SpherePack()
# generate randomly spheres with uniform radius distribution
sp.makeCloud((0,0,0),(box_size,box_size,height*2),rMean=rayon,rRelFuzz=.0005,seed=init_seed)

# add the sphere pack to the simulation
sp.toSimulation(color=(1,0,0),mask=mask1,material=SphereMat)
nb_particles=len(O.bodies)
O.save('0st-step_'+filename_yade)

# save file text mode; beginning of the run
filename1='./output_'+nom_file+'_'+str(i_pas)+'.txt'
export.textExt(filename1,format='x_y_z_r_attrs',attrs=['b.state.vel[0]','b.state.vel[1]','b.state.vel[2]'],comment='V_x V_y V_z')

# simulation loop 
O.engines=[
 		ForceResetter(),
		InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
		InteractionLoop(
			# handle sphere+sphere and facet+sphere collisions
			[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
			[Ip2_FrictMat_FrictMat_FrictPhys(frictAngle=0.0)],
                        [Law2_ScGeom_FrictPhys_CundallStrack()]
#			[Law2_L3Geom_FrictPhys_ElPerfPl()]
		),
		NewtonIntegrator(damping=0.75,gravity=(0,-gravity_y,-gravity_z),label='Newton_integrator'),
		# call the checkUnbalanced function (defined below) every 2 seconds
		PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
]


# timesteps
O.dt=.5*PWaveTimeStep()


traitement_file = nom_file+'.traitement'
traitement = open(traitement_file, 'w')
traitement.close()
# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
 global converg_min,theta_max,step0,step1,step2, pas_box,i_pas,step_precedent,rayon,init_seed,ratio,nom_file,filename_yade, nb_cycles,numero_cycle,num_Angle
 #print("iteration= ",utils.unbalancedForce(),O.iter,step_precedent)# numéro de l'étape
 if (O.iter)<(5000): return
 if (O.iter-step0)<(5000): return
 #print ("bal iter steps :",utils.unbalancedForce(),O.iter,step0,step1)
#########################################################################
# the rest will be run only if unbalanced is < <converg_min (stabilized packing)
 if utils.unbalancedForce()>converg_min: return 
 #print ("bal iter steps :",utils.unbalancedForce(),O.iter,step0,step1)
#########################################################################
 if step0==0:
  step0 = O.iter
  i_pas=1
  print ("len-O.bodies = ",len(O.bodies))
  O.save('1st-step_'+filename_yade)
# save file text mode; beginning of the run
  filename1='./output_'+nom_file+'_'+str(i_pas)+'.txt'
  export.textExt(filename1,format='x_y_z_r_attrs',attrs=['b.state.vel[0]','b.state.vel[1]','b.state.vel[2]'],comment='V_x V_y V_z')
  print ("fin de step 1: ",step0)
  step1=1
 #  export.textExt('./output_creation_2packs.txt',format='x_y_z_r_attrs',attrs=['b.state.vel[0]','b.state.vel[1]','b.state.vel[2]'],comment='V_x V_y V_z')
  print ("len-O.bodies (2 packs) = ",len(O.bodies))
  O.pause()


O.run()
# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()

and the table for batch command:
model_type packing_fraction pas_box theta_max nb_cycles converg_min nb_layers init_seed friction ratio
1 0.6 0.15 25.0 4  0.00005 10 440 0.15 1
1 0.6 0.15 22.0 4  0.00005 10 440 0.15 1

and the first lines of the output for the two runs:
#format x_y_z_r_attrs
# x y z r V_x V_y V_z
0.297853	0.453123	0.284843	0.0250031	5.65698e-06	-2.8491e-06	6.80765e-06
0.342439	0.740925	0.585506	0.0249955	-5.30826e-05	5.30518e-06	-1.07065e-05
and
#format x_y_z_r_attrs
# x y z r V_x V_y V_z
0.297853	0.453123	0.284843	0.0250031	4.84423e-06	-4.98896e-06	6.7201e-06
0.342438	0.740925	0.585506	0.0249955	-3.45893e-05	9.91756e-06	-2.91766e-06

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.