\( \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<} \)

# Data mining

## Read

### Local data

All data of the simulation are accessible from python; when you open the Inspector, blue labels of various data can be clicked – left button for getting to the documentation, middle click to copy the name of the object (use `Ctrl-V` or middle-click to paste elsewhere). The interesting objects are among others (see Omega for a full list):

1. O.engines

   Engines are accessed by their index (position) in the simulation loop:

       O.engines[0]      # first engine
       O.engines[-1]     # last engine

   Note

   The index can change if O.engines is modified. Labeling introduced in the section below is a better solution for reliable access to a particular engine.

2. O.bodies

   Bodies are identified by their id, which is guaranteed to not change during the whole simulation:

       O.bodies[0]                                                         # first body
       [b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]    # list of radii of all spherical bodies
       sum([b.state.mass for b in O.bodies])                               # sum of masses of all bodies
       numpy.average([b.state.vel[0] for b in O.bodies])                   # average velocity in x direction

   Note

   Uniqueness of Body.id is not guaranteed, since newly created bodies might recycle ids of deleted ones.

3. O.forces

   Generalized forces (forces, torques) acting on each particle. They are (usually) reset at the beginning of each step with ForceResetter, subsequently forces from individual interactions are accumulated in InteractionLoop. To access the data, use:

       O.forces.f(0)     # force on #0
       O.forces.t(1)     # torque on #1

4. O.interactions

   Interactions are identified by ids of the respective interacting particles (they are created and deleted automatically during the simulation):

       O.interactions[0,1]   # interactions of #0 with #1
       O.interactions[1,0]   # the same object
       O.bodies[0].intrs()   # all interactions of body #0
       for i in O.bodies[12].intrs(): print (i.isReal,i.id1,i.id2)    # get some info about interactions of body #12
       [(i.isReal,i.id1,i.id2) for i in O.bodies[12].intrs()]         # same thing, but make a list

#### Labels

Engines and functors can be labeled, which means that python variable of that name is automatically created.

    Yade [1]: O.engines=[
       ...:    NewtonIntegrator(damping=.2,label='newtonCustomLabel')
       ...: ]
       ...:

    Yade [2]: newtonCustomLabel.damping=.4

    Yade [3]: O.engines[0].damping              # O.engines[0] and newtonCustomLabel are the same objects
    Out[3]: 0.4

    Yade [4]: newtonCustomLabel==O.engines[0]   # O.engines[0] and newtonCustomLabel are the same objects
    Out[4]: True

Exercises

1. Find meaning of this expression:

       max([b.state.vel.norm() for b in O.bodies])

2. Run the Gravity deposition script, pause after a few seconds of simulation. Write expressions that compute

   1. kinetic energy \(\sum \frac{1}{2} m_i |v_i| ^2\)

   2. average mass (hint: use numpy.average)

   3. maximum \(z\)-coordinate of all particles

   4. number of interactions of body #1

### Global data

Useful measures of what happens in the simulation globally:

unbalanced force

  ratio of maximum contact force and maximum per-body force; measure of staticity, computed with unbalancedForce.

porosity

  ratio of void volume and total volume; computed with porosity.

coordination number

  average number of interactions per particle, avgNumInteractions

stress tensor (periodic boundary conditions)

  averaged force in interactions, computed with normalShearStressTensors

fabric tensor

  distribution of contacts in space (not yet implemented); can be visualized with plotDirections

#### Energies

Evaluating energy data for all components in the simulation (such as gravity work, kinetic energy, plastic dissipation, damping dissipation) can be enabled with

    O.trackEnergy=True

Subsequently, energy values are accessible in the O.energy; it is a dictionary where its entries can be retrived with `keys()` and their values with `O.energy[key]`.

## Save

### PyRunner

To save data that we just learned to access, we need to call Python from within the simulation loop. PyRunner is created just for that; it inherits periodicy control from PeriodicEngine and takes the code to run as text (must be quoted, i.e. inside `'...'`) attribute called command. For instance, adding this to O.engines will print the current step number every one second wall clock time:

    O.engines=O.engines+[ PyRunner(command='print(O.iter)',realPeriod=1) ]

Writing complicated code inside command is awkward; in such case, we define a function that will be called:

    def myFunction():
            '''Print step number, and pause the simulation is unbalanced force is smaller than 0.05.'''
            print(O.iter)
            if unbalancedForce()<0.05:
                    print('Unbalanced force is smaller than 0.05, pausing.')
                    O.pause()

Now this function can be added to O.engines:

    O.engines+=[PyRunner(command='myFunction()',iterPeriod=100)]

or, in general, like that:

    O.engines=[
            # ...
            PyRunner(command='myFunction()',iterPeriod=100) # call myFunction every 100 steps
    ]

Warning

If a function was declared inside a live yade session (ipython) and PyRunner attribute updateGlobals is set to False then an error `NameError: name 'myFunction' is not defined` will occur unless python globals() are updated with command

    globals().update(locals())

Exercises

1. Run the Gravity deposition simulation, but change it such that:

   1. utils.unbalancedForce is printed every 2 seconds.

   2. check every 1000 steps the value of unbalanced force

      - if smaller than 0.2, set damping to 0.8 (hint: use labels)

      - if smaller than 0.1, pause the simulation

### Keeping history

Yade provides the plot module used for storing and plotting variables (plotting itself will be discussed later). Let us start by importing this module and declare variable names that will be plotted:

    from yade import plot
    plot.plots={'t':('coordNum','unForce',None,'Ek')}                # kinetic energy will have legend on the right as indicated by None separator.

Periodic storing of data is done with PyRunner and the plot.addData function. Also let’s enable energy tracking:

    O.trackEnergy=True
    def addPlotData():
            # this function adds current values to the history of data, under the names specified
            plot.addData(t=O.time,Ek=kineticEnergy(),coordNum=avgNumInteractions(),unForce=unbalancedForce())

Now this function can be added to O.engines:

    O.engines+=[PyRunner(command='addPlotData()',iterPeriod=20)]

or, in general, like that:

    O.engines=[  # ...,
            PyRunner(command='addPlotData()',iterPeriod=20)         # call the addPlotData function every 20 iterations
    ]

History is stored in plot.data, and can be accessed using the variable name, e.g. `plot.data['Ek']`, and saved to text file (for post-processing outside yade) with plot.saveDataTxt.

## Plot

plot provides facilities for plotting history saved with plot.addData as 2d plots. Data to be plotted are specified using dictionary plot.plots

    plot.plots={'t':('coordNum','unForce',None,'Ek')}

History of all values is given as the name used for plot.addData; keys of the dictionary are \(x\)-axis values, and values are sequence of data on the \(y\) axis; the `None` separates data on the left and right axes (they are scaled independently). The plot itself is created with

    plot.plot()         # on the command line, F8 can be used as shorthand

While the plot is open, it will be updated periodically, so that simulation evolution can be seen in real-time.

### Energy plots

Plotting all energy contributions would be difficult, since names of all energies might not be known in advance. Fortunately, there is a way to handle that in Yade. It consists in two parts:

1. plot.addData is given all the energies that are currently defined:

       plot.addData(i=O.iter,total=O.energy.total(),**O.energy)

   The O.energy.total functions, which sums all energies together. The `**O.energy` is special python syntax for converting dictionary (remember that O.energy is a dictionary) to named functions arguments, so that the following two commands are identical:

       function(a=3,b=34)              # give arguments as arguments
       function(**{'a':3,'b':34})      # create arguments from dictionary

2. Data to plot are specified using a function that gives names of data to plot, rather than providing the data names directly:

       plot.plots={'i':['total']+O.energy.keys()}

   where `total` is the name we gave to `O.energy.total()` above, while `O.energy.keys()` will always return list of currently defined energies.

#### Energy plot example

Plotting energies inside a live yade session, for example by launching examples/test/triax-basic-without-plots.py would look following:

    from yade import plot
    O.trackEnergy=True
    O.step()                          # performing a single simulation step is necessary to populate O.energy.keys()
    plot.plots={'t':O.energy.keys()+['total']}

    def addPlotData():
            # this function adds current values to the history of data, under the names specified
            plot.addData( t=O.time , total=O.energy.total() , **O.energy )

    O.engines+=[PyRunner(command='addPlotData()',iterPeriod=20)]

    globals().update(locals())        # do this only because this is an example of a live yade session

Press F8 to show plot window and F11 to show 3D view, then press ▶ to start simulation.

#### Using multiple plots

It is also possible to make several separate plots, for example like this:

    plot.plots={ 't':('total','kinetic') , 't ':['elastPotential','gravWork'] , 't  ':('nonviscDamp') }

Warning

There cannot be duplicate names declared in separate plots. This is why spaces were used above to indicate the same variable `t`.

With the caveat above, a following example inside a live yade session launched on examples/test/triax-basic-without-plots.py would look following:

    from yade import plot
    O.trackEnergy=True
    plot.plots={ 't':('total','kinetic') , 't ':['elastPotential','gravWork'] , 't  ':('nonviscDamp') }

    def addPlotData():
            # assign value to all three: 't', 't ' and 't  ' with single t=... assignment
            plot.addData( t=O.time , total=O.energy.total() , **O.energy )

    O.engines+=[PyRunner(command='addPlotData()',iterPeriod=20)]

    globals().update(locals())        # do this only because this is an example of a live yade session

    plot.plot(subPlots=False)         # show plots in separate windows

    plot.plot(subPlots=True)          # same as pressing F8: close current plot windows and reopen a single new one

Press F8 to show plot window and F11 to show 3D view, then press ▶ to start simulation, see video below:

Exercises

1. Calculate average momentum in y direction.

2. Run the Gravity deposition script, plotting unbalanced force and kinetic energy.

3. While the script is running, try changing the NewtonIntegrator.damping parameter (do it from both Inspector and from the command-line). What influence does it have on the evolution of unbalanced force and kinetic energy?

4. Think about and write down all energy sources (input); write down also all energy sinks (dissipation).

5. Simulate Gravity deposition and plot all energies as they evolve during the simulation.

See also

Most Examples with tutorial use plotting facilities of Yade, some of them also track energy of the simulation.