\( \usepackage{euler} % must be loaded before fontspec for the whole doc (below); this must be kept for pngmath, however \usepackage{hyperref} \usepackage{amsmath} \usepackage{amsbsy} %\usepackage{mathabx} \usepackage{underscore} \usepackage[all]{xy} % Metadata of the pdf output \hypersetup{pdftitle={Yade Documentation}} \hypersetup{pdfauthor={Smilauer, V., et al.}} \hypersetup{pdfkeywords={Discrete element method, granular materials, contact law, hydromechanical coupling}} % symbols \renewcommand{\mat}{\boldsymbol} \renewcommand{\vec}{\boldsymbol} \renewcommand{\tens}{\boldsymbol} %normalized and locfram here \def\normalized#1{\widehat{#1}} \def\locframe#1{\widetilde{#1}} % timestep \def\Dt{\Delta t} \def\Dtcr{\Dt_{\rm cr}} % algorithm complexity \def\bigO#1{{\mathcal{O}(#1)}} % variants for greek symbols \let\epsilon\varepsilon \let\theta\vartheta \let\phi\varphi % shorthands \let\sig\sigma \let\eps\epsilon % variables at different points of time \def\prev#1{#1^-} \def\pprev#1{#1^\ominus} \def\curr#1{#1^{\circ}} \def\nnext#1{#1^\oplus} \def\next#1{#1^+} % shorthands for geometry \def\currn{\curr{\vec{n}}} \def\currC{\curr{\vec{C}}} \def\uT{\vec{u}_T} \def\curruT{\curr{\vec{u}}_T} \def\prevuT{\prev{\vec{u}}_T} \def\currn{\curr{\vec{n}}} \def\prevn{\prev{\vec{n}}} % motion \def\pprevvel{\pprev{\dot{\vec{u}}}} \def\nnextvel{\nnext{\dot{\vec{u}}}} \def\curraccel{\curr{\ddot{\vec{u}}}} \def\prevpos{\prev{\vec{u}}} \def\currpos{\curr{\vec{u}}} \def\nextpos{\next{\vec{u}}} \def\curraaccel{\curr{\dot{\vec{\omega}}}} \def\pprevangvel{\pprev{\vec{\omega}}} \def\nnextangvel{\nnext{\vec{\omega}}} \def\loccurr#1{\curr{\locframe{#1}}} \def\numCPU{n_{\rm cpu}} \DeclareMathOperator{\Align}{Align} \DeclareMathOperator{\sign}{sgn} % sorting algorithms \def\isleq#1{\currelem{#1}\ar@/^/[ll]^{\leq}} \def\isnleq#1{\currelem{#1}\ar@/^/[ll]^{\not\leq}} \def\currelem#1{\fbox{$#1$}} \def\sortSep{||} \def\sortInv{\hbox{\phantom{||}}} \def\sortlines#1{\xymatrix@=3pt{#1}} \def\crossBound{||\mkern-18mu<} \)

# Hands-on

## Shell basics

### Directory tree

Directory tree is hierarchical way to organize files in operating systems. A typical (reduced) tree in linux looks like this:

    /            Root
    ├──boot        System startup
    ├──bin         Low-level programs
    ├──lib         Low-level libraries
    ├──dev         Hardware access
    ├──sbin        Administration programs
    ├──proc        System information
    ├──var         Files modified by system services
    ├──root        Root (administrator) home directory
    ├──etc         Configuration files
    ├──media       External drives
    ├──tmp         Temporary files
    ├──usr         Everything for normal operation (usr = UNIX system resources)
    │    ├──bin       User programs
    │    ├──sbin      Administration programs
    │    ├──include   Header files for c/c++
    │    ├──lib       Libraries
    │    ├──local     Locally installed software
    │    └──doc       Documentation
    └──home        Contains the user's home directories
         ├──user      Home directory for user
         └──user1     Home directory for user1

Note that there is a single root `/`; all other disks (such as USB sticks) attach to some point in the tree (e.g. in `/media`).

### Shell navigation

Shell is the UNIX command-line, interface for conversation with the machine. Don’t be afraid.

#### Moving around

The shell is always operated by some `user`, at some concrete `machine`; these two are constant. We can move in the directory structure, and the current place where we are is current directory. By default, it is the home directory which contains all files belonging to the respective user:

    user@machine:~$                         # user operating at machine, in the directory ~ (= user's home directory)
    user@machine:~$ ls .                    # list contents of the current directory
    user@machine:~$ ls foo                  # list contents of directory foo, relative to the dcurrent directory ~ (= ls ~/foo = ls /home/user/foo)
    user@machine:~$ ls /tmp                 # list contents of /tmp
    user@machine:~$ cd foo                  # change directory to foo
    user@machine:~/foo$ ls ~                # list home directory (= ls /home/user)
    user@machine:~/foo$ cd bar              # change to bar (= cd ~/foo/bar)
    user@machine:~/foo/bar$ cd ../../foo2   # go to the parent directory twice, then to foo2 (cd ~/foo/bar/../../foo2 = cd ~/foo2 = cd /home/user/foo2)
    user@machine:~/foo2$ cd                 # go to the home directory (= ls ~ = ls /home/user)
    user@machine:~$

Users typically have only permissions to write (i.e. modify files) only in their home directory (abbreviated `~`, usually is `/home/user`) and `/tmp`, and permissions to read files in most other parts of the system:

    user@machine:~$ ls /root    # see what files the administrator has
    ls: cannot open directory /root: Permission denied

#### Keys

Useful keys on the command-line are:

<tab> | show possible completions of what is being typed (use abundantly!)
^C (=Ctrl+C) | delete current line
^D | exit the shell
↑↓ | move up and down in the command history
^C | interrupt currently running program
^\ | kill currently running program
Shift-PgUp | scroll the screen up (show past output)
Shift-PgDown | scroll the screen down (show future output; works only on quantum computers)

#### Running programs

When a program is being run (without giving its full path), several directories are searched for program of that name; those directories are given by `$PATH`:

    user@machine:~$ echo $PATH     # show the value of $PATH
    /usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games
    user@machine:~$ which ls       # say what is the real path of ls

The first part of the command-line is the program to be run (`which`), the remaining parts are arguments (`ls` in this case). It is up to the program which arguments it understands. Many programs can take special arguments called options starting with `-` (followed by a single letter) or `--` (followed by words); one of the common options is `-h` or `--help`, which displays how to use the program (try `ls --help`).

Full documentation for each program usually exists as manual page (or man page), which can be shown using e.g. `man ls` (`q` to exit)

### Starting yade

If yade is installed on the machine, it can be (roughly speaking) run as any other program; without any arguments, it runs in the “dialog mode”, where a command-line is presented:

    user@machine:~$ yade
    Welcome to Yade 2019.01a
    TCP python prompt on localhost:9002, auth cookie `adcusk'
    XMLRPC info provider on http://localhost:21002
    [[ ^L clears screen, ^U kills line. F12 controller, F11 3d view, F10 both, F9 generator, F8 plot. ]]
    Yade [1]:                                            #### hit ^D to exit
    Do you really want to exit ([y]/n)?
    Yade: normal exit.

The command-line is in fact `python`, enriched with some yade-specific features. (Pure python interpreter can be run with `python` or `ipython` commands).

Instead of typing commands on-by-one on the command line, they can be be written in a file (with the .py extension) and given as argument to Yade:

    user@machine:~$ yade simulation.py

For a complete help, see `man yade`

Exercises

1. Open the terminal, navigate to your home directory

2. Create a new empty file and save it in `~/first.py`

3. Change directory to `/tmp`; delete the file `~/first.py`

4. Run program `xeyes`

5. Look at the help of Yade.

6. Look at the manual page of Yade

7. Run Yade, exit and run it again.

## Python basics

We assume the reader is familar with Python tutorial and only briefly review some of the basic capabilities. The following will run in pure-python interpreter (`python` or `ipython`), but also inside Yade, which is a super-set of Python.

Numerical operations and modules:

    Yade [1]: (1+3*4)**2        # usual rules for operator precedence, ** is exponentiation
    Out[1]: 169

    Yade [2]: import math       # gain access to "module" of functions

    Yade [3]: math.sqrt(2)      # use a function from that module
    Out[3]: 1.4142135623730951

    Yade [4]: import math as m  # use the module under a different name

    Yade [5]: m.cos(m.pi)
    Out[5]: -1.0

    Yade [6]: from math import *  # import everything so that it can be used without module name

    Yade [7]: cos(pi)
    Out[7]: -1.0

Variables:

    Yade [8]: a=1; b,c=2,3       # multiple commands separated with ;, multiple assignment

    Yade [9]: a+b+c
    Out[9]: 6

### Sequences

#### Lists

Lists are variable-length sequences, which can be modified; they are written with braces `[...]`, and their elements are accessed with numerical indices:

    Yade [10]: a=[1,2,3]          # list of numbers

    Yade [11]: a[0]               # first element has index 0
    Out[11]: 1

    Yade [12]: a[-1]              # negative counts from the end
    Out[12]: 3

    Yade [13]: a[3]               # error
    [31m---------------------------------------------------------------------------[39m
    [31mIndexError[39m                                Traceback (most recent call last)
    [36mCell[39m[36m [39m[32mIn[13][39m[32m, line 1[39m
    [32m----> [39m[32m1[39m [43ma[49m[43m[[49m[32;43m3[39;49m[43m][49m               [38;5;66;03m# error[39;00m

    [31mIndexError[39m: list index out of range

    Yade [14]: len(a)             # number of elements
    Out[14]: 3

    Yade [15]: a[1:]              # from second element to the end
    Out[15]: [2, 3]

    Yade [16]: a+=[4,5]           # extend the list

    Yade [17]: a+=[6]; a.append(7) # extend with single value, both have the same effect

    Yade [18]: 9 in a             # test presence of an element
    Out[18]: False

Lists can be created in various ways:

    Yade [19]: range(10)
    Out[19]: range(0, 10)

    Yade [20]: range(10)[-1]
    Out[20]: 9

List of squares of even number smaller than 20, i.e. \(\left\{a^2\;\forall a\in \{0,\cdots,19\} \;\middle|\; 2 \| a\right\}\) (note the similarity):

    Yade [21]: [a**2 for a in range(20) if a%2==0]
    Out[21]: [0, 4, 16, 36, 64, 100, 144, 196, 256, 324]

#### Tuples

Tuples are constant sequences:

    Yade [22]: b=(1,2,3)

    Yade [23]: b[0]
    Out[23]: 1

    Yade [24]: b[0]=4              # error
    [31m---------------------------------------------------------------------------[39m
    [31mTypeError[39m                                 Traceback (most recent call last)
    [36mCell[39m[36m [39m[32mIn[24][39m[32m, line 1[39m
    [32m----> [39m[32m1[39m [43mb[49m[43m[[49m[32;43m0[39;49m[43m][49m=[32m4[39m              [38;5;66;03m# error[39;00m

    [31mTypeError[39m: 'tuple' object does not support item assignment

#### Dictionaries

Mapping from keys to values:

    Yade [25]: ende={'one':'ein' , 'two':'zwei' , 'three':'drei'}

    Yade [26]: de={1:'ein' , 2:'zwei' , 3:'drei'}; en={1:'one' , 2:'two' , 3:'three'}

    Yade [27]: ende['one']         ## access values
    Out[27]: 'ein'

    Yade [28]: de[1], en[2]
    Out[28]: ('ein', 'two')

### Functions, conditionals

    Yade [29]: 4==5
    Out[29]: False

    Yade [30]: a=3.1

    Yade [31]: if a<10:
       ....:     b=-2          # conditional statement
       ....: else:
       ....:     b=3
       ....:

    Yade [32]: c=0 if a<1 else 1      # trenary conditional expression

    Yade [33]: b,c
    Out[33]: (-2, 1)

    Yade [34]: def square(x): return x**2    # define a new function
       ....:

    Yade [35]: square(2)                     # and call that function
    Out[35]: 4

Exercises

1. Read the following code and say what wil be the values of `a` and `b`:

       a=range(5)
       b=[(aa**2 if aa%2==0 else -aa**2) for aa in a]

## Yade basics

Yade objects are constructed in the following manner (this process is also called “instantiation”, since we create concrete instances of abstract classes: one individual sphere is an instance of the abstract Sphere, like Socrates is an instance of “man”):

    Yade [36]: Sphere           # try also Sphere?
    Out[36]: yade.wrapper.Sphere

    Yade [37]: s=Sphere()       # create a Sphere, without specifying any attributes

    Yade [38]: s.radius         # 'nan' is a special value meaning "not a number" (i.e. not defined)
    Out[38]: nan

    Yade [39]: s.radius=2       # set radius of an existing object

    Yade [40]: s.radius
    Out[40]: 2.0

    Yade [41]: ss=Sphere(radius=3)   # create Sphere, giving radius directly

    Yade [42]: s.radius, ss.radius     # also try typing s.<tab> to see defined attributes
    Out[42]: (2.0, 3.0)

### Particles

Particles are the “data” component of simulation; they are the objects that will undergo some processes, though do not define those processes yet.

#### Singles

There is a number of pre-defined functions to create particles of certain type; in order to create a sphere, one has to (see the source of sphere for instance):

1. Create Body

2. Set Body.shape to be an instance of Sphere with some given radius

3. Set Body.material (last-defined material is used, otherwise a default material is created)

4. Set position and orientation in Body.state, compute mass and moment of inertia based on Material and Shape

In order to avoid such tasks, shorthand functions are defined in the utils module; to mention a few of them, they are utils.sphere, utils.facet, utils.wall. The utils module is imported at startup in such a way that the `utils.` prefix is not necessary for accessing them.

    Yade [43]: s=sphere((0,0,0),radius=1)    # create sphere particle centered at (0,0,0) with radius=1

    Yade [44]: s.shape                       # s.shape describes the geometry of the particle
    Out[44]: <Sphere instance at 0x24285170>

    Yade [45]: s.shape.radius                # we already know the Sphere class
    Out[45]: 1.0

    Yade [46]: s.state.mass, s.state.inertia # inertia is computed from density and geometry
    Out[46]:
    (4188.790204786391,
     Vector3(1675.516081914556253,1675.516081914556253,1675.516081914556253))

    Yade [47]: s.state.pos                   # position is the one we prescribed
    Out[47]: Vector3(0,0,0)

    Yade [48]: s2=sphere((-2,0,0),radius=1,fixed=True)     # explanation below

In the last example, the particle was fixed in space by the `fixed=True` parameter to sphere; such a particle will not move, creating a primitive boundary condition.

A particle object is not yet part of the simulation; in order to do so, a special function O.bodies.append (also see Omega::bodies and Scene) is called:

    Yade [49]: O.bodies.append(s)            # adds particle s to the simulation; returns id of the particle(s) added
    Out[49]: 26

#### Packs

There are functions to generate a specific arrangement of particles in the pack module; for instance, cloud (random loose packing) of spheres can be generated with the pack.SpherePack class:

    Yade [50]: from yade import pack

    Yade [51]: sp=pack.SpherePack()                   # create an empty cloud; SpherePack contains only geometrical information

    Yade [52]: sp.makeCloud((1,1,1),(2,2,2),rMean=.2) # put spheres with defined radius inside box given by corners (1,1,1) and (2,2,2)
    Out[52]: 7

    Yade [53]: for c,r in sp: print(c,r)               # print center and radius of all particles (SpherePack is a sequence which can be iterated over)
       ....:
    Vector3(1.713652662502263668,1.66425785124973169,1.774021581779061885) 0.2
    Vector3(1.674034048872393488,1.627750978698175421,1.341140339843993612) 0.2
    Vector3(1.280648149400611668,1.449855563106149603,1.718767618078457016) 0.2
    Vector3(1.27740625110324002,1.747045396026098318,1.223225846831202723) 0.2
    Vector3(1.642699430388902204,1.225579590216237413,1.221307881119183758) 0.2
    Vector3(1.630461830625127284,1.253598927475262137,1.740564718622163998) 0.2
    Vector3(1.212332950739578585,1.302293173532684811,1.325494639911212191) 0.2

    Yade [54]: sp.toSimulation()                      # create particles and add them to the simulation
    Out[54]: [27, 28, 29, 30, 31, 32, 33]

#### Boundaries

facet (triangle Facet), wall (infinite axes-aligned plane Wall) and box (finite axes-aligned cuboids Box) geometries are typically used to define boundaries. For instance, a “floor” for the simulation can be created like this:

    Yade [55]: O.bodies.append(wall(-1,axis=2))
    Out[55]: 34

There are other convenience functions (like aabbWall for creating closed or open rectangular box, or family of ymport functions)

### Look inside

The simulation can be inspected in several ways. All data can be accessed from python directly:

    Yade [56]: len(O.bodies)
    Out[56]: 35

    Yade [57]: O.bodies[10].shape.radius   # radius of body #10 (will give error if not sphere, since only spheres have radius defined)
    Out[57]: 0.16

    Yade [58]: O.bodies[12].state.pos      # position of body #12
    Out[58]: Vector3(1.521292810892004077,1.058202337666694248,1.831570188729839455)

Besides that, Yade says this at startup (the line preceding the command-line):

    [[ ^L clears screen, ^U kills line. F12 controller, F11 3d view, F10 both, F9 generator, F8 plot. ]]

Controller

  Pressing `F12` brings up a window for controlling the simulation. Although typically no human intervention is done in large simulations (which run “headless”, without any graphical interaction), it can be handy in small examples. There are basic information on the simulation (will be used later).

3d view

  The 3d view can be opened with F11 (or by clicking on button in the Controller – see below). There is a number of keyboard shortcuts to manipulate it (press `h` to get basic help), and it can be moved, rotated and zoomed using mouse. Display-related settings can be set in the “Display” tab of the controller (such as whether particles are drawn).

Inspector

  Inspector is opened by clicking on the appropriate button in the Controller. It shows (and updates) internal data of the current simulation. In particular, one can have a look at engines, particles (Bodies) and interactions (Interactions). Clicking at each of the attribute names links to the appropriate section in the documentation.

Exercises

1. What is this code going to do?

       Yade [59]: O.bodies.append([sphere((2*i,0,0),1) for i in range(1,20)])
       Out[59]: [35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53]

2. Create a simple simulation with cloud of spheres enclosed in the box `(0,0,0)` and `(1,1,1)` with mean radius .1. (hint: pack.SpherePack.makeCloud)

3. Enclose the cloud created above in box with corners `(0,0,0)` and `(1,1,1)`; keep the top of the box open. (hint: aabbWall; type `aabbWall?` or `aabbWall??` to get help on the command line)

4. Open the 3D view, try zooming in/out; position axes so that \(z\) is upwards, \(y\) goes to the right and \(x\) towards you.

### Engines

Engines define processes undertaken by particles. As we know from the theoretical introduction, the sequence of engines is called simulation loop. Let us define a simple interaction loop:

    Yade [60]: O.engines=[                   # newlines and indentations are not important until the brace is closed
       ....:    ForceResetter(),
       ....:    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
       ....:    InteractionLoop(           # dtto for the parenthesis here
       ....:        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
       ....:        [Ip2_FrictMat_FrictMat_FrictPhys()],
       ....:        [Law2_ScGeom_FrictPhys_CundallStrack()]
       ....:    ),
       ....:    NewtonIntegrator(damping=.2,label='newtonCustomLabel')      # define a label newtonCustomLabel under which we can access this engine easily
       ....: ]
       ....:

    Yade [61]: O.engines
    Out[61]:
    [<ForceResetter instance at 0x24047a60>,
     <InsertionSortCollider instance at 0x240f95b0>,
     <InteractionLoop instance at 0x242ebf30>,
     <NewtonIntegrator instance at 0x23f4dfd0>]

    Yade [62]: O.engines[-1]==newtonCustomLabel    # is it the same object?
    Out[62]: True

    Yade [63]: newtonCustomLabel.damping
    Out[63]: 0.2

Instead of typing everything into the command-line, one can describe simulation in a file (script) and then run yade with that file as an argument. We will therefore no longer show the command-line unless necessary; instead, only the script part will be shown. Like this:

    O.engines=[                   # newlines and indentations are not important until the brace is closed
            ForceResetter(),
            InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
            InteractionLoop(           # dtto for the parenthesis here
                     [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                     [Ip2_FrictMat_FrictMat_FrictPhys()],
                     [Law2_ScGeom_FrictPhys_CundallStrack()]
            ),
            GravityEngine(gravity=(0,0,-9.81)),                    # 9.81 is the gravity acceleration, and we say that
            NewtonIntegrator(damping=.2,label='newtonCustomLabel') # define a label under which we can access this engine easily
    ]

Besides engines being run, it is likewise important to define how often they will run. Some engines can run only sometimes (we will see this later), while most of them will run always; the time between two successive runs of engines is timestep (\(\Dt\)). There is a mathematical limit on the timestep value, called critical timestep, which is computed from properties of particles. Since there is a function for that, we can just set timestep using PWaveTimeStep:

    O.dt=PWaveTimeStep()

Each time when the simulation loop finishes, time `O.time` is advanced by the timestep `O.dt`:

    Yade [64]: O.dt=0.01

    Yade [65]: O.time
    Out[65]: 0.0

    Yade [66]: O.step()

    Yade [67]: O.time
    Out[67]: 0.01

For experimenting with a single simulations, it is handy to save it to memory; this can be achieved, once everything is defined, with:

    O.saveTmp()

Exercises

1. Define engines as in the above example, run the Inspector and click through the engines to see their sequence.

2. Write a simple script which will

   1. define particles as in the previous exercise (cloud of spheres inside a box open from the top) but with a smaller radius (*rMean*=.04)

   2. define a simple simulation loop, as the one given above

   3. set \(\Dt\) equal to the critical P-Wave \(\Dt\)

   4. save the initial simulation state to memory

3. Run the previously-defined simulation multiple times, while changing the value of timestep (use the ⟳ button to reload the initial configuration).

   1. Try changing the gravity parameter, before running the simulation.

   2. See what happens as you increase \(\Dt\) above the P-Wave value.

   3. Try changing damping

4. Reload the simulation, open the 3d view, open the Inspector, select a particle in the 3d view (shift-click). Then run the simulation and watch how forces on that particle change; pause the simulation somewhere in the middle, look at interactions of this particle.

5. At which point can we say that the deposition is done, so that the simulation can be stopped?

See also

The Bouncing sphere example shows a basic simulation.