CHOLMOD warning: matrix not positive definite

Asked by Zhicheng Gao

Hello,
I reinstalled Yade from the source code,the version is yade-2021-06-08.git-4852d56, and the version of Ubuntu is 18.04. I tried to delete the particles. There is no problem with the FlowEngine, but the PeriodicFlowEngine will give the following warning.
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
Does it matter to the simulation? Here is my code:
##______________ First section, generate sample_________

from __future__ import print_function
from yade import pack, qt, plot
from math import *

nRead=readParamsFromTable(
        ## model parameters
        num_spheres=1000,
        targetPorosity= .4,
        confiningPressure=-100000,
        ## material parameters
        compFricDegree=15,#contact friction during the confining phase
        finalFricDegree=30,#contact friction during the deviatoric loading
        young=2e8,
        poisson=.2,
        density=2600,
        alphaKr=7.5,
        alphaKtw=0,
        competaRoll=.22,
        finaletaRoll=.22,
        etaTwist=0,
        normalCohesion=0,
        shearCohesion=0,
        ## fluid parameters
        fluidDensity=1000,
        dynamicViscosity=.001,
        ## control parameters
        damp=0,
        stabilityThreshold=.001,
        ## output specifications
        filename='suffusion',
        unknowOk=True
)

from yade.params.table import *

O.periodic=True
O.cell.hSize=Matrix3(1,0,0, 0,1,0, 0,0,1)
# create materials for spheres
#shear strength is the sum of friction and adhesion, so the momentRotationLaw=True
O.materials.append(CohFrictMat(alphaKr=alphaKr,alphaKtw=alphaKtw,density=density,etaRoll=competaRoll,etaTwist=etaTwist,frictionAngle=radians(compFricDegree),momentRotationLaw=True,normalCohesion=normalCohesion,poisson=poisson,shearCohesion=shearCohesion,young=young,label='spheres'))

# generate particles packing
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),-1,0.3333,num_spheres,False, 0.95,seed=1)
sp.toSimulation(material='spheres')

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom6D()],
            [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='contact',setCohesionNow=False,setCohesionOnNewContacts=False)],
            [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True)],
 ),
        PeriodicFlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        PeriTriaxController(label='triax',
            # specify target values and whether they are strains or stresses
            goal=(confiningPressure,confiningPressure,confiningPressure), stressMask=7,
            # type of servo-control, the strain rate isn't determined, it shloud check the unbalanced force
            dynCell=True,maxStrainRate=(10,10,10),
            # wait until the unbalanced force goes below this value
            maxUnbalanced=stabilityThreshold,relStressTol=1e-3,
            doneHook='compactionFinished()'
            ),
        NewtonIntegrator(damping=0)
]

import sys
def compactionFinished():
 # after sample preparation, save the state
        O.save('compactedState'+filename+'.yade.gz')
        print('Compacted state saved', 'porosity', porosity())
        # next time, called python command
        triax.doneHook=''
        O.pause()
O.run()
O.wait()

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

O.run(1,1)
csdList=flow.getConstrictionsFull()
print(len(O.bodies),len(csdList),'finished')
flow.dead=1
print(O.bodies[10].shape.radius)

for i in range(10,50):
    O.bodies.erase(i)
print(len(O.bodies))
O.run(1000,True)

flow.dead=0
flow.updateTriangulation=True

O.run(1,1)
csdList1=flow.getConstrictionsFull()
print(len(csdList1))

Question information

Language:
English Edit question
Status:
Expired
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Zhicheng Gao (zhichenggao) said :
#1

Can anyone provide some help, thank you very much!

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

Hi,
It can happen that LDL fails for numerical reason, then the solver switches to LDLt, which is more stable.
I _could_ be a safe warning.

That said, I'm struggling to understand why you would use periodic FlowEngine with (it seems) a non-periodic problem.
Also this makes no sense in periodic conditions since there are no boundaries:
flow.bndCondIsPressure=[0,0,1,1,0,0]

I hope it helps.

Bruno

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#3

Hello, Bruno. Thank you very much for your reply. I'm doing research similar to the paper " A discrete numerical model involving partial fluid-solid coupling to describe suffusion effects in soils "[1]. In this paper, the periodic boundary condition is chosen and it says "In order to avoid any boundary effect we use periodic boundary conditions as if the same grain arrangement and the same flow field were replicated periodically in space", so I constructed scene in periodic boundary conditions and used the PeriodicFlowEngine to simulate flow. Is that right?
I want to make the fluid flow from top to bottom under a certain pressure gradient, and no fluid flows out from the side. So I set flow.bndCondIsPressure=[0,0,1,1,0,0], flow.bndCondValue=[0,0,1,0,0,0]. If my approach is wrong, then how should I set up to achieve what I want?
And how to switch solver to LDlt?

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#4
Revision history for this message
Zhicheng Gao (zhichenggao) said :
#5

Hi,
I have searched for the answers about the PeriodicFlowEngine. I concluded that since there is no boundary condition for periodic boundary conditions, the PeriodicFlowEngine does not need to set boundary conditions,that is, it does not need to set the bndCondIsPressure, bndCondValue, boundaryUseMaxMin. As I want to achieve a similar simulation like the paper " A discrete numerical model involving partial fluid-solid coupling to describe suffusion effects in soils "[1], In short, it is to achieve the fluid flow from top to bottom under a certain pressure gradient, and no fluid flows out from the side. So I just need to apply a macroscopic pressure gradient by flow.gradP=Vector3(0, i,0), where i is the macroscopic pressure gradient. Is that right?
Thank you for any help!
[1]https://www.sciencedirect.com/science/article/pii/S0266352X17303087

Revision history for this message
Launchpad Janitor (janitor) said :
#6

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