New question #697642 on Yade:
https://answers.launchpad.net/yade/+question/697642

After deleting some particles, I encountered a lot of problems in the 
PeriodicFlowEngine. Such as:
(1) GS did not converge in 20k iterations (maybe because the reference pressure 
is 0?);
(2) segmentation fault (core dumped);
(3) Periodicity is broken;
(4) 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

For errors (3) and (4),  I have solved it by changing the value of 
PeriodicFlowEngine.duplicateThreshold and the PeriodicFlowEngine.useSolver.

For errors (1) and (2), I spent a lot of time and tried a lot, but I still 
can't solve them. Please help me!

 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, and some 
particles are removed during the calculation. 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?
The errors (1) and (2) still exist, this is my simplified 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=100,
        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 
flow.dead=0
flow.defTolerance=-1
flow.meshUpdateInterval=-1
flow.useSolver=0
flow.permeabilityFactor=1
flow.viscosity=10
flow.gradP=Vector3(0,1,0)
flow.duplicateThreshold=0.3 
flow.updateTriangulation=True
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,90):
    O.bodies.erase(i)
print(len(O.bodies))
O.run(1000,True)

flow.dead=0
flow.defTolerance=-1
flow.meshUpdateInterval=-1
flow.useSolver=0
flow.permeabilityFactor=1
flow.viscosity=10
flow.gradP=Vector3(0,1,0)
flow.duplicateThreshold=0.3
flow.updateTriangulation=True

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

[1]https://www.sciencedirect.com/science/article/pii/S0266352X17303087



-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to     : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Reply via email to