yade-users team mailing list archive

[Aoxi Zhang <question701438@xxxxxxxxxxxxxxxxxxxxx>]

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.