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

# Examples with tutorial

The online version of this tutorial contains embedded videos.

## Bouncing sphere

Following example is in file doc/sphinx/tutorial/01-bouncing-sphere.py.

    # basic simulation showing sphere falling ball gravity,
    # bouncing against another sphere representing the support

    # DATA COMPONENTS

    # add 2 particles to the simulation
    # they the default material (utils.defaultMat)
    O.bodies.append(
            [
                    # fixed: particle's position in space will not change (support)
                    sphere(center=(0, 0, 0), radius=.5, fixed=True),
                    # this particles is free, subject to dynamics
                    sphere((0, 0, 2), .5)
            ]
    )

    # FUNCTIONAL COMPONENTS

    # simulation loop -- see presentation for the explanation
    O.engines = [
            ForceResetter(),
            InsertionSortCollider([Bo1_Sphere_Aabb()]),
            InteractionLoop(
                    [Ig2_Sphere_Sphere_ScGeom()],  # collision geometry
                    [Ip2_FrictMat_FrictMat_FrictPhys()],  # collision "physics"
                    [Law2_ScGeom_FrictPhys_CundallStrack()]  # contact law -- apply forces
            ),
            # Apply gravity force to particles. damping: numerical dissipation of energy.
            NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.1)
    ]

    # set timestep to a fraction of the critical timestep
    # the fraction is very small, so that the simulation is not too fast
    # and the motion can be observed
    O.dt = .5e-4 * PWaveTimeStep()

    # save the simulation, so that it can be reloaded later, for experimentation
    O.saveTmp()

## Gravity deposition

Following example is in file doc/sphinx/tutorial/02-gravity-deposition.py.

    # gravity deposition in box, showing how to plot and save history of data,
    # and how to control the simulation while it is running by calling
    # python functions from within the simulation loop

    # import yade modules that we will use below
    from yade import pack, plot

    # create rectangular box from facets
    O.bodies.append(geom.facetBox((.5, .5, .5), (.5, .5, .5), wallMask=31))

    # create empty sphere packing
    # sphere packing is not equivalent to particles in simulation, it contains only the pure geometry
    sp = pack.SpherePack()
    # generate randomly spheres with uniform radius distribution
    sp.makeCloud((0, 0, 0), (1, 1, 1), rMean=.05, rRelFuzz=.5)
    # add the sphere pack to the simulation
    sp.toSimulation()

    O.engines = [
            ForceResetter(),
            InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
            InteractionLoop(
                    # handle sphere+sphere and facet+sphere collisions
                    [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                    [Ip2_FrictMat_FrictMat_FrictPhys()],
                    [Law2_ScGeom_FrictPhys_CundallStrack()]
            ),
            NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
            # call the checkUnbalanced function (defined below) every 2 seconds
            PyRunner(command='checkUnbalanced()', realPeriod=2),
            # call the addPlotData function every 200 steps
            PyRunner(command='addPlotData()', iterPeriod=100)
    ]
    O.dt = .5 * 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

    # if the unbalanced forces goes below .05, the packing
    # is considered stabilized, therefore we stop collected
    # data history and stop
    def checkUnbalanced():
       if unbalancedForce() < .05:
          O.pause()
          plot.saveDataTxt('bbb.txt.bz2')
          # plot.saveGnuplot('bbb') is also possible

    # collect history of data which will be plotted
    def addPlotData():
       # each item is given a names, by which it can be the unsed in plot.plots
       # the **O.energy converts dictionary-like O.energy to plot.addData arguments
       plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

    # define how to plot data: 'i' (step number) on the x-axis, unbalanced force
    # on the left y-axis, all energies on the right y-axis
    # (O.energy.keys is function which will be called to get all defined energies)
    # None separates left and right y-axis
    plot.plots = {'i': ('unbalanced', None, O.energy.keys)}

    # show the plot on the screen, and update while the simulation runs
    plot.plot()

    O.saveTmp()

## Oedometric test

Following example is in file doc/sphinx/tutorial/03-oedometric-test.py.

    # gravity deposition, continuing with oedometric test after stabilization
    # shows also how to run parametric studies with yade-batch

    # The components of the batch are:
    # 1. table with parameters, one set of parameters per line (ccc.table)
    # 2. readParamsFromTable which reads respective line from the parameter file
    # 3. the simulation muse be run using yade-batch, not yade
    #
    # $ yade-batch --job-threads=1 03-oedometric-test.table 03-oedometric-test.py
    #

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

    # create box with free top, and ceate loose packing inside the box
    from yade import pack, plot

    O.bodies.append(geom.facetBox((.5, .5, .5), (.5, .5, .5), wallMask=31))
    sp = pack.SpherePack()
    sp.makeCloud((0, 0, 0), (1, 1, 1), rMean=rMean, rRelFuzz=rRelFuzz)
    sp.toSimulation()

    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, -9.81), damping=0.5),
            # 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 = .5 * PWaveTimeStep()

    # 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 < 5000:
          return
       # the rest will be run only if unbalanced is < .1 (stabilized packing)
       if unbalancedForce() > .1:
          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, -.1)
       # 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]) > maxLoad:
          plate.state.vel *= -1
          # 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]) < minLoad:
          # 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()

    O.run()
    # when running with yade-batch, the script must not finish until the simulation is done fully
    # this command will wait for that (has no influence in the non-batch mode)
    waitIfBatch()

### Batch table

To run the same script doc/sphinx/tutorial/03-oedometric-test.py in batch mode to test different parameters, execute command `yade-batch 03-oedometric-test.table 03-oedometric-test.py`, also visit page http://localhost:9080 to see the batch simulation progress.

    rMean rRelFuzz maxLoad
    .05 .1 1e6
    .05 .2 1e6
    .05 .3 1e6

## Periodic simple shear

Following example is in file doc/sphinx/tutorial/04-periodic-simple-shear.py.

    # encoding: utf-8

    # script for periodic simple shear test, with periodic boundary
    # first compresses to attain some isotropic stress (checkStress),
    # then loads in shear (checkDistorsion)
    #
    # the initial packing is either regular (hexagonal), with empty bands along the boundary,
    # or periodic random cloud of spheres
    #
    # material friction angle is initially set to zero, so that the resulting packing is dense
    # (sphere rearrangement is easier if there is no friction)
    #

    # setup the periodic boundary
    O.periodic = True
    O.cell.hSize = Matrix3(2, 0, 0, 0, 2, 0, 0, 0, 2)

    from yade import pack, plot

    # the "if 0:" block will be never executed, therefore the "else:" block will be
    # to use cloud instead of regular packing, change to "if 1:" or something similar
    if 0:
       # create cloud of spheres and insert them into the simulation
       # we give corners, mean radius, radius variation
       sp = pack.SpherePack()
       sp.makeCloud((0, 0, 0), (2, 2, 2), rMean=.1, rRelFuzz=.6, periodic=True)
       # insert the packing into the simulation
       sp.toSimulation(color=(0, 0, 1))  # pure blue
    else:
       # in this case, add dense packing
       O.bodies.append(pack.regularHexa(pack.inAlignedBox((0, 0, 0), (2, 2, 2)), radius=.1, gap=0, color=(0, 0, 1)))

    # create "dense" packing by setting friction to zero initially
    O.materials[0].frictionAngle = 0

    # simulation loop (will be run at every step)
    O.engines = [
            ForceResetter(),
            InsertionSortCollider([Bo1_Sphere_Aabb()]),
            InteractionLoop(
                    # interaction loop
                    [Ig2_Sphere_Sphere_ScGeom()],
                    [Ip2_FrictMat_FrictMat_FrictPhys()],
                    [Law2_ScGeom_FrictPhys_CundallStrack()]
            ),
            NewtonIntegrator(damping=.4),
            # run checkStress function (defined below) every second
            # the label is arbitrary, and is used later to refer to this engine
            PyRunner(command='checkStress()', realPeriod=1, label='checker'),
            # record data for plotting every 100 steps; addData function is defined below
            PyRunner(command='addData()', iterPeriod=100)
    ]

    # set the integration timestep to be 1/2 of the "critical" timestep
    O.dt = .5 * PWaveTimeStep()

    # prescribe isotropic normal deformation (constant strain rate)
    # of the periodic cell
    O.cell.velGrad = Matrix3(-.1, 0, 0, 0, -.1, 0, 0, 0, -.1)

    # when to stop the isotropic compression (used inside checkStress)
    limitMeanStress = -5e5

    # called every second by the PyRunner engine
    def checkStress():
       # stress tensor as the sum of normal and shear contributions
       # Matrix3.Zero is the intial value for sum(...)
       stress = getStress().trace() / 3.
       print('mean stress', stress)
       # if mean stress is below (bigger in absolute value) limitMeanStress, start shearing
       if stress < limitMeanStress:
          # apply constant-rate distorsion on the periodic cell
          O.cell.velGrad = Matrix3(0, 0, .1, 0, 0, 0, 0, 0, 0)
          # change the function called by the checker engine
          # (checkStress will not be called anymore)
          checker.command = 'checkDistorsion()'
          # block rotations of particles to increase tanPhi, if desired
          # disabled by default
          if 0:
             for b in O.bodies:
                # block X,Y,Z rotations, translations are free
                b.state.blockedDOFs = 'XYZ'
                # stop rotations if any, as blockedDOFs block accelerations really
                b.state.angVel = (0, 0, 0)
          # set friction angle back to non-zero value
          # tangensOfFrictionAngle is computed by the Ip2_* functor from material
          # for future contacts change material (there is only one material for all particles)
          O.materials[0].frictionAngle = .5  # radians
          # for existing contacts, set contact friction directly
          for i in O.interactions:
             i.phys.tangensOfFrictionAngle = tan(.5)

    # called from the 'checker' engine periodically, during the shear phase
    def checkDistorsion():
       # if the distorsion value is >.3, exit; otherwise do nothing
       if abs(O.cell.trsf[0, 2]) > .5:
          # save data from addData(...) before exiting into file
          # use O.tags['id'] to distinguish individual runs of the same simulation
          plot.saveDataTxt(O.tags['id'] + '.txt')
          # exit the program
          #import sys
          #sys.exit(0) # no error (0)
          O.pause()

    # called periodically to store data history
    def addData():
       # get the stress tensor (as 3x3 matrix)
       stress = sum(normalShearStressTensors(), Matrix3.Zero)
       # give names to values we are interested in and save them
       plot.addData(exz=O.cell.trsf[0, 2], szz=stress[2, 2], sxz=stress[0, 2], tanPhi=(stress[0, 2] / stress[2, 2]) if stress[2, 2] != 0 else 0, i=O.iter)
       # color particles based on rotation amount
       for b in O.bodies:
          # rot() gives rotation vector between reference and current position
          b.shape.color = scalarOnColorScale(b.state.rot().norm(), 0, pi / 2.)

    # define what to plot (3 plots in total)
    ## exz(i), [left y axis, separate by None:] szz(i), sxz(i)
    ## szz(exz), sxz(exz)
    ## tanPhi(i)
    # note the space in 'i ' so that it does not overwrite the 'i' entry
    plot.plots = {'i': ('exz', None, 'szz', 'sxz'), 'exz': ('szz', 'sxz'), 'i ': ('tanPhi',)}

    # better show rotation of particles
    Gl1_Sphere.stripes = True

    # open the plot on the screen
    plot.plot()

    O.saveTmp()

## 3d postprocessing

Following example is in file doc/sphinx/tutorial/05-3d-postprocessing.py. This example will run for 20000 iterations, saving `*.png` snapshots, then it will make a video `3d.mpeg` out of those snapshots.

    # demonstrate 3d postprocessing with yade
    #
    # 1. qt.SnapshotEngine saves images of the 3d view as it appears on the screen periodically
    #    makeVideo is then used to make real movie from those images
    # 2. VTKRecorder saves data in files which can be opened with Paraview
    #    see the User's manual for an intro to Paraview

    # generate loose packing
    from yade import pack, qt

    sp = pack.SpherePack()
    sp.makeCloud((0, 0, 0), (2, 2, 2), rMean=.1, rRelFuzz=.6, periodic=True)
    # add to scene, make it periodic
    sp.toSimulation()

    O.engines = [
            ForceResetter(),
            InsertionSortCollider([Bo1_Sphere_Aabb()]),
            InteractionLoop(
                    # interaction loop
                    [Ig2_Sphere_Sphere_ScGeom()],
                    [Ip2_FrictMat_FrictMat_FrictPhys()],
                    [Law2_ScGeom_FrictPhys_CundallStrack()]
            ),
            NewtonIntegrator(damping=.4),
            # save data for Paraview
            VTKRecorder(fileName='3d-vtk-', recorders=['all'], iterPeriod=1000),
            # save data from Yade's own 3d view
            qt.SnapshotEngine(fileBase='3d-', iterPeriod=200, label='snapshot'),
            # this engine will be called after 20000 steps, only once
            PyRunner(command='finish()', iterPeriod=20000)
    ]
    O.dt = .5 * PWaveTimeStep()

    # prescribe constant-strain deformation of the cell
    O.cell.velGrad = Matrix3(-.1, 0, 0, 0, -.1, 0, 0, 0, -.1)

    # we must open the view explicitly (limitation of the qt.SnapshotEngine)
    qt.View()

    # this function is called when the simulation is finished
    def finish():
       # snapshot is label of qt.SnapshotEngine
       # the 'snapshots' attribute contains list of all saved files
       makeVideo(snapshot.snapshots, '3d.mpeg', fps=10, bps=10000)
       O.pause()

    # set parameters of the renderer, to show network chains rather than particles
    # these settings are accessible from the Controller window, on the second tab ("Display") as well
    rr = yade.qt.Renderer()
    rr.shape = False
    rr.intrPhys = True

## Periodic triaxial test

Following example is in file doc/sphinx/tutorial/06-periodic-triaxial-test.py. A variant of this exemple includes capillary forces, see doc/sphinx/tutorial/06-periodic-triaxial-test-capillarity.py

    # encoding: utf-8

    # periodic triaxial test simulation
    #
    # The initial packing is either
    #
    # 1. random cloud with uniform distribution, or
    # 2. cloud with specified granulometry (radii and percentages), or
    # 3. cloud of clumps, i.e. rigid aggregates of several particles
    #
    # The triaxial consists of 2 stages:
    #
    # 1. isotropic compaction, until sigmaIso is reached in all directions;
    #    this stage is ended by calling compactionFinished()
    # 2. constant-strain deformation along the z-axis, while maintaining
    #    constant stress (sigmaIso) laterally; this stage is ended by calling
    #    triaxFinished()
    #
    # Controlling of strain and stresses is performed via PeriTriaxController,
    # of which parameters determine type of control and also stability
    # condition (maxUnbalanced) so that the packing is considered stabilized
    # and the stage is done.
    #

    sigmaIso = -1e5

    #import matplotlib
    #matplotlib.use('Agg')

    # generate loose packing
    from yade import pack, qt, plot

    O.periodic = True
    sp = pack.SpherePack()
    if 0:
       ## uniform distribution
       sp.makeCloud((0, 0, 0), (2, 2, 2), rMean=.1, rRelFuzz=.3, periodic=True)
    else:
       ## create packing from clumps
       # configuration of one clump
       c1 = pack.SpherePack([((0, 0, 0), .03333), ((.03, 0, 0), .017), ((0, .03, 0), .017)])
       # make cloud using the configuration c1 (there could c2, c3, ...; selection between them would be random)
       sp.makeClumpCloud((0, 0, 0), (2, 2, 2), [c1], periodic=True, num=500)

    # setup periodic boundary, insert the packing
    sp.toSimulation()

    O.engines = [
            ForceResetter(),
            InsertionSortCollider([Bo1_Sphere_Aabb()]),
            InteractionLoop([Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
            PeriTriaxController(
                    label='triax',
                    # specify target values and whether they are strains or stresses
                    goal=(sigmaIso, sigmaIso, sigmaIso),
                    stressMask=7,
                    # type of servo-control
                    dynCell=True,
                    maxStrainRate=(10, 10, 10),
                    # wait until the unbalanced force goes below this value
                    maxUnbalanced=.1,
                    relStressTol=1e-3,
                    # call this function when goal is reached and the packing is stable
                    doneHook='compactionFinished()'
            ),
            NewtonIntegrator(damping=.2),
            PyRunner(command='addPlotData()', iterPeriod=100),
    ]
    O.dt = .5 * PWaveTimeStep()

    def addPlotData():
       plot.addData(
               unbalanced=unbalancedForce(),
               i=O.iter,
               sxx=triax.stress[0],
               syy=triax.stress[1],
               szz=triax.stress[2],
               exx=triax.strain[0],
               eyy=triax.strain[1],
               ezz=triax.strain[2],
               # save all available energy data
               Etot=O.energy.total(),
               **O.energy
       )

    # enable energy tracking in the code
    O.trackEnergy = True

    # define what to plot
    plot.plots = {
            'i': ('unbalanced',),
            'i ': ('sxx', 'syy', 'szz'),
            ' i': ('exx', 'eyy', 'ezz'),
            # energy plot
            ' i ': (O.energy.keys, None, 'Etot'),
    }
    # show the plot
    plot.plot()

    def compactionFinished():
       # set the current cell configuration to be the reference one
       O.cell.trsf = Matrix3.Identity
       # change control type: keep constant confinement in x,y, 20% compression in z
       triax.goal = (sigmaIso, sigmaIso, -.2)
       triax.stressMask = 3
       # allow faster deformation along x,y to better maintain stresses
       triax.maxStrainRate = (1., 1., .1)
       # next time, call triaxFinished instead of compactionFinished
       triax.doneHook = 'triaxFinished()'
       # do not wait for stabilization before calling triaxFinished
       triax.maxUnbalanced = 10

    def triaxFinished():
       print('Finished')
       O.pause()

## Fluid injection

Following example is in file doc/sphinx/tutorial/07-fluid-injection.py. The video below results from post-processing with paraview

    # This script simulates the injection of a fluid in a localized region below immersed particles
    # The simulation is periodic along z-axis.
    # at first execution, the simulation starts by depositing the particles in the container then saves the scene before proceeding to injection
    # further execution will reload the deposited layer and start injection directly to gain time
    # WARNING: changes in some input parameters like dimensions of the box or number of particles will not be reflected as long as the saved state is present on disk,
    #          remember to erase it to force a new generation, or set newSample=True below

    from yade import pack, export
    import yade.timing
    from math import *
    from pylab import rand
    import os.path
    import numpy as np
    import matplotlib.pyplot as plt

    O.periodic = True

    # Dimensions of the box and of the injection zone
    width = 0.8  #
    height = 1
    depth = 0.4
    aperture = 0.05 * width

    # number of spheres
    numSpheres = 2000
    # contact friction during deposition
    compFricDegree = 10
    # porosity of the initial cloud
    porosity = 0.8
    # Cundall's damping (zero recommanded)
    damp = 0.0
    # fluid viscosity
    mu = 0.01
    # flow rate at the inlet
    injectedFlux = -0.001
    # name of output folder
    key = 'output0'

    newSample = False  #turn this true if you want to generate new sample by pluviation each time you run the script

    # Deduced mean size for generating the cloud and consistency check
    filename = "init" + key + str(numSpheres)
    volume = width * height * depth
    meanRad = pow(volume * (1 - porosity) / (pi * (4 / 3.) * numSpheres), 1 / 3.)
    if (meanRad * 6 > depth):
       print("INCOMPATIBLE SIZES. INCREASE DEPTH OR INCREASE NUM_SPHERES")

    # if no deposited layer has been found (first execution), generate bodies
    if not os.path.isfile(filename + ".yade") or newSample:  #we create new sample if it does not exist...
       O.cell.hSize = Matrix3(width, 0, 0, 0, 3. * height, 0, 0, 0, depth)
       O.materials.append(FrictMat(young=400000.0, poisson=0.5, frictionAngle=compFricDegree / 180.0 * pi, density=1600, label='spheres'))
       O.materials.append(FrictMat(young=400000.0, poisson=0.5, frictionAngle=radians(15), density=1000, label='walls'))
       lowBox = box(center=(width / 2.0, height, width / 2.0), extents=(width * 1000.0, 0, width * 1000.0), fixed=True, wire=False)
       O.bodies.append(lowBox)
       topBox = box(center=(width / 2.0, 2 * height + 4 * meanRad, width / 2.0), extents=(width * 1000.0, 0, width * 1000.0), fixed=True, wire=False)
       O.bodies.append(topBox)

       sp = pack.SpherePack()
       sp.makeCloud((0, height + 2 * meanRad, 0), (width, 2 * height + 2 * meanRad, depth), -1, .002, numSpheres, periodic=True, porosity=porosity, seed=2)
       O.bodies.append([sphere(s[0], s[1], color=(0.6 + 0.15 * rand(), 0.5 + 0.15 * rand(), 0.15 + 0.15 * rand()), material='spheres') for s in sp])
    else:  #... else we re-use the previous one
       O.load(filename + ".yade")

    # look better
    Gl1_Sphere.stripes = True

    # define the fluid solver with appropriate parameter (see PeriodicFlowEngine's documentation)
    flow = PeriodicFlowEngine(
            dead=1,
            meshUpdateInterval=40,
            defTolerance=-1,
            permeabilityFactor=1.0,
            useSolver=3,
            duplicateThreshold=depth,
            wallIds=[-1, -1, 0, 1, -1, -1],
            bndCondIsPressure=[0, 0, 0, 1, 0, 0],
            viscosity=mu,
            label="flow"
    )

    newton = NewtonIntegrator(damping=damp, gravity=(0, -9.81, 0))

    O.engines = [
            ForceResetter(),
            InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Box_Aabb()], allowBiggerThanPeriod=1),
            InteractionLoop([Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
            GlobalStiffnessTimeStepper(active=1, timeStepUpdateInterval=1000, timestepSafetyCoefficient=0.3, defaultDt=utils.PWaveTimeStep(), label='timestepper'),
            flow,  #  <====== the solver is inserted here, for the moment it is "dead" (doing nothing)
            newton,
            # some recorders for post-processing
            PyRunner(command="flow.saveVtk(key)", iterPeriod=25, dead=1, label="vtkE"),
            VTKRecorder(recorders=["spheres", "velocity", "stress"], iterPeriod=25, dead=1, fileName=key + '/', label="vtkR")
    ]

    ########## if this is fresh execution, get static equilibrium and save the result for later use  #############
    if not os.path.isfile(filename + ".yade") or newSample:
       O.run(1000, 1)
       while unbalancedForce() > 0.01:
          O.run(100, 1)
       # turn the recorders and the solver on
       vtkE.dead = vtkR.dead = flow.dead = 0
       #add a recorder and define what to plot
       O.engines = O.engines + [PyRunner(command=('myAddPlotData()'), iterPeriod=50, label="recorder")]
       O.save(filename + ".yade")

    O.saveTmp()

    ########## define what to plot ##############
    from yade import plot
    import pylab

    # First, find particle at the top of the sample (for evaluating initial height of the layer)
    maxY = 0
    for s in O.bodies:
       if isinstance(s.shape, Sphere):
          pos = s.state.pos
          if pos[1] > maxY:
             maxY = pos[1]

    def myAddPlotData():
       index = flow.getCell(0.5 * width, height, 0.5 * depth)
       if index > 0:
          ########## find particle at the top of the sample ##############
          simpleH = 0
          for s in O.bodies:
             if isinstance(s.shape, Sphere):
                pos = s.state.pos
                if pos[1] > simpleH:
                   simpleH = pos[1]
          ########## function to compute hf ##############
          cavityh = height
          for s in O.bodies:
             v = s.state.vel
             magvel = pow((v[0] * v[0] + v[1] * v[1] + v[2] * v[2]), 0.5)
             if magvel > 0.14:
                pos = s.state.pos
                if pos[1] > cavityh:
                   cavityh = pos[1]

          #########################
          plot.addData(
                  t=O.time,
                  i=O.iter,
                  p=flow.getPorePressure((0.5 * width, height, 0.5 * depth)),
                  q=-injectedFlux,
                  Ho=maxY - 1,
                  hf=cavityh - 1,
                  H=simpleH - 1
          )

          H = simpleH - 1
          hf = cavityh - 1

    ################ impose the costant flux  ###############

    # In this function we find the elements of the mesh which have a face in the injection region, to distribute the inlet flux
    # it is inserted in the solver below
    def imposeFlux(valF):
       found = 0
       listF1 = []
       for k in range(flow.nCells()):
          if 0 in flow.getVertices(k) and flow.getCellCenter(k)[0] > ((width - aperture) / 2.) and flow.getCellCenter(k)[0] < ((width + aperture) / 2.):
             listF1.append(k)
       flow.clearImposedFlux()
       if len(listF1) == 0:
          flow.imposeFlux((0.5 * width, height, 0.5 * depth), valF)
       else:
          for k in listF1:
             flow.imposeFlux(flow.getCellBarycenter(k), valF / float(len(listF1)))

    injectedFlux = -1  #this one is large and it will fluidize violently, see below for smoother evolutions

    # very important: the flux needs to be imposed with this "hook", which plugs our custom function in the right place in the solving sequence
    flow.blockHook = "imposeFlux(injectedFlux)"

    ###############  EXERCISE ###################
    # 1- use trial and error in the shell (see the note below) to find approximately an upper-bound of "injectedFlux", below which there is only limited movements of the particles
    # 2- implement in this script a progressive increase of the flux, starting from a small fraction of max value above, and exceeding it by a factor 2 at the end
    # 3- choose an approriate plot to show that the pressure response is initially linear, then sublinear. Is the upper-bound from question 1 also an upper-bound for the linear response?
    # 4- compare with fig. 10 from https://doi.org/10.1103/PhysRevE.94.052905 (also available at https://arxiv.org/pdf/1703.02319)
    # 5- use paraview to generate a video similar to /watch?v=gH585XaQEcY (it can be with a constant flux)

    #NOTE:
    # to change the injected flux interactively, change it in the global scope:
    # globals()["injectedFlux"]=-0.03
    # else we have two variables with the same name in different scopes and