Commit e5b994fb authored by François Kneib's avatar François Kneib Committed by Janek Kozicki

[Python3] Remove the "old_div" from the futurize --stage2 as it is buggy.

parent c77d408f
......@@ -3,7 +3,6 @@
# syntax:python
from __future__ import print_function
from past.builtins import execfile
from builtins import str
import sys,os,os.path,time
......
from __future__ import print_function
#!${pyExecutable}
# encoding: utf-8
#
# vim: syntax=python
# portions © 2008 Václav Šmilauer <eudoxos@arcig.cz>
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from builtins import str
from builtins import range
from builtins import object
from past.utils import old_div
import os, sys, _thread, time, logging, pipes, socket, xmlrpc.client, re, shutil, random
# Add search path for yade Python-modules
......@@ -115,8 +112,8 @@ finished: %s
if info:
ret+='<td>'
if info['stopAtIter']>0:
ret+='<nobr>%2.2f%% done</nobr><br/><nobr>step %d/%d</nobr>'%(old_div(info['iter']*100.,info['stopAtIter']),info['iter'],info['stopAtIter'])
finishTime = str(time.ctime(time.time()+int((old_div(round(info['stopAtIter'] - info['iter']),info['speed'])))))
ret+='<nobr>%2.2f%% done</nobr><br/><nobr>step %d/%d</nobr>'%(info['iter']*100./info['stopAtIter'],info['iter'],info['stopAtIter'])
finishTime = str(time.ctime(time.time()+int((round(info['stopAtIter'] - info['iter'])/info['speed']))))
ret+='<br/><font size="1"><nobr>%s finishes</nobr></font><br/>'%finishTime
else: ret+='<nobr>step %d</nobr>'%(info['iter'])
if info['realtime']!=0: ret+='<br/><nobr>speed %g/sec</nobr>'%(info['speed'])
......@@ -566,10 +563,10 @@ if opts.timing>0:
import math
for i,l in enumerate(useLines):
jobTimes=[j.durationSec for j in jobs if j.lineNo==l and j.durationSec!=None]
tSum=sum(jobTimes); tAvg=old_div(tSum,len(jobTimes))
tSum=sum(jobTimes); tAvg=tSum/len(jobTimes)
tMin,tMax=min(jobTimes),max(jobTimes)
tDev=math.sqrt(old_div(sum((t-tAvg)**2 for t in jobTimes),len(jobTimes)))
tRelDev=old_div(tDev,tAvg)
tDev=math.sqrt(sum((t-tAvg)**2 for t in jobTimes)/len(jobTimes))
tRelDev=tDev/tAvg
out.write('%d\t%d\t%.2f\t%.2f\t%.3g\t%.2f\t%.2f\t|\t'%(l,len(jobTimes),tAvg,tDev,tRelDev,tMin,tMax)+'\t'.join([params[l][p] for p in paramNames])+'\n')
if not gnuplotOut:
......
from __future__ import print_function
#!${pyExecutable}
# encoding: utf-8
#
......@@ -7,8 +8,6 @@
# This script is to be used with OAR task scheduler. May be an example to use use with other task scheduler for clusters
# Adapted from yade-batch
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import str
......
......@@ -14,8 +14,6 @@
# setup the periodic boundary
from __future__ import print_function
from __future__ import division
from past.utils import old_div
O.periodic=True
O.cell.refSize=(2,2,2)
......@@ -111,7 +109,7 @@ def addData():
# get the stress tensor (as 3x3 matrix)
stress=sum(normalShearStressTensors(),Matrix3.Zero)
# give names to values we are interested in and save them
plot.addData(exz=O.cell.trsf[0,2],szz=stress[2,2],sxz=stress[0,2],tanPhi=old_div(stress[0,2],stress[2,2]),i=O.iter)
plot.addData(exz=O.cell.trsf[0,2],szz=stress[2,2],sxz=stress[0,2],tanPhi=stress[0,2]/stress[2,2],i=O.iter)
# color particles based on rotation amount
for b in O.bodies:
# rot() gives rotation vector between reference and current position
......
......@@ -3,10 +3,8 @@
#
# module documentation
#
from __future__ import division
from builtins import str
from builtins import range
from past.utils import old_div
import sys,os,os.path
reload(sys)
sys.setdefaultencoding('utf8')
......@@ -136,7 +134,7 @@ def inheritanceDiagram(klass,willBeLater):
# margin size is in inches. The text area on page in .pdf is 6.3in by 9.8in. I'll use a default that each class uses one fourth of page width. If maxDepth>=5 then the image just gets smaller.
pageWidth=6.3
pageFraction=4
fixPdfMargin=(old_div(pageWidth,pageFraction))*max(0,pageFraction-maxDepth)
fixPdfMargin=(pageWidth/pageFraction)*max(0,pageFraction-maxDepth)
ret=""
extraCaption=["",0]
extraPdfCaptionSet=set()
......
from __future__ import division
#!/usr/local/bin/yade-trunk -x
# -*- coding: utf-8 -*-
#Author : Kneib François, francois.kneib@gmail.com
......@@ -6,7 +5,6 @@ from __future__ import division
#/!\ this is just a DISPLAY FEATURE, computed particles still are SPHERICAL.
from builtins import range
from past.utils import old_div
from yade import qt
X=1
Y=1
......@@ -41,7 +39,7 @@ for i in O.bodies:
i.shape.color=(0.,0.,1.)
#Fix some spheres to make the simulation "interesting".
for i in range(0,int(old_div(len(O.bodies),10))):
for i in range(0,int(len(O.bodies)/10)):
O.bodies[i].state.blockedDOFs="xyzXYZ"
O.bodies[i].shape.color=(1.,0.,0.)
......
......@@ -7,8 +7,6 @@
from __future__ import print_function
from __future__ import division
from past.utils import old_div
from yade import utils, plot, timing
from yade import pack
......@@ -31,9 +29,9 @@ o.dt = 1.0e-5
particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho
Vi1 = math.sqrt(old_div(k1,particleMass))*DeltaPMax*Chi1
Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1 = old_div(DeltaPMax*(kp-k1)*(r1+r2),(kp*2*r1*r2))
PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)
#*************************************
......
......@@ -11,9 +11,7 @@ How to run this script:
/path/to/yade ./biaxialSmooth.py
Please amend these instructions if you find that they do not work.
"""
from __future__ import division
from builtins import str
from past.utils import old_div
from esys.escript import *
from esys.finley import Rectangle
from esys.weipa import saveVTK
......@@ -88,14 +86,14 @@ while t < 100: # apply 100 load steps
tractTop = traction*topSurf
forceTop = integrate(tractTop,where=FunctionOnBoundary(dom))
lengthTop = integrate(topSurf,where=FunctionOnBoundary(dom))
fout.write(str(old_div(t*vel,ly))+' '+str(forceTop[1])+' '+str(lengthTop)+'\n')
fout.write(str(t*vel/ly)+' '+str(forceTop[1])+' '+str(lengthTop)+'\n')
vR=prob.getLocalVoidRatio()
fabric=prob.getLocalFabric()
strain = prob.getCurrentStrain()
saveGauss2D(name='./result/gauss/time_'+str(t)+'.dat',strain=strain,stress=stress,fabric=fabric)
volume_strain = trace(strain)
dev_strain = symmetric(strain) - old_div(volume_strain*k,dim)
dev_strain = symmetric(strain) - volume_strain*k/dim
shear = sqrt(2*inner(dev_strain,dev_strain))
saveVTK("./result/vtk/biaxialSmooth_%d.vtu"%t,disp=disp,shear=shear,e=vR)
......
......@@ -11,10 +11,8 @@ How to run this script:
/path/to/yade ./footing.py
Please amend these instructions if you find that they do not work.
"""
from __future__ import division
from builtins import str
from builtins import range
from past.utils import old_div
from esys.escript import *
from esys.finley import ReadGmsh
from esys.weipa import saveVTK
......@@ -96,7 +94,7 @@ while t < 58: # apply 58 loading step; further loading would abort the program d
strain = prob.getCurrentStrain()
saveGauss2D(name='./result/gauss/time_'+str(t)+'.dat',strain=strain,stress=stress,fabric=fabric)
volume_strain = trace(strain)
dev_strain = symmetric(strain) - old_div(volume_strain*k,dim)
dev_strain = symmetric(strain) - volume_strain*k/dim
shear = sqrt(2*inner(dev_strain,dev_strain))
saveVTK("./result/vtk/footing_%d.vtu"%t,disp=disp,stress=stress,shear=shear,e=vR,rot=rotation)
......
......@@ -11,10 +11,8 @@ How to run this script:
/path/to/yade ./retainingSmooth.py
Please amend these instructions if you find that they do not work.
"""
from __future__ import division
from builtins import str
from builtins import range
from past.utils import old_div
from esys.escript import *
from esys.finley import Rectangle
from esys.weipa import saveVTK
......@@ -103,7 +101,7 @@ while t < 100:
strain = prob.getCurrentStrain()
saveGauss2D(name='./result/gauss/time_'+str(t)+'.dat',strain=strain,stress=stress,fabric=fabric)
volume_strain = trace(strain)
dev_strain = symmetric(strain) - old_div(volume_strain*k,dim)
dev_strain = symmetric(strain) - volume_strain*k/dim
shear = sqrt(2*inner(dev_strain,dev_strain))
saveVTK("./result/vtk/retainingSmooth_%d.vtu"%t,disp=disp,stress=stress,shear=shear,e=vR,rot=rotation)
......
......@@ -11,10 +11,8 @@ How to run this script:
/path/to/yade ./triaxialRough.py
Please amend these instructions if you find that they do not work.
"""
from __future__ import division
from builtins import str
from builtins import range
from past.utils import old_div
from esys.escript import *
from esys.finley import Brick
from esys.weipa import saveVTK
......@@ -87,7 +85,7 @@ while t < 100:
tractTop = traction*topSurf
forceTop = integrate(tractTop,where=FunctionOnBoundary(dom))
areaTop = integrate(topSurf,where=FunctionOnBoundary(dom))
fout.write(str(old_div(t*vel,lz))+' '+str(forceTop[2])+' '+str(areaTop)+'\n')
fout.write(str(t*vel/lz)+' '+str(forceTop[2])+' '+str(areaTop)+'\n')
vR=prob.getLocalVoidRatio()
rotation=prob.getLocalAvgRotation()
......@@ -95,7 +93,7 @@ while t < 100:
strain = prob.getCurrentStrain()
saveGauss3D(name='./result/gauss/time_'+str(t)+'.dat',strain=strain,stress=stress,fabric=fabric)
volume_strain = trace(strain)
dev_strain = symmetric(strain) - old_div(volume_strain*k,dim)
dev_strain = symmetric(strain) - volume_strain*k/dim
shear = sqrt(2./3.*inner(dev_strain,dev_strain))
saveVTK("./result/vtk/triaxialRough_%d.vtu"%t,disp=disp,shear=shear,e=vR,rot=rotation)
......
......@@ -11,9 +11,7 @@ How to run this script:
/path/to/yade ./undrain.py
Please amend these instructions if you find that they do not work.
"""
from __future__ import division
from builtins import str
from past.utils import old_div
from esys.escript import *
from esys.finley import Rectangle
from esys.weipa import saveVTK
......@@ -36,7 +34,7 @@ except OSError as exc:
confining = -2.e5; pore = 1.e5 # initial pore pressure
perm = old_div(0.001**2,(180.*8.9e-4)); # unscaled permeability, using KC equation
perm = 0.001**2/(180.*8.9e-4); # unscaled permeability, using KC equation
kf = 2.2e9 # fluid bulk modulus
dt = .1; vel = -0.0001 # time step and loading speed
lx = 0.05; ly = 0.1 # sample dimension
......@@ -79,10 +77,10 @@ while t < 400:
stress = prob.getCurrentStress() # effective stress at GP
strain = prob.getCurrentStrain() # disp. grad at GP
volume_strain = trace(strain) # volumetric strain
dev_strain = symmetric(strain) - old_div(volume_strain*k,dim) # deviatoric strain
dev_strain = symmetric(strain) - volume_strain*k/dim # deviatoric strain
shear = sqrt(2.*inner(dev_strain,dev_strain)) # shear strain
fab = prob.getLocalFabric() # fabric tensor at GP
dev_fab = 4.*(fab-old_div(trace(fab),dim*k))
dev_fab = 4.*(fab-trace(fab)/dim*k)
anis = sqrt(.5*inner(dev_fab,dev_fab))
p = prob.getEquivalentPorosity() # porosity at GP
rot = prob.getLocalAvgRotation() # average rotation at GP
......@@ -101,7 +99,7 @@ while t < 400:
tractTop = traction*topSurf
forceTop = integrate(tractTop,where=FunctionOnBoundary(dom))
lengthTop = integrate(topSurf,where=FunctionOnBoundary(dom))
fout.write(str(old_div(t*vel*dt,ly))+' '+str(forceTop[1])+' '+str(lengthTop)+'\n')
fout.write(str(t*vel*dt/ly)+' '+str(forceTop[1])+' '+str(lengthTop)+'\n')
prob.getCurrentPacking(time=t,prefix='./result/packing/')
time_elapse = time.time() - time_start
......
from __future__ import print_function
from __future__ import division
##############################################################################################################################
#Authors: Luc Sibille luc.sibille@3sr-grenoble.fr
# Franck Lomine franck.lomine@insa-rennes.fr
......@@ -10,7 +9,6 @@ from __future__ import division
##############################################################################################################################
from past.utils import old_div
from yade import pack,timing,utils
......@@ -91,34 +89,34 @@ O.materials.append(FrictMat(young=50e6,poisson=.5,frictionAngle=0.0,density=3000
wallOversizeFactor=1.001
thickness=0.00001;
#bottom box
center= (old_div((lowerCornerW[0]+upperCornerW[0]),2),lowerCornerW[1]-thickness/2.0,old_div((lowerCornerW[2]+upperCornerW[2]),2))
halfSize= (old_div(wallOversizeFactor*fabs(lowerCornerW[0]-upperCornerW[0]),2)+thickness,thickness/2.0,old_div(wallOversizeFactor*fabs(lowerCornerW[2]-upperCornerW[2]),2)+thickness)
center= ((lowerCornerW[0]+upperCornerW[0])/2,lowerCornerW[1]-thickness/2.0,(lowerCornerW[2]+upperCornerW[2])/2)
halfSize= (wallOversizeFactor*fabs(lowerCornerW[0]-upperCornerW[0])/2+thickness,thickness/2.0,wallOversizeFactor*fabs(lowerCornerW[2]-upperCornerW[2])/2+thickness)
b1=utils.box(center=[center[0],center[1],center[2]],extents=[halfSize[0],halfSize[1],halfSize[2]],color=[0,1,0],fixed=True,wire=True,material='walls')
O.bodies.append(b1)
#--
#Top box
center=(old_div((lowerCornerW[0]+upperCornerW[0]),2),upperCornerW[1]+thickness/2.0,old_div((lowerCornerW[2]+upperCornerW[2]),2))
halfSize =(old_div(wallOversizeFactor*fabs(lowerCornerW[0]-upperCornerW[0]),2)+thickness,thickness/2.0,old_div(wallOversizeFactor*fabs(lowerCornerW[2]-upperCornerW[2]),2)+thickness)
center=((lowerCornerW[0]+upperCornerW[0])/2,upperCornerW[1]+thickness/2.0,(lowerCornerW[2]+upperCornerW[2])/2)
halfSize =(wallOversizeFactor*fabs(lowerCornerW[0]-upperCornerW[0])/2+thickness,thickness/2.0,wallOversizeFactor*fabs(lowerCornerW[2]-upperCornerW[2])/2+thickness)
b2=utils.box(center=[center[0],center[1],center[2]],extents=[halfSize[0],halfSize[1],halfSize[2]],color=[0,1,0],fixed=True,wire=True,material='walls')
O.bodies.append(b2)
#--
center=(lowerCornerW[0]-thickness/2.0,old_div((lowerCornerW[1]+upperCornerW[1]),2),old_div((lowerCornerW[2]+upperCornerW[2]),2))
halfSize=(thickness/2.0,old_div(wallOversizeFactor*fabs(lowerCornerW[1]-upperCornerW[1]),2)+thickness,old_div(wallOversizeFactor*fabs(lowerCornerW[2]-upperCornerW[2]),2)+thickness)
center=(lowerCornerW[0]-thickness/2.0,(lowerCornerW[1]+upperCornerW[1])/2,(lowerCornerW[2]+upperCornerW[2])/2)
halfSize=(thickness/2.0,wallOversizeFactor*fabs(lowerCornerW[1]-upperCornerW[1])/2+thickness,wallOversizeFactor*fabs(lowerCornerW[2]-upperCornerW[2])/2+thickness)
b3=utils.box(center=[center[0],center[1],center[2]],extents=[halfSize[0],halfSize[1],halfSize[2]],color=[0,1,0],fixed=True,wire=True,material='walls')
O.bodies.append(b3)
#--
center=(upperCornerW[0]+thickness/2.0,old_div((lowerCornerW[1]+upperCornerW[1]),2),old_div((lowerCornerW[2]+upperCornerW[2]),2))
halfSize=(thickness/2.0,old_div(wallOversizeFactor*fabs(lowerCornerW[1]-upperCornerW[1]),2)+thickness,old_div(wallOversizeFactor*fabs(lowerCornerW[2]-upperCornerW[2]),2)+thickness)
center=(upperCornerW[0]+thickness/2.0,(lowerCornerW[1]+upperCornerW[1])/2,(lowerCornerW[2]+upperCornerW[2])/2)
halfSize=(thickness/2.0,wallOversizeFactor*fabs(lowerCornerW[1]-upperCornerW[1])/2+thickness,wallOversizeFactor*fabs(lowerCornerW[2]-upperCornerW[2])/2+thickness)
b4=utils.box(center=[center[0],center[1],center[2]],extents=[halfSize[0],halfSize[1],halfSize[2]],color=[0,1,0],fixed=True,wire=True,material='walls')
O.bodies.append(b4)
#--
center=(old_div((lowerCornerW[0]+upperCornerW[0]),2),old_div((lowerCornerW[1]+upperCornerW[1]),2),lowerCornerW[2]-thickness/2.0)
halfSize=(old_div(wallOversizeFactor*fabs(lowerCornerW[0]-upperCornerW[0]),2)+thickness,old_div(wallOversizeFactor*fabs(lowerCornerW[1]-upperCornerW[1]),2)+thickness,thickness/2.0)
center=((lowerCornerW[0]+upperCornerW[0])/2,(lowerCornerW[1]+upperCornerW[1])/2,lowerCornerW[2]-thickness/2.0)
halfSize=(wallOversizeFactor*fabs(lowerCornerW[0]-upperCornerW[0])/2+thickness,wallOversizeFactor*fabs(lowerCornerW[1]-upperCornerW[1])/2+thickness,thickness/2.0)
b5=utils.box(center=[center[0],center[1],center[2]],extents=[halfSize[0],halfSize[1],halfSize[2]],color=[0,1,0],fixed=True,wire=True,material='walls')
O.bodies.append(b5)
#--
center=(old_div((lowerCornerW[0]+upperCornerW[0]),2),old_div((lowerCornerW[1]+upperCornerW[1]),2),upperCornerW[2]+thickness/2.0)
halfSize=(old_div(wallOversizeFactor*fabs(lowerCornerW[0]-upperCornerW[0]),2)+thickness,old_div(wallOversizeFactor*fabs(lowerCornerW[1]-upperCornerW[1]),2)+thickness,thickness/2.0);
center=((lowerCornerW[0]+upperCornerW[0])/2,(lowerCornerW[1]+upperCornerW[1])/2,upperCornerW[2]+thickness/2.0)
halfSize=(wallOversizeFactor*fabs(lowerCornerW[0]-upperCornerW[0])/2+thickness,wallOversizeFactor*fabs(lowerCornerW[1]-upperCornerW[1])/2+thickness,thickness/2.0);
b6=utils.box(center=[center[0],center[1],center[2]],extents=[halfSize[0],halfSize[1],halfSize[2]],color=[0,1,0],fixed=True,wire=True,material='walls')
O.bodies.append(b6)
......@@ -127,7 +125,7 @@ O.bodies.append(b6)
# Creation of the assembly of particles
# defintiion of the material for particles
O.materials.append(FrictMat(young=50e6,poisson=.5,frictionAngle=old_div(35*3.1415926535,180),density=3000,label='spheres'))
O.materials.append(FrictMat(young=50e6,poisson=.5,frictionAngle=35*3.1415926535/180,density=3000,label='spheres'))
## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
sp.makeCloud(lowerCornerS,upperCornerS,-1,0.33,nbSpheres,False, 0.5,seed=1) #"seed" make the "random" generation always the same
......
......@@ -3,9 +3,7 @@
# The script was used to generate the result and supplementary material (video) of [1]. The only difference is the problem size (40k particles in the paper vs. 1k (default) in this version)
# [1] Yuan, C., & Chareyre, B. (2017). A pore-scale method for hydromechanical coupling in deformable granular media. Computer Methods in Applied Mechanics and Engineering, 318, 1066-1079. (http://www.sciencedirect.com/science/article/pii/S0045782516307216)
from __future__ import division
from builtins import str
from past.utils import old_div
import matplotlib; matplotlib.rc('axes',grid=True)
from yade import pack
import pylab
......@@ -62,7 +60,7 @@ O.engines=[
while 1:
O.run(1000,True)
unb=unbalancedForce()
if unb<0.01 and old_div(abs(triax.goal1-triax.meanStress),abs(triax.goal1))<0.001:
if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
break
#############################
......@@ -78,7 +76,7 @@ setContactFriction(radians(finalFricDegree))
while 1:
O.run(1000,True)
unb=unbalancedForce()
if unb<0.001 and old_div(abs(triax.goal1-triax.meanStress),abs(triax.goal1))<0.001:
if unb<0.001 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
break
triax.depth0=triax.depth
......@@ -132,7 +130,7 @@ unsat.surfaceTension = 10
file=open('pcSwStrain.txt',"w")
for pg in arange(1.0e-5,15.0,0.1):
#increase gaz pressure at the top boundary
unsat.bndCondValue=[0,0,old_div((-1.0)*pg*unsat.surfaceTension,meanDiameter),0,0,0]
unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]
#compute the evolution of interfaces
unsat.invasion()
#save the phases distribution in vtk format, to be displayed by paraview
......
......@@ -23,9 +23,7 @@
## ______________ First section, similar to triax-tutorial/script-session1.py _________________
from __future__ import print_function
from __future__ import division
from builtins import range
from past.utils import old_div
from yade import pack
num_spheres=1000# number of spheres
......@@ -44,8 +42,8 @@ sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "r
sp.toSimulation(material='spheres')
triax=TriaxialStressController(
maxMultiplier=1.+old_div(2e4,young), # spheres growing factor (fast growth)
finalMaxMultiplier=1.+old_div(2e3,young), # spheres growing factor (slow growth)
maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
thickness = 0,
stressMask = 7,
max_vel = 0.005,
......@@ -73,7 +71,7 @@ triax.goal1=triax.goal2=triax.goal3=-10000
while 1:
O.run(1000, True)
unb=unbalancedForce()
if unb<0.001 and old_div(abs(-10000-triax.meanStress),10000)<0.001:
if unb<0.001 and abs(-10000-triax.meanStress)/10000<0.001:
break
setContactFriction(radians(finalFricDegree))
......@@ -115,7 +113,7 @@ O.dynDt=False
O.run(1,1)
Qin = flow.getBoundaryFlux(2)
Qout = flow.getBoundaryFlux(3)
permeability = old_div(abs(Qin),1.e-4) #size is one, we compute K=V/∇H
permeability = abs(Qin)/1.e-4 #size is one, we compute K=V/∇H
print("Qin=",Qin," Qout=",Qout," permeability=",permeability)
#C. now the oedometer test, drained at the top, impermeable at the bottom plate
......@@ -128,18 +126,18 @@ newton.damping=0
#keep in mind that we are not in an homogeneous material and the small strain
#assumption is not verified => we don't expect perfect match
#there can be also an overshoot of pressure in the very beginning due to dynamic effects
Cv=old_div(permeability*modulus,1e4)
Cv=permeability*modulus/1e4
zeroTime=O.time
zeroe22 = - triax.strain[1]
dryFraction=0.05 #the top layer is affected by drainage on a certain depth, we account for it here
drye22 = old_div(1000,modulus*dryFraction)
drye22 = 1000/modulus*dryFraction
wetHeight=1*(1-dryFraction)
def consolidation(Tv): #see your soil mechanics handbook...
U=1
for k in range(50):
M=old_div(pi,2*(2*k+1))
U=U-old_div(2,M**2*exp(-M**2*Tv))
M=pi/2*(2*k+1)
U=U-2/M**2*exp(-M**2*Tv)
return U
triax.goal2=-11000
......@@ -150,7 +148,7 @@ from yade import plot
## a function saving variables
def history():
plot.addData(e22=-triax.strain[1]-zeroe22,e22_theory=drye22+old_div((1-dryFraction)*consolidation(old_div((O.time-zeroTime)*Cv,wetHeight**2))*1000.,modulus),t=O.time,p=flow.getPorePressure((0.5,0.1,0.5)),s22=-triax.stress(3)[1]-10000)
plot.addData(e22=-triax.strain[1]-zeroe22,e22_theory=drye22+(1-dryFraction)*consolidation((O.time-zeroTime)*Cv/wetHeight**2)*1000./modulus,t=O.time,p=flow.getPorePressure((0.5,0.1,0.5)),s22=-triax.stress(3)[1]-10000)
#plot.addData(e22=-triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5)))
O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]
......
from __future__ import print_function
from __future__ import division
#########################################################################################################################################################################
# Author: Raphael Maurin, raphael.maurin@imft.fr
# 07/07/2016
......@@ -12,7 +11,6 @@ from __future__ import division
############################################################################################################################################################################
#Import libraries
from past.utils import old_div
from yade import pack, plot
import math
import random as rand
......@@ -41,7 +39,7 @@ gravityVector = Vector3(0,0.0,-9.81) #Gravity vector
#Particles contact law/material parameters
normalStiffness = 1e4
youngMod = old_div(normalStiffness,diameterPart) #Young modulus of the particles from the stiffness wanted.
youngMod = normalStiffness/diameterPart #Young modulus of the particles from the stiffness wanted.
poissonRatio = 0.5 #poisson's ratio of the particles. Classical values, does not have much influence
#Material of particle 1, 2, and 3, with different density densPart1, densPart2 and densPart3 defined above (1500, 1000 and 500kg/m3 by default)
......@@ -85,7 +83,7 @@ O.engines = [
[Law2_ScGeom_ViscElPhys_Basic()]
,label = 'interactionLoop'),
#Apply an hydrodynamic force to the particles
HydroForceEngine(densFluid = densFluidPY,viscoDyn = kinematicViscoFluid*densFluidPY,zRef = groundPosition,gravity = gravityVector,deltaZ = old_div(fluidHeight,ndimz),expoRZ = 0.,lift = False,nCell = ndimz,vCell = 1.,vxFluid = np.zeros(ndimz),phiPart = np.zeros(ndimz),vxPart = np.zeros(ndimz),vFluctX = np.zeros(len(O.bodies)),vFluctY = np.zeros(len(O.bodies)),vFluctZ = np.zeros(len(O.bodies)),ids = idApplyForce, label = 'hydroEngine'),
HydroForceEngine(densFluid = densFluidPY,viscoDyn = kinematicViscoFluid*densFluidPY,zRef = groundPosition,gravity = gravityVector,deltaZ = fluidHeight/ndimz,expoRZ = 0.,lift = False,nCell = ndimz,vCell = 1.,vxFluid = np.zeros(ndimz),phiPart = np.zeros(ndimz),vxPart = np.zeros(ndimz),vFluctX = np.zeros(len(O.bodies)),vFluctY = np.zeros(len(O.bodies)),vFluctZ = np.zeros(len(O.bodies)),ids = idApplyForce, label = 'hydroEngine'),
#To plot the wall normal position of the spheres with time
PyRunner(command = 'plotPos()', virtPeriod = 0.01, label = 'plot'),
# Integrate the equation and calculate the new position/velocities...
......@@ -98,7 +96,7 @@ O.dt = 5e-7 #Low no purpose, in order to observe the sedimentation
#Plot the wall normal position of the spheres with time
def plotPos():
plot.addData(z1 = old_div(O.bodies[2].state.pos[2],fluidHeight),z2 = old_div(O.bodies[3].state.pos[2],fluidHeight),z3 = old_div(O.bodies[4].state.pos[2],fluidHeight), time = O.time)
plot.addData(z1 = O.bodies[2].state.pos[2]/fluidHeight,z2 = O.bodies[3].state.pos[2]/fluidHeight,z3 = O.bodies[4].state.pos[2]/fluidHeight, time = O.time)
if O.time>endTime:
print('\nEnd of the simulation, {0}s simulated as asked!\n'.format(endTime))
O.pause()
......
from __future__ import print_function
from __future__ import division
#########################################################################################################################################################################
# Author: Raphael Maurin, raphael.maurin@imft.fr
# 08/07/2016
......@@ -22,7 +21,6 @@ from __future__ import division
#Import libraries
from past.utils import old_div
from yade import pack, plot
import math
import random as rand
......@@ -58,7 +56,7 @@ endTime = 10 #Time simulated (in seconds)
expoDrag_PY = 3.1 # Richardson Zaki exponent for the hindrance function of the drag force applied to the particles
ndimz = 20 #Number of cells in the height
dz = old_div(widthCell*diameterPart,ndimz) # Fluid discretization step in the wall-normal direction
dz = widthCell*diameterPart/ndimz # Fluid discretization step in the wall-normal direction
# Initialization of the main vectors
vxFluidPY = np.ones(ndimz)*18.5 # Vertical fluid velocity profile: u^f = u_x^f(z) e_x, with x the streamwise direction and z the wall-normal
......@@ -75,7 +73,7 @@ gravityVector = Vector3(-9.81,0.0,0.) #Gravity vector. Inclined along x to have
#Particles contact law/material parameters
maxPressure = (densPart-densFluidPY)*phiPartMax*Nlayer*diameterPart*9.81#Estimated max particle pressure from the static load
normalStiffness = maxPressure*diameterPart*1e4 #Evaluate the minimal normal stiffness to be in the rigid particle limit (cf Roux and Combe 2002)