CHOLMOD warning: matrix not positive definite with PeriodicFlowEngine

Asked by Katia Boschi on 2021-03-02

Dear All,

I have tried to run the following code (simple shear test upon a saturated granular sample under constant pressure conditions):

from yade import pack,qt,utils,ymport
from yade import export
from yade import timing
from yade import plot
from builtins import range
import numpy
import time
from math import *

O.load('prova.xml')

num_spheres=5000
young=750.e6
finalFricDegree = 0.35
graindensity=2600.0
poisson=0.17
en=0.9
sigmaIso=-150.e3
radiusR=0.001/2.
D50=radiusR*2.0

#O.materials.append(FrictViscoMat(young=young,poisson=poisson,frictionAngle=compFricDegree,density=graindensity,betan=0.0336,label='spheres'))

triax=PeriTriaxController(
 dynCell=True,
 mass=0.2,
 goal=[sigmaIso,sigmaIso,sigmaIso],
 stressMask=7,
 maxStrainRate=[0.0001,0.0001,0.0001],
 label='triax')

newton=NewtonIntegrator(damping=0.1)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
         [Ig2_Sphere_Sphere_ScGeom()],
         [Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys()],
         [Law2_ScGeom_FrictViscoPhys_CundallStrackVisco(traceEnergy=True)]
 ),
 PeriodicFlowEngine(dead=1,label="flow"),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
 triax,
 newton
]

while 1:
 O.run(1000, True)
 unb=unbalancedForce()
 if unb<0.00001:
  break

triax.dead=0
flow.dead=0

flow.pressureForce=True
flow.viscousShear=True
flow.viscousNormalBodyStress=True
flow.viscousShearBodyStress=True
flow.normalLubrication=True
flow.shearLubrication=True

flow.defTolerance=0.1
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=0.001
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,0]
GlobalStiffnessTimeStepper.dead=True
O.dt=min(0.5*PWaveTimeStep(),0.5*1./1200.*pi/flow.viscosity*graindensity*(D50/2.)**2)
O.dynDt=False

O.run(10,1)

while 1:
 O.run(5000, True)
 unb=unbalancedForce()
 if unb<0.0001:
  if O.cell.velGrad[0,1]<0.01:
   O.cell.velGrad=O.cell.velGrad+Matrix3(0,0.001,0, 0,0,0, 0,0,0)
   O.dt=min(0.5*PWaveTimeStep(),0.5*1./1200.*pi/flow.viscosity*graindensity*(D50/2.)**2)
  if O.cell.velGrad[0,1]>0.01:
   break

with "prova.xml" a dry granular sample properly generated in a periodic cell.

Even if I switch off both the PeriTriaxialController and the cell strain rate imposition, after some cycles, the following error appears:

CHOLMOD warning: matrix not positive definite. file: ../Supernodal/t_cholmod_super_numeric.c line: 911
something went wrong in Cholesky factorization, use LDLt as fallback this time1

Would you mind helping me with this error?

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Katia Boschi
Solved:
2021-03-22
Last query:
2021-03-22
Last reply:
2021-03-05

Hi Katia, I didn't check the details but did you impose a pressure somewhere? Else with an incompressible fluid the problem is ill-posed (and then singular matrix is expected).
B

Katia Boschi (katia.boschi) said : #2

Dear Bruno,

yes! I have imposed a zero pressure in the z direction for this reason.

Katia

Robert Caulk (rcaulk) said : #3

Have you tried decreasing the timestep? Your definition of O.dt looks suspicious to me.

Cheers,

Robert

Katia Boschi (katia.boschi) said : #4

Dear Robert,

yes, I have tried to change it: I have both decreased it several orders of magnitude and kept the GlobalStiffnessTimeStepper active.

Katia

Katia Boschi (katia.boschi) said : #5

Dear All,

I attach the same script with few particles' sample generation:
(Even by running this script with the latest yadedaily package stable version (yade-2021.01a), I get the same error. )

from yade import pack,qt,utils,ymport
from yade import export
from yade import timing
from yade import plot
from builtins import range
import numpy
import time
from math import *

num_spheres=200
young=750.e6
finalFricDegree = 0.35
graindensity=2600.0
poisson=0.17
en=0.9
sigmaIso=-150.e3
radiusR=0.001/2.
D50=radiusR*2.0

O.periodic=True
O.cell.setBox(0.01,0.01,0.01)
mn,mx=Vector3(0,0,0),Vector3(0.01,0.01,0.01)

O.materials.append(FrictViscoMat(young=young,poisson=poisson,frictionAngle=finalFricDegree,density=graindensity,betan=0.0336,label='spheres'))
sp=pack.SpherePack()
sp.makeCloud(mn, mx, rMean=radiusR, rRelFuzz=0.0, num=num_spheres, periodic=True, seed=1)
sp.toSimulation(material='spheres')

triax=PeriTriaxController(
 dynCell=True,
 mass=0.2,
 goal=[sigmaIso,sigmaIso,sigmaIso],
 stressMask=7,
 globUpdate=5,
 maxStrainRate=[1,1,1],
 label='triax')

newton=NewtonIntegrator(damping=0.1)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
         [Ig2_Sphere_Sphere_ScGeom()],
         [Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys()],
         [Law2_ScGeom_FrictViscoPhys_CundallStrackVisco(traceEnergy=True)]
 ),
 PeriodicFlowEngine(dead=1,label="flow"),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
 triax,
 newton
]

while 1:
 O.run(1000, True)
 unb=unbalancedForce()
 if unb<0.0001:
  if triax.relStressTol<0.1:
   break

flow.dead=0
flow.pressureForce=True
flow.viscousShear=True
flow.viscousNormalBodyStress=True
flow.viscousShearBodyStress=True
flow.normalLubrication=True
flow.shearLubrication=True

flow.defTolerance=0.1
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=0.0001
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,0]
GlobalStiffnessTimeStepper.dead=True
O.dt=min(0.5*PWaveTimeStep(),0.5*1./1200.*pi/flow.viscosity*graindensity*(D50/2.)**2)
#O.dynDt=True

O.run(10,1)

while 1:
 O.run(5000, True)
 unb=unbalancedForce()
 if unb<0.0001:
  if O.cell.velGrad[0,1]<0.01:
   O.cell.velGrad=O.cell.velGrad+Matrix3(0,0.005,0, 0,0,0, 0,0,0)
   O.dt=min(0.5*PWaveTimeStep(),0.5*1./1200.*pi/flow.viscosity*graindensity*(D50/2.)**2)
  if O.cell.velGrad[0,1]>0.01:
   break

Katia Boschi (katia.boschi) said : #6

I have defined boundary conditions on hypothetical walls, but in fact there are no walls...