Different stress values from different Yade versions when running PeriodicFlowEngine

Asked by Hien N.G. NGUYEN on 2020-07-02

I’m working on a Fluid-DEM coupling simulation script and trying to export two things : (1) the total stress acting at the boundary and (2) the effective stress with changing pore water pressure.

What the script will do :
State I : Execute an iso-compression using PeriTriaxController engine
State II : The PeriTriaxController engine is switched off, the PeriodicFlowEngine is then applied with fluid pressure increasing by stepping value Delta_P = 100 kPa after every 1000 iterations. Once the fluid pressure reaches 500 kPa , the simulation stops and a stress plot will be given.

The script runs fine without any technical error ; however, when running with YADE 2018.02b (installed from Ubuntu repo), the output stress is calculated as it’s supposed to be : the boundary stress is constant while the effective stress decreases according to the imposed fluid pressure; but when running with Yadedaily (version Yade 20200701-4022~5cf8771~bionic1 at spoken time) , I receive a different behaviour: the boundary stress increases gradually.

Can you please tell me what changes have been made in YADE between the two versions that produces such different results on the same script? What modification should be done in latest Yade version to obtain the result given in YADE 2018.2b?

More info:
1) I’m running Ubuntu 18.04.04
2) useSolver was set to zero because Yade 2018.2b gives me error about not-compiled CHOLMOD when I used 3, but I think I should ask in another question about this issue.

The MWE script is below:

from __future__ import print_function
from yade import pack,plot,export,utils
import math

sp=pack.SpherePack()
O.periodic=True

# dimensions of sample (fixed by particle size such as L/D~X)
DIAMETER =1.e-1
RADIUS =0.5*DIAMETER
length =10*(DIAMETER)
height =length
width =length
thickness =length/100.

# microproperties
DENS =2600
E =1.e9
P =0.25
compFRIC =1.
FRIC =30.
TENS =0.
COH =0.

# boundary conditions
PI =1.e5 # for sample preparation: isotropic compaction up to PI pressure
PC =PI # confining pressure for assembly preparation
SN =5.e6 # normal stress

# simulation control
DAMPSHEAR =0.
ITER =1e6
OUT ='test'

#### create sample and loading boxes

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

O.materials.append(CohFrictMat(isCohesive=True,density=DENS,young=E,poisson=P,frictionAngle=radians(0.),normalCohesion=1e100,shearCohesion=1e100,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=DENS,young=E,poisson=P,frictionAngle=radians(compFRIC),normalCohesion=TENS,shearCohesion=COH,label='sphereMat'))

upBox = utils.box(center=(length/2.0,2*height+thickness/2.0,width/2.0),orientation=Quaternion(1,0,0,0),extents=(length,thickness/2.,width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2.0,height-thickness/2.0,width/2.0),orientation=Quaternion(1,0,0,0),extents=(length,thickness/2.,width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
O.bodies.append([upBox,lowBox])

sp.makeCloud((0,height+1.5*RADIUS,0),(length,2*height-1.5*RADIUS,width),rMean=RADIUS,rRelFuzz=0.2,periodic=True)
O.bodies.append([utils.sphere(s[0],s[1],color=(0,0.5,1),material='sphereMat') for s in sp])

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

#### engines

flow=PeriodicFlowEngine(
    isActivated=0 # no flow calculation if this is False
    ,useSolver=0 # value 3 not work with YADE 2018.2b package version
    ,defTolerance=+1 # with this value, triangulation is done according to meshUpdateInterval
    ,meshUpdateInterval=1000
    ,duplicateThreshold=0.5
    ,boundaryUseMaxMin=[0,0,1,1,0,0]
    ,wallIds=[-1,-1,1,0,-1,-1]
    ,wallThickness=thickness
    ,bndCondIsPressure=[0,0,0,0,0,0]
        ,bndCondValue=[0,0,0,0,0,0]
    ,permeabilityFactor=1
    ,viscosity=1e-3
    ,fluidBulkModulus=2.2e9
    ,label='flowEng'
    )

O.engines=[
    ForceResetter()
    ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
    ,InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
    )
    ,flow
    ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
    ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep(),label='timeStepper')
    ,NewtonIntegrator(damping=0.3,label='newton')
    ,PyRunner(command='dataRecorder()',iterPeriod=1,label='recData',dead=True)
]

inputP=0
def dataRecorder():
    global inputP
    h=vol=vol_s=nb_s=0.
    h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
    vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
    contactStress=getStress(vol)
    for o in O.bodies:
        if isinstance(o.shape,Sphere):
            nb_s += 1
            vol_s += 4.*pi/3.*(o.shape.radius)**3
        n = 1-vol_s/vol
        nbFrictCont=0.
    for i in O.interactions:
        if i.isReal and i.phys.cohesionBroken:
            nbFrictCont+=1
    plot.addData(
        iter=O.iter
        ,stress_upWall0=abs(O.forces.f(0)[0]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall1=abs(O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall2=abs(O.forces.f(0)[2]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
        ,contactStress00=(contactStress[0,0]),contactStress01=(contactStress[0,1]),contactStress02=(contactStress[0,2]),contactStress10=(contactStress[1,0]),contactStress11=abs(contactStress[1,1]),contactStress12=(contactStress[1,2]),contactStress20=(contactStress[2,0]),contactStress21=(contactStress[2,1]),contactStress22=(contactStress[2,2])
        ,x=O.bodies[0].state.pos[0]
        ,pf=inputP
        ,height=h
        ,volume=vol
        ,porosity=n
        ,k=2.0*nbFrictCont/nb_s
        ,Ek=kineticEnergy()
        ,unbF=unbalancedForce()
    )

phase=0
def triaxDone():
    global phase
    volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])
    if phase==0:
        h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
        vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
        contactStress=getStress(vol)
        vol_s=Rmean=Rmax=nbSph=0
        Rmin=1e6
        for o in O.bodies:
            if isinstance(o.shape,Sphere):
                nbSph+=1
                Rmean+=o.shape.radius
                if o.shape.radius>Rmax: Rmax=o.shape.radius
                if o.shape.radius<Rmin: Rmin=o.shape.radius
                vol_s += 4.*pi/3.*(o.shape.radius)**3
        Rmean=Rmean/nbSph
        n = 1-vol_s/vol
        print('DONE! iter =',O.iter,'| sample generated: nb spheres =',nbSph,', Rmean =',Rmean,', Rratio =',Rmax/Rmin,', porosity =',n)
        print('Changing contact properties now')
        utils.setContactFriction(radians(FRIC))
        print('APPLYING CONFINING PRESSURE: sx,sy and sz will go to PC =',PC)
        triax.goal=(-PC/volRatio,-PC/volRatio,-PC/volRatio) # - - -
        phase+=1
    elif phase==1:
        print('DONE! iter =',O.iter,'| isotropic confinement done: stresses =',volRatio*triax.stress)
        triax.dead=True
        O.pause()

#### Initialization
print('SAMPLE PREPARATION!')

O.run(1000000,1)
O.step()

#### Applying normal stress
print('NORMAL LOADING! iter =',O.iter)

stage=0
stiff=fnPlaten=currentSN=0.
def servo():
    global stage,stiff,fnPlaten,currentSN
    if stage==0:
        currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
        unbF=unbalancedForce()
        boundaryVel=copysign(min(0.1,abs(0.5*(currentSN-SN))),currentSN-SN)
        O.bodies[0].state.vel[1]=boundaryVel
        if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
            stage+=1
            fnPlaten=O.forces.f(0)[1]
            print('Normal stress =',currentSN,' | unbF =',unbF)
            for i in O.interactions.withBody(O.bodies[0].id):
                stiff+=i.phys.kn
            print('Normal stiffness =',stiff)
            print('DONE! iter =',O.iter)
            O.pause()
    if stage==1:
        fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
        boundaryVel=copysign(min(0.1,abs(0.35*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired) # why 0.35?
        O.bodies[0].state.vel[1]=boundaryVel

O.engines = O.engines[:4]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[4:]
O.run(1000000,1)

recData.dead=False
O.run(1,1)

#### injecting fluid
print('FLUID! iter =',O.iter)
flowEng.isActivated=1
flowEng.boundaryUseMaxMin=[1,1,0,0,1,1]
flowEng.bndCondIsPressure=[0,0,1,0,0,0]
newton.damping=0.

O.run(1,1)

# automatic
iter0=O.iter
iterMax=5e3
deltaP=1e5
while O.iter<(iter0+iterMax):
  O.run(1000,True)
  inputP+=deltaP
  flowEng.bndCondValue=[0.,0.,+inputP,0.,0.,0.]
  flowEng.updateBCs()
  print('updateBCs! inputP=',inputP)
  O.run(1,1)

plot.plots={'iter':('pf','stress_upWall1','contactStress11')}

plot.plot()

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-07-06
Last reply:
2020-07-21
Jérôme Duriez (jduriez) said : #1

Hi,

To inspect changes, git diff [*] source code is the way to go. After identifying which commit 2018.02b corresponds to.

[*] https://git-scm.com/docs/git-diff

Robert Caulk (rcaulk) said : #2

>>2) useSolver was set to zero because Yade 2018.2b gives me error about not-compiled CHOLMOD when I used 3, but I think I should ask in another question about this issue.

Please recompile correctly with CHOLMOD and redo your comparison using useSolver=3 or 4. useSolver = 0 is really not supported or check tested anymore.

Hien N.G. NGUYEN (rin) said : #3

Many thanks for the comments.

@Jerome: I've never done that before, can you please show me a simple example?

@Robert, as you suggested, I downloaded the source code of Yade 2018.2b and Yade 2019.01a and compile with CHOLMOD flag ON, the script runs fine with useSolver=3 and gives me the correct stress values.

However, for Yadedaily 2020-ish, the stresses still act strangely no matter how I configure the useSolver value, I tried both values 0 and 3.

Robert Caulk (rcaulk) said : #4

>and compile with CHOLMOD flag ON

What does that mean? We don't have an explicit cholmod flag during compilation. Please indicate your cmake command and cmake command output.

>Yadedaily 2020-ish

What does that mean? I dont think we have a "Yadedaily 2020-ish" branch.

Hien N.G. NGUYEN (rin) said : #5

Just to clarify, the flag I turned on for the compilation is CHOLMOD_GPU, which is not related to this issue, my bad.

Below is the cmake command and the cmake command output of YADE-2019.01a on my PC

cmake -DCMAKE_INSTALL_PREFIX=/home/AD/nguyen39/Programmes/YADE2019/install /home/AD/nguyen39/Programmes/YADE2019/trunk -DCHOLMOD_GPU=ON -DSUITESPARSEPATH=/usr/local/SuiteSparse-5.7.2
-- The C compiler identification is GNU 7.5.0
-- The CXX compiler identification is GNU 7.5.0
-- Check for working C compiler: /usr/bin/cc
-- Check for working C compiler: /usr/bin/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Check for working CXX compiler: /usr/bin/c++
-- Check for working CXX compiler: /usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Found PythonInterp: /usr/bin/python (found version "2.7.17")
-- Found OpenMP_C: -fopenmp (found version "4.5")
-- Found OpenMP_CXX: -fopenmp (found version "4.5")
-- Found OpenMP: TRUE (found version "4.5")
qmake: could not exec '/usr/lib/x86_64-linux-gnu/qt4/bin/qmake': No such file or directory
-- Found unsuitable Qt version "" from NOTFOUND
-- Found PkgConfig: /usr/bin/pkg-config (found version "0.29.1")
-- Version is set to 2019.01a
-- Found OpenGL: /usr/lib/x86_64-linux-gnu/libOpenGL.so
-- GTS using gts-config /usr/bin/gts-config
-- Using GTS from /usr
-- Found GL2PS: /usr/lib/x86_64-linux-gnu/libgl2ps.so
-- Found CGAL: /usr/include, /usr/lib/x86_64-linux-gnu/libCGAL.so
-- Found NumPy: version "1.13.3" /usr/lib/python2.7/dist-packages/numpy/core/include
-- Found Loki: /usr/include
-- GCC Version >= 4.8. Adding -ftrack-macro-expansion=0
-- GCC Version >= 4.8. Adding -save-temps
-- GCC Version >= 4.9. Adding -fstack-protector-strong
-- Looking for pthread.h
-- Looking for pthread.h - found
-- Looking for pthread_create
-- Looking for pthread_create - not found
-- Looking for pthread_create in pthreads
-- Looking for pthread_create in pthreads - not found
-- Looking for pthread_create in pthread
-- Looking for pthread_create in pthread - found
-- Found Threads: TRUE
-- Boost version: 1.65.1
-- Found the following Boost libraries:
-- python
-- thread
-- filesystem
-- iostreams
-- regex
-- serialization
-- system
-- date_time
-- chrono
-- atomic
-- Boost_VERSION: 106501
-- Boost_LIB_VERSION: 1_65_1
-- Boost_INCLUDE_DIRS: /usr/include
-- Boost_LIBRARIES: /usr/lib/x86_64-linux-gnu/libboost_python.so/usr/lib/x86_64-linux-gnu/libboost_thread.so/usr/lib/x86_64-linux-gnu/libboost_filesystem.so/usr/lib/x86_64-linux-gnu/libboost_iostreams.so/usr/lib/x86_64-linux-gnu/libboost_regex.so/usr/lib/x86_64-linux-gnu/libboost_serialization.so/usr/lib/x86_64-linux-gnu/libboost_system.so/usr/lib/x86_64-linux-gnu/libboost_date_time.so/usr/lib/x86_64-linux-gnu/libboost_chrono.so/usr/lib/x86_64-linux-gnu/libboost_atomic.so/usr/lib/x86_64-linux-gnu/libpthread.so
-- Found Eigen3: /usr/include/eigen3 (Required is at least version "2.91.0")
-- Found BZip2: /usr/lib/x86_64-linux-gnu/libbz2.so (found version "1.0.6")
-- Looking for BZ2_bzCompressInit
-- Looking for BZ2_bzCompressInit - found
-- Found ZLIB: /usr/lib/x86_64-linux-gnu/libz.so (found version "1.2.11")
-- Found PythonLibs: /usr/lib/x86_64-linux-gnu/libpython2.7.so (found version "2.7.17")
-- Found Eigen3, version: 3.3.4
-- Disable vectorization
-- The imported target "vtkRenderingPythonTkWidgets" references the file
   "/usr/lib/x86_64-linux-gnu/libvtkRenderingPythonTkWidgets.so"
but this file does not exist. Possible reasons include:
* The file was deleted, renamed, or moved to another location.
* An install or uninstall procedure did not complete successfully.
* The installation package was faulty and contained
   "/usr/lib/cmake/vtk-6.3/VTKTargets.cmake"
but not all the files it references.

-- The imported target "vtk" references the file
   "/usr/bin/vtk"
but this file does not exist. Possible reasons include:
* The file was deleted, renamed, or moved to another location.
* An install or uninstall procedure did not complete successfully.
* The installation package was faulty and contained
   "/usr/lib/cmake/vtk-6.3/VTKTargets.cmake"
but not all the files it references.

-- Found VTK
-- Found OpenMP_C: -fopenmp (found version "4.5")
-- Found OpenMP_CXX: -fopenmp (found version "4.5")
-- Found OpenMP
-- GTS using gts-config /usr/bin/gts-config
-- Using GTS from /usr
-- Checking for one of the modules 'glib-2.0'
-- Found GLib2: glib-2.0 /usr/include/glib-2.0;/usr/lib/x86_64-linux-gnu/glib-2.0/include
-- Looking for include file glib/gregex.h
-- Looking for include file glib/gregex.h - not found
-- Found GTS
-- Found GLUT: /usr/lib/x86_64-linux-gnu/libglut.so
-- USE QT5
-- Found QGLVIEWER-qt5: /usr/include/QGLViewer
-- Found GUI-Qt5-LIBS
-- GMP libs: /usr/lib/x86_64-linux-gnu/libgmp.so /usr/lib/x86_64-linux-gnu/libgmpxx.so
-- Found GMP: /usr/include/x86_64-linux-gnu
-- GMP libs: /usr/lib/x86_64-linux-gnu/libgmp.so /usr/lib/x86_64-linux-gnu/libgmpxx.so
-- Found CGAL
-- SuiteSparse version 5.1.2 found, CHOLMOD direct solver for CPU activated.
-- Found Cholmod: /usr/lib/x86_64-linux-gnu/libcholmod.so
-- Found OpenBlas: /usr/lib/x86_64-linux-gnu/libopenblas.so
-- Found Metis: /usr/include
-- Found Cholmod in /usr/lib/x86_64-linux-gnu/libcholmod.so
-- Found OpenBlas in /usr/lib/x86_64-linux-gnu/libopenblas.so
-- Found Metis in /usr/lib/x86_64-linux-gnu/libmetis.so
-- Found CuBlas: /usr/lib/x86_64-linux-gnu/libcublas.so
-- Found Lapack: /usr/lib/x86_64-linux-gnu/liblapack.so
-- Found CuBlas in /usr/lib/x86_64-linux-gnu/libcublas.so
-- Found Lapack in /usr/lib/x86_64-linux-gnu/liblapack.so
-- Found GL2PS
LBMFLOW is still experimental, building and running LBM engine are at your own risk!
Yade will be installed to /home/AD/nguyen39/Programmes/YADE2019/install
-- Suffix is set to -2019.01a
-- LIBRARY_OUTPUT_PATH is set to lib/x86_64-linux-gnu
-- runtimePREFIX is set to /home/AD/nguyen39/Programmes/YADE2019/install
-- Found gts: /usr/lib/python2.7/dist-packages/gts
-- Use system gts version
-- Found minieigen: /usr/lib/python2.7/dist-packages/minieigen.so
-- Found Tkinter: /usr/lib/python2.7/lib-tk/Tkinter.pyc
-- VTK version >5 and <8 is found
-- ===========================================================
-- Yade configured with following features: Odeint VTK OpenMP GTS GUI-Qt5 CGAL PFVFLOW LINSOLV CHOLMOD_GPU TWOPHASEFLOW GL2PS LBMFLOW
-- Disabled features: SPH DEFORM LIQMIGRATION MASK_ARBITRARY THERMAL PROFILING PotentialParticles PotentialBlocks
-- Optimized build
-- ===========================================================
-- Configuring done
-- Generating done
CMake Warning:
  Manually-specified variables were not used by the project:

    SUITESPARSEPATH

-- Build files have been written to: /home/AD/nguyen39/Programmes/YADE2019/bin

We have tried with earlier Yade versions (compiled from source) including v1.20, v2018.02b, v2019.01a and the stress output was correct.

However, both Yadedaily binary install from `sudo apt install yadedaily` and compiled from source (version Yade -2020.01a) provided incorrect stress values.

Anyway, currently the script works fine with Yade-2019.01a, so I will stick with this version. This question is considered solved. About what changes have been exactly made, I might make another query when I get a hand on `git diff` as Jerome suggested.
Thank you.

Launchpad Janitor (janitor) said : #6

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Robert Caulk (rcaulk) said : #7

Changing to answered to avoid the automatic launchpad janitor cleanup.

Can you help with this problem?

Provide an answer of your own, or ask Hien N.G. NGUYEN for more information if necessary.

To post a message you must log in.