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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp