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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp