Yade

The speed of run is too slow ! What are the ways to run the code faster?

Asked by ehsan benabbas

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

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
ehsan benabbas
Solved:
Last query:
Last reply:
Revision history for this message
Robert Caulk (rcaulk) said : 		#1

Sure, use the -j argument while launching Yade [1]. Fun fact yade -h lets you see the available arguments when launching yade.

[1]https://yade-dem.org/doc/introduction.html#starting-yade

Revision history for this message
ehsan benabbas (ehsanben) said : 		#2

Thank you for the answer. By doing -j4 the iter/sec will increase.
But, when I use the following command:
cat /proc/cpuinfo | grep processor -c
I get "8" as an output which means my PC has 8 processor. (is that right?)
At the same time, if I use the follwoing command:
cat /proc/cpuinfo
This will be my output

processor : 0
vendor_id : GenuineIntel
cpu family : 6
model : 60
model name : Intel(R) Core(TM) i7-4770 CPU @ 3.40GHz
stepping : 3
microcode : 0x27
cpu MHz : 3780.467
cache size : 8192 KB
physical id : 0
siblings : 8
core id : 0
cpu cores : 4
apicid : 0
initial apicid : 0
fpu : yes
fpu_exception : yes
cpuid level : 13
wp : yes
flags : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm cpuid_fault epb invpcid_single pti ssbd ibrs ibpb stibp tpr_shadow vnmi flexpriority ept vpid ept_ad fsgsbase tsc_adjust bmi1 avx2 smep bmi2 erms invpcid xsaveopt dtherm ida arat pln pts md_clear flush_l1d
bugs : cpu_meltdown spectre_v1 spectre_v2 spec_store_bypass l1tf mds swapgs itlb_multihit
bogomips : 6784.72
clflush size : 64
cache_alignment : 64
address sizes : 39 bits physical, 48 bits virtual
power management:

Which means my PC has 4 cores. Watching the "System Monitor" of my Ubuntu, I see 8 CUPs (CUP1, CUP2,...,CUP8). So, it means that these 4 cores have 8 processors overall (is that correct?)

When I run my code by the following command
/home/ehsan/yade/install/bin/yade-2019-08-08.git-775ae74 -j4 packing.py
Only 4 CUPs of those 8 CUPs in "System Monitor" will be engaged and the other 4 will not.

When I run my code with the following command
/home/ehsan/yade/install/bin/yade-2019-08-08.git-775ae74 -j8 packing.py
I get lower rate for iter/sec in comparison with -j4 while in "System Monitor" all 8 CPUs are engaged.

Revision history for this message
ehsan benabbas (ehsanben) said : 		#3

Additional information about my PC and its CPU:

ehsan@ehsan:~$ sudo dmidecode --type processor
[sudo] password for ehsan:
# dmidecode 3.1
Getting SMBIOS data from sysfs.
SMBIOS 2.7 present.

Handle 0x0013, DMI type 4, 42 bytes
Processor Information
 Socket Designation: SOCKET 0
 Type: Central Processor
 Family: Core i7
 Manufacturer: Intel
 ID: C3 06 03 00 FF FB EB BF
 Signature: Type 0, Family 6, Model 60, Stepping 3
 Flags:
  FPU (Floating-point unit on-chip)
  VME (Virtual mode extension)
  DE (Debugging extension)
  PSE (Page size extension)
  TSC (Time stamp counter)
  MSR (Model specific registers)
  PAE (Physical address extension)
  MCE (Machine check exception)
  CX8 (CMPXCHG8 instruction supported)
  APIC (On-chip APIC hardware supported)
  SEP (Fast system call)
  MTRR (Memory type range registers)
  PGE (Page global enable)
  MCA (Machine check architecture)
  CMOV (Conditional move instruction supported)
  PAT (Page attribute table)
  PSE-36 (36-bit page size extension)
  CLFSH (CLFLUSH instruction supported)
  DS (Debug store)
  ACPI (ACPI supported)
  MMX (MMX technology supported)
  FXSR (FXSAVE and FXSTOR instructions supported)
  SSE (Streaming SIMD extensions)
  SSE2 (Streaming SIMD extensions 2)
  SS (Self-snoop)
  HTT (Multi-threading)
  TM (Thermal monitor supported)
  PBE (Pending break enabled)
 Version: Intel(R) Core(TM) i7-4770 CPU @ 3.40GHz
 Voltage: 1.2 V
 External Clock: 100 MHz
 Max Speed: 3400 MHz
 Current Speed: 3400 MHz
 Status: Populated, Enabled
 Upgrade: Socket BGA1155
 L1 Cache Handle: 0x0015
 L2 Cache Handle: 0x0014
 L3 Cache Handle: 0x0016
 Serial Number: Not Specified
 Asset Tag: Not Specified
 Part Number: Not Specified
 Core Count: 4
 Core Enabled: 4
 Thread Count: 8
 Characteristics:
  64-bit capable

Revision history for this message
Robert Caulk (rcaulk) said : 		#4

You have 4 cores 8 threads. Threads!=cores.

Revision history for this message
Chareyre (bruno-chareyre-9) said : 		#5

Hi,
Hardware is factor. Yade implementation is another. Your script is another.

The three ways to run faster are to improve one of the three factors. It's
hard to say more based on your question.
Cheers
Bruno

Le ven. 7 févr. 2020 22:38, ehsan benabbas <
<email address hidden>> a écrit :

> Question #688575 on Yade changed:
> https://answers.launchpad.net/yade/+question/688575
>
> ehsan benabbas gave more information on the question:
> Additional information about my PC and its CPU:
>
> ehsan@ehsan:~$ sudo dmidecode --type processor
> [sudo] password for ehsan:
> # dmidecode 3.1
> Getting SMBIOS data from sysfs.
> SMBIOS 2.7 present.
>
> Handle 0x0013, DMI type 4, 42 bytes
> Processor Information
> Socket Designation: SOCKET 0
> Type: Central Processor
> Family: Core i7
> Manufacturer: Intel
> ID: C3 06 03 00 FF FB EB BF
> Signature: Type 0, Family 6, Model 60, Stepping 3
> Flags:
> FPU (Floating-point unit on-chip)
> VME (Virtual mode extension)
> DE (Debugging extension)
> PSE (Page size extension)
> TSC (Time stamp counter)
> MSR (Model specific registers)
> PAE (Physical address extension)
> MCE (Machine check exception)
> CX8 (CMPXCHG8 instruction supported)
> APIC (On-chip APIC hardware supported)
> SEP (Fast system call)
> MTRR (Memory type range registers)
> PGE (Page global enable)
> MCA (Machine check architecture)
> CMOV (Conditional move instruction supported)
> PAT (Page attribute table)
> PSE-36 (36-bit page size extension)
> CLFSH (CLFLUSH instruction supported)
> DS (Debug store)
> ACPI (ACPI supported)
> MMX (MMX technology supported)
> FXSR (FXSAVE and FXSTOR instructions supported)
> SSE (Streaming SIMD extensions)
> SSE2 (Streaming SIMD extensions 2)
> SS (Self-snoop)
> HTT (Multi-threading)
> TM (Thermal monitor supported)
> PBE (Pending break enabled)
> Version: Intel(R) Core(TM) i7-4770 CPU @ 3.40GHz
> Voltage: 1.2 V
> External Clock: 100 MHz
> Max Speed: 3400 MHz
> Current Speed: 3400 MHz
> Status: Populated, Enabled
> Upgrade: Socket BGA1155
> L1 Cache Handle: 0x0015
> L2 Cache Handle: 0x0014
> L3 Cache Handle: 0x0016
> Serial Number: Not Specified
> Asset Tag: Not Specified
> Part Number: Not Specified
> Core Count: 4
> Core Enabled: 4
> Thread Count: 8
> Characteristics:
> 64-bit capable
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>
>
>

Revision history for this message
Chareyre (bruno-chareyre-9) said : 		#6

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 address hidden>> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
ehsan benabbas (ehsanben) said : 		#7

Thank you Bruno for you answer and a special thanks for [1]. As I already run the original code (the one that I postet above with 20000 particles) and my PC was engaging with that, I just defined I new one with the parameters as [1] and run this code while the previous one was running at the same time.

[1] https://gitlab.com/yade-dev/trunk/blob/master/examples/triax-tutorial/script-session1.py

Name Count Time Rel. time
-------------------------------------------------------------------------------------------------------
ForceResetter 71000 2434809.955us 0.55%
InsertionSortCollider 48793 302790858.018us 68.68%
InteractionLoop 71000 65393251.802us 14.83%
GlobalStiffnessTimeStepper 890 195757.733us 0.04%
TriaxialStressController 71000 6225534.471us 1.41%
"macro_recorder" 709 366551.038us 0.08%
TriaxialStateRecorder 710 303926.099us 0.07%
NewtonIntegrator 71000 63175600.748us 14.33%
TOTAL 440886289.864us 100.00%

Base on this table, "InsertionSortCollider" takes the most time with almost 69%, then "InteractionLoop" and "NewtonIntegrator" each of those with almost 15%.

I will be really appreciated if you let me know your opinion. Is this related to my hardware, Yade implementation, or Your script? or a combination of those. What do I do to improve the running speed?

Thank you,
Ehsan

Revision history for this message
ehsan benabbas (ehsanben) said : 		#8

any comment with the new information I provided?

Revision history for this message
Robert Caulk (rcaulk) said : 		#9

Using -j4 is the fastest you will get with a 4 core CPU.

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#10

There is no good reason for the collider to run so often (48793 / 71000), and then you are waisting most of the cpu time therein.
I would guess the default verletDist goes wrong for some reason (you can assign it explicitely to solve the problem), or your scene is completely unstable and some particles keep jumping around at high speed.

In the end, most of the example you posted is useless since post-processing is not significant, and your hardware is irrelevant too.
Please post minimal scripts if you want someone to test them.

Cheers

Bruno

Revision history for this message
ehsan benabbas (ehsanben) said : 		#11

Minimal scripts:

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
poisson=Kn/Kt # 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
target_strain=0.25 ### FIXME %
print ("The target strain has been set to:", target_strain)
sigmaIso=-5e5 # Pa ### FIXME
particleDensity=2000 #kg/m3

##################################################
################# 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)
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])

################################################
################# DEFINING ENGINES #################

print ('============ DEFINING ENGINES ============')
newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
    PyRunner(iterPeriod=IP,command='history()',label='macro_recorder'),
 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 ##################')

##################################################
################# DEVIATORIC LOADING #################

print ('============ APPLYING DEVIATORIC LOADING ============')
triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))

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

triax.stressControl_1=True
triax.stressControl_2=False
triax.stressControl_3=True

triax.strainRate2=strainRate

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')

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)

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')

O.pause()

Revision history for this message
Chareyre (bruno-chareyre-9) said : 		#12

> os.chdir('/home/ehsan/Desktop/3axresults')

Seriously?
B

Le lun. 10 févr. 2020 19:53, ehsan benabbas <
<email address hidden>> a écrit :

> Question #688575 on Yade changed:
> https://answers.launchpad.net/yade/+question/688575
>
> Status: Answered => Open
>
> ehsan benabbas is still having a problem:
> Minimal scripts:
>
>
> 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
> poisson=Kn/Kt # 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
> target_strain=0.25 ### FIXME %
> print ("The target strain has been set to:", target_strain)
> sigmaIso=-5e5 # Pa ### FIXME
> particleDensity=2000 #kg/m3
>
> ##################################################
> ################# 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)
> 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])
>
> ################################################
> ################# DEFINING ENGINES #################
>
> print ('============ DEFINING ENGINES ============')
> newton=NewtonIntegrator(damping=damp)
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_FrictPhys()],
> [Law2_ScGeom_FrictPhys_CundallStrack()]
> ),
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
> triax,
> PyRunner(iterPeriod=IP,command='history()',label='macro_recorder'),
> 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 ##################')
>
> ##################################################
> ################# DEVIATORIC LOADING #################
>
> print ('============ APPLYING DEVIATORIC LOADING ============')
> triax.internalCompaction=False
> setContactFriction(radians(finalFricDegree))
>
> 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
>
> triax.stressControl_1=True
> triax.stressControl_2=False
> triax.stressControl_3=True
>
> triax.strainRate2=strainRate
>
> 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')
>
> 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)
>
> 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')
>
> O.pause()
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>
>
>

Revision history for this message
ehsan benabbas (ehsanben) said : 		#13

I realized the most important factor in the speed of running is the "young" parameter. In the way I define this parameter (or change it to an arbitrary value) it influences the run time significantly. Knowing this, I changed "young" to the value you defined already in [1] and the code works really fast for 20000 particles (about 15 minutes). However, there is another problem. As I want to define Kn and Ks directly, I need to change FrictMat to ViscoElMat. But as this is another question, I will post a new question that anyone can use. Thank you for your helps.