Yade

simulation explodes as soon as flow engine is activated

Asked by Luc Scholtès

Hello there,

I have been struggling with this for a while now and could not figure out what is wrong. My simulation (shear test with periodic boundary conditions) explodes as soon as the flow engine is activated (see the last 2 lines) and I don't know if it is due to the technical problems that I am facing with the flow engine or if there is something wrong with my script. Could someone please try it and see if it works on his side?

Many thanks

Luc

###

from yade import pack

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

# dimensions of sample
RADIUS=0.05
length=10*(2*RADIUS)
height=length
width=length
thickness=RADIUS/10.

# microproperties
DENS=2600
E=1e8
P=0.5
compFRIC=10.
FRIC=30.
TENS=0.
COH=0.

# loading conditions
SN=5e6
RATE=0.1

####

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=(0,2*height+thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(0,height-thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,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,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
 ,useSolver=3
 ### what should we define in the next 2 lines?
 ,defTolerance=-1 # with this value, triangulation is done according to meshUpdateInterval
 ,meshUpdateInterval=1000 # this value by default
 ,duplicateThreshold=0.5 # how do we define this?
 ,boundaryUseMaxMin=[0,0,0,0,0,0] # what is this for? should we have [0,0,1,1,0,0] because of the top and bottom walls?
 ,wallIds=[-1,-1,1,0,-1,-1] # top wall id=0, bottom wall id=0
 ,wallThickness=thickness
 ,bndCondIsPressure=[1,1,0,0,1,1] # drained in the horizontal plane
        ,bndCondValue=[0,0,0,0,0,0]
 ,permeabilityFactor=1.e-3
 ,viscosity=1.
 ,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()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep())
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-1.e5/volRatio,-1.e5/volRatio,-1.e5/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='compactionDone()',label='triax')
 ,flow
 ,NewtonIntegrator(damping=0.2,label='newton')
]

def compactionDone():
 print 'Changing contact properties now'
 tt=TriaxialCompressionEngine()
 tt.setContactProperties(FRIC)
 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)

#### 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 'DONE! iter=',O.iter
   O.pause()
 if stage==1:
  fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  boundaryVel=copysign(min(RATE,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)

print 'Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
print 'Normal stress (contacts) = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1]

#### preparing for shearing
print 'Gluing spheres to boundary walls'

## gluing particles in contact with the walls
for i in O.interactions:
  if i.isReal:
   if isinstance(O.bodies[i.id1].shape,Box):
    O.bodies[i.id2].shape.color[0]=1
   if isinstance(O.bodies[i.id2].shape,Box):
    O.bodies[i.id1].shape.color[0]=1

for i in O.interactions:
  if i.isReal and ( O.bodies[i.id1].shape.color[0]==1 and O.bodies[i.id2].shape.color[0]==1 ):
   O.bodies[i.id1].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
   O.bodies[i.id2].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
   O.bodies[i.id1].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
   O.bodies[i.id2].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
   i.phys.initCohesion=True

#### shearing
print 'SHEARING! iter=',O.iter

flow.isActivated=1 # sample explodes as soon as fluid is activated???
O.bodies[0].state.vel[0]=RATE

O.run(10)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Luc Scholtès
Solved:
Last query:
Last reply:
Revision history for this message
Robert Caulk (rcaulk) said : 		#1

Hey Luc,

I tried your script, and I can confirm the explosion.

I managed to stabilized the simulation by reducing the timestep:

flow.isActivated=1 # sample explodes as soon as fluid is activated???
O.dynDt=False
O.dt =1e-8

Let me know if that works.

Best,

Robert

Revision history for this message
lingran.zhang (lingran-zhang) said : 		#2

Hello Luc,

When running your script with reduced time step, my simulation explodes shortly (in several time steps) after the flow engine is activated (I guess yade on both of our computers are compiled from the same source code).

Regards,
Lingran

Revision history for this message
Luc Scholtès (luc) said : 		#3

Hi Robert (and Lingran),

Unfortunately, even with a timestep equal to 1e-9, I get the same behavior (explosion)... In addition, I now get this at the first iteration and at different iteration before the explosion:

CHOLMOD warning: matrix not positive definite
something went wrong in Cholesky factorization, use LDLt as fallback this time1

I get the same behavior with yadedaily and yade2017.

I must be doing something wrong.

Thanks anyway for your interest.

Luc

Revision history for this message
Robert Caulk (rcaulk) said : 		#4

Hey Luc,

I receive the same error as you after running the simulation out :-(

I think this:

CHOLMOD warning: matrix not positive definite
something went wrong in Cholesky factorization, use LDLt as fallback this time1

suggests improper boundary conditions.

So I think you have a problem with one of these:

 ,boundaryUseMaxMin=[0,0,0,0,0,0] # what is this for? should we have [0,0,1,1,0,0] because of the top and bottom walls?
 ,wallIds=[-1,-1,1,0,-1,-1] # top wall id=0, bottom wall id=0

I am unsure how boundary conditions should be properly assigned for periodic simulations that do not use 6 walls, but my quick glance at the source code makes me think the -1 for wall ids is preventing yade from adding boundary conditions to those sides of the cube.

I will try to look a bit deeper at it this weekend.

Cheers,

Robert

Revision history for this message
Luc Scholtès (luc) said : 		#5

Yes, I figured it could be a problem with the BC...

My model is a cube set up with (xmin, xmax, ymin, ymax) as periodic boundaries and (zmin,zmax) as non periodic (loading walls).

It seem that the only configurations that do not produce an explosion is when I apply a pressure gradient in the z direction and a no flux condition in the x and y directions (the simulation does not crash). This case corresponds to boundaryUseMaxMin=[0,0,1,1,0,0] or when boundaryUseMaxMin=[0,0,0,1,0,0] (I used that one in a previous study).

Unfortunately, I am interested in a case where a pressure gradient is applied along the x direction (boundaryUseMaxMin=[1,1,0,0,0,0]) and this case always crashes.

I have conceptual difficulties with periodic BC and I am probably doing something wrong... I have seen results (another DEM code) with a similar set up though...

I don't know if I have to rethink my simulation or if it is a limitation of the code (applying a pressure gradient across periodic BC)...

Luc

Revision history for this message
Jérôme Duriez (jduriez) said : 		#6

Hi,

Maybe it's worth for you to give a look to http://people.3sr-grenoble.fr/users/bchareyre/pubs/Marzougui2015.pdf and its § "Periodic Boundary Conditions" ?

This suggests applying a pressure gradient along a periodic boundary (if this is what you're interested in ?) has been done at one point ?

Jérôme

Revision history for this message
Luc Scholtès (luc) said : 		#7

Thank you Jerome. The link you sent me is broken (I get a picture of a squirrel instead...) but I could find the article anyway.

I'll have a look and do some testing (there is a dedicated function in the periodic flow engine to apply a pressure gradient but I don't know yet how to combine it with the BCs).

I'll close this question and may open another one specifically asking for

how to apply a pressure gradient along periodic boundary conditions?

later on.

Best,

Luc

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#8

Hi Luc,
First, the pressure gradient is not "along periodic boundary conditions". It is just the pressure gradient, a vector.
Second, what made you think that boundaryUseMaxMin=[0,0,1,1,0,0] was somehow related to defining the pressure gradient?? Is there something in the doc along this line (I shall then fix)?
boundaryUseMaxMin defines if the problem is (=0) bounded by a stateful yade.body (a plane typically) which displacements are reflected in the flow problem, or (=1) by an imaginary static plane placed at the min/max particles coordinates (its position is adjusted in time but its velocity is always zero from the fluid point of view).
If boundaryUseMaxMin=1 of course the corresponding id of the wall need to be defined.
If wallIds=[-1,-1,1,0,-1,-1] and boundaryUseMaxMin=[1,1,0,0,0,0] it means the position of the boundaries at x- and x+ are retrieved by body ids -1 and -1.
I guess you can imagine what happens after that (well, I can't imagine really, but nothing good for sure).
Cheers
Bruno

Revision history for this message
Luc Scholtès (luc) said : 		#9

Hi Bruno,

Thank you for your message. I figured out the pressure gradient thanks to the reference given byJerome. Also, I made a mistake in my previous message regarding the names of the attributes. Sorry for that.

The problem I have is not related to boundaryUseMaxMin but to bndCondIsPressure (at least, that's my guess).

Actually, I set up boundaryUseMaxMin=[0,0,0,0,0,0] even though I have 2 boundary walls at ymin and ymax and it seems to work as long as wallIds=[-1,-1,1,0,-1,-1] (if wallIds is not specified, the simulation crashes with "periodicity is broken....").

What I would like to impose is bndCondIsPressure=[1,1,0,0,1,1] (with bndCondValue=[0,0,0,0,0,0])
because this is what I thought would be representative of a test with no flux conditions along ymin and ymax (my loading walls) and drainage in the (x,z) plane but it seems that I am doing something wrong.

Luc

Revision history for this message
Luc Scholtès (luc) said : 		#10

Based on my latest searches and testing, here are some conclusions (please correct me if I am wrong):

1) Due to the periodic boundaries in x and z, assigning bndCondIsPressure=[0,0,X,X,0,0], [1,1,X,X,1,1] , [0,0,X,X,1,1] or [1,1,X,X,0,0] does not make any difference in the problem solved

2) due to 1), assigning bndCondIsPressure=[X,X,1,1,X,X] or [X,X,0,0,X,X] actually defines if the problem is respectively drained or undrained (drainage can also be defined with [X,X,1,0,X,X] or [X,X,0,1,X,X]).

3) due to 2), the case bndCondIsPressure=[X,X,0,0,X,X] cannot work when using an incompressible fluid but works OK when using a compressible one.

4) assigning a pressure gradient with gradP works independently of the bndCondIsPressure

Luc

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#11

Hi,
You are right and the reason for that is that a periodic problem has no boundary, hence no boundary condition.
bndCondIsPressure=[X,X,0,0,X,X] is ok but you need to impose the pressure in one internal pore to fix a reference value (arbitrary).
Adding a constant pressure to a periodic pressure field still results in a valid solution. No unique solution <=> the linear system reflecting the Stokes problem is not inversible if no pressure is imposed anywhere, then crash.
Bruno