How to do undraind test with PFV (FlowEngine)?

Hi,
I've been trying to include water in the undrained test wih FlowEngine, but failed.

The strain will not increase as it is told to. In fact, it will not change at all.

Everyting is ok if the flowengine is not included.

Can someone tell me why? Or can someone give me some example code? Thanks a lot!

Here is my code:

# -*- coding: utf-8 -*-
#*************************************************************************
# Copyright (C) 2010 by Bruno Chareyre *
# bruno.chareyre_at_grenoble-inp.fr *
# *
# This program is free software; it is licensed under the terms of the *
# GNU General Public License v2 or later. See file LICENSE for details. *
#*************************************************************************/

## This script details the simulation of a triaxial test on sphere packings using Yade
## See the associated pdf file for detailed exercises
## the algorithms presented here have been used in published papers, namely:
## * Chareyre et al. 2002 (http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)
## * Chareyre and Villard 2005 (https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
## * Scholtès et al. 2009 (http://dx.doi.org/10.1016/j.ijengsci.2008.07.002)
## * Tong et al.2012 (http://dx.doi.org/10.2516/ogst/2012032)
##
## Most of the ideas were actually developped during my PhD.
## If you want to know more on micro-macro relations evaluated by triaxial simulations
## AND if you can read some french, it is here: http://tel.archives-ouvertes.fr/docs/00/48/68/07/PDF/Thesis.pdf

from yade import pack

## control to open the file
Not_Open_flag = 1
if Not_Open_flag:
   stress_info = open("/home/caowei/Documents/myYade/trunk/examples/stress_info.txt",'a')
   contact_info = open("/home/caowei/Documents/myYade/trunk/examples/contactOri.txt",'a')
   particle_info = open("/home/caowei/Documents/myYade/trunk/examples/particle.txt",'a')
   Not_Open_flag = 0
############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
 num_spheres=700,# number of spheres
 compFricDegree = 35, # contact friction during the confining phase
 key='_triax_base_', # put you simulation's name here
 unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity = 0.425 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 35 # contact friction during the deviatoric loading
rate=-0.01 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=1e8 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing

## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()

clumps=False #turn this true for the same example with clumps
if clumps:
 ## approximate mean rad of the futur dense packing for latter use
 volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
 mean_rad = pow(0.09*volume/num_spheres,0.3333)
 ## define a unique clump type (we could have many, see clumpCloud documentation)
 c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
 ## generate positions and input them in the simulation
 sp.makeClumpCloud(mn,mx,[c1],periodic=False)
 sp.toSimulation(material='spheres')
 O.bodies.updateClumpProperties()#get more accurate clump masses/volumes/inertia
else:
 sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
 O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
 #or alternatively (higher level function doing exactly the same):
 #sp.toSimulation(material='spheres')

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
 ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
 ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
 ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
 ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
 ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
 ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
 stressMask = 7,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton,
 #PyRunner(command='checkStress0()',iterPeriod=500,label='checker')
]
O.dt=.5*utils.PWaveTimeStep()
# for the samples of this size (diametre is around 1.8cm), the O.dt is 8.8e-5 s, which is about a hundred times of
# the small sample, and this leads to much faster calculation speed.

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

## UNCOMMENT THE FOLLOWING SECTIONS ONE BY ONE
## DEPENDING ON YOUR EDITOR, IT COULD BE DONE
## BY SELECTING THE CODE BLOCKS BETWEEN THE SUBTITLES
## AND PRESSING CTRL+SHIFT+D

#######################################
### APPLYING CONFINING PRESSURE ###
#######################################

#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal2=triax.goal3=-100000

while 1:
  O.run(1000, True)
  ##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-100000-triax.meanStress)/100000<0.001:
    break

O.save('confinedState'+key+'.yade.gz')
print "### Isotropic state saved ###"
print triax.porosity
print triax.meanStress
print len(O.bodies)

###################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
###################################################

### We will reach a prescribed value of porosity with the REFD algorithm
### (see http://dx.doi.org/10.2516/ogst/2012032 and
### http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)

import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
 ## we decrease friction value and apply it to all the bodies and contacts
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 ## while we run steps, triax will tend to grow particles as the packing
 ## keeps shrinking as a consequence of decreasing friction. Consequently
 ## porosity will decrease
 O.run(500,1)

O.save('compactedState'+key+'.yade.gz')
print "### Compacted state saved ###"

while 1:
  O.run(100, True)
  ##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-100000-triax.meanStress)/100000<0.001:
    break

##############################
### DEVIATORIC LOADING ###
##############################

##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

##set stress control on x and z, we will impose strain rate on y
triax.stressMask = 0
##now goal2 is the target strain rate
triax.goal2=rate
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-rate / 2.0
triax.goal3=-rate / 2.0

##we can change damping here. What is the effect in your opinion?
newton.damping=0.2

##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
O.saveTmp()

#####################################################
### Example of how to record and plot data ###
#####################################################

from yade import plot

### a function saving variables
def history():
   plot.addData(unbalanced=unbalancedForce(),e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
        ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
      s11=-triax.stress(triax.wall_right_id)[0],
      s22=-triax.stress(triax.wall_top_id)[1],
      s33=-triax.stress(triax.wall_front_id)[2],
                    devi = -triax.stress(triax.wall_top_id)[1] - (-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2]) / 2.0,
                    p = triax.meanStress,
      i=O.iter)
        stress_info.write(" %d " %(O.iter))
        stress_info.write(" %6.4f " %(triax.stress(triax.wall_right_id)[0]))
        stress_info.write(" %6.4f " %(triax.stress(triax.wall_top_id)[1]))
        stress_info.write(" %6.4f " %(triax.stress(triax.wall_front_id)[2]))
        #print("stress:")
        #print(triax.stress(triax.wall_right_id)[0])
        #print(triax.stress(triax.wall_top_id)[1])
        #print(triax.stress(triax.wall_front_id)[2])
        #print("strain:")
        #print(triax.strain[0])
        #print(triax.strain[1])
        #print(triax.strain[2])
        for i in range(3):
           stress_info.write("%6.6f " %(triax.strain[i]))
        stress_info.write("\n")

        if O.iter % 2000 == 0:

           for i in O.interactions:
              contact_info.write("%d " %(i.id1))
              contact_info.write("%d " %(i.id2))
              contact_info.write("%.4f " %(i.phys.normalForce[0]))
              contact_info.write("%.4f " %(i.phys.normalForce[1]))
              contact_info.write("%.4f " %(i.phys.normalForce[2]))
              contact_info.write("%.4f " %(i.phys.shearForce[0]))
              contact_info.write("%.4f " %(i.phys.shearForce[1]))
              contact_info.write("%.4f " %(i.phys.shearForce[2]))
              contact_info.write("\n")

           for i in range(len(O.bodies)):
              if i > 5:
                 particle_info.write("%d " %(i))
                 if O.bodies[i].isClump:
                    particle_info.write("\n")
                 else:
                    particle_info.write("%f " %(O.bodies[i].state.pos[0]))
                    particle_info.write("%f " %(O.bodies[i].state.pos[1]))
                    particle_info.write("%f " %(O.bodies[i].state.pos[2]))
                    particle_info.write("%f " %(O.bodies[i].shape.radius))
                    particle_info.write("%d " %(O.bodies[i].clumpId))
                    particle_info.write("\n")

if 1:
  ## include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=200,command='history()',label='recorder')]+O.engines[5:7]
  ##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
else:
  ## With the line above, we are recording some variables twice. We could in fact replace the previous
  ## TriaxialRecorder
  ## by our periodic engine. Uncomment the following line:
  O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')

### declare what is to plot. "None" is for separating y and y2 axis
plot.plots={'i':('unbalanced',),'i ':('s11','s22','s33'),' i':('e11','e22','e33'),'e22':('ev'),'p':('devi')}
### the traditional triaxial curves would be more like this:
##plot.plots={'e22':('s11','s22','s33',None,'ev')}

## display on the screen (doesn't work on VMware image it seems)
plot.plot()

devi_test = -triax.stress(triax.wall_top_id)[1] - (-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2]) / 2.0

def stress_control1():
  global devi_test
  while devi_test < 50000:
     O.run(100)
     devi_test = -triax.stress(triax.wall_top_id)[1] - (-triax.stress(triax.wall_right_id)[0]\
     -triax.stress(triax.wall_front_id)[2]) / 2.0

def stress_control2():
  global devi_test
  while devi_test > -50000:
     O.run(100)
     devi_test = -triax.stress(triax.wall_top_id)[1] - (-triax.stress(triax.wall_right_id)[0]\
     -triax.stress(triax.wall_front_id)[2]) / 2.0

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.dead=0
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=10
flow.bndCondIsPressure=[0,0,0,0,0,0]
flow.bndCondValue=[0,0,0,0,0,0]
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
O.dt=0.1e-3
O.dynDt=False

O.run(1,1)
Qin = flow.getBoundaryFlux(2)
Qout = flow.getBoundaryFlux(3)
tic_toc = 0
print(O.iter)
print(triax.stressMask)
print(triax.goal1)
print(triax.goal2)
print(triax.goal3)