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
