Question #701438 on Yade changed:
https://answers.launchpad.net/yade/+question/701438
Status: Answered => Open
Aoxi Zhang is still having a problem:
Hi Karol and Jan,
Thanks for your feedback!
Below are the scripts of the MWE:
#### phase1.py####
from __future__ import division
from yade import pack, plot
import math
import numpy as np
import random
from random import gauss
import timeit
import pickle
num_spheres=5000
targetPorosity = 0.405
compFricDegree = 30
finalFricDegree = 19
rate=-0.01
damp=0.6
stabilityThreshold=0.001
confinement=100e3
mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)
MatWall=O.materials.append(FrictMat(young=2e9,poisson=0.3,frictionAngle=0,density=0,label='walls'))
MatSand =
O.materials.append(CohFrictMat(isCohesive=True,young=2e8,alphaKr=0.15,alphaKtw=0.15,\
poisson=0.3,frictionAngle=radians(30),etaRoll=0.15,etaTwist=0.15,\
density=2650.0,normalCohesion=0, shearCohesion=0,\
momentRotationLaw=True,label='sand'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='sand') for center,rad in sp])
Gl1_Sphere.quality=3
newton=NewtonIntegrator(damping=damp)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="Ip2Coh"),
Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,label="TimeStepper"),
TriaxialStressController(
maxMultiplier=1.4,
finalMaxMultiplier=1.004,
thickness = 0,
stressMask = 7,
internalCompaction=True,
label='triax'),
newton,
]
Gl1_Sphere.stripes=0
triax.goal1=triax.goal2=triax.goal3=-confinement
while 1:
O.run(1000, True)
unb=unbalancedForce()
print 'unbF:',unb,' meanStress:
',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1]
if unb<stabilityThreshold and
abs(-confinement-triax.meanStress)/confinement<0.0001:
break
print "### state 1 Isotropic completed ###"
triax.internalCompaction=False
import sys
while triax.porosity - targetPorosity>0.0001:
compFricDegree = 0.95*compFricDegree
setContactFriction(radians(compFricDegree))
print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
sys.stdout.flush()
O.run(500,1)
print "### state 2 Reach target porosity completed ###"
print "### state 3 - click run to start stress relaxation ###"
setContactFriction(radians(finalFricDegree))
NewtonIntegrator.damping=0.9
O.engines=O.engines+[PyRunner(command='stopStressRelaxation()',
iterPeriod=1,label="StopRelaxation")]
for i in range(6):
O.bodies[i].state.vel = Vector3(0,0,0)
O.bodies[i].state.blockedDOFs='xyzXYZ'
triax.dead=True
def stopStressRelaxation():
print('unb',unbalancedForce())
if unbalancedForce()<1e-7:
O.pause()
NewtonIntegrator.damping=0.0#newton.damping=0
print "stress relaxation finished, damping has been set as 0"
print "Please go to phase2"
O.save('sample_generatedFromPhase1_afterRelax.yade.gz')
#####
##### phase2.py#####
from __future__ import division
from yade import pack, plot
import math
import numpy as np
import random
import pickle
from yade import ymport,export
utils.readParamsFromTable(dt=1e-7)
from yade.params import table
finalFricDegree=19
O.load('sample_generatedFromPhase1_afterRelax.yade.gz')
RunTimeLength=3e-4 # how long for the time of data recording
Odt=table.dt
O.dt=Odt
O.engines[3].dead=True## turn of GlobalStiffnessTimeStepper
# TimeStepper.dead=True## Note, use label cannot successfully close time stepper
StopRelaxation.dead=True## turn of StopRelaxation in previous phase
NewtonIntegrator.damping=0.0
monitoredSand=[1000,2001,3002,4003,5004]
def avgVel(idList):
vel=0.0
avg=0.0
for i in idList:
vel+=O.bodies[i].state.vel[1]
avg=vel/len(idList)
return avg
print("O.dt=",O.dt)
O.engines=O.engines+[PyRunner(command='outputData()',
iterPeriod=1,label="outputData")]
O.engines=O.engines+[PyRunner(command='stopSimulation()',
iterPeriod=1,label="stopSimu")]
def stopSimulation():
print ('Running, O.time,',O.time)
if O.time > RunTimeLength:
O.pause()
print ('Running finished, O.realtime,',O.realtime)
def outputData():
f = open("ResultsOfPhase2_dt{}.txt".format(Odt), "a+")
f.write('{} {} {}\n'.format(O.time,O.dt,avgVel(monitoredSand)))
f.close()
O.resetTime()
O.run()
utils.waitIfBatch()
######
###### phase3.py#####
from yade import pack, plot
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator
import matplotlib.ticker as ticker
import numpy as np
import pandas as pd
font1 = {'family': 'Arial', 'fontweight': 'normal', 'size': 14,
'style':'normal'}
linewidth=2
fig , ax = plt.subplots(figsize=(6.25,3))
bwith = 0.5# here control the width of axis
ax.spines['bottom'].set_linewidth(bwith)
# ax.spines['bottom'].set_color('black')
ax.spines['left'].set_linewidth(bwith)
# ax.spines['left'].set_color('black')
ax.spines['top'].set_linewidth(bwith)
# ax.spines['top'].set_color('black')
ax.spines['right'].set_linewidth(bwith)
# ax.spines['right'].set_color('black')
def plotResults(dt):
f=pd.read_csv("ResultsOfPhase2_dt{}.txt".format(dt),sep=" ",header=None)
f.columns=["Otime","Odt","avgVel"]
Otime,Odt,avgVel=
np.array(f["Otime"]),np.array(f["Odt"]),np.array(f["avgVel"])
ax.plot(Otime,avgVel,label="dt={}".format(dt))
plotResults(1e-7)
plotResults(1e-8)
plotResults(1e-9)
plotResults(5e-9)
# =============================================================================
# Figure setting
#1 set legend
ax.legend(loc=0,frameon=False,prop={'family': 'Arial','size': 11})
ax.set_xlabel('Time (s)',font1)
ax.set_ylabel('Velocity (m/s)',font1)
ax.tick_params(which='both',labelsize=10,direction='in')
plt.xticks(fontproperties = 'Arial', size = 10)
plt.yticks(fontproperties = 'Arial', size = 10)
# plt.ylim(-3e-5,-2e-5)
plt.xlim(0,3e-4)
plt.tight_layout(pad=0.4)
plt.savefig('Fig_compare_dt.jpg', dpi=600)
plt.show()
#####
#### For run phase2.py in batch mode, another txt file could be as below:
dt
1e-7
5e-8
1e-8
5e-9
1e-9
Thanks
Aoxi
--
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@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
