Batch file question

Asked by Huan

Hello,

Since my particles is almost 5 million, so my professor suggest that I use batch file to run it. But, when I read the file of Setting up a simulation from this "https://yade-dem.org/doc/tutorial-geo.html?highlight=batch." I'm couldn't understand what exactly is readParamsFromTable(damping =.2) and where should I create it.

##This is my code##
import random
import math
from yade import geom, pack, utils, plot, ymport

# Define material properties
youngModulus = 1e7
poissonRatio = 0.25
density = 2700
frictionAngle =0.7

# Create material
material = O.materials.append(FrictMat(young=youngModulus, poisson=poissonRatio, density=density))
material1 = O.materials.append(FrictMat(young=youngModulus, poisson=poissonRatio, density=density, frictionAngle=frictionAngle))

# Define cylinder with funnel parameters
center = (0, 0, 0)
diameter = 0.102
height = 0.18

# create cylindrical body with radius 0.102 m and height 0.064 m
cylinder = geom.facetCylinder(center=center, radius=diameter/2, height=height, segmentsNumber=80, wallMask=6)

# assign material to each body in the cylinder
for body in cylinder:
    body.bodyMat = material

# add cylinder to simulation
O.bodies.append(cylinder)

# Define cylinder with funnel parameters
center1 = (0,0,height/2)
dBunker = 0.4
dOutput = 0.102
hBunker = 0
hOutput = 0.15
hPipe = 0

# create funnel as a bunker with diameter 0.102 m, height 0.064 m
funnel = geom.facetBunker(center=center1, dBunker=dBunker, dOutput=dOutput, hBunker=hBunker,hOutput=hOutput, hPipe=hPipe, segmentsNumber=80, wallMask=4)

# assign material to each body in the funnel
for body in funnel:
    body.bodyMat = material

# add funnel to simulation
O.bodies.append(funnel)

# define sphere parameters and number of spheres
rMean1 = (0.0125+0.019)/4
rRelFuzz1 = (0.019-0.0125)/4/rMean1
num1 = 13
rMean2 = (0.0095+0.0125)/4
rRelFuzz2 = (0.0125-0.0095)/4/rMean2
num2 = 51
rMean3 = (0.00475+0.0095)/4
rRelFuzz3 = (0.0095-0.00475)/4/rMean3
num3 = 563
rMean4 = (0.00236+0.00475)/4
rRelFuzz4 = (0.00475-0.00236)/4/rMean4
num4 = 5101
rMean5 = (0.00118+0.00236)/4
rRelFuzz5 = (0.00236-0.00118)/4/rMean4
num5 = 18369
rMean6 = 0.00118/4
num6 = 950584
rMean7 = 0.0003/4
num7 = 4003898

#create empty sphere packing
sp = pack.SpherePack()

# generate randomly sphere
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean1, rRelFuzz = rRelFuzz1, num = num1)
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean2, rRelFuzz = rRelFuzz2, num = num2)
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean3, rRelFuzz = rRelFuzz3, num = num3)
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean4, rRelFuzz = rRelFuzz4, num = num4)
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean5, rRelFuzz = rRelFuzz5, num = num5)
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean6, num = num6)
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean7, num = num7)

# add the sphere pack to the simulation
sp.toSimulation(material = material1)

for body in O.bodies:
   if not isinstance(body.shape, Sphere):
       continue
   if body.shape.radius >= rMean1 :
       body.shape.color = (0,0,1) #blue
   if body.shape.radius <= rMean1 and body.shape.radius > rMean2:
       body.shape.color = (0,0,1) #blue
   if body.shape.radius <= rMean2 and body.shape.radius > rMean3:
       body.shape.color = (1,0,0) #red
   if body.shape.radius <= rMean3 and body.shape.radius > rMean4:
       body.shape.color = (0,1,0) #green
   if body.shape.radius <= rMean4 and body.shape.radius > rMean5:
       body.shape.color = (1,1,0) #yellow
   if body.shape.radius <= rMean5 :
       body.shape.color = (1,0,1) #magenta
   if body.shape.radius <= rMean6 :
       body.shape.color = (1,1,1) #white
   if body.shape.radius <= rMean7 :
       body.shape.color = (0,0,0) #black

O.engines = [
        ForceResetter(),
        # sphere, facet, wall
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(), Bo1_Wall_Aabb()]),
        InteractionLoop(
                # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -1000), damping=0.3),
        # the label creates an automatic variable referring to this engine
        # we use it below to change its attributes from the functions called
        PyRunner(command='checkUnbalanced()', realPeriod=2, label='checker'),
]
O.dt = PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy = True

# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time
# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
 # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
 if O.iter < 25000:
  return
 # add plate at the position on the top of the packing
 # the maximum finds the z-coordinate of the top of the topmost particle
 O.bodies.append(wall(max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]), axis=2, sense=-1))
 global plate # without this line, the plate variable would only exist inside this function
 plate = O.bodies[-1] # the last particles is the plate
 # Wall objects are "fixed" by default, i.e. not subject to forces
 # prescribing a velocity will therefore make it move at constant velocity (downwards)
 plate.state.vel = (0, 0, -0.8)
 # start plotting the data now, it was not interesting before
 O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=200)]
 # next time, do not call this function anymore, but the next one (unloadPlate) instead
 checker.command = 'unloadPlate()'

def unloadPlate():
 # if the force on plate exceeds maximum load, start unloading
 if abs(O.forces.f(plate.id)[2]) > 1e3:
  plate.state.vel *= -0.8
  # next time, do not call this function anymore, but the next one (stopUnloading) instead
  checker.command = 'stopUnloading()'

def stopUnloading():
    if abs(O.forces.f(plate.id)[2]) < 1e2:
        # calculate the volume of the packing
        volume_packing = sum(4/3 * math.pi * b.shape.radius**3 for b in O.bodies if isinstance(b.shape, Sphere))
        num_spheres = len([b for b in O.bodies if isinstance(b.shape,Sphere)])
        # print the number of spheres and volume of packking
        print("Number of spheres:", "{:d}".format(num_spheres))
        print("V Packing:", "{:e}".format(volume_packing))
        # print the height of the plate
        top = abs(plate.state.pos[2] + height/2)
        print("Plate height:", top)
        # calculate the volume of the cylinder
        new_volume_cylinder = math.pi * (diameter/2)**2 * top
        print("V Cylinder:", "{:e}".format(new_volume_cylinder))
        # calculate the porosity and porosity percentage
        new_porosity = (new_volume_cylinder - volume_packing) / new_volume_cylinder
        new_porosity_percent = new_porosity * 100
        print("Porosity:", "{:.2f}".format(new_porosity))
        print("Porosity:", "{:.2f}%".format(new_porosity_percent))
        # O.tags can be used to retrieve unique identifiers of the simulation
        # if running in batch, subsequent simulation would overwrite each other's output files otherwise
        # d (or description) is simulation description (composed of parameter values)
        # while the id is composed of time and process number
        plot.saveDataTxt(O.tags['d.id'] + '.txt')
        O.pause()

def addPlotData():
 if not isinstance(O.bodies[-1].shape, Wall):
  plot.addData()
  return
 Fz = O.forces.f(plate.id)[2]
 plot.addData(Fz=Fz, w=plate.state.pos[2] - plate.state.refPos[2], unbalanced=unbalancedForce(), i=O.iter)

# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots = {'i': ('unbalanced',), 'w': ('Fz',)}
plot.plot()

Does using batch file mean I can divide the simulation into 2 parts?
1. Gravity deposition and calculate the value of porosity (loose pack)
2. Oedometric test and calculate the value of porosity (compacted pack)

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said :
#1

Hello,

see also [1]

> Does using batch file mean I can divide the simulation into 2 parts?

No.
Batch run in Yade context means that you run the same script multiple times with given (usually varying) parameters.
E.g. you could run your simulation with varying value of gravity (artificial example) with values 9.81, 20, 100 and 500.

Batch run can be run all jobs sequentially or all jobs in parallel or a combination of both.
The main advantage is simplicity of the setting.

> I'm couldn't understand what exactly is readParamsFromTable(damping =.2)

readParamsFromTable reads parameters from table if given, otherwise it uses default value (damping=0.2 in this case)

For the gravity example above, you can create the "grav.table" file:
### grav.table
gravity
9.81
20
100
500
###

### script.py
readParamsFromTable(
    gravity = 9.81,
)
from yade.params.table import *
print(gravity)
###

and then run it as
yade-batch grav.table script.py

Cheers
Jan

[1] https://yade-dem.org/doc/user.html#batch-queuing-and-execution-yade-batch

Revision history for this message
Huan (huan-liu) said :
#2

Hi,

I tried your example and it created 4 log file.

But I try to do mine by creating a oedometric.table with the following:
### oedometric.table
rMean rRelFuzz Num maxLoad
(0.0125+0.019)/4 (0.019-0.0125)/4/(0.0125+0.019)/4 13 1e3
(0.0095+0.0125)/4 (0.0125-0.0095)/4/(0.0095+0.0125)/4 51 1e3
(0.00475+0.0095)/4 (0.0095-0.00475)/4/(0.00475+0.0095)/4 563 1e3
(0.00236+0.00475)/4 (0.00475-0.00236)/4/(0.00236+0.00475)/4 5101 1e3
(0.00118+0.00236)/4 (0.00236-0.00118)/4/(0.00118+0.00236)/4 18369 1e3
0.00118/4 0 950584 1e3
0.0003/4 0 40038981e3
###

for the following code name cm2b.py:
import random
import math
from yade import geom, pack, utils, plot, ymport

# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(rMean=rMean, rRelFuzz=rRelFuzz,num=Num maxLoad=1e3, minLoad=1e2)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *
print(rMean, rRelFuzz, Num, maxLoad)

# Define material properties
youngModulus = 1e7
poissonRatio = 0.25
density = 2700
frictionAngle =0.7

# Create material
material = O.materials.append(FrictMat(young=youngModulus, poisson=poissonRatio, density=density))
material1 = O.materials.append(FrictMat(young=youngModulus, poisson=poissonRatio, density=density, frictionAngle=frictionAngle))

# Define cylinder with funnel parameters
center = (0, 0, 0)
diameter = 0.102
height = 0.18

# create cylindrical body with radius 0.102 m and height 0.064 m
cylinder = geom.facetCylinder(center=center, radius=diameter/2, height=height, segmentsNumber=80, wallMask=6)

# assign material to each body in the cylinder
for body in cylinder:
    body.bodyMat = material

# add cylinder to simulation
O.bodies.append(cylinder)

# Define cylinder with funnel parameters
center1 = (0,0,height/2)
dBunker = 0.4
dOutput = 0.102
hBunker = 0
hOutput = 0.15
hPipe = 0

# create funnel as a bunker with diameter 0.102 m, height 0.064 m
funnel = geom.facetBunker(center=center1, dBunker=dBunker, dOutput=dOutput, hBunker=hBunker,hOutput=hOutput, hPipe=hPipe, segmentsNumber=80, wallMask=4)

# assign material to each body in the funnel
for body in funnel:
    body.bodyMat = material

# add funnel to simulation
O.bodies.append(funnel)

#create empty sphere packing
sp = pack.SpherePack()

# generate randomly sphere
sp.makeCloud((-dBunker/4,-dBunker/4,1.3*height),(dBunker/4,dBunker/4,2*height), rMean = rMean, rRelFuzz = rRelFuzz, num = num)

# add the sphere pack to the simulation
sp.toSimulation(material = material1)

O.engines = [
        ForceResetter(),
        # sphere, facet, wall
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(), Bo1_Wall_Aabb()]),
        InteractionLoop(
                # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -1000), damping=0.3),
        # the label creates an automatic variable referring to this engine
        # we use it below to change its attributes from the functions called
        PyRunner(command='checkUnbalanced()', realPeriod=2, label='checker'),
]
O.dt = PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.rMean rRelFuzz maxLoadenergy['energyName'] subsequently
O.trackEnergy = True

# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time
# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
 # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
 if O.iter < 25000:
  return
 # add plate at the position on the top of the packing
 # the maximum finds the z-coordinate of the top of the topmost particle
 O.bodies.append(wall(max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]), axis=2, sense=-1))
 global plate # without this line, the plate variable would only exist inside this function
 plate = O.bodies[-1] # the last particles is the plate
 # Wall objects are "fixed" by default, i.e. not subject to forces
 # prescribing a velocity will therefore make it move at constant velocity (downwards)
 plate.state.vel = (0, 0, -0.8)
 # start plotting the data now, it was not interesting before
 O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=200)]
 # next time, do not call this function anymore, but the next one (unloadPlate) instead
 checker.command = 'unloadPlate()'

def unloadPlate():
 # if the force on plate exceeds maximum load, start unloading
 if abs(O.forces.f(plate.id)[2]) > 1e3:
  plate.state.vel *= -0.8
  # next time, do not call this function anymore, but the next one (stopUnloading) instead
  checker.command = 'stopUnloading()'

def stopUnloading():
    if abs(O.forces.f(plate.id)[2]) < 1e2:
        # calculate the volume of the packing
        volume_packing = sum(4/3 * math.pi * b.shape.radius**3 for b in O.bodies if isinstance(b.shape, Sphere))
        num_spheres = len([b for b in O.bodies if isinstance(b.shape,Sphere)])
        # print the number of spheres and volume of packking
        print("Number of spheres:", "{:d}".format(num_spheres))
        print("V Packing:", "{:e}".format(volume_packing))
        # print the height of the plate
        top = abs(plate.state.pos[2] + height/2)
        print("Plate height:", top)
        # calculate the volume of the cylinder
        new_volume_cylinder = math.pi * (diameter/2)**2 * top
        print("V Cylinder:", "{:e}".format(new_volume_cylinder))
        # calculate the porosity and porosity percentage
        new_porosity = (new_volume_cylinder - volume_packing) / new_volume_cylinder
        new_porosity_percent = new_porosity * 100
        print("Porosity:", "{:.2f}".format(new_porosity))
        print("Porosity:", "{:.2f}%".format(new_porosity_percent))
        # O.tags can be used to retrieve unique identifiers of the simulation
        # if running in batch, subsequent simulation would overwrite each other's output files otherwise
        # d (or description) is simulation description (composed of parameter values)
        # while the id is composed of time and process number
        plot.saveDataTxt(O.tags['d.id'] + '.txt')
        O.pause()

def addPlotData():
 if not isinstance(O.bodies[-1].shape, Wall):
  plot.addData()
  return
 Fz = O.forces.f(plate.id)[2]
 plot.addData(Fz=Fz, w=plate.state.pos[2] - plate.state.refPos[2], unbalanced=unbalancedForce(), i=O.iter)

# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots = {'i': ('unbalanced',), 'w': ('Fz',)}
plot.plot()

Run the following command in terminal:
yade-batch oedometric.table cm2b.py
And get error:
FileNotFoundError: [Errno 2] No such file or directory: 'cm2b.py.rMean=(0.00118+0.00236)/4,rRelFuzz=(0.00236-0.00118)/4/(0.00118+0.00236)/4,Num=18369,maxLoad=1e6.log'

So I don't understand which part did I do wrong. But if I do batch would I still get the value of porosity?

Revision history for this message
Jan Stránský (honzik) said :
#3

> (0.0125+0.019)/4 (0.019-0.0125)/4/(0.0125+0.019)/4 13 1e3
> FileNotFoundError: [Errno 2] No such file or directory: 'cm2b.py.rMean=(0.00118+0.00236)/4,rRelFuzz=(0.00236-0.00118)/4/(0.00118+0.00236)/4,Num=18369,maxLoad=1e6.log'

Do not use slash "/" (linux path separator) in table file.
Log files are named after values, therefore pointing to non-existent file location.
instead of /4 use *0.25

Alternatively use "description" column to name the log files [1]

> 0.0003/4 0 40038981e3

are you sure with this last one? Shouldn't there be a space before 1e3?

> But if I do batch would I still get the value of porosity?

Yes, porosity will be computed in each run.
If you print it, it will be recorded in the log file.
Another option is to save the values to a different file.

Cheers
Jan

[1] https://yade-dem.org/doc/user.html#batch-queuing-and-execution-yade-batch

Revision history for this message
Huan (huan-liu) said :
#4

Ah the slash was the problem. So, I modify the table to:

### oedometric.table
rMean rRelFuzz num maxLoad minLoad
0.007875 0.2063 13 1e3 1e2
0.0055 0.1364 51 1e3 1e2
0.0035625 0.3333 563 1e3 1e2
0.0017775 0.3361 5101 1e3 1e2
0.000085 0.3333 18369 1e3 1e2
0.0000295 0 950584 1e3 1e2
0.0000075 0 4003898 1e3 1e2
###

the log files are built correctly, but I get this error in the log files:
>readParamsFromTable(rMean=rMean, rRelFuzz=rRelFuzz,num=num, maxLoad=1e3, minLoad=1e2)
>NameError: name 'rMean' is not defined

I did a typo and I did fix the spacing part

If I understand correctly, I assume that since I print the value in the stopUnloading part if it run properly it would print the porosity into the log files

Revision history for this message
Jan Stránský (honzik) said :
#5

> readParamsFromTable(rMean=rMean, rRelFuzz=rRelFuzz,num=num, maxLoad=1e3, minLoad=1e2)
> NameError: name 'rMean' is not defined

Please provide a complete script.
But as the error says, rMean is not defined. I.e. the variable rMean is first used at this line, but there is nothing about it before in the script.
I guess similar situation would be with rRelFuzz and num.

You have to either use a literal value:
###
readParamsFromTable(rMean=123,...)
###

or use defined variable, e.g.:
###
rMean = 123
readParamsFromTable(rMean=rMean,...)
from yade.params.table import *
###
but then rMean value is that from readParamsFromTable and the first line does not make sense.. e.g.
###
rMean = 123
readParamsFromTable(rMean=234,...)
from yade.params.table import *
print(rMean) # 234
###

The way I used to use it is:
###
readParamsFromTable(
    # default values
    rMean = 123,
    ...
)
from yade.params.table import *
###
now you can use rMean variable (imported from yade.params.table) which is either the value given from table in batch run or the specified default value in case of plain yade run or if the value in table is not given.

> If I understand correctly, I assume that since I print the value in the stopUnloading part if it run properly it would print the porosity into the log files

Yes

Cheers
Jan

Can you help with this problem?

Provide an answer of your own, or ask Huan for more information if necessary.

To post a message you must log in.