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