CU triaxial (PFV)

Hello Sam,

CHOLMOD usually throws this error when the boundary conditions are improperly defined (hence the not positive definite). But it might also throw this error if the time step is too large or if your packing has exploded.

One way or the other, there is no way to know which problem is yours since you have simply posted the error and you have not posted a minimal working example [1].

[1]https://yade-dem.org/wiki/Howtoask

Cheers,

Robert

Dear Robert,

Sorry for the lack of information.
BTW, you can find herewith my script, as you requested:

############################
### 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()]
 ),
 ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
 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
]

#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=-10000

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(-10000-triax.meanStress)/10000<0.001:
    break

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

#####################################################
### 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)

if 1:
  # include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,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')

O.run(100,True)

## declare what is to plot. "None" is for separating y and y2 axis
plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')}
## 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=1.3e-3
flow.fluidBulkModulus = 2e9
flow.imposePressure(Vector3(triax.width/2,triax.height/2,triax.depth/2),0.001)
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)

Regards
Sam

Hello Sam,

Your script is not an MWE.py. It does not have a packing associated with it and it contains syntax errors e.g. NameError: name 'young' is not defined. Consider reading this webpage -->> https://yade-dem.org/wiki/Howtoask <<-- to expedite this process.

Looking at the segment of code that you posted - you are trying to compress an incompressible fluid. The problem is ill-posed. You also appear to be imposing pressure at the center of the system, but your boundary conditions are all no flow. Another ill-posed system.

-rc

Dear Robert,

Thank you for your valuable comment.

BTW, I have just two more questions about your last comment (#5). You mentioned I am trying to compress an incompressible fluid, whereas Water is incompressible fluid and so what do you mean exactly about it??
Also, you said I imposed the pressure at the centre and my boundary conditions are not flow. While I want to model Triaxial CU (Undrained Condition) and water cannot flow out from the assembly??

Regards
Sam

You can find herewith MWE.py, as you requested:

from yade import pack

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

nRead=readParamsFromTable(
 num_spheres=700,
 compFricDegree = 35,
 key='_Triax_PFV_CU_',
 unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.425
compFricDegree = table.compFricDegree
finalFricDegree = 35
rate=-0.01
damp=0.2
stabilityThreshold=0.01
young=1e8
mn,mx=Vector3(0,0,0),Vector3(1,1,1)

O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2650,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')

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

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)

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"),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton,
]
O.dt=.5*utils.PWaveTimeStep()

Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

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

triax.goal1=triax.goal2=triax.goal3=-100000

while 1:
  O.run(1000, True)
  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 ###
###################################################

import sys
while triax.porosity>targetPorosity:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)

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

while 1:
  O.run(100, True)
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-100000-triax.meanStress)/100000<0.001:
    break

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

triax.internalCompaction=False

setContactFriction(radians(finalFricDegree))

triax.stressMask = 0
triax.goal2=rate
triax.goal1=-rate / 2.0
triax.goal3=-rate / 2.0

newton.damping=0.1

O.saveTmp()

Hello Sam,

The script you posted on Comment #6 does not contain a fluid model (it is inactive when you set flow.dead=1). Actually, the script from Comment #6 also does not contain a packing, so it is pretty far from yielding any kind of FlowEngine error like the one you posted in the OP. You are welcome to post the script that causes the FlowEngine error shown in the OP and we can try to help fix it. Or if you are having problems related to a different topic, please open a new thread.

My responses below are w.r.t to your other partial script posted in Comment #4, since that is the script pertaining to your questions.

>>
You mentioned I am trying to compress an incompressible fluid, whereas Water is incompressible fluid and so what do you mean exactly about it??
>>

Actually, looking at your script from Comment #4, I see you are choosing to model the water as a compressible fluid when you set flow.fluidBulkModulus=2e9 [1]. So I was mistaken when I said you were trying to compress an incompressible fluid. Instead, you are compressing a compressible fluid so your problem may actually be well-posed. But you still have not posted an MWE that yields the error shown in the OP (is that what we are trying to fix?), so I am still speculating.

>>
Also, you said I imposed the pressure at the centre and my boundary conditions are not flow. While I want to model Triaxial CU (Undrained Condition) and water cannot flow out from the assembly??
>>

I'm sorry, what is the question? I will just give information in an attempt to answer: you are choosing to model the boundaries as no-flux boundaries when you set flow.bndCondIsPressure=[0,0,0,0,0,0] [2]. Similarly, you are imposing a constant pressure in the center of your system when you set flow.imposePressure(Vector3(triax.width/2,triax.height/2,triax.depth/2),0.001) [3].

Cheers,

Robert

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.fluidBulkModulus
[2https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.bndCondIsPressure
[3]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.imposePressure

Thank you, Robert