yade install problem

Welcome to yade!

Is there a reason you are compiling from sources? Since you are new to yade, I recommend installing yadedaily from prebuilt packages to become familiar with the software before moving on to compiling from sources. This can be achieved simply with the following commands:

sudo bash -c 'echo "deb http://www.yade-dem.org/packages/ xenial/" >> /etc/apt/sources.list'
wget -O - http://www.yade-dem.org/packages/yadedev_pub.gpg | sudo apt-key add -
sudo apt-get update
sudo apt-get install yadedaily

(note, you will want to replace xenial with your distribution).

That error is telling you that it wants to find the file CMakeLists.txt in /home/tanvir/myyade/trunk, but it is unable to. You should navigate to that folder and check to see if CMakeLists.txt exists there. We already know it doesn't, and this means your path/to/sources is incorrect or you did not clone trunk properly. So find out the path to the folder that contains CMakeLists.txt, and use that as your path/to/sources.

I am guessing you created a folder called "trunk" and then cloned yade inside of it. Am I right? This would mean your path to sources is actually:

/home/tanvir/myyade/trunk/trunk

Or maybe you did not clone yade at all? In which case you will want to start the installation guide [1] from the beginning and follow each step.

[1]https://yade-dem.org/doc/installation.html

Hi, thank your for your reply. I want to install yadedaily.

sudo bash -c 'echo "deb http://www.yade-dem.org/packages/ xenial/" >> /etc/apt/sources.list'
wget -O - http://www.yade-dem.org/packages/yadedev_pub.gpg | sudo apt-key add -
sudo apt-get update
sudo apt-get install yadedaily

is this command sufficient for all the procedure or do I need follow prerequisites, compilation and all other steps?

Yes, that series of commands is sufficient, assuming you are using a debian or ubuntu distribution of linux.

Just a comment, for Ubuntu (preferably a newer version like 16.04), you can install yade just by

sudo apt-get install yade

Jan

my version is ubuntu16.04.. I have install yade using ''sudo apt-get install yade''. now I am trying to run a python script of but getting following error,

Welcome to Yade 1.20.0
TCP python prompt on localhost:9000, auth cookie `dksesy'
XMLRPC info provider on http://localhost:21000
Running script 20160108PeriodicStreamVersionCOUPLING3D.py
Traceback (most recent call last):
  File "/usr/bin/yade", line 182, in runScript
    execfile(script,globals())
  File "20160108PeriodicStreamVersionCOUPLING3D.py", line 38, in <module>
    from nsmp1d_yade3D525Phi061 import * #tanvir
ImportError: No module named nsmp1d_yade3D525Phi061
/home/tanvir/.local/lib/python2.7/site-packages/IPython/config.py:13: ShimWarning: The `IPython.config` package has been deprecated since IPython 4.0. You should import from traitlets.config instead.
  "You should import from traitlets.config instead.", ShimWarning)
/home/tanvir/.local/lib/python2.7/site-packages/IPython/core/interactiveshell.py:448: UserWarning: As of IPython 5.0 `PromptManager` config will have no effect and has been replaced by TerminalInteractiveShell.prompts_class
  warn('As of IPython 5.0 `PromptManager` config will have no effect'
[[ ^L clears screen, ^U kills line. F12 controller, F11 3d view (use h-key for showing help), F10 both, F9 generator, F8 plot. ]]

Is there any possibility that yade is not installed properly? Or do I need to install something else along with yade?

Hi Tanvir,

it seems that Yade is installed correctly. The problem seems to be inside your script, namely:

 File "20160108PeriodicStreamVersionCOUPLING3D.py", line 38, in <module>
    from nsmp1d_yade3D525Phi061 import * #tanvir
ImportError: No module named nsmp1d_yade3D525Phi061

could you please include the content of your script? Also please read [1]

Sometimes this error is present even if the module imported is inside the directory where yade is started. In this case
import sys
sys.path.append('.')
might help..

cheers
Jan

[1] https://yade-dem.org/wiki/Howtoask

Hi Jan, here is the script.
#New version, amelioration wrt 20150715 :
#- IMPORTANT Modify the length of the vector evaluated from ndimz to ndimz. Works with Yade6 and is not adapted to the other cases
#- Add a limiter to avoid negative tau transmitted to the fluid (unphysical and source of numerical instabilities
#- Add a new configuration for the turbulentFluctuation function. Correspond to a new formulation of the C++ HydroForceEngine. Possibility to use a new turbulent fluctuation model where each particle fluctuation lifetime depends on its wall-normal position.
#- Add an option for fluidProfile imposed case.
#- Register the averaged velocity of the part1 and 2 in twoSize mode (Yade6)

#If activated, adapt the script to the associated computer path
Cluster = 0
philippeComputer = 0
julienComputer = 0
lamiaComputer = 0
raphaelComputer = 0

#Import libraries
from yade import pack, plot
import math
import random as rand
import numpy as np
np.set_printoptions(threshold=numpy.nan)
import yade.timing
import pylab as pyl
from matplotlib import pyplot as pyp
import matplotlib.gridspec as gridspec
#Register the path for later use (tests + save)
scriptPath = os.path.abspath(os.path.dirname(sys.argv[-1]))

##
## Main parameters of the simulation
##
diameterPart = 6e-3 #Diameter of the particles, in meter
densPart = 2500 #density of the particles, in kg/m3
waterDepth = 7. #Water depth in diameter
ndimz = 525 #Number of cells in the height
from nsmp1d_yade3D525Phi061 import *

lengthCell = 15 #Streamwise length of the periodic cell
widthCell = 15 #Spanwise length of the periodic cell
Nlayer = 10. #nb of layer of particle, in diameter

fluidHeight = (Nlayer+waterDepth)*diameterPart #Height of the flow from the bottom of the sample

quasi2D = 0
twoSize = 0

#If twoSize option activated, overwrite the previous input #tanvir
if twoSize:
 Nlayer1 = 10 #nb of layer of particle type 1
 Nlayer2 = 1 #nb of layer of particle type 2
 diameterPart1 = 6e-3 #Diameter of the particles type 1, in meter
 diameterPart2 = 4e-3 #Diameter of the particles type 2, in meter
 fluidHeight = Nlayer1*diameterPart1 + Nlayer2*diameterPart2 + 7*diameterPart1 #free surface position
 diameterPart = diameterPart1 #Size 1 taken as reference

#Fluid library importation and precisions
phiPartMax = 0.61 #Value of the maximum volume fraction used for the thresholding

##
## Option of the simulation
##

#Imposed fluid profile? If yes, require to add a file fluid.py in the folder with inside the value of vxFluid vector + turbStress vector if you want to add turbulent fluctuations also.
imposedFluidPro = 0

#Temporary imposed profile: impose a fluid profile for 50 seconds (considered as equilibrium) and switch to the fluid resolution
initImposedFluidPro = 0

#Activation and choice of the fluid turbulent fluid velocity fluctuation. Cannot choose both. If none, will not apply fluct.
turbModel1 = 1
turbModel2 = 0

#Stress averaging: if activated, compute and save the stress tensor depth profile at each execution of measure()
stressSave = 0
pythonStressCalculation = 0 #use the python function instead of C++

#Option to put lateral walls, to break the biperiodicity and study the influence of walls.
lateralWalls = 0

#Option for time averaging of the DEM data before fluid resolution
timeAveraging = 0 #1 activated, 0 not activated
nbAv = 5. #number of averaging between the fluid resolution

#Averaging parameter : if activated, evaluate the average using python function partConcVelocity() instead of the C++ one
pythonAverage = 0

#debug stress mode, record all the quantities necessary to study the momentum balance in detail.
debugStress = 0

#Video
video = 1 #1, make a movie
if video:
 snapshotEngine = 0 #To make the movie later using mencoder (require the qt.viewer open, obtain what you see)
 paraview = 1 #To make the movie later using paraview (do not require the qt.viewer open, faster)
 fps = 1. #frame per second #tanvir

##
## Secondary parameters of the simulation (less changed)
##

#Particles parameters
restitCoef = 0.5 #Restitution coefficient of the particles
partFrictAngle = atan(0.4) #friction angle of the particles, in radian
slope = 0.1 #Angle of the slope in radian #tanvir

# Fluid parameters
expoRZ = 3.1 # Richardson Zaki exponent for the hindrance function of the drag
dpdx = 0.0 #pressure gradient along streamwise direction
densFluid = 1000. #Density of the fluid
kinematicViscoFluid = 1e-6 #kinematic viscosity of the fluid

# Fluid simulation parameters (RANS), see fluidResolution function
dsig = np.zeros(ndimz)
dz = fluidHeight/(1.0*(ndimz-1))
sig=pyl.linspace(0e0,1e0,ndimz)
for i in range(ndimz-1): # calcul du pas d'espace
    dsig[i] = sig[i+1] - sig[i]

#In order to fasten the script, some correspondance are supposed. Should be respected otherwise it does not work well
dtFluid = 1e-5 #Time step for the fluid resolution
fluidResolPeriod = 1e-2 #1/fluidResolPeriod = frequency of the calculation of the fluid profile
measurePeriod = (1e-1)*5 # We perform measurement every measurePeriod seconds #tanvir
savePeriod = 1e1 # We save the simulation (to be able to reload) every savePeriod seconds.

#Quasi 2D option
if quasi2D ==1:
 print '\nQuasi 2D option activated\n'
 lateralWalls = 1
 widthCell = 6.5/6.
 phiPartMax = 0.51
 print 'phiPartMax = 0.51\n'

#
## Tests on the script
#
#Width check and Value of the fluid library and size of the system
if widthCell<2 and phiPartMax>0.6:
 if widthCell<=1:
  print 'Impossible to run simulation with width exactly equal to the diameter particle size !\n'
  exit()
 else:
  print '\nCareful !! Uncoherence between the 2D character and the value of phiPartMax !! \n'
#Value of ndimz and dz
if ndimz != int(round(30*fluidHeight/diameterPart)):
 print('\n!!Problem dz = d/{0} !!\nUse ndimz = {1} instead to have dz = d/30 (and change accordingly the fluid library import)\n'.format(diameterPart/fluidHeight*ndimz,int(round(30*fluidHeight/diameterPart))))
# if quasi2D!=1:
# exit()

#Option compatibility
if measurePeriod<fluidResolPeriod:
 print '\nNot possible to have measurePeriod<fluidResolPeriod in the present configuration, modify the script to make measurement inside measure() instead of fluidResolution() !\n'
 exit()

turbModel = 0
if turbModel1 ==1 or turbModel2==1:
 turbModel = 1
 if turbModel1 ==1 and turbModel2==1:
  print "It is not possible to have the two fluctuations model in the same time ! \n Modify turbModel1 or turbModel2"
  exit()

#Reloading simulation
reloadingSimulation = 0
if os.path.exists(scriptPath +'/paramSave.py'):
 reloadingSimulation = 0#1 Out for now, useless, it is not working

##
## Initialization of the parameters
##

# Averaging/Save
vxPart = np.zeros(ndimz) # streamwise particle velocity depth profile
vyPart = np.zeros(ndimz) # spanmwise particle velocity depth profile
vzPart = np.zeros(ndimz) # normal particle velocity depth profile
qsMean = 0 #Mean bead flux
dragTerm = np.zeros(ndimz) # drag term transmitted to the fluid vector
averageDrag = np.zeros(ndimz) # drag field vector
stressProfile = getStressProfile(1.0,ndimz,dz,0.0,np.zeros(ndimz),np.zeros(ndimz),np.zeros(ndimz)) #Matrix full of zero
positionPartX = np.zeros(len(O.bodies)) #Defined global in order to use
positionPartZ = np.zeros(len(O.bodies)) #Defined global in order to use
velocityPartX = np.zeros(len(O.bodies)) # the saving function
velocityPartZ = np.zeros(len(O.bodies)) # the saving function
if twoSize:
 phiPart1 = np.zeros(ndimz) # Concentration vector
 vxPart1 = np.zeros(ndimz) # streamwise particle velocity depth profile
 vyPart1 = np.zeros(ndimz) # spanwise particle velocity depth profile
 vzPart1 = np.zeros(ndimz) # wall-normal particle velocity depth profile
 phiPart2 = np.zeros(ndimz) # Concentration vector
 vxPart2 = np.zeros(ndimz) # streamwise particle velocity depth profile
 vyPart2 = np.zeros(ndimz) # spanwise particle velocity depth profile
 vzPart2 = np.zeros(ndimz) # wall-normal particle velocity depth profile

#Initialization of the variables necessary for timeAveraging option
if timeAveraging:
 n_dt = 1.
 vxPart_dt = np.zeros(ndimz)
 phiPart_dt = np.zeros(ndimz)
 averageDrag_dt = np.zeros(ndimz)
 if twoSize:
  phiPart1_dt = np.zeros(ndimz)
  averageDrag1_dt = np.zeros(ndimz)
  vxPart1_dt = np.zeros(ndimz)
  vyPart1_dt = np.zeros(ndimz)
  vzPart1_dt = np.zeros(ndimz)
  phiPart2_dt = np.zeros(ndimz)
  averageDrag2_dt = np.zeros(ndimz)
  vxPart2_dt = np.zeros(ndimz)
  vyPart2_dt = np.zeros(ndimz)
  vzPart2_dt = np.zeros(ndimz)

#Parameter initialized only if it is the initial simulation, otherwise they are charged
if reloadingSimulation !=1:
 #Parameters which should be loaded and not initiallized if we recover the simulation
 phiPart = np.zeros(ndimz) # Concentration vector
 turbulentViscosity = np.zeros(ndimz) # Turbulent viscocity
 waterDepth = 0.0 #Water depth : height of the "clear" water
 vxFluid = np.zeros(ndimz) # Fluid velocity at time step n
 turbStress = np.zeros(ndimz) #Turbulent stress vector
 fileNumber = 0 #counter to save the files

###########################
## PARAMETERS DECLARATION##
###########################

#Geometrical configuration, define useful quantities
height = 5*fluidHeight #heigth of the periodic cell, bigger than the real height to take into account jumps
length = lengthCell*diameterPart #length of the stream
width = widthCell*diameterPart #width of the stream
gravityVector = Vector3(9.81*sin(slope),0.0,-9.81*cos(slope))
gravityVectorApplied = Vector3(0.0,0.0,-9.81*cos(slope))
groundPosition = height/4.0
VLayer = dz*length*width #Volume of a layer
partVolume = pi/6.*pow(diameterPart,3)
if twoSize:
 partVolume1 = partVolume
 partVolume2 = pi/6.*pow(diameterPart2,3)

#Particles contact law parameters
maxPressure = (densPart-densFluid)*phiPartMax*Nlayer*diameterPart*abs(gravityVector[2]) #Estimated max particle pressure ("hydrostatic")
normalStiffness = maxPressure*diameterPart*1e4 #Evaluate the minimal normal stiffness to be in the rigid particle limit (cf Roux and Combe 2002)
youngMod = normalStiffness/diameterPart #Young modulus of the particles from the stiffness wanted. Need that in order to use the restitution coefficient definition + more practical when polydisperse sample.
poissonRatio = 0.5 #poisson's ratio of the particles. Classical values, not much influence (see ref in Da Cruz et al 2005)
if twoSize:
 youngMod1 = normalStiffness/diameterPart1
 youngMod2 = normalStiffness/diameterPart2
 O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod1, poisson=poissonRatio, density=densPart, frictionAngle=partFrictAngle, label='Mat1'))
 O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod2, poisson=poissonRatio, density=densPart, frictionAngle=partFrictAngle, label='Mat2'))
#Material
O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod, poisson=poissonRatio, density=densPart, frictionAngle=partFrictAngle, label='Mat'))

# if it is the initial simulation create the framework and the simulation loop, otherwise load it
if reloadingSimulation !=1:
 ##########################################
 ## FRAMEWORK CREATION (sphere, wall,...)##
 ##########################################

 #Definition of the semi-periodic cell
 O.periodic = True
 O.cell.setBox(length,width,height)

 #To be compatible with lateral walls
 leftLimitY = 0.0
 rightLimitY = width
 centerLimitY = width/2.0
 #If lateralWalls, add the walls and increase the cell size to avoid particle touching the walls to appear on the other side
 if lateralWalls:
  #Warn the user
  print '\nlateralWalls option activated: mono-periodic boundary condition only !\n'
  #Increase the cell size
  O.cell.setBox(length,width+2*diameterPart,height)
  #Modify accordingly the position of the center of the cell and the wall right and left position
  leftLimitY = diameterPart
  rightLimitY = width+diameterPart
   centerLimitY = diameterPart + width/2.0
  #Define the wall and add to the simulation
  sidePlaneL = box(center= (length/2.0,leftLimitY,height/2.0),extents=(2000,0,height*10),fixed=True,wire = True,color = (1,0,0), material = 'Mat')
  sidePlaneR = box(center= (length/2.0,rightLimitY,height/2.0),extents=(2000,0,height*10.0),fixed=True,wire=True, material = 'Mat',color = (0,0,1))
  O.bodies.append([sidePlaneR,sidePlaneL])

 if quasi2D == 1:
  # Regular arrangement of spheres sticked at the bottom with random height
  a = range(0,int(length/(diameterPart))) #The length is divided in diameter of particle
  for b in a: #loop creating a set of sphere sticked at the bottom with a (uniform) random altitude comprised between 0.5 (diameter/12) and 5.5mm (11diameter/12) with steps of 0.5mm. The repartition is made around groundPosition. Exactly the same set up as the one from the experiment.
   n = rand.randrange(0,12,1)/12.0*diameterPart #Define a number between 0 and 11/12 diameter with steps of 1/12 diameter (0.5mm in the experiment)
   O.bodies.append(sphere((b*diameterPart,centerLimitY,groundPosition - 11*diameterPart/12.0/2.0 + n),diameterPart/2.,color=(0,0,0),fixed = True,material = 'Mat'))

 else:
  #Execute the file containing the database which is a 2D list containing the positions and radius of the particles composing the random fixed bed. Each position x,y (integer) contain a list of particle with position and radius at this place
  if Cluster==1:
   pathFixedBot = '/home/maurin/'
  elif philippeComputer==1:
   pathFixedBot = '/home/philippe/fixedBot/'
  elif tanvir==1:
   pathFixedBot = '/home/tanvir/output/'
  elif julienComputer==1:
   pathFixedBot = '/home/julien/Yade/scriptsRaphael/'
  elif lamiaComputer==1:
   pathFixedBot = '/home/lamia/periodicStream/'
  elif raphaelComputer:
   pathFixedBot = '/home/raphael/these/scripts/twoWayCoupling/'
  else:
   print '\nNo precised computer used ! Cannot resolve the fixed bottom file path accordingly ! Precise the computer used at line 12\n'
   exit()

  if diameterPart==6e-3:
   if lengthCell>= 10:
     execfile(pathFixedBot+'fixedBotDataLongMONO.py')
   else:
     execfile(pathFixedBot+'fixedBotDataShort.py')
     print '\n Use the short database for random fixed bottom generation ! \n Be careful, polydisperse bottom !\n'
  elif diameterPart==3e-3:
   execfile(pathFixedBot+'fixedBotDataMONOd3.py')
  elif diameterPart==12e-3:
   execfile(pathFixedBot+'fixedBotDataMONOd12.py')
  else:
   print '\nno fixed bottom of the right size !! quit the simulation ! \n'
   exit()

  xMax = len(fixedBottomData) # Maximum size of the data base along x axis
  yMax = len(fixedBottomData[0]) # Maximum size of the data base along y axis
  x = rand.randrange(0,xMax,1) # Generate a uniform random number between 0 and the size of the data base along x (in diameterPart)
  y = rand.randrange(0,yMax,1) # Generate a uniform random number between 0 and the size of the data base along y (in diameterPart)
  #Starting from the random x and y position, loop over the length and width of the require fixed bed increasing x and y position (in diameterPart). Create then a sphere starting from 0 for x and y.
  for i in range(0,lengthCell):
   for j in range(0,int(widthCell)):
    for k in fixedBottomData[(i+x)%xMax][(j+y)%yMax]:
     O.bodies.append(sphere(center = ((i+k[0]/diameterPart-int(k[0]/diameterPart))*diameterPart,(j+k[1]/diameterPart-int(k[1]/diameterPart))*diameterPart+leftLimitY,groundPosition+k[2]),radius = k[3],material = 'Mat',fixed = True,color = (0,0,0)))

 # Ground reference and Wall on the side
 lowPlane = box(center= (length/2.0,centerLimitY,groundPosition),extents=(200,200,0),fixed=True,wire=False,color = (0.,1.,0.),material = 'Mat') #Build a plane to have a reference for the eyes
 WaterSurface = box(center= (length/2.0,centerLimitY,groundPosition+fluidHeight),extents=(2000,width/2.0,0),fixed=True,wire=False,color = (0,0,1),material = 'Mat',mask = 0) #Build a plane to have a reference for the eyes of the position of the water surface. No interaction with the sphere
 O.bodies.append([lowPlane,WaterSurface]) #add to simulation

 #Count the initial number of particles (for later check)
 initNumber = len(O.bodies)

 #Create a loose cloud of particle inside the cell
 partCloud = pack.SpherePack()

 if twoSize:
  diameterPart = diameterPart1
  partVolume = partVolume1
  Nlayer = Nlayer1

 #Define the number of particle from the density of particle per unit length
 if quasi2D==1:#to keep the correspondance with previous work where I putter NbPart = Nl*L/d
  partNumber = int(Nlayer*lengthCell)
 else:
  partNumber = int(Nlayer*phiPartMax*diameterPart*length*width/partVolume) #Volume of beads
 #Define the deposition height considering that the packing realised by make cloud is 0.2
 depositionHeight = Nlayer*phiPartMax/0.2*diameterPart #Consider that the packing realised by make cloud is 0.2
 #Create a cloud of partNumber particles
 partCloud.makeCloud(minCorner=(0,leftLimitY,groundPosition+diameterPart+1e-4),maxCorner=(length,rightLimitY,groundPosition+depositionHeight),rRelFuzz=0., rMean=diameterPart/2.0, num = partNumber)
 #Send this packing to simulation with material Mat
 partCloud.toSimulation(material='Mat')
 partNumber += initNumber

 #Same for a second size of particles
 if twoSize:
  partCloud2 = pack.SpherePack()
  if quasi2D==1:#to keep the correspondance with previous work where I putter NbPart = Nl*L/d
   partNumber2 = int(Nlayer2*lengthCell*diameterPart1/diameterPart2)
  else:
   partNumber2 = int(Nlayer2*phiPartMax*diameterPart2*length*width/partVolume2) #Volume of beads
  #Define the deposition height considering that the packing realised by make cloud is 0.2
  depositionHeight2 =Nlayer2*diameterPart2*phiPartMax/0.2
  partCloud2.makeCloud(minCorner=(0,leftLimitY,groundPosition+depositionHeight),maxCorner=(length,rightLimitY,groundPosition+depositionHeight+depositionHeight2),rRelFuzz=0., rMean=diameterPart2/2.0, num = partNumber2)
  partCloud2.toSimulation(material='Mat2',color = (0,0,1)) #send this packing to simulation with material Mat2
  #Increment the counters for check and deposition time input
  depositionHeight+=depositionHeight2
  partNumber+=partNumber2
 #Evaluate the deposition time considering the free-fall time of the highest particle to the ground (overestimation)
 depoTime = sqrt(depositionHeight*2/abs(gravityVector[2]))

 #Check
 if len(O.bodies)<partNumber:
  print '\nProblem : too low number of beads in the sample compared with what was asked !\n'
  exit()

 # Collect the ids of the spheres which are dynamic and to which we should add the hydroforce
 dv =1e-2 #Scale of the velocity imposed to the particle, no influence if not extreme case en = 1
 idApplyForce = []
 for b in O.bodies:
  if isinstance(b.shape,Sphere) and b.dynamic:
   if widthCell<2: #For the quasi 2D case, need to push slightly the spheres randomly in order to break the 2D configuration
    deltaV = - dv + rand.uniform(0,2*dv)
    b.state.vel[1]+=deltaV
   idApplyForce+=[b.id]

 #########################
 #### SIMULATION LOOP#####
 #########################

 O.engines = [
  # Reset the forces
  ForceResetter(),
  # Detect the potential contacts
  InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod = True),
  # Calculate the different interactions
  InteractionLoop(
     [Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
     [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
     [Law2_ScGeom_ViscElPhys_Basic()]
  ,label = 'interactionLoop'),
  #Apply an hydrodynamic force to the particles
  HydroForceEngine(densFluid = densFluid,viscoDyn = kinematicViscoFluid*densFluid,zRef = groundPosition,gravity = gravityVectorApplied,deltaZ = dz,expoRZ = expoRZ,lift = False,nCell = ndimz,vCell = length*width*dz ,vxFluid = vxFluid,phiPart = phiPart,vxPart = vxPart,ids = idApplyForce, label = 'addFluidForce', dead = True),
  #Time averaging of the DEM data for fluid resolution
  PyRunner(command = 'timeAveragingDEMdata()', virtPeriod = fluidResolPeriod/nbAv, label = 'timeAve', dead = True),
  #Calculate the fluid profile using the routine
  PyRunner(command = 'fluidResolution(fluidResolPeriod)', virtPeriod = fluidResolPeriod, label = 'fluidResol', dead = True),
  #Calculate the turbulent fluctuations from a DRW random walk model
  PyRunner(command = 'turbulentFluctuation()', virtPeriod = 0.1, label = 'turbFluct', dead = True),
  #Measurement, output files
  PyRunner(command = 'measure()', virtPeriod = measurePeriod, label = 'measurement', dead = True),
  # Check if the packing is stabilized, if yes activate the hydro force on the grains and the slope.
  PyRunner(command='gravityDeposition(depoTime)',virtPeriod = 0.01,label = 'packing'),
  #GlobalStiffnessTimeStepper, determine the time step
  GlobalStiffnessTimeStepper(defaultDt = 1e-4, viscEl = True,timestepSafetyCoefficient = 0.7, label = 'GSTS'),
  # Integrate the equation and calculate the new position/velocities...
  NewtonIntegrator(damping=0.2, gravity=gravityVector, label='newton')
  ]
 #save the initial configuration to be able to recharge the simulation starting configuration easily
 O.saveTmp()
 #run
 O.run()

 #Add engines to make a video
 if video:
  if snapshotEngine:
   print '\n video snapshot activated : takes snapshot of the qt.view as you see it\n'
   O.pause()
   O.engines += [SnapshotEngine(virtPeriod=1.0/fps,fileBase=scriptPath+'/figure/',label='snapshooter')]
  elif paraview:
   print '\n video paraview activated : register vtk file for later use of paraview (video...)\n'
   O.engines += [VTKRecorder(virtPeriod=1/fps,recorders=['spheres','colors'],fileName=scriptPath+'/Paraview/p1-')]
  else:
   print '\nBe careful the video is not activated : need to specify snapshotEngine or paraview value.\n'
   O.pause()

 #Add a pyRunner to re-activate the fluctuations at 50s
 if initImposedFluidPro:
  O.engines += [PyRunner(command='activateFluidResol()',virtPeriod = 10,label = 'actFlResol')]

# else, it is a reloading of simulation, we load the necessary parameters and the simulation
else :
 execfile(scriptPath +'/paramSave.py') #Execute the file containing the parameter registered previously
 O.load(scriptPath + '/save.xml') #load the simulation where it stopped
 print '\n RELOAD THE PREVIOUS SIMULATION\n'

if imposedFluidPro==1 or initImposedFluidPro==1:
 execfile('fluid.py') #Import the fluid profile (+turbStress if necessary)

###################################################################################################################################
################################################### FUNCTION DEFINITION #########################################################
###################################################################################################################################

###### ######
### LET THE TIME FOR THE GRAVITY DEPOSITION AND ACTIVATE THE FLUID AT THE END ###
###### ######
def gravityDeposition(lim):
 #At the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
 #The rest will be run only if O.time>lim = estimation of the deposition time
 if O.time<lim : return
 else :
  newton.damping = 0.0
     packing.dead = True # Remove the present engine for the following

  addFluidForce.dead = False # Activate the fluid force
  fluidResol.dead = False # Activate the fluid profile resolution
# O.timingEnabled = True # Activate a function to check the time taken by the different function

  # If option activated, Activate the turbulent fluctuation calculation
  if turbModel==1:
   turbFluct.dead = False
  # If option activated, activate the function to perform cumulated time averaging on DEM for fluidResol
  if timeAveraging:
   timeAve.dead = False
  #If twoSize C++ option activated, express the size of the two particles in the engine for later averaging
  if twoSize and pythonAverage!=1:
   addFluidForce.radiusPart1 = diameterPart1/2.
   addFluidForce.radiusPart2 = diameterPart2/2.
   addFluidForce.twoSize = True

  if imposedFluidPro==1 or initImposedFluidPro==1:
   fluidResol.dead = True #Force no fluid resolution
   addFluidForce.vxFluid = vxFluid #pass the fluid velocity vector to add drag force value
   addFluidForce.simplifiedReynoldStresses = turbStress #and reynolds stress tensor for fluctuations
  else:
   #Initialization of the fluid resolution
   fluidResolution(fluidResolPeriod)

  measurement.dead = False # Activate the Measurement
 return
###############
#########################################

############# ###################
# Activate the fluid resolution. Used only when initImposedFluidPro option activated ###
############# ###################
def activateFluidResol():
 global initImposedFluidPro
 if O.time>50:
  fluidResol.dead = False #Activate the fluid resolution
  initImposedFluidPro = 0 #Put the option to zero, to avoid confusion in the script
  actFlResol.dead = True #Kill the engine

#######################################################

###################################################################################################
############## Function to time-average the DEM data transmitted to fluid resolution ###########
###################################################################################################
def timeAveragingDEMdata():
 global phiPart_dt
 global vxPart_dt
 global averageDrag_dt
 global phiPart1_dt
 global averageDrag1_dt
 global vxPart1_dt
 global vyPart1_dt
 global vzPart1_dt
 global phiPart2_dt
 global averageDrag2_dt
 global vxPart2_dt
 global vyPart2_dt
 global vzPart2_dt
 global n_dt

 #Re-initialized at each new averaging step
 if n_dt==0:
  vxPart_dt = np.zeros(ndimz)
  phiPart_dt = np.zeros(ndimz)
  averageDrag_dt = np.zeros(ndimz)
  if twoSize:
   phiPart1_dt = np.zeros(ndimz)
   averageDrag1_dt = np.zeros(ndimz)
   vxPart1_dt = np.zeros(ndimz)
   vyPart1_dt = np.zeros(ndimz)
   vzPart1_dt = np.zeros(ndimz)
   phiPart2_dt = np.zeros(ndimz)
   averageDrag2_dt = np.zeros(ndimz)
   vxPart2_dt = np.zeros(ndimz)
   vyPart2_dt = np.zeros(ndimz)
   vzPart2_dt = np.zeros(ndimz)

 #Evaluate the average vx,vy,vz,phi,drag profiles and store it in addFluidForce
 addFluidForce.averageProfile()

 #Store the values in
 phiPart_dt += np.array(addFluidForce.phiPart)
 vxPart_dt += np.array(addFluidForce.vxPart)
 averageDrag_dt += np.array(addFluidForce.averageDrag)
 if twoSize:
  phiPart1_dt += np.array(addFluidForce.phiPart1)
  averageDrag1_dt += np.array(addFluidForce.averageDrag1)
  vxPart1_dt += np.array(addFluidForce.vxPart1)
  vyPart1_dt += np.array(addFluidForce.vyPart1)
  vzPart1_dt += np.array(addFluidForce.vzPart1)
  phiPart2_dt += np.array(addFluidForce.phiPart2)
  averageDrag2_dt += np.array(addFluidForce.averageDrag2)
  vxPart2_dt += np.array(addFluidForce.vxPart2)
  vyPart2_dt += np.array(addFluidForce.vyPart2)
  vzPart2_dt += np.array(addFluidForce.vzPart2)

 n_dt+=1.

###############
#########################################

################ #########################
########## EVALUATE THE TURBULENT FLUCTUATION ASSOCIATED TO EACH PARTICLE##########
################# #########################
# Argument: ndimz, turbulentViscosity, vxFluid, phiPart,
# Modification : turbFluct.virtPeriod, addFluidForce.bedElevation, addFluidForce.simplifiedReynoldStresses, addFluidForce.velFluct
addFluidForce.simplifiedReynoldStresses = np.zeros(ndimz)
addFluidForce.bedElevation = 0.
addFluidForce.turbulentFluctuation()
addFluidForce.fluctTime = np.zeros(len(O.bodies))
def turbulentFluctuation():
 global waterDepth
 global turbStress

 #For stability requirement at the initialization stage
 if O.time<depoTime+0.5:
  print 'No turbulent fluctuation in the initialization process for stability reasons!'
  turbFluct.virtPeriod = 0.5
 else:
  # Evaluate nBed, the position of the bed which is assumed to be located around the first maximum of concentration when considering decreasing z.
  nBed = 0.
  bedElevation = 0.
  for n in range(1,ndimz):
   # if there is a peak and its value is superior to 0.5, we consider it to be the position of the bed
   if phiPart[ndimz - n] < phiPart[ndimz - n-1] and phiPart[ndimz - n] > 0.5 :
    nBed = ndimz - n
    waterDepth = (ndimz-1 - nBed)*dz
    bedElevation = fluidHeight - waterDepth #Evaluate the bed elevation for the following
    break
  #(Re)Define the bed elevation over which fluid turbulent fluctuations will be applied.
  addFluidForce.bedElevation = bedElevation

  #First turbulent fluctuation model: a unique constant lifetime for the turbulent fluctuation, flucTimeScale
  if turbModel1:
   vMeanAboveBed = sum(vxFluid[nBed:])/(ndimz-nBed) # fluid elocity scale in the water depth
   flucTimeScale = waterDepth/vMeanAboveBed # time scale of the fluctuation w_d/v, eddy turn over time
   # New evaluation of the random fluid velocity fluctuation for each particle.
   addFluidForce.turbulentFluctuation()
   turbFluct.virtPeriod = flucTimeScale #Actualize when will be calculated the next fluctuations.

  elif turbModel2:
   # New evaluation of the random fluid velocity fluctuation for each particle.
   addFluidForce.turbulentFluctuationBIS()
   turbFluct.virtPeriod = min(addFluidForce.fluctTime)#Actualize when will be calculated the next fluct.
   addFluidForce.dtFluct = min(addFluidForce.fluctTime)#Actualize the time step of the function

  if O.time<depoTime+1: #To avoid important virtPeriod at the start, impose it for 1s. Fasten the transient.
   turbFluct.virtPeriod = 0.1
   if turbModel2: addFluidForce.dtFluct=0.1

###############
#########################################

###### ######
### CALL THE FLUID PROFILE RESOLUTION ROUTINE###
###### ######
#Arguments : ndimz, vxFluid, vxPart, dragTerm, densFluid, fluidHeight, sig, dsig, diameterPart, phiPartFluid, kinematicViscoFluid, densPart, dpdx, slope, tfin,d tFluid
#Modify : vxFluid, addFluidForce.vxFluid
testFluidResol = 0 # Test variable for the fluid resolution
def fluidResolution(tfin):
 global phiPart
 global vxPart
 global vyPart
 global vzPart
 global vxFluid
 global averageDrag
 global turbulentViscosity
 global testFluidResol
 global dragTerm
 global taufsi
 global n_dt
 global turbStress
 if twoSize:
  global averageDrag1
  global phiPart1
  global vxPart1
  global vyPart1
  global vzPart1
  global averageDrag2
  global phiPart2
  global vxPart2
  global vyPart2
  global vzPart2

 #Two different way of evaluating the average profiles: phiPart, vxPart, averageDrag, (+ same for 1 and 2 with twoSize)
 #From python
 if pythonAverage==1 or debugStress==1:
  averageProfile()
 #Or from C++
 elif timeAveraging:#with multiple time averaging made previously in timeAveragingDEMdata
  #Normalize the time averaged data
  vxPart = vxPart_dt/n_dt
  phiPart = phiPart_dt/n_dt
  averageDrag = averageDrag_dt/n_dt
  if twoSize:
   phiPart1 = phiPart1_dt/n_dt
   averageDrag1 = averageDrag1_dt/n_dt
   vxPart1 = vxPart1_dt/n_dt
   vyPart1 = vyPart1_dt/n_dt
   vzPart1 = vzPart1_dt/n_dt
   phiPart2 = phiPart2_dt/n_dt
   averageDrag2 = averageDrag2_dt/n_dt
   vxPart2 = vxPart2_dt/n_dt
   vyPart2 = vyPart2_dt/n_dt
   vzPart2 = vzPart2_dt/n_dt

  #And report the other one (for usual measurement)
  vyPart = np.array(addFluidForce.vyPart)
  vzPart = np.array(addFluidForce.vzPart)
  n_dt = 0

 else:#with simple time averaging
  if raphaelComputer!=1:
   #Evaluate the average vx,vy,vz,phi,drag profiles and store it in addFluidForce
   addFluidForce.averageProfile()
  else:#On old ubuntu/Yade version, needs to do it differently:
   #At next time step, Evaluate the average vx,vy,vz,phi,drag profiles and store it in addFluidForce
   if testFluidResol == 0:
    addFluidForce.activateAverage = True #Activate the average for next time step
    fluidResol.iterPeriod = 1 #For it to be re-executed at next time step
    testFluidResol = 1
    return

  ## Import the necessary average variables for the fluid resolution from addFluidForce : phiPart,vxPart, averageDrag
  vxPart = np.array(addFluidForce.vxPart)
  vyPart = np.array(addFluidForce.vyPart)
  vzPart = np.array(addFluidForce.vzPart)
  phiPart = np.array(addFluidForce.phiPart)
  averageDrag = np.array(addFluidForce.averageDrag)

  #Back to the usual resolution period for next step
  fluidResol.iterPeriod = -1 #For it not to be re-executed next time
  fluidResol.virtPeriod = fluidResolPeriod #Back to the normal period
  testFluidResol = 0 #Re-initialize the test

  #If two particle size in the simulation, store complementary necessary average
  if twoSize:
   vCell = dz*width*length
   #C++ function evaluating the depth profiles of solid volume fraction and solid velocities
   phiPart1=np.array(addFluidForce.phiPart1)
   phiPart2=np.array(addFluidForce.phiPart2)
   averageDrag1=np.array(addFluidForce.averageDrag1)
   averageDrag2=np.array(addFluidForce.averageDrag2)
   vxPart1=np.array(addFluidForce.vxPart1)
   vxPart2=np.array(addFluidForce.vxPart2)
   vyPart1=np.array(addFluidForce.vyPart1)
   vyPart2=np.array(addFluidForce.vyPart2)
   vzPart1=np.array(addFluidForce.vzPart1)
   vzPart2=np.array(addFluidForce.vzPart2)

 #
 ## Prepare the variables for the fluid resolution
 #

 # Threshold the solid vol. fraction in the bed to remove the oscillations due to layering and avoid turbulence in the bed
 phiPartFluid = np.zeros(ndimz)
 for n in range(1,ndimz+1):
  # If the concentration inside the bed reach the maximum concentration (imposed in the fluid solver). Then all the values of concentration under this one are going to be thresholded at phiPartMax, in order to damp the turbulence inside the bed.
  if phiPart[ndimz - n] > phiPartMax:
   for i in range(n,ndimz+1):
     phiPartFluid[ndimz - i] = phiPartMax
   break
  #Otherwise we report the normal value.
  phiPartFluid[ndimz - n] = phiPart[ndimz - n]

 #Density of particle per unit volume N/Vcell:
 nPart = phiPart[:ndimz]/partVolume

 #Calculate the momentum transfer from the drag force interaction : N<fd>/Vcell
 dragTerm = nPart*averageDrag[:ndimz]
 if twoSize: #if two size, sum the two contributions n1<f1> + n2 <f2>, formulation necessary to keep energy conservation.
  dragTerm = phiPart1[:ndimz]/partVolume1*averageDrag1[:ndimz] + phiPart2[:ndimz]/partVolume2*averageDrag2[:ndimz]

 #Create Taufsi/rhof = dragTerm/(rhof(vf-vxp)) to transmit to the fluid code
 taufsi = np.zeros(ndimz)
 for i in range(1,ndimz):
  urel = max(abs(vxFluid[i] - vxPart[i]),1e-5) #Limit the value to avoid division by 0
  #Limit the max value of taufsi to the fluid resolution limit, i.e. 1/(10dt) and required positive (charact. time)
  taufsi[i] = max(0,min(dragTerm[i]/urel/densFluid, pow(10*dtFluid,-1)))

 #
 ## Call the routine for fluid resolution
 #
 [vxFluid_np,turbulentViscosity] = nsmp1d_yade(fluidHeight,sig,dsig,diameterPart,vxFluid,1e0-phiPartFluid[0:ndimz],densFluid,kinematicViscoFluid,vxPart[0:ndimz],phiPartFluid[0:ndimz],densPart,dpdx,slope,tfin,dtFluid,taufsi)

 #Debugstress option
 if debugStress==1:
  global taufsiBIS
  taufsiBIS = np.zeros(ndimz)
  for i in range(1,ndimz):
   urel = max(abs(vxFluid[i] - vxPart[i]),1e-5)
   taufsiBIS[i] = min(dragTermBIS[i]/urel/densFluid, pow(10*dtFluid,-1))
   global vxFluid_before
   vxFluid_before = vxFluid
   vxFluid = vxFluid_np
   path = scriptPath + '/'+ str(fileNumber)+'_debug.py'
   globalParam = ['vxPart','vxPartBIS','phiPart','phiPartBIS','vxFluid','vxFluid_before', 'averageDrag','averageDragBIS','dragTerm','dragTermBIS','taufsi','taufsiBIS']
   Save(path, globalParam)

 #update
 vxFluid = vxFluid_np

 #Modify the velocity in the HydroForceEngine which apply the force
 addFluidForce.vxFluid = vxFluid

 #Evaluate the turbulent shear stress for the turbulent fluctuations model and Shields number evaluation
 turbStress = np.zeros(ndimz)
 n = 0
 for nu_t in turbulentViscosity[:-2]: # Turbulent viscosity is global, calculated in fluidResolution
  turbStress[n] = nu_t*(vxFluid[n+1]-vxFluid[n])/dz #tau/rho
  n+=1
 addFluidForce.simplifiedReynoldStresses = turbStress #Actualize the Reynolds stress tensor (/rho)value

###############
#########################################

####### ########
### OUTPUT ###
####### ########
#Function which register the different quantities asked for post processing
#Arguments : len(O.bodies), O.bodies, scriptPath, fileNumber, measurePeriod, savePeriod, qsMean, phiPart, vxPart, vxFluid, turbulentViscosity, dragTerm, waterDepth
#Modify : fileNumber, 'paramSave.py', 'save.xml'
def measure():
 global fileNumber
 global qsMean
 global stressProfile
 global vxPart
 global vyPart
 global vzPart
 global phiPart
 global averageDrag
 global phiPart1
 global phiPart2
 global averageDrag1
 global averageDrag2
 global vxPart1
 global vxPart2
 global vyPart1
 global vyPart2
 global vzPart1
 global vzPart2

 if imposedFluidPro==1 or initImposedFluidPro==1:#Then you should evaluate the average (as not done in fluidResolution)
  if pythonAverage==1 or debugStress==1:
   averageProfile()
  #Or from C++
  else:
   #Evaluate the average vx,vy,vz,phi,drag profiles and store it in addFluidForce
   addFluidForce.averageProfile()

   #Replace it in the global python variables
   vxPart = np.array(addFluidForce.vxPart)
   vyPart = np.array(addFluidForce.vyPart)
   vzPart = np.array(addFluidForce.vzPart)
   phiPart = np.array(addFluidForce.phiPart)
   averageDrag = np.array(addFluidForce.averageDrag)

   #If two particle size in the simulation, store complementary necessary average
   if twoSize:
    vCell = dz*width*length
    #C++ function evaluating the depth profiles of solid volume fraction and solid velocities
    phiPart1=np.array(addFluidForce.phiPart1)
    phiPart2=np.array(addFluidForce.phiPart2)
    averageDrag1=np.array(addFluidForce.averageDrag1)
    averageDrag2=np.array(addFluidForce.averageDrag2)
    vxPart1=np.array(addFluidForce.vxPart1)
    vxPart2=np.array(addFluidForce.vxPart2)
    vyPart1=np.array(addFluidForce.vyPart1)
    vyPart2=np.array(addFluidForce.vyPart2)
    vzPart1=np.array(addFluidForce.vzPart1)
    vzPart2=np.array(addFluidForce.vzPart2)

 # Save all the following variables/split in two file for faster post processing

 if stressSave == 1:
  if pythonStressCalculation:
   stressProfile = stressFieldCalculation(ndimz,dz*width*length,vxPart,dz)
  else:
   stressProfile = getStressProfile(dz*width*length,ndimz,dz,groundPosition,vxPart,vyPart,vzPart)
  for n in range(0,ndimz):
   stressProfile[0][n] = np.array(stressProfile[0][n]).reshape(1,9)[0]
   stressProfile[1][n] = np.array(stressProfile[1][n]).reshape(1,9)[0]
  globalParam = ['stressProfile']
  path = scriptPath + '/'+ str(fileNumber)+'_stress.py'
  Save(path, globalParam)

 #if necessary output file with the positions and velocities of the particles
 if 0:
  global positionPartX
  global positionPartZ
  global velocityPartX
  global velocityPartZ
  positionPartX = np.zeros(len(O.bodies)) #Defined global in order to use
  positionPartZ = np.zeros(len(O.bodies)) #Defined global in order to use
  velocityPartX = np.zeros(len(O.bodies)) # the saving function
  velocityPartZ = np.zeros(len(O.bodies)) # the saving function
  for b in O.bodies:
   if isinstance(b.shape,Sphere): # If it is a sphere
    positionPartX[b.id] = b.state.pos[0]
    positionPartZ[b.id] = b.state.pos[2]
    velocityPartX[b.id] = b.state.vel[0]
    velocityPartZ[b.id] = b.state.vel[2]
  path = scriptPath + '/'+ str(fileNumber)+'_2.py'
  globalParam = ['positionPartX', 'positionPartZ','velocityPartX', 'velocityPartZ']
  Save(path, globalParam)

 #Evaluate the dimensionless sediment rate to save
 qsMean = sum(phiPart*vxPart)*dz/sqrt((densPart/densFluid - 1)*abs(gravityVector[2])*pow(diameterPart,3))

 path = scriptPath + '/'+ str(fileNumber)+'.py'
 globalParam = ['qsMean','phiPart','vxPart','vyPart','vzPart','vxFluid','turbulentViscosity','dragTerm','turbStress']
 if twoSize:
  globalParam =['qsMean','phiPart','vxPart','vyPart','vzPart','vxFluid','turbulentViscosity','dragTerm','turbStress', 'phiPart1','vxPart1','vyPart1','vzPart1','phiPart2','vxPart2','vyPart2','vzPart2']
 Save(path, globalParam)
 fileNumber +=1

#PHF if twoSize!=1 and O.time>=150:
# if O.time>=150:
# print '\n End of the simulation, simulated 150s as required !\n '
# O.pause()
# O.exitNoBacktrace()

 #For reload, save sometimes the variables necessary to reload. (does not work anymore...)
 if (fileNumber*measurePeriod)%savePeriod == 0:
  #Save the simulation of yade the global python parameter in order to be able to recover the simulation
  O.save(scriptPath + '/'+'save.xml')
  # Save the global parameters which shouldn't be initialized to zero at the reloading.
  globalParam = ['phiPart','fileNumber','vxFluid','turbulentViscosity','waterDepth','turbStress']
  Save(scriptPath + '/'+'paramSave.py', globalParam)

 ##
 # Evaluation and printing of the Shields number or qsMean.

 if imposedFluidPro==1 or initImposedFluidPro==1:
  print 'qsMean', qsMean
 else:
  if twoSize:
   #Consider the average diameter in the layer in motion (which is the one to consider for Shields number)
   dummy = np.where(phiPart[:ndimz]<0.8*phiPartMax)
   phi1Sum = sum(phiPart1[dummy])
   phi2Sum = sum(phiPart2[dummy])
   #Volume weighted average diameter: N1Vp1 d1 + N2Vp2 d2/(N1Vp1 + N2Vp2) (*VCell/VCell)
   diameterMix = (diameterPart1*phi1Sum + diameterPart2*phi2Sum)/(phi1Sum + phi2Sum)
  else:
   #Only one size in the simulation...
   diameterMix = diameterPart
  shieldsNumber = densFluid*max(turbStress)/((densPart-densFluid)*diameterMix*abs(gravityVector[2]))
  print 'Shields number', shieldsNumber

###############
#########################################

########## ##########
### EVALUATE THE MEAN VOLUME FRACTION, DRAG AND STREAMWISE VELOCITY OF THE PARTICLES IN FUNCTION OF THE DEPTH###
########## ##########
#Arguments : ndimz, dz, length, Vlayer, groundPosition, radius, phiPart
#Modify : phiPart, vxPart, averageDrag
def averageProfile():
 global phiPart
 global averageDrag
 global vxPart
 global vyPart
 global vzPart
 global qsMean

 global phiPart1
 global averageDrag1
 global phiPart2
 global averageDrag2

 #Initialization of the concentration and the streamwise velocity. (two times the water height in order to take into account the particles which are also above the fluid in the concentration and velocity profile
 phiPart_loc = np.zeros(ndimz) #Mean particle concentration in function of depth
 vxPart_loc = np.zeros(ndimz) #Mean particle velocity in function of depth
 vyPart_loc = np.zeros(ndimz) #Mean particle velocity in function of depth
 vzPart_loc = np.zeros(ndimz) #Mean particle velocity in function of depth
 averageDrag_loc = np.zeros(ndimz) #Mean drag

 phiPart1_loc = np.zeros(ndimz)
 averageDrag1_loc = np.zeros(ndimz)
 phiPart2_loc = np.zeros(ndimz)
 averageDrag2_loc = np.zeros(ndimz)
 if debugStress==1:
  global vxPartBIS
  global phiPartBIS
  global averageDragBIS
  global dragTermBIS
  averageDragBIS_loc = np.zeros(ndimz)
  phiPartBIS_loc = np.zeros(ndimz)
  vxPartBIS_loc = np.zeros(ndimz)
  dragTermBIS_loc = np.zeros(ndimz)

 #Import the vector of velocity fluctuation associated to each particle in C++ code
 VxFluct = addFluidForce.vFluctX
 VzFluct = addFluidForce.vFluctZ

 shiftPosition = 1.2*diameterPart #To shift the position in order to avoid to have negative position (not compatible with the volume calculation used later)

 ## Loop over the file to determine the concentration and the mean velocity
 for b in O.bodies:
  if isinstance(b.shape,Sphere): # If it is a sphere
   z = b.state.pos[2]-groundPosition # Altitude of the particle redefined with the position of the bottom wall as the
   radius = b.shape.radius
   Np = int(z/dz) #Layer corresponding to the position of the center of the particle

   #Evaluate the drag force undergone by the particle (Dallavalle formulation + Richardson Zaki correction)
   uRel = Vector3(0.0,.0,0.0)
   if Np<ndimz and Np>=0 : uRel = Vector3(vxFluid[Np]+VxFluct[b.id],0.0,VzFluct[b.id]) - b.state.vel
   fDrag = 0.5*pi*pow(radius,2)*densFluid*(0.44*uRel.norm()+24.4*kinematicViscoFluid/(2*radius))*pow(1-phiPart[Np],-expoRZ)*uRel

   if debugStress==1:
    dragTermBIS_loc[Np]+=fDrag[0]
    phiPartBIS_loc[Np]+=1.0
    vxPartBIS_loc[Np]+=b.state.vel[0]

   #Shift the position to avoid z<0 which induces some problems in the volume calculation. This is done here in order not to change the value of the drag. Evaluate all the parameters after this shift. (Np in particular is reevaluated
   z+=shiftPosition
   Np = int(z/dz) #new evaluation of Np after the shift. DO NOT REMOVE, IMPORTANT
   minZ = int((z-radius)/dz) #Lowest layer occupied by the particle
   maxZ = int((z+radius)/dz) #Highest layer occupied by the particle
   deltaCenter = z - Np*dz #length between the real center of the particle and the position of the layer in which it is considered to be (Np).

   #loop from the bottom to the top layer inside which the particle is contained
   numLayer = minZ
   while numLayer <= maxZ:
    if numLayer>=0 and numLayer<ndimz:
     zInf = (numLayer-Np-1)*dz + deltaCenter
     zSup = (numLayer-Np)*dz + deltaCenter

     if numLayer==minZ: zInf = -radius
     if numLayer==maxZ: zSup = radius

     #evaluate the volume of the part of the particle contained in the layer considered (the formula results from the analytical resolution of the integral of a layer of sphere in cylindrical coordinate)
     V = math.pi*radius*radius*(zSup-zInf +(pow(zInf,3)-pow(zSup,3))/(3*radius*radius))

     #Add the quantities to the (not normalized for now) mean particle velocity and concentration
     phiPart_loc[numLayer]+=V
     vxPart_loc[numLayer]+=V*b.state.vel[0]
     vyPart_loc[numLayer]+=V*b.state.vel[1]
     vzPart_loc[numLayer]+=V*b.state.vel[2]
     averageDrag_loc[numLayer]+=V*fDrag[0]
     if twoSize:
      if radius==diameterPart1/2.0:
       phiPart1_loc[numLayer]+=V
       averageDrag1_loc[numLayer]+=V*fDrag[0]
      elif radius==diameterPart2/2.0:
       phiPart2_loc[numLayer]+=V
       averageDrag2_loc[numLayer]+=V*fDrag[0]
     #Increment for the next loop
    numLayer+=1

 #Normalize the weighted velocity/different mean quantities of each layer by the total particle volume in this layer
 n = 0
 for i in phiPart_loc:
  if i: #if non zero, to avoid division by zero
   vxPart_loc[n]/=i
   vyPart_loc[n]/=i
   vzPart_loc[n]/=i
   averageDrag_loc[n]/=i
   if twoSize:
    if phiPart1_loc[n]:
     averageDrag1_loc[n]/=phiPart1_loc[n]
    if phiPart2_loc[n]:
     averageDrag2_loc[n]/=phiPart2_loc[n]

  n += 1
 #Normalization of the concentration
 phiPart_loc/=VLayer
 phiPart1_loc/=VLayer
 phiPart2_loc/=VLayer

 if debugStress==1:
  #Normalization
  n = 0
  for i in phiPartBIS_loc:#Loop over the number of part contained in each layer
   #Normalize the average drag and particle velocity
   if i:
    averageDragBIS_loc[n] = dragTermBIS_loc[n]/i
    vxPartBIS_loc[n]/=i
   n+=1
  phiPartBIS_loc*=(partVolume/VLayer)
  dragTermBIS_loc/=VLayer
  #Report the value in global variables
  vxPartBIS = vxPartBIS_loc
  phiPartBIS = phiPartBIS_loc
  averageDragBIS = averageDragBIS_loc
  dragTermBIS = dragTermBIS_loc

 #Resize the vector to remove the effect of the shift previously introduced
 phiPart_loc = Resize(phiPart_loc,int(shiftPosition/dz))
 vxPart_loc = Resize(vxPart_loc,int(shiftPosition/dz))
 vyPart_loc = Resize(vyPart_loc,int(shiftPosition/dz))
 vzPart_loc = Resize(vzPart_loc,int(shiftPosition/dz))
 averageDrag_loc = Resize(averageDrag_loc,int(shiftPosition/dz))

 phiPart1_loc = Resize(phiPart1_loc,int(shiftPosition/dz))
 averageDrag1_loc = Resize(averageDrag1_loc,int(shiftPosition/dz))
 phiPart2_loc = Resize(phiPart2_loc,int(shiftPosition/dz))
 averageDrag2_loc = Resize(averageDrag2_loc,int(shiftPosition/dz))
 #Impose variables at the bottom to be consistent with the fluid formulation
 averageDrag_loc[0] = 0.0
 vxPart_loc[0] = 0.0

 ###Actualize the global variables
 phiPart = phiPart_loc
 vxPart = vxPart_loc
 vyPart = vyPart_loc
 vzPart = vzPart_loc
 averageDrag = averageDrag_loc

 phiPart1 = phiPart1_loc
 averageDrag1 = averageDrag1_loc
 phiPart2 = phiPart2_loc
 averageDrag2 = averageDrag2_loc
 #And the one associated to the fluid force application
 addFluidForce.phiPart = phiPart_loc

 return
###############
#########################################

###### ######
## Function to resize the vector removing the n first case##
###### ######
def Resize(a,n):
 import numpy as np
 l = len(a)
 b = np.zeros(l)
 k = 0
 for i in range(n,l):
  b[k] = a[i]
  k+=1
 return b
###############
#########################################

###### ######
### Save global variables ###
###### ######

#Function to save global variables in a python file which can be re-executed for later
def Save(filePathName, globalVariables):
 f = open(filePathName,'w')
 f.write('from numpy import *\n')
 for i in globalVariables:
  f.write(i + ' = '+repr(globals()[i]) + '\n')
 f.close()
###############
#########################################

###### ######
### python Stress field calculation ###
###### ######
def stressFieldCalculation(ndimz,vCell,vxPart,dz):
 #Initialization
 granularTemp_loc = np.zeros(ndimz)
 N_loc = np.zeros(ndimz)
 stressTensorField_loc = []
 kineticStressTensorField_loc = []
 for i in range(0,ndimz):
  stressTensorField_loc += [Matrix3(0,0,0,0,0,0,0,0,0)]
  kineticStressTensorField_loc += [Matrix3(0,0,0,0,0,0,0,0,0)]
 ###
 ## Dynamical stress tensor part
 ###
 for b in O.bodies:
  Np = int((b.state.pos[2] - groundPosition)/dz)
  if Np>=0 and Np<ndimz:
   vFluct = b.state.vel - Vector3(vxPart[Np],vyPart[Np],vzPart[Np])
   stressTensorField_loc[Np]+= -1/vCell*b.state.mass*vFluct.outer(vFluct)
   kineticStressTensorField_loc[Np]+= -1/vCell*b.state.mass*vFluct.outer(vFluct)
   granularTemp_loc[Np] += 1/3.0*(pow(vFluct[0],2) + pow(vFluct[1],2) + pow(vFluct[2],2))
   N_loc[Np]+=1.0
 ###
 ## Love-Weber stress tensor part (contact
 ###
 for Int in O.interactions:
  if O.bodies[Int.id1].dynamic or O.bodies[Int.id2].dynamic:
   Np1 = int((O.bodies[Int.id1].state.pos[2] - groundPosition)/dz)
   Np2 = int((O.bodies[Int.id2].state.pos[2] - groundPosition)/dz)
   #Vector joining the two particles (from id2 to id1)
   x12 = O.bodies[Int.id1].state.pos - O.bodies[Int.id2].state.pos
   #Contact vector of the two particles (defined as applied by id1 on id2)
   fContact = Int.phys.normalForce + Int.phys.shearForce
   #Remove the artificial effect of the periodic BC if one particle goes to the other side
   x12 -= O.cell.hSize*Int.cellDist
   #No reason to have one upon another, so check and define the min and max correpondingly
   if Np2==Np1: #If the contact vector is entirely inside the layer, add it directly.
    if Np1>=0 and Np1<ndimz:
     stressTensorField_loc[Np1]+= 1/vCell*fContact.outer(x12) #outer product of contact force and branch vector
   else:
    if Np1>Np2: #Particle 1 above 2 (z component)
     minZ = Np2
     minPosZ = O.bodies[Int.id2].state.pos[2] - groundPosition
     maxZ = Np1
     maxPosZ = O.bodies[Int.id1].state.pos[2] - groundPosition
    elif Np2>Np1: #Particle 2 above 1 (z component)
     minZ = Np1
     minPosZ = O.bodies[Int.id1].state.pos[2] - groundPosition
     maxZ = Np2
     maxPosZ = O.bodies[Int.id2].state.pos[2] - groundPosition

    #Loop over the layers containing the vector connecting the two particles
    numLayer = minZ
    while numLayer<=maxZ:
     if numLayer>=0 and numLayer<ndimz:
      deltaZ = dz
      if numLayer==minZ:
                                                       deltaZ = dz - (minPosZ - minZ*dz)
      elif numLayer==maxZ:
                                                       deltaZ = maxPosZ - maxZ*dz

      branchVectCell=deltaZ/abs(x12[2])*Vector3(x12[0],x12[1],x12[2])
      #Add to stress tensor of the considered layer the part considered
                                                stressTensorField_loc[numLayer]+= 1/vCell*fContact.outer(branchVectCell) #outer product
     #Increment the layer
     numLayer+=1
 #Normalization
 for n in range(0,len(N_loc)):
  if N_loc[n]>0:
   granularTemp_loc[n]/=N_loc[n]

 return [stressTensorField_loc,kineticStressTensorField_loc,granularTemp_loc]

###############
#########################################

###DEBUG/COMPARISON WITH getStress
#To obtain exactly the same as getStress, you need to shift groundPosition -> groundPosition - diameterPart + to remove the condition on the dynamics of the particles, to avoid the non contribution of the bottom wall and particles.
#This function might be also useful, it is reproducing exactly getStress in python
def getStressBIS(vCell):
 stressTensorField = Matrix3(0,0,0,0,0,0,0,0,0)
 for Int in O.interactions:
  x12 = O.bodies[Int.id1].state.pos - O.bodies[Int.id2].state.pos
  x12 -= O.cell.hSize*Int.cellDist
  fContact = Int.phys.normalForce + Int.phys.shearForce
  stressTensorField+= fContact.outer(x12) #outer product of contact force a
 return stressTensorField/vCell
###############
#########################################

Hi,
There is a chance that the problem is just not related to yade.
Could you run try this code as a python script and tell us?
B

# testImport.py
import math
import random as rand
import numpy as np
np.set_printoptions(threshold=numpy.nan)
import yade.timing
import pylab as pyl
from matplotlib import pyplot as pyp
import matplotlib.gridspec as gridspec
from nsmp1d_yade3D525Phi061 import *

hi I have ran the code as python script and got following error,

# testImport.py
import math
import random as rand
import numpy as np
np.set_printoptions(threshold=numpy.nan)
import yade.timing
import pylab as pyl
from matplotlib import pyplot as pyp
import matplotlib.gridspec as gridspec
from nsmp1d_yade3D525Phi061 import *

Which is the error?
ps. "import yade.timing" should not be in the list, sorry.

B

hi the error says,

Welcome to Yade 2017.01a-11-11c276f~xenial
TCP python prompt on localhost:9000, auth cookie `kceydu'
XMLRPC info provider on http://localhost:21000
Running script test.py
Traceback (most recent call last):
  File "/usr/bin/yadedaily", line 182, in runScript
    execfile(script,globals())
  File "test.py", line 10, in <module>
    from nsmp1d_yade3D525Phi061 import *
ImportError: No module named nsmp1d_yade3D525Phi061

Please run it with python, not with Yade.
B

hi, It says,

Traceback (most recent call last):
  File "test.py", line 5, in <module>
    np.set_printoptions(threshold=numpy.nan)
NameError: name 'numpy' is not defined

Consider spending some time reading through the python tutorial [1].

The line should be corrected to:

np.set_printoptions(threshold=np.nan)

since we imported numpy as np.

[1]https://docs.python.org/2.7/tutorial/