yade-users team mailing list archive

[Luc OGER <question693194@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Luc OGER is still having a problem:
Hi Jérôme,
here the full code and the table for the batch test with only the capillray pressure changes : 1000 and 5000
exactly the same particle positions (611 moving) for all of them after 1913000 steps !!
if I am testing type = 3 it is alaso the same thing : no difference??
of course the results are different between the three cases (1 dry, 2 hertz,3 linear) 

So I don't catch why these engine loops are not playing well?

------------------------------------------------ 
table :
model_type packing_fraction pas_box theta_max converg_min cover_pack_fraction init_seed friction ratio Cap_Pres_Imposed
2 .55 .1 80.0 .002 0.3 625 0.5 2. 1000
2 .55 .1 80.0 .002 0.3 625 0.5 2. 5000

and the code :
# avec directory : +'/'+
# 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
#3 .55 .1 80.0 .002 0.3 625 0.5 2.5 5000
readParamsFromTable(model_type=1, packing_fraction=0.6,pas_box=0.1, theta_max=45.0, converg_min=0.00005,cover_pack_fraction=.2,init_seed=10, friction=0.5,ratio=2.5,Cap_Pres_Imposed=1000)
# 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,friction # size ratio between the glued spheres and the moving ones
global cover_pack_fraction ,converg_min, init_seed,nombre_moving # 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
global angle_actif,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
position_init	= numpy.arange(1200*4, dtype=float)
position_init	= position_init.reshape((1200,4))
numpy.zeros_like(position_init)
Cap_Pressure = 0.0          #valeur initiale 
angle			= math.atan(friction)	# value in radians for friction Angle
angle_init		= math.atan(0.10)
angle_actif		= angle_init
# initialisation coordonnees initiales
#newTable("position_init",1200,4) # Create a new table with 5 rows and 3 column
#creating a material (FrictMat):
#frictionAngle(=.5)
# Contact friction angle (in radians). Hint : use ‘radians(degreesValue)’ in python scripts.

id_SphereMat=O.materials.append(FrictMat(young=young_modulus,poisson=poisson_ratio,density=2500,frictionAngle=angle_actif,label="glass_beads"))
SphereMat=O.materials[id_SphereMat]

   
#1 = >hertz model sans capillary forces
if model_type == 1:#hertz model with capillary forces
	converg_min = converg_min
elif model_type == 2: #hertz model with capillary forces
	converg_min = 10.0*converg_min
else:   #for linear model only
	converg_min = 10*converg_min
box_size = 4.0
i_pas = 0 

step0 = 0
step1 = 0
step2 = 0
step_precedent = 0
mask1=0b01
mask2=0b10
mask3=0b11

gravity_y = 9.81*sin( pas_box*i_pas*3.14/180.0)
gravity_z = 9.81*cos( pas_box*i_pas*3.14/180.0)

O.reset

# create rectangular box from facets
O.bodies.append(geom.facetBox((box_size/2.0,box_size/2.0,box_size*0.62),(box_size/2.0,box_size/2.0,box_size*0.62),wallMask=31))

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

rayon = 0.025
height = 5. * rayon
nombre= int(packing_fraction*(box_size*box_size)/(rayon*rayon*3.14159) )
print("nombre = ",nombre)
print("friction angle = ",angle_actif)
print("Model type (1 sec- 2 hertz, 3 lineaire) = ",model_type)
sp=pack.SpherePack()
# generate randomly spheres with uniform radius distribution
sp.makeCloud((0,0,0),(box_size,box_size,height),rMean=rayon,rRelFuzz=.0005,num= nombre,seed=init_seed)

# add the sphere pack to the simulation
sp.toSimulation(color=(1,0,0),mask=mask1,material=SphereMat)


# simulation loop 
#define engines:
if model_type == 1:
	O.engines=[
 		ForceResetter(),
		InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()], avoidSelfInteractionMask = mask2),
		InteractionLoop(
			# handle sphere+sphere and facet+sphere collisions
			[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
			[Ip2_FrictMat_FrictMat_FrictPhys(frictAngle=angle_actif)],
                        [Law2_ScGeom_FrictPhys_CundallStrack()]
#			[Law2_L3Geom_FrictPhys_ElPerfPl()]
		),
		NewtonIntegrator(damping=local_damping ,gravity=(0,-gravity_y,-gravity_z),label='Newton_integrator'),
		# call the checkUnbalanced function (defined below) every 2 seconds
		PyRunner(command='checkUnbalanced()',iterPeriod=1000,label='checker'),
	]
elif model_type == 2: #hertz model with capillary forces
	O.engines=[
		ForceResetter(),
		InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()], avoidSelfInteractionMask = mask2),
		InteractionLoop(
			[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
			[Ip2_FrictMat_FrictMat_MindlinCapillaryPhys(label='ContactModel')],#for hertz model only
			[Law2_ScGeom_MindlinPhys_Mindlin()]                                 #for hertz model only
		),
		Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=Cap_Pressure ),
                #It is also necessary to disable interactions removal by the constitutive law (Law2). 
                #The only combinations of laws supported are currently capillary law + Law2_ScGeom_FrictPhys_CundallStrack
                #and capillary law + Law2_ScGeom_MindlinPhys_Mindlin (and the other variants of Hertz-Mindlin).
		#for hertz model only Value of the capillary pressure Uc defined as Uc=Ugas-Uliquid 10000 initial
		#Capillary pressure uc, to be defined equal to Law2_ScGeom_CapillaryPhys_Capillarity.capillaryPressure. 
		NewtonIntegrator(damping=local_damping ,gravity=(0,-gravity_y,-gravity_z),label='Newton_integrator'),
		# call the checkUnbalanced function (defined below) every 2 seconds
		PyRunner(command='checkUnbalanced()',iterPeriod=1000,label='checker'),	
	]
	ContactModel.betan=viscous_normal
	ContactModel.betas=viscous_shear
	ContactModel.useDamping=True
else:   #for linear model only
	O.engines=[
		ForceResetter(),
		InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()], avoidSelfInteractionMask = mask2),
		InteractionLoop(
			[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
			[Ip2_FrictMat_FrictMat_CapillaryPhys()],	#for linear model only
			[Law2_ScGeom_FrictPhys_CundallStrack()],	#for linear model only
		),
		Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=Cap_Pressure),#for linear model only
		NewtonIntegrator(damping=local_damping,gravity=(0,-gravity_y,-gravity_z),label='Newton_integrator'),
   		# call the checkUnbalanced function (defined below) every 2 seconds
		PyRunner(command='checkUnbalanced()',iterPeriod=1000,label='checker'),
	]


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

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy=True

nom_file=str(model_type)+"_"+str(packing_fraction)+"_"+str(cover_pack_fraction)+"_"+str(ratio)+'_'+str(friction)+'_'+str(pas_box)+'_'+str(converg_min)+'_'+str(init_seed)+'_'+str(Cap_Pres_Imposed)
filename_yade=nom_file+'.yade'
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,friction,theta_max,step0,step1,step2, pas_box,i_pas,step_precedent,rayon,init_seed,ratio,cover_pack_fraction,nom_file,filename_yade,nombre_moving
# print("iteration= ",utils.unbalancedForce(),O.iter,step_precedent)# numéro de l'étape
 if (O.iter)<(5000): return
 if (O.iter-step0)<(5000): return
 if (O.iter-step_precedent)>100000 and step_precedent>0 and step2==1:
  O.save('3rd-step_interrompu_'+filename_yade)
  O.pause()
#########################################################################
# 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,step2)
#########################################################################
 if step0==0:
  step0 = O.iter
  print ("len-O.bodies = ",len(O.bodies))
  for b in O.bodies: 
			#b1.state.vel = Vector3(0,0,0)
			b.state.vel[0] = 0		# stop	
			b.state.vel[1] = 0		# tous les
			b.state.vel[2] = 0		# deplacements translations
			b.state.angVel = Vector3(0,0,0)# stop les rotations 
			#b.angVel[0]=0
			#b.angVel[1]=0
			#b.angVel[2]=0
			b.state.blockedDOFs='xyzXYZ'# réinisialise les accelerations
			b.state.dynamic=False #réinisialise les vitesses
  O.save('1st-step_'+filename_yade)
# save file text mode; beginning of the run
# export.textExt('./output_fin.txt',format='x_y_z_r_attrs',attrs=('b.state.pos.norm()','b.state.pos'),comment='dstN dstV_x dstV_y dstV_z')
#  export.textExt('./output_fin.txt',format='x_y_z_r_attrs',attrs=['b.state.vel[0]','b.state.vel[1]','b.state.vel[2]'],comment='Vx Vy Vz')
  print ("fin de step 1: ",step0)
##################creation deuxieme couche #######################################################
 if step0 ==O.iter and step1==0:
  step1=1
  nombre_moving = int(cover_pack_fraction*(box_size*box_size)/(ratio*rayon*ratio*rayon*3.14159))
  print ("nombre_moving =", nombre_moving)
#generate particles:
  sp=pack.SpherePack()
# generate randomly spheres with uniform radius distribution
  sp.makeCloud((rayon,rayon,2.*rayon),(box_size-rayon,box_size-rayon,(3.2+.1)*rayon*ratio),rMean=rayon*ratio,rRelFuzz=.005*ratio,num= nombre_moving,seed=2*init_seed)
# O.bodies.append([sphere(c,r,material=SphereMat,color=(0,2,0)) for c,r in sp])
# add the sphere pack to the simulation
  sp.toSimulation(color=(0,1,0),mask = mask3,material=SphereMat)#,material=SphereMat
#  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))
 

#########################################################################
 if O.iter>step0 and step1==1:
   step1=O.iter
   step2=1
   step_precedent=O.iter
   O.save('2nd-step_'+filename_yade)
   angle_actif = angle
   O.materials[0].frictionAngle=angle_actif # radians
   # for existing contacts, set contact friction directly
   for i in O.interactions: i.phys.tangensOfFrictionAngle=tan(angle_actif)
   print("friction angle (step2) = ",angle_actif)
   pos_init()
# save file text mode; beginning of the run
# export.textExt('./output_fin.txt',format='x_y_z_r_attrs',attrs=('b.state.pos.norm()','b.state.pos'),comment='dstN dstV_x dstV_y dstV_z')
#   export.textExt('./output_fin-2ndstep.txt',format='x_y_z_r_attrs',attrs=['b.state.vel[0]','b.state.vel[1]','b.state.vel[2]'],comment='Vx Vy Vz')
   print ("fin de step2 : ", step1)
   Cap_Pressure=Cap_Pres_Imposed
#########################################################################
 if O.iter>step1 and step2==1:
  if O.iter<step_precedent+10000 : return
  step_precedent=O.iter
  num_Angle=pas_box*i_pas
  str_Angle= f"{(num_Angle):.3f}"
  pos_moved()
#  filename1='./output_'+nom_file+'_'+str(pas_box*i_pas)+'_'+str(i_pas)+'.txt'
  filename1='./output_'+nom_file+'_'+str_Angle+'_'+str(i_pas)+'.txt'
  filenamevtk='./output_'+nom_file+'_'+str(i_pas)+'.vtk'
  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')
#  export.text('/tmp/test.txt')
# # now open paraview, click "Tools" menu -> "Python Shell", click "Run Script", choose "pv_section.py" from this directiory
# or just run python pv_section.py
# see /home/oger/yade_full_install_local/trunk/examples/test/paraview-spheres-solid-section/pv_section.py
#  text2vtk and text2vtkSection function can be copy-pasted from yade/py/export.py into separate py file to avoid executing yade or to use pure python
  export.text2vtk(filename1,filenamevtk)
  i_pas = i_pas+ 1
  gravity_y = 9.81*sin( pas_box*i_pas*3.14/180.0)
  gravity_z = 9.81*cos( pas_box*i_pas*3.14/180.0)
  Newton_integrator.gravity=Vector3(0,-gravity_y,-gravity_z)
  print ("info : ", (pas_box*i_pas), gravity_y, gravity_z, O.iter)
  ForceResetter()
  #print("iteration= ", O.iter)# numéro de l'étape
  if (pas_box*i_pas)>=theta_max :
   # save output file txt after  end of second sedimentation
####batch :     export.textExt('./output_fin_3rd-step.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')
   O.save('3rd-step_'+filename_yade)
   O.pause()
   

def pos_init():
 global position_init,nombre_mobile,nombre_moving
 nombre_mobile=0
 for b in O.bodies:
  if b.mask==mask3:
   position_init[nombre_mobile][0]=b.state.pos[0]
   position_init[nombre_mobile][1]=b.state.pos[1]
   position_init[nombre_mobile][2]=b.state.pos[2]
   position_init[nombre_mobile][3]=b.shape.radius
#   print( position_init[nombre_mobile][0], position_init[nombre_mobile][1], position_init[nombre_mobile][2], position_init[nombre_mobile][3])
   nombre_mobile = nombre_mobile+1

def pos_moved():
 global position_init,str_Angle,nombre_moving,traitement_file,i_pas,ratio,pas_box
 nombre_mobile = 0
 nombre_moved = 0
 for b in O.bodies:
  if b.mask==mask3:
   dist=((position_init[nombre_mobile][0]-b.state.pos[0])*(position_init[nombre_mobile][0]-b.state.pos[0]) + \
   (position_init[nombre_mobile][1]-b.state.pos[1])*(position_init[nombre_mobile][1]-b.state.pos[1]) + \
   (position_init[nombre_mobile][2]-b.state.pos[2])*(position_init[nombre_mobile][2]-b.state.pos[2]))
   if dist>( position_init[nombre_mobile][3]* position_init[nombre_mobile][3]):
    nombre_moved=nombre_moved + 1
#    print("====>",nombre_mobile,nombre_moved)
   nombre_mobile = nombre_mobile+1
 traitement = open(traitement_file, 'a')
 print (i_pas,nombre_moved,nombre_moving,i_pas*pas_box,nombre_moved/nombre_moving, file=traitement) 
 traitement.close()
# save tmp pour rerun convergence
O.saveTmp()


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()

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