yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #21997
Re: [Question #688575]: The speed of run is too slow ! What are the ways to run the code faster?
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 <
question688575@xxxxxxxxxxxxxxxxxxxxx> 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 : yade-users@xxxxxxxxxxxxxxxxxxx
> 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.
