yade-users team mailing list archive

[Luc Scholtès <question659557@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi folks,

I face a very strange behavior with the flow engine. With useSolver=3, the same simulation can either work on not. More precisely, the flow can be computed or not depending on... well, I don't know...

For instance, running a simple Darcy like test (constant head through a non dynamic packing), I can either obtain the expected result (permeability of the packing) or nothing... Here is a MWE:

from yade import pack, ymport

### material
def sphereMat(): return JCFpmMat(type=1,density=3000,young=9e9,poisson=0.2,frictionAngle=radians(18))
def wallMat(): return JCFpmMat(type=0,density=3000,young=9e9,poisson=0.2,frictionAngle=radians(0))

### bodies
mn,mx=Vector3(0,0,0),Vector3(1,1,1)
O.bodies.append(aabbWalls([mn,mx],oversizeFactor=1.5,thickness=0.1,material=wallMat))
sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((0.,0.,0.),(1,1.,1.)),spheresInCell=2000,radius=1/40.,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation(material=sphereMat)

dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

### engines
flow=FlowEngine(
        isActivated=1,
        ### choose solver to use (0: Gauss Seidel, 1: Taucs, 2: Pardiso, 3: CHOLMOD)
        useSolver=3, # 3 should be used by default but does not work everytime... 
         boundaryUseMaxMin = [0,0,0,0,0,0],
        bndCondIsPressure = [1,1,0,0,0,0],
        bndCondValue = [1,0,0,0,0,0],
        permeabilityFactor=1,
        viscosity=0.001,
        #debug=1
)

intR=1.1
O.engines=[
	ForceResetter()
	,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')])
	,InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
	)
	,flow
	,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8)
	,NewtonIntegrator()
]

### simulation starts here

## to block the particles (better for permeability test)
for o in O.bodies:
 o.state.blockedDOFs+='xyzXYZ'
 
O.run(1,True)
## getBoundaryFlux get the total discharge [m3/s]
Qin = flow.getBoundaryFlux(0) 
Qout = flow.getBoundaryFlux(1)
# if Qout is the total discharge, we can compute k=Q*nu*Length/(Area*(Pout-Pin))
# if Qout is the flux, we can compute k=Q*nu*Length/(Pout-Pin) -> getFlux gives total discharge -> Qout (m3/s)!
permeability = abs(Qout)*flow.viscosity*X/(Y*Z) # !!! if Pout=1, Pin=0
permeability2 = flow.averageVelocity()*flow.viscosity*X # !!! if Pout=1, Pin=0
conductivity = permeability*1000*9.82/flow.viscosity # K=rho*g*k/nu
print "Qin=",Qin," Qout=",Qout," ARE THEY EQUAL? IF NOT-> NO FLOW!"
print "Permeability [m2]=",permeability," || Hydraulic conductivity [m/s]=",conductivity
flow.saveVtk() #  to see the result in Paraview

If I launch it 10 times, it fails 2 or 3 times to give the correct result which is:

Qin= -3.93561858335  Qout= 3.93561858335  ARE THEY EQUAL? IF NOT-> NO FLOW!
Permeability [m2]= 0.00392697615833  || Hydraulic conductivity [m/s]= 38562.9058748

With useSolver=0, I get it 10 times out of 10. 

Also, it seems that the percentage of success is higher (100%) when less particles are involved...

I have yade installed on Ubuntu xenial 16.04.2 LTS with the latest git version (but I checked and I can reproduce the same behavior with yade or yadedaily).

Any suggestion would be much appreciated.

Thanks

Luc


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