Question #688575 on Yade changed: https://answers.launchpad.net/yade/+question/688575
Chareyre proposed the following answer: Use yade.timing to discover where runtime is spent and report it here. That's a first step. Bruno Le ven. 7 févr. 2020 06:29, ehsan benabbas < [email protected]> a écrit : > New question #688575 on Yade: > https://answers.launchpad.net/yade/+question/688575 > > Hello everyone, > > I will be really appreciated it if anyone can help me with this question. > > I am using Ubuntu 18.04, and Yade 2019-08-08.git-775ae74 > > I coded a Triaxial test based on [1] and developed it to store the micro + > macro variables. The problem is it takes sooo long to be run when I define > my main sample for the triaxial test. From last week, I run the code and it > has been still running and I believe it would be running for more 2 weeks > with respect to the finishing criterion. Do you know any way to make the > running faster? because it is not useful in the way it is. > > [1] > https://gitlab.com/yade-dev/trunk/blob/master/examples/triax-tutorial/script-session1.py > > My code: > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > > ###################################################################################################### > ######### TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z > IS THE FRONT AXIS ######### > > ###################################################################################################### > > # > ============================================================================= > # Please copy the file to the Desktop and then run the code. > # - to record the micro in the same iterations as the history() function, > only modify the pyRunner iterPeriod : "PyRunner(iterPeriod=IP, > command='record_micro_data()',label='micro_recorder',dead=True)"; > # - to change the micro iterPeriod, change the "micro_record_iterPeriod" > variable to what you wish; > # - to record the micro at each "ORN" iteration, basically remove the > pyRunner (for the record micro function), and call "record_micro_data()" > just below "export.text(str(O.iter/ORN))" call in the deviatoric while loop. > # > ============================================================================= > > print ('************** START **************') > import numpy as np > import time > import datetime, os > start_time=time.time() > from datetime import datetime > import math > from yade import qt, export, utils > from yade import pack > > ###################################### > ######### DEFINING VARIABLES ######### > > print ('============ DEFINING VARIABLES ============') > > nRead=readParamsFromTable( > num_spheres=20000, > compFricDegree = 29, > key='_triax_', > unknownOk=True > ) > > from yade.params import table > > num_spheres=table.num_spheres > key=table.key > targetPorosity = 0.4 > compFricDegree = table.compFricDegree > finalFricDegree = 29 > IP=100 # iteration period to record data and stuff > micro_record_iterPeriod=IP > ORN=3000 # O.Run Number of iterations > micro_record_enable_normal_branch=True > micro_record_float_type = np.float32 > damp=0.2 > thick=0 > stabilityThreshold=0.01 > PCPC=0.0001 # Precision of Confining Pressure Convergence > r_min=0.1*1e-3 # m > d_min=2*r_min # m > r_max=0.3*1e-3 # m > d_max=2*r_max # m > r_avr=(r_min+r_max)/2 # m > d_avr=2*r_avr # m > r_fuz=(r_max/r_avr)-1 # m > Kn=10e8*(d_avr) ### FIXME > Kt=10e8*(d_avr) ### FIXME > young=Kn/r_avr # 2 (E r1 E r2 / E r1 + E r2) >>> E = Kn/r_avr > #young=5e6 # contact stiffness > poisson=Kn/Kt # Kt/Kn > #poisson=0.3 # Kt/Kn > Ls=0.02 # m length of specimen ### FIXME > L_REV=7*(d_avr) # m > if Ls < L_REV: > sys.exit("*** ERROR! The specimen's dimension is too samll! ***") > elif Ls==L_REV: > print ("*** This is the minimum specimen's dimension you can take! > ***") > else: > print ("*** The specimen's dimension is good enough! ***") > mn,mx=Vector3(0,0,0),Vector3(Ls,Ls,Ls) > Vt=-1*1e-3 # m/s # negative sign describes the compression direction > strainRate=Vt/Ls # %/sec > #strainRate=-0.01 # %/sec ### FIXME > target_strain=0.25 ### FIXME % > print ("The target strain has been set to:", target_strain) > sigmaIso=-5e5 # Pa ### FIXME > particleDensity=2000 #kg/m3 > #particleDensity=2600 > > ###################################################### > ################# DEFINING MATERIALS ################# > > print ('============ DEFINING MATERIALS ============') > > O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=particleDensity,label='spheres')) > > O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls')) > > #################################################### > ################# DEFINING PACKING ################# > > print ('============ DEFINING PACKING ============') > walls=aabbWalls([mn,mx],thickness=thick,material='walls') > for w in walls:w.shape.radius=0 > wallIds=O.bodies.append(walls) > sp=pack.SpherePack() > clumps=False > sp.makeCloud(mn,mx,r_avr,r_fuz,num_spheres,False, 0.95,seed=1) > #sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make > the "random" generation always the same > O.bodies.append([sphere(center,rad,material='spheres') for center,rad in > sp]) > > from yade import export > os.mkdir('3axresults') > os.chdir('/home/ehsan/Desktop/3axresults') > export.text('InitialPackingData') > > ########################################################## > ################# DEFINING TRIAXIAL TEST ################# > > print ('============ DEFINING TRIAXIAL TEST ============') > triax=TriaxialStressController( > maxMultiplier=1.+2e4/young, > finalMaxMultiplier=1.+2e3/young, > thickness = thick, > stressMask = 7, > internalCompaction=True, > ) > > ###################################################### > ################# DEFINING FUNCTIONS ################# > > print ('============ DEFINING FUNCTIONS ============') > from yade import plot > def history(): > plot.addData( > e11 = -triax.strain[0], > e22 = -triax.strain[1], > e33 = -triax.strain[2], > ev = -triax.strain[0]-triax.strain[1]-triax.strain[2], > s11 = -triax.stress(triax.wall_right_id)[0], > s22 = -triax.stress(triax.wall_top_id)[1], > s33 = -triax.stress(triax.wall_front_id)[2], > i = O.iter, > t = O.time, # virtual (yade) time --- time of simulation > fab = utils.fabricTensor()[0]) > > #FK: Micro data recording > #NOTE: Data structure in hdf5 is: > # hdf5-file/ > # ├── bodies/ > # | ├── positions(3) > # | ├── radius(1) > # └── interactions/ > # ├── ids(2) > # ├── normalForce(3) > # └── shearForce(3) > #IMPORTANT NOTE: the array indexation in bodies/positions and > bodies/radius is the same as in interactions/ids : it is the Yade > indexation (O.bodies[my_index]) > #It means that the values in interactions/ids hdf5 array can be used as > index in bodies/positions and bodies/radius. > > import h5py # works with ORN and macro period of recording > def record_micro_data(name): > if not os.path.exists("./"+key+"micro/"): os.mkdir("./"+key+"micro/") > filename="./"+key+"micro/"+name > with h5py.File(filename+".h5", "w") as f: > #Record time: > f.create_dataset("time",data=O.time) > #Create a "folder" for bodies: > f.create_group("bodies") > #Add bodies positions > f.create_dataset("bodies/positions",data = np.asarray([b.state.pos > for b in O.bodies]).astype(micro_record_float_type) ) > #Add bodies radius. Radius=0 means its a wall > f.create_dataset("bodies/radius",data = np.asarray([b.shape.radius > for b in O.bodies]).astype(micro_record_float_type) ) > #Create a "folder" for interactions: > f.create_group("interactions") > #Prepare the output array, performance is way better when you set > its size a-priori. > out_array = > np.empty([O.interactions.countReal(),8],dtype=micro_record_float_type) > #Fill the output_array > for i,I in enumerate(O.interactions): > out_array[i,0]=I.id1 > out_array[i,1]=I.id2 > out_array[i,2:5]=I.phys.normalForce > out_array[i,5:8]=I.phys.shearForce > f.create_dataset("interactions/ids",data = > out_array[:,0:2].astype(int) ) > f.create_dataset("interactions/normalForce",data = > out_array[:,2:5] ) > f.create_dataset("interactions/shearForce",data = out_array[:,5:8] > ) > if micro_record_enable_normal_branch: > with h5py.File(filename+"_branch_normal.h5", "w") as f: > out_array = > np.empty([O.interactions.countReal(),8],dtype=micro_record_float_type) > #prepare the output array, performance is way better when you set its size > a-priori. > for i,I in enumerate(O.interactions): > out_array[i,0]=I.id1 > out_array[i,1]=I.id2 > out_array[i,2:5]=I.geom.normal > > out_array[i,5:8]=O.bodies[I.id2].state.pos-O.bodies[I.id1].state.pos > #branch vector computation > f.create_dataset("ids",data = out_array[:,0:2].astype(int)) > f.create_dataset("normal",data = out_array[:,2:5]) > f.create_dataset("branch",data = out_array[:,5:8]) > > #################################################### > ################# DEFINING ENGINES ################# > > print ('============ DEFINING ENGINES ============') > newton=NewtonIntegrator(damping=damp) > > O.engines=[ > ForceResetter(), > #GravityEngine(gravity=(0,-9.806,0),warnOnce=False), > InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), > InteractionLoop( > [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], > [Ip2_FrictMat_FrictMat_FrictPhys()], > [Law2_ScGeom_FrictPhys_CundallStrack()] > ), > ## We will use the global stiffness of each body to determine an > optimal timestep (see > https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf) > > GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), > triax, > PyRunner(iterPeriod=IP,command='history()',label='macro_recorder'), > > #PyRunner(iterPeriod=micro_record_iterPeriod,command='record_micro_data()',label='micro_recorder',dead=True), > TriaxialStateRecorder(iterPeriod=IP,file='WallStresses'+table.key), > newton > ] > > Gl1_Sphere.stripes=True > if nRead==0: yade.qt.Controller(), yade.qt.View() > > ############################################################### > ################# APPLYING CONFINING PRESSURE ################# > > print ('============ APPLYING CONFINING PRESSURE ============') > > triax.goal1 = sigmaIso > triax.goal2 = sigmaIso > triax.goal3 = sigmaIso > > while 1: > O.run(ORN,True) > unb = unbalancedForce() > > meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3 > ConfStressRatio=abs(sigmaIso-triax.meanStress)/abs(sigmaIso) > print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, > 50kPa ~~~~~~~~~~~~~') > print ('mean stress engine:',triax.meanStress,' mean stress > (Calculated):',meanS, ' ConfSratio:',ConfStressRatio,' step:', O.iter/ORN, > ' Time:',O.time, ' TimeStep',O.dt) > print ('porosity:',triax.porosity, ' void > ratio:',triax.porosity/(1-triax.porosity)) > if unb<stabilityThreshold and ConfStressRatio<PCPC: > break > > export.text('FinalPhase01PackingData') > e22Check=-triax.strain[1] ###%%%*** > print ('Axial Strain',e22Check) > print ('Mean stress engine: ',triax.meanStress) > print ('Mean stress (Calculated): ',meanS) > print ('################## Isotropic phase is finished and saved > successfully ##################') > > ############################################################# > ################# REACHING TARGET POROSITY ################## > > print ('============ REACHING TARGET POROSITY ============') > import sys > while triax.porosity>targetPorosity: > compFricDegree = 0.95*compFricDegree > setContactFriction(radians(compFricDegree)) > print ('\r Friction: ',compFricDegree,' > porosity:',triax.porosity,'step= ',O.iter/ORN,' Time:',O.time, ' > TimeStep',O.dt) > sys.stdout.flush() > O.run(ORN,True) > > print ('################## Target porosity is reached and compacted state > saved successfully ##################') > > ############################################################# > ################ PRINT SOME CHECK VARIABLES ################# > > RRmax=max([b.shape.radius for b in O.bodies]) > RRmin=min([b.shape.radius for b in O.bodies]) > print('Maximum Radius:',RRmax) > print('Minimum Radius:',RRmin) > print ('Number of elements:', len(O.bodies)) > print ('Box Volume engine:', triax.boxVolume) > Vt = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2]) # total volume of the > specimen (box) > print ('Box Volume calculated:', Vt) > if Vt == triax.boxVolume: > print ("*** Volume calculation is Correct. ***") > else: > sys.exit("*** ERROR! Volume calculation is WRONG. ***") > Vs=triax.particlesVolume > print('Total volume of particles (Vs):',Vs) > Vv=Vt-Vs > print('Total volume of voids Calculated (Vv):',Vv) > print ('porosity:',triax.porosity) > n=Vv/Vt > print ('porosity Calculated (n):',n) > if n == triax.porosity: > print ("*** Porosity calculation is Correct. ***") > e=n/(1-n) > print ('Void ratio Calculated (e):',e) > else: > sys.exit("*** ERROR! Porosity calculation is WRONG. ***") > print ('step that starts the deviatoric loading ', O.iter/ORN) > > ###################################################### > ################# DEVIATORIC LOADING ################# > > print ('============ APPLYING DEVIATORIC LOADING ============') > triax.internalCompaction=False > setContactFriction(radians(finalFricDegree)) > > # Allow walls to move. >>> False means fixed , True means free to move > triax.wall_top_activated=True > triax.wall_bottom_activated=False > triax.wall_left_activated=True > triax.wall_right_activated=True > triax.wall_back_activated=True > triax.wall_front_activated=True > > # Make Stress Control on X,Y,Z directions by setting True, False means > Strain Control > triax.stressControl_1=True > triax.stressControl_2=False > triax.stressControl_3=True > > # set the raite of strain to the desirable axis > #triax.strainRate1=100*rate > triax.strainRate2=strainRate > #triax.strainRate3=100*rate > > # set the deviatoric conditions > triax.stressMask = 5 > triax.goal1 = sigmaIso > triax.goal2 = strainRate > triax.goal3 = sigmaIso > > newton.damping=0.1 ### FIXME > os.mkdir('Phase02PackingData') > os.chdir('/home/ehsan/Desktop/3axresults/Phase02PackingData') > > #micro_recorder.dead=False > > while 1: > O.run(ORN,True) > export.text(str(O.iter/ORN)) > record_micro_data(str(O.iter/ORN)) > unb=unbalancedForce() > axialS=triax.stress(triax.wall_top_id)[1] > eps2=triax.strain[1] > eps2cal=(triax.height-Ls)/Ls > print ('~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, > Strain Rate ~~~~~~~~~~~~~') > print ('sigma2: ',axialS, ' q: ', axialS-sigmaIso,' step= ', > O.iter/ORN,' Time:',O.time, ' TimeStep',O.dt) > print ('Strain2 Calculated: ',eps2cal, ' Strain2 engine: ',eps2,'axial > deformation (%)', (eps2-e22Check)*100) > print('A:', abs(eps2-e22Check), ' B:',target_strain) > if abs(eps2-e22Check)>=target_strain: #triax.sigma2 > break > > > print ('################## Deviatoric phase is finished and saved > successfully ##################') > > ####################################### > ################# END ################# > > print ('************** END **************') > O.realtime > print ('Analysis has been taken for',O.realtime, 'seconds or', > O.realtime/60, 'minutes') > print('Real time of run:',(time.time() - start_time), 'seconds > or',(time.time() - start_time)/60, 'minutes') > > ############################################################## > ################# RECORD Macro DATA and Plot ################# > > print ('============ RECORD AND PLOT DATA ============') > > O.run(ORN,True) # Show the deformed shape of specimen, otherwise the > specimen shows no deformation schematically > # > ============================================================================= > # plot.plots={'e22':('ev',),'j':('s22',),'e22 ':('s22',),' > j':('s11','s33',)} > # plot.labels={'s11':'$\sigma_{11}$' , 's22':'$\sigma_{22}$' , > 's33':'$\sigma_{33}$' , 'e22':'$\epsilon_{22}$' , 'ev':'$\epsilon_{V}$' , > 'i':'Time Step'} > # plot.plot() > # > ============================================================================= > os.chdir('/home/ehsan/Desktop/3axresults') > plot.saveDataTxt('Macro_results') > > from matplotlib import pyplot as plt > > e22 = plot.data["e22"] > ev = plot.data["ev"] > s22 = plot.data["s22"] > s11 = plot.data["s11"] > s33 = plot.data["s33"] > i = plot.data["i"] > t = plot.data["t"] > > fig1, ax1=plt.subplots(figsize=(15, 10)) > fig2, ax2=plt.subplots(figsize=(15, 10)) > fig3, ax3=plt.subplots(figsize=(15, 10)) > fig4, ax4=plt.subplots(figsize=(15, 10)) > > ax1.plot(e22,ev,color='k',linewidth=1,label='Volumetric Strain') > ax1.legend() > #ax1.set_title('sample graph') > ax1.set_xlabel('Axial Strain') > ax1.set_ylabel('Volumetric Strain') > ax1.grid(True) > fig1.tight_layout() > fig1.savefig('plot1.eps', format='eps') > > ax2.plot(t,s22,color='k',linewidth=1,label='Axial Stress') > ax2.legend() > #ax2.set_title('sample graph') > ax2.set_xlabel('Time') > ax2.set_ylabel('Axial Stress') > ax2.grid(True) > fig2.tight_layout() > fig2.savefig('plot2.eps', format='eps') > > ax3.plot(e22,s22,color='k',linewidth=1,label='Stress-Strain') > ax3.legend() > #ax3.set_title('sample graph') > ax3.set_xlabel('Axial Strain') > ax3.set_ylabel('Axial Stress') > ax3.grid(True) > fig3.tight_layout() > fig3.savefig('plot3.eps', format='eps') > > ax4.plot(t,s11,color='b',linewidth=1,label='Confining Stress - 1') > ax4.plot(t,s33,color='r',linewidth=1,label='Confining Stress - 3') > ax4.legend() > #ax4.set_title('sample graph') > ax4.set_xlabel('Time') > ax4.set_ylabel('Confining Pressure') > ax4.grid(True) > fig4.tight_layout() > fig4.savefig('plot4.eps', format='eps') > > plt.show() > > O.pause() > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > Thank you, > Ehsan > > -- > 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 : [email protected] > Unsubscribe : https://launchpad.net/~yade-users > More help : https://help.launchpad.net/ListHelp > -- 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 : [email protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
