Flow engine working randomly with CHOLMOD

Asked by Luc Scholtès

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

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Luc Scholtès
Solved:
Last query:
Last reply:

This question was reopened

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

Hi Luc,
I've heard the same story from different sides in the last months. It seems the linear solver has issues.
Can you show here the versions of the libraries used with "ldd /usr/lib/x86_64-linux-gnu/yadedaily/libyade.so"?
Bruno

Revision history for this message
Luc Scholtès (luc) said :
#2

Hi Bruno,

Thank you for the prompt answer.

Here is what the ldd command gives:

 linux-vdso.so.1 => (0x00007ffeb87d0000)
 libboost_python-py27.so.1.58.0 => /usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.58.0 (0x00007f54512dd000)
 libboost_thread.so.1.58.0 => /usr/lib/x86_64-linux-gnu/libboost_thread.so.1.58.0 (0x00007f54510b6000)
 libboost_filesystem.so.1.58.0 => /usr/lib/x86_64-linux-gnu/libboost_filesystem.so.1.58.0 (0x00007f5450e9e000)
 libboost_iostreams.so.1.58.0 => /usr/lib/x86_64-linux-gnu/libboost_iostreams.so.1.58.0 (0x00007f5450c85000)
 libboost_regex.so.1.58.0 => /usr/lib/x86_64-linux-gnu/libboost_regex.so.1.58.0 (0x00007f545097c000)
 libboost_serialization.so.1.58.0 => /usr/lib/x86_64-linux-gnu/libboost_serialization.so.1.58.0 (0x00007f545071b000)
 libboost_system.so.1.58.0 => /usr/lib/x86_64-linux-gnu/libboost_system.so.1.58.0 (0x00007f5450517000)
 libpthread.so.0 => /lib/x86_64-linux-gnu/libpthread.so.0 (0x00007f54502f9000)
 libpython2.7.so.1.0 => /usr/lib/x86_64-linux-gnu/libpython2.7.so.1.0 (0x00007f544fd6b000)
 libCGAL.so.11 => /usr/lib/x86_64-linux-gnu/libCGAL.so.11 (0x00007f544fb44000)
 libgmp.so.10 => /usr/lib/x86_64-linux-gnu/libgmp.so.10 (0x00007f544f8c3000)
 libcholmod.so.3.0.6 => /usr/lib/x86_64-linux-gnu/libcholmod.so.3.0.6 (0x00007f544f5ef000)
 libopenblas.so.0 => /usr/lib/libopenblas.so.0 (0x00007f544d55b000)
 librt.so.1 => /lib/x86_64-linux-gnu/librt.so.1 (0x00007f544d352000)
 libvtkIOXML-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtkIOXML-6.2.so.6.2 (0x00007f544d05e000)
 _GLViewer.so => /usr/lib/x86_64-linux-gnu/yadedaily/py/yade/qt/_GLViewer.so (0x00007f544cbbe000)
 libglut.so.3 => /usr/lib/x86_64-linux-gnu/libglut.so.3 (0x00007f544c975000)
 libvtkIOCore-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtkIOCore-6.2.so.6.2 (0x00007f544c6fe000)
 libvtkCommonDataModel-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtkCommonDataModel-6.2.so.6.2 (0x00007f544c18f000)
 libdl.so.2 => /lib/x86_64-linux-gnu/libdl.so.2 (0x00007f544bf8a000)
 libvtkCommonCore-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtkCommonCore-6.2.so.6.2 (0x00007f544ba5c000)
 libQt5Core.so.5 => /usr/lib/x86_64-linux-gnu/libQt5Core.so.5 (0x00007f544b586000)
 libGLU.so.1 => /usr/lib/x86_64-linux-gnu/libGLU.so.1 (0x00007f544b316000)
 libGL.so.1 => /usr/lib/x86_64-linux-gnu/mesa/libGL.so.1 (0x00007f544b0a2000)
 libstdc++.so.6 => /usr/lib/x86_64-linux-gnu/libstdc++.so.6 (0x00007f544ad20000)
 libm.so.6 => /lib/x86_64-linux-gnu/libm.so.6 (0x00007f544aa16000)
 libgomp.so.1 => /usr/lib/x86_64-linux-gnu/libgomp.so.1 (0x00007f544a7f4000)
 libgcc_s.so.1 => /lib/x86_64-linux-gnu/libgcc_s.so.1 (0x00007f544a5de000)
 libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f544a213000)
 libbz2.so.1.0 => /lib/x86_64-linux-gnu/libbz2.so.1.0 (0x00007f544a003000)
 libz.so.1 => /lib/x86_64-linux-gnu/libz.so.1 (0x00007f5449de9000)
 libicui18n.so.55 => /usr/lib/x86_64-linux-gnu/libicui18n.so.55 (0x00007f5449986000)
 libicuuc.so.55 => /usr/lib/x86_64-linux-gnu/libicuuc.so.55 (0x00007f54495f2000)
 /lib64/ld-linux-x86-64.so.2 (0x00005617c8ea4000)
 libutil.so.1 => /lib/x86_64-linux-gnu/libutil.so.1 (0x00007f54493ee000)
 libblas.so.3 => /usr/lib/libblas.so.3 (0x00007f544918d000)
 liblapack.so.3 => /usr/lib/liblapack.so.3 (0x00007f54489aa000)
 libamd.so.2.4.1 => /usr/lib/x86_64-linux-gnu/libamd.so.2.4.1 (0x00007f54487a0000)
 libcolamd.so.2.9.1 => /usr/lib/x86_64-linux-gnu/libcolamd.so.2.9.1 (0x00007f5448599000)
 libcamd.so.2.4.1 => /usr/lib/x86_64-linux-gnu/libcamd.so.2.4.1 (0x00007f544838f000)
 libccolamd.so.2.9.1 => /usr/lib/x86_64-linux-gnu/libccolamd.so.2.9.1 (0x00007f5448183000)
 libsuitesparseconfig.so.4.4.6 => /usr/lib/x86_64-linux-gnu/libsuitesparseconfig.so.4.4.6 (0x00007f5447f80000)
 libgfortran.so.3 => /usr/lib/x86_64-linux-gnu/libgfortran.so.3 (0x00007f5447c55000)
 libvtkIOXMLParser-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtkIOXMLParser-6.2.so.6.2 (0x00007f5447a38000)
 libvtkCommonExecutionModel-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtkCommonExecutionModel-6.2.so.6.2 (0x00007f544778a000)
 libvtkCommonMisc-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtkCommonMisc-6.2.so.6.2 (0x00007f5447573000)
 libvtkCommonSystem-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtkCommonSystem-6.2.so.6.2 (0x00007f544735f000)
 libvtksys-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtksys-6.2.so.6.2 (0x00007f544711a000)
 libQGLViewer.so.2 => /usr/lib/x86_64-linux-gnu/libQGLViewer.so.2 (0x00007f5446e74000)
 libgl2ps.so.0 => /usr/lib/libgl2ps.so.0 (0x00007f5446c61000)
 libQt5Xml.so.5 => /usr/lib/x86_64-linux-gnu/libQt5Xml.so.5 (0x00007f5446c25000)
 libQt5OpenGL.so.5 => /usr/lib/x86_64-linux-gnu/libQt5OpenGL.so.5 (0x00007f5446bca000)
 libQt5Widgets.so.5 => /usr/lib/x86_64-linux-gnu/libQt5Widgets.so.5 (0x00007f544653d000)
 libQt5Gui.so.5 => /usr/lib/x86_64-linux-gnu/libQt5Gui.so.5 (0x00007f5445ff5000)
 libX11.so.6 => /usr/lib/x86_64-linux-gnu/libX11.so.6 (0x00007f5445cba000)
 libXi.so.6 => /usr/lib/x86_64-linux-gnu/libXi.so.6 (0x00007f5445aaa000)
 libXxf86vm.so.1 => /usr/lib/x86_64-linux-gnu/libXxf86vm.so.1 (0x00007f54458a3000)
 libvtkCommonTransforms-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtkCommonTransforms-6.2.so.6.2 (0x00007f5445674000)
 libvtkCommonMath-6.2.so.6.2 => /usr/lib/x86_64-linux-gnu/libvtkCommonMath-6.2.so.6.2 (0x00007f5445452000)
 libpcre16.so.3 => /usr/lib/x86_64-linux-gnu/libpcre16.so.3 (0x00007f54451eb000)
 libglib-2.0.so.0 => /lib/x86_64-linux-gnu/libglib-2.0.so.0 (0x00007f5444eda000)
 libexpat.so.1 => /lib/x86_64-linux-gnu/libexpat.so.1 (0x00007f5444cb0000)
 libxcb-dri3.so.0 => /usr/lib/x86_64-linux-gnu/libxcb-dri3.so.0 (0x00007f5444aad000)
 libxcb-present.so.0 => /usr/lib/x86_64-linux-gnu/libxcb-present.so.0 (0x00007f54448aa000)
 libxcb-sync.so.1 => /usr/lib/x86_64-linux-gnu/libxcb-sync.so.1 (0x00007f54446a3000)
 libxshmfence.so.1 => /usr/lib/x86_64-linux-gnu/libxshmfence.so.1 (0x00007f544449f000)
 libglapi.so.0 => /usr/lib/x86_64-linux-gnu/libglapi.so.0 (0x00007f5444270000)
 libXext.so.6 => /usr/lib/x86_64-linux-gnu/libXext.so.6 (0x00007f544405e000)
 libXdamage.so.1 => /usr/lib/x86_64-linux-gnu/libXdamage.so.1 (0x00007f5443e5a000)
 libXfixes.so.3 => /usr/lib/x86_64-linux-gnu/libXfixes.so.3 (0x00007f5443c54000)
 libX11-xcb.so.1 => /usr/lib/x86_64-linux-gnu/libX11-xcb.so.1 (0x00007f5443a52000)
 libxcb-glx.so.0 => /usr/lib/x86_64-linux-gnu/libxcb-glx.so.0 (0x00007f5443838000)
 libxcb-dri2.so.0 => /usr/lib/x86_64-linux-gnu/libxcb-dri2.so.0 (0x00007f5443633000)
 libxcb.so.1 => /usr/lib/x86_64-linux-gnu/libxcb.so.1 (0x00007f5443411000)
 libdrm.so.2 => /usr/lib/x86_64-linux-gnu/libdrm.so.2 (0x00007f5443201000)
 libicudata.so.55 => /usr/lib/x86_64-linux-gnu/libicudata.so.55 (0x00007f544174a000)
 libquadmath.so.0 => /usr/lib/x86_64-linux-gnu/libquadmath.so.0 (0x00007f544150a000)
 libgobject-2.0.so.0 => /usr/lib/x86_64-linux-gnu/libgobject-2.0.so.0 (0x00007f54412b7000)
 libpng12.so.0 => /lib/x86_64-linux-gnu/libpng12.so.0 (0x00007f5441091000)
 libharfbuzz.so.0 => /usr/lib/x86_64-linux-gnu/libharfbuzz.so.0 (0x00007f5440e33000)
 libpcre.so.3 => /lib/x86_64-linux-gnu/libpcre.so.3 (0x00007f5440bc2000)
 libXau.so.6 => /usr/lib/x86_64-linux-gnu/libXau.so.6 (0x00007f54409be000)
 libXdmcp.so.6 => /usr/lib/x86_64-linux-gnu/libXdmcp.so.6 (0x00007f54407b7000)
 libffi.so.6 => /usr/lib/x86_64-linux-gnu/libffi.so.6 (0x00007f54405af000)
 libfreetype.so.6 => /usr/lib/x86_64-linux-gnu/libfreetype.so.6 (0x00007f5440304000)
 libgraphite2.so.3 => /usr/lib/x86_64-linux-gnu/libgraphite2.so.3 (0x00007f54400de000)

Revision history for this message
Robert Caulk (rcaulk) said :
#3

Hello Luc,

I ran your script 10 times with each of useSolver=0 and 3 (also 4 to check the GPU consistency). I was unable to reproduce the failure you describe, I simply get the following results 100% of the time:

Qin= -4.12495019931 Qout= 4.12493778356 ARE THEY EQUAL? IF NOT-> NO FLOW!
Permeability [m2]= 0.0041254746407 || Hydraulic conductivity [m/s]= 40512.1609716

I am running Yade from sources (git version as of Sep 28, 2017, which is the latest w.r.t basic FlowEngine and FlowBoundingSphere) on Ubuntu 16.04 LTS.

If I were you, I would clean my build directory and re make/install.

Best,

Robert

Revision history for this message
Luc Scholtès (luc) said :
#4

Hi Robert,

Thank you for your advice but I already did that... twice... I even removed yade and yadedaily and reinstalled them. Always the same problem... This and the qt issue, I am really thinking about moving back to Ubuntu 14.04. This is really annoying.

Just a stupid question: I got this message adding -DDFNFLOW=ON at the cmake step:

CMake Warning:
  Manually-specified variables were not used by the project:

    DFNFLOW

I thought it would not work so I uncommented the ifdef DFNFlow in the DFNFlow.cpp but was it necessary?

Luc

Revision history for this message
Robert Caulk (rcaulk) said :
#5

Hmmmm, when you mention that it fails 2 or 3 times to give the correct result. What result is the failed result exactly? You mention that the flow is simply not computed 20-30% of the time? So you get 0 values for everything?

>I thought it would not work so I uncommented the ifdef DFNFlow in the DFNFlow.cpp but was it necessary?

I believe it was not necessary. Since the cmake step doesn't actually parse any of the cpp files, it doesn't realize we will use DFNFLOW later. That's my guess.

>This and the qt issue,

Side question: what kind of graphics card are you using? I can't believe how well this qt problem has eluded us. I've tried a few times to really dig into it but I can't come up with anything. For some reason, I do not have the qt problem on my new computer with 16.04, but I do have the problem on my old comp with 16.04. I honestly have no idea why.

Revision history for this message
Luc Scholtès (luc) said :
#6

Typically, here is what I get:

Qin= -308.240063332 Qout= 0.0
Permeability [m2]= 0.0 || Hydraulic conductivity [m/s]= 0.0

and when I open the vtk file, I can clearly see 0 pressure values in every cells (1e-16) and high velocity values on several cells located on the boundary where the pressure head is applied (1e5).

Regarding my laptop configuration (it is quite old but still running OK, bought it in 2012):

Dell Precision M4700
Memory: 77 GIB
Processor: Intel® Core™ i7-3840QM CPU @ 2.80GHz × 8
Graphics: Gallium 0.4 on AMD CAPE VERDE (DRM 2.49.0 / 4.10.0-37-generic, LLVM 3.8.0)
OS-type: 64 bits

Another interesting fact: I always get a wrong result (or no result if you prefer) when I am launching this script while another simulation (different, no flow involved, working alright) is running at the same time... Is it a sign that my laptop is at the end of its carreer :)

Luc

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

>CMake Warning:
  Manually-specified variables were not used by the project:
    DFNFLOW

Syntax problem. DFNFLOW is a flag for gcc, not for cmake. Consequently you have to pass it as a gcc flag:
-DCMAKE_CXX_FLAGS="-DDFNFLOW"

Revision history for this message
Luc Scholtès (luc) said :
#8

Before reinstalling ubuntu 14, I would liek to ask if this problem coul dbe related to the libsuitesparse-metis-dev that needs to be install (cf. yade installation page).

When I tried to install it

sudo apt-get install libsuitesparse-metis-dev

Here is what I get:

Reading package lists... Done
Building dependency tree
Reading state information... Done
Package libsuitesparse-metis-dev is not available, but is referred to by another package.
This may mean that the package is missing, has been obsoleted, or
is only available from another source
However the following packages replace it: libsuitesparse-dev:i386 libsuitesparse-dev
E: Package 'libsuitesparse-metis-dev' has no installation candidate

Following this message, I tried to install libsuitesparse-dev:i386 and libsuitesparse-dev but, and it seems that I missed it the first time, libsuitesparse-dev:i386 cannot be installed.

Here is what I get when I sudo apt-get install libsuitesparse-dev:i386:

Reading package lists... Done
Building dependency tree
Reading state information... Done
Some packages could not be installed. This may mean that you have
requested an impossible situation or if you are using the unstable
distribution that some required packages have not yet been created
or been moved out of Incoming.
The following information may help to resolve the situation:
The following packages have unmet dependencies:
 libsuitesparse-dev:i386 : Depends: libcholmod3.0.6:i386 (= 1:4.4.6-1) but it is not going to be installed
                           Depends: libumfpack5.7.1:i386 (= 1:4.4.6-1) but it is not going to be installed
                           Depends: libspqr2.0.2:i386 (= 1:4.4.6-1) but it is not going to be installed
                           Depends: libblas-dev:i386 but it is not going to be installed or
                                    libblas.so:i386
                           Depends: liblapack-dev:i386 but it is not going to be installed or
                                    liblapack.so:i386
E: Unable to correct problems, you have held broken packages.

What should I do? Any suggestions about this? What do you guys have done on Ubuntu 16.04?

Luc

Revision history for this message
Robert Caulk (rcaulk) said :
#9

Hello Luc,

Actually, you should be able to use these packages as a replacement:

sudo apt-get -y install libopenblas-dev libsuitesparse-dev libmetis-dev

Best,

Robert

Revision history for this message
Luc Scholtès (luc) said :
#10

Alright...

I reinstalled ubuntu (16.04.3) and reinstalled yade (git version):

It seems to work correctly now.

Why?

I don't know... Probably a problem in the dependencies somewhere.

But, I now have another, not critical but still annoying, problem:

FlowEngine requires the walls to be appended in the simulation before the spheres, right?

Up to now, I could simply import the sample, get its dimensions (with utils.aabbExtrema()) and then reset the scene with O.reset() so as to append the walls with the right dimensions and then append the packing again. It was a relatively simple way to automatize the simulations and I could run the same simulation with different packings without bothering about their dimensions (which was pretty convenient).

The problem is that, now, when I use the same process, the FlowEngine does not work (again the same behavior with no flow computed at all)... I suppose this is due to the use of the O.reset() in the script but why would that be like this?

If you feel that I should open a bug report or another question for this issue, I'll do it and then I'll close this one.

Luc

Revision history for this message
Robert Caulk (rcaulk) said :
#11

We might a swell sort it out here since FlowEngine's behavior doesn't change.

Hmmm, I use the same method you describe for getting packing dimensions:

sp = O.bodies.append(ymport.textExt('4cmEdge_1mm.spheres', 'x_y_z_r',color=(0.1,0.1,0.9)))

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

mn,mx=Vector3(xinf, yinf,zinf),Vector3(xsup, ysup, zsup)
print "mn ", mn, "mx ", mx

O.reset()

followed by the actual simulation (including FlowEngine). And it works for me.

Can you post the MWE here so I can test it with my setup?

Best,

Robert

Revision history for this message
Luc Scholtès (luc) said :
#12

Here you go. It fails every time. Of course, you'll have to create your own packing (cf PACK parameter in the first lines).

Also, I tried to simulate an injection and it does not look good at all compared to what I used to obtain. Any chance you would have a MWE that I could try on my side?

Thanks for your help.

Luc

from yade import pack, ymport

#### packing you want to test
PACK='111_5k'
intR=1.262

#### if you want to adjust the macroscopic permeability
pFactor=8.e-18
# for k=1e-20 m2: 8.e-18 for 111_5k

#### material definition
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))

#### preprocessing to get dimensions
O.bodies.append(ymport.text(PACK+'.spheres',scale=1.,shift=Vector3(0,0,0),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

R=0
Rmax=0
Rmin=1000
numSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   numSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
   if o.shape.radius<Rmin:
     Rmin=o.shape.radius
Rmean=R/numSpheres

#print "nb spheres=",numSpheres, " | mean Diameter=", 2*Rmean, ",X,Y,Z,", X,Y,Z

#### reset the scene now
O.reset()

#### simulation is built up now
mn,mx=Vector3(xinf,yinf,zinf),Vector3(xsup,ysup,zsup)
walls=aabbWalls([mn,mx],oversizeFactor=1.5,thickness=X/10.,material=wallMat)
wallIds=O.bodies.append(walls)

O.bodies.append(ymport.text(PACK+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

flow=DFNFlowEngine(
        isActivated=1,
        useSolver=3,
        boundaryUseMaxMin = [0,0,0,0,0,0],
        permeabilityFactor=pFactor,
        viscosity=0.001
)

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(smoothJoint=True)]
 )
 ,flow
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=0.1*utils.PWaveTimeStep())
 ,NewtonIntegrator(damping=0.4)
]

## if you want to block the particles (better for permeability test)
for o in O.bodies:
 o.state.blockedDOFs+='xyzXYZ'

## flow along X direction
print 'Flow along X!'
flow.bndCondIsPressure = [1,1,0,0,0,0]
flow.bndCondValue = [1,0,0,0,0,0]
O.step()
## 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() # if you want to see the result in Paraview

Revision history for this message
Robert Caulk (rcaulk) said :
#13

The script runs with success and produces:

Qin= -5.06465867124e-20 Qout= 5.06465867124e-20 ARE THEY EQUAL? IF NOT-> NO FLOW!
Permeability [m2]= 6.33158067289e-22 || Hydraulic conductivity [m/s]= 6.21761222077e-15

But I don't have the most recent version of DFNFlow. Does this script (with the automatic dimension detection) work with the regular FlowEngine for you?

>Also, I tried to simulate an injection and it does not look good at all compared to what I used to obtain. Any chance you would have a MWE that I could try on my side?

What happens? What were you obtaining before? You can apply flux with:

flow.imposeFlux(flow.getCellCenter(flow.getCell(0, 0, 0)), -flowRate)

if this doesn't work I can try to put an MWE together for you

Revision history for this message
Luc Scholtès (luc) said :
#14

Yes, same problem with either FlowEngine or DFNFlowEngine (which is kind of reassuring and "logical"): the flow is not computed.

I guess the initial problem persists... The linear solver just does not work all the time (and it seems to be affected with the O.reset())... It would explain why my injection test is also not working... That's totally crazy and just incomprehensible...

Should I burn my laptop at this stage...

ANW, I'll try to install ubuntu 14.04 to see if the problem persists. If it does, well... I don't know.

Thanks for your efforts ANW.

Luc

Revision history for this message
Luc Scholtès (luc) said :
#15

FYI, the O.reset() script works with yade installed from sudo-apt get installed yade (yade 1.20.0).

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

>FlowEngine requires the walls to be appended in the simulation before the spheres, right?

Not at all. But you have to specify their ids if they are not 0,...,5. There's a FlowEngine's attribute for this (can't remember the name now).
As always I would suggest to not impose importing a pack for scripts in general and for MWE specifically. There is not rational in doing that (*), and it can be an additional source of problems (ref. your problem with appending walls, unrelated to the solver issue).
I have really no "pack" ready on my hard drive, not even an idea of how to get one (I'm serious), so I am unable to test your script.
B

(*) No it is not saving cpu time (in case that would be the next argument) since a self-generating script can re-use it's own pack multiple times.

Revision history for this message
Luc Scholtès (luc) said :
#17

Hi guys,

I reinstalled everything from scratch (ubuntu 16.04 LTS and yade from git). The problem persists and I cannot figure out what is wrong.

Several things:

- flow without surrounding walls: boundaryUseMaxMin=[1,1,1,1,1,1] should let us run a simulation without walls around the packing (at least, this is what I understood from the doc). I don't know if you tried but it does not work. It is certainly my fault but if you could check on your side, I would be curious to hear about it (I have a script for that if needed).

- flow with walls appended after the packing: I tried to use the wallIds option in Flow Engine but, again, it seems that it does not work. Here again, it is certainly my fault but if you could check on your side, I would be curious to hear about it (I also have a script for that if needed).

- Regarding the script with the predefined packing, as mentioned in a previous message, it happened to trigger the problem at every try as opposed to the case where the packing is built up inside the simulation. That was the only reason why I mentioned it. I understand that a MWE should not be built up with such option but Robert seems to use the same procedure so I wanted to make sure I did not introduce any misleading line in my script (and I did not since Robert could run it successfully on his side). Besides this, I must say that the average coordination number is a parameter of my simulations with bonded particles. I need to calibrate the interaction range so that the number of bonds per particle corresponds to the value identified during the calibration of my models. Like this, I ensure that the simulated behavior is not affected by, among other things, the number of particles. Thus, I need to run my simulations on packings for which I predetermine the interaction range which gives the right coordination number. This is why I need to import packings in my simulations. Up to now, it was working fine.

- Coming back to the flow solver, I must say that I am totally clueless at the moment... Bruno, you mentioned in a previous message that you've heard stories about similar issues. How did these stories end? Could you find a solution (apart from reading incantations while burning the computer on any sort of dedicated altar)? I was not able to find anything on the mailing list...

Thank you for keeping up on this. I was supposed to work on HM coupled problems 100% of my time during the coming months...

Luc

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

> boundaryUseMaxMin=[1,1,1,1,1,1]

It means that instead of walls positions it will use static positions defined by max/min (the static positions can change with time but in terms of boundary condition it is always as if the wall was immobile). It does not mean "no wall", at some point some fluid forces on the walls are calculated and they need to be assigned to a body (that's at leas the current code). So, you still need the corresponding bodies (not necessarily of the Box class).

For wallIds, yes please show a script.

>How did these stories end?

They did not end unfortunately, but for the moment I am not able to tell if it is a problem on yade side or on user's side.
Too vague statements/crash reports to know if it is wrong timestep, wrong allocation of ressources on a cluster, or really unexpected behavior by the fluid solver.
I understand the reason for the reloading exemple, ok.
The need to predetermine interaction range does not make reloading a packing mandatory (imagine you insert the python script doing this calibration inside the script doing the actual simulation?), but that's another story. We can discuss it separately if you like.
Cheers
Bruno

Revision history for this message
Robert Caulk (rcaulk) said :
#19

If you post the two scripts that you have created for testing boundaryUseMaxMin() and wallIds(), I will gladly test them on my setup. I have never used either of those functions so I do not have scripts that readily employ them.

Also, don't forget to turn flow.debug on so that you can get some more detailed information!

Revision history for this message
Luc Scholtès (luc) said :
#20

- OK for the boundaryUseMaxMin=[1,1,1,1,1,1] . I misunderstood the doc.

- For the calibration of the coordination number, yes, we could do that with a dedicated python script but I must say that I am not a big fan of the randomDensePack() function, even though my way to generate dense packing is probably not better.

- Regarding the problem with the solver, I guess I'll need to figure it out by myself... Not an easy task...

- Regarding the walls appended after the packing, no error message shows up really (even with debug=1). It is just that the computation is wrong (Qin and Qout are not equal) and that the VTK file cannot be open with paraview (error reading ascii cell data). Could this be related to the tesselation?

Here is the script:

from yade import pack, ymport

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

#### bodies
X=Y=Z=1

sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((0.,0.,0.),(X,Y,Z)),spheresInCell=1000,radius=X/20.,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

Rmean=0
nbSph=0
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
Rmean=Rmean/nbSph

print "nb spheres=",nbSph, " | mean Diameter=", 2*Rmean, ",X,Y,Z=", X,Y,Z

mn,mx=Vector3(xinf,yinf,zinf),Vector3(xsup,ysup,zsup)
#mn,mx=Vector3(0,0,0),Vector3(X,Y,Z)
walls=utils.aabbWalls(extrema=[mn,mx],oversizeFactor=1.5,thickness=min(X,Y,Z)/100.,material=sphereMat)
ids=O.bodies.append(walls)
print 'walls ids=',ids

#### engines
flow=FlowEngine(
        isActivated=1,
        useSolver=3,
        wallIds=ids,
        boundaryUseMaxMin=[0,0,0,0,0,0],
        permeabilityFactor=1,
        viscosity=0.001,
        debug=1
)

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()])
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
 )
 ,flow
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=0.1*utils.PWaveTimeStep())
 ,NewtonIntegrator(damping=0.)
]

#### simulation
## to block the particles (better for permeability test)
for o in O.bodies:
 o.state.blockedDOFs+='xyzXYZ'

## flow along X direction
print 'Flow along X!'
flow.bndCondIsPressure = [1,1,0,0,0,0]
flow.bndCondValue = [1,0,0,0,0,0]
O.step()
## 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?"
print "Permeability [m2]=",permeability," || Hydraulic conductivity [m/s]=",conductivity
flow.saveVtk() # if you want to see the result in Paraview

If you want to compare with the script where the walls are appended before the packing (which gives a reasonable result with Qin=Qout and a vtk file that can be open in Paraview), here it is:

from yade import pack, ymport

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

#### bodies
X=Y=Z=1

mn,mx=Vector3(0,0,0),Vector3(X,Y,Z)
walls=utils.aabbWalls(extrema=[mn,mx],oversizeFactor=1.5,thickness=min(X,Y,Z)/100.,material=sphereMat)
ids=O.bodies.append(walls)
print 'wall ids=',ids

sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((0.,0.,0.),(X,Y,Z)),spheresInCell=1000,radius=X/20.,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

Rmean=0
nbSph=0
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
Rmean=Rmean/nbSph

print "nb spheres=",nbSph, " | mean Diameter=", 2*Rmean, ",X,Y,Z=", X,Y,Z

#### engines
flow=FlowEngine(
        isActivated=1,
        useSolver=3,
        wallIds=ids,
        boundaryUseMaxMin=[0,0,0,0,0,0],
        permeabilityFactor=1,
        viscosity=0.001,
        debug=1
)

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()])
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
 )
 ,flow
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=0.1*utils.PWaveTimeStep())
 ,NewtonIntegrator(damping=0.)
]

#### simulation
## to block the particles (better for permeability test)
for o in O.bodies:
 o.state.blockedDOFs+='xyzXYZ'

## flow along X direction
print 'Flow along X!'
flow.bndCondIsPressure=[1,1,0,0,0,0]
flow.bndCondValue=[1,0,0,0,0,0]
O.step()
## 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?"
print "Permeability [m2]=",permeability," || Hydraulic conductivity [m/s]=",conductivity
flow.saveVtk() # if you want to see the result in Paraview

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

> for the boundaryUseMaxMin=[1,1,1,1,1,1] . I misunderstood the doc.

Or the doc is not clear... feel free to improve when you spot missing pieces of information.

> For the calibration of the coordination number, yes, we could do that with a dedicated python script but I must say that I am not a big fan of the randomDensePack() function

I never thought about using randomDensePack.

> Regarding the problem with the solver, I guess I'll need to figure it out by myself... Not an easy task...

I can't reproduce your problem, sorry.

> Regarding the walls appended after the packing, no error message shows up really (even with debug=1). It is just that the computation is wrong (Qin and Qout are not equal)

Your script reads:
Qin = flow.getBoundaryFlux(0)
Qout = flow.getBoundaryFlux(1)
and you get:
Qin= -1.62684535981e-16 Qout= 1.17717250878e-16

Yade [1]: flow.getBoundaryFlux?
Docstring:
getBoundaryFlux( (FlowEngineT)arg1, (int)boundary) -> float :
    Get total flux through boundary defined by its body id.

Are bodies 0 and 1 the walls when you append the walls after the sĥeres?

> and that the VTK file cannot be open with paraview (error reading ascii cell data). Could this be related to the tesselation?

It could be that someone assumed wall ids 0,...,5 for some vtk export. You could check that in the code. It should not interfere in the actual resolution (but it is still annoying).

Bruno

Revision history for this message
lingran.zhang (lingran-zhang) said :
#22

Dear all,

I confirm that I am facing exactly the same problem as Luc Scholtès, that is, Flow engine works randomly! Not even to mention that strange results were obtained during the injection test. My yade was installed on Ubuntu 16.04 LTS (see computer details below) with the latest git version (trunck-master.zip file downloaded on 24/10/2017, uncomment the ifdef DFNFlow).

Particularly,

1. For script No.1 provided by Luc, with useSolver=3, the simulations give Qin>>Qout=0 for the first run and give Qin=Qout!=0 for the other runs. However, with useSolver=0, the simulations always give Q_in=Q_out!=0.
Same things occur for script No.4.

2. For script No.2, runing in a whole script gives Qin>>Qout=0, while running line by line of the same script gives Qin=Qout!=0. In fact, depending on the script, this could happen the other way around.

3. For script No.3, the simuluations give different values of Qin and Qout for each run, but always give Qin!=Qout!=0

Maybe there is something wrong in my computer configuration or in the code?

Here is my desktop detail:
Ubuntu 16.04 LTS
Memory : 62,8 GiB
processor : Intel® Xeon(R) CPU E5-2690 v3 @ 2.60GHz x 48
Graphics : Quadro FX570/PCIe/SSE2
OS type: 64-bit
Disk: 3,9 TB

Thank you!
Regards,
Lingran

Revision history for this message
Robert Caulk (rcaulk) said :
#23

- Regarding the problem with the solver, I guess I'll need to figure it out by myself... Not an easy task...

I will gladly help you down that path if there really is a problem with the solver, but the good news is I don't think either of us need to go down it. I think Bruno solved this...

I've run both of the most recent scripts posted by Luc. The first one (walls appended before spheres) runs just fine. And according to Luc's comment, the first script works just fine for him too(?) So the only issue is the second script.

The second script (walls appended after spheres) simply had the problem of using improper wall Ids for the getBoundaryFlux() functions (as Bruno suggested). So the fix for that second script would be:

Qin = flow.getBoundaryFlux(ids[0])
Qout = flow.getBoundaryFlux(ids[1])

So just to be clear, both scripts run successfully on my setup (16.04 LTS and trunk compiled Yade).

Revision history for this message
Luc Scholtès (luc) said :
#24

Well, you are right, I made a mistake in my script when asking for the Qin and Qout. It gives the correct result if, of course, you ask for the right ids (...).

Nonethless, this is not what caused the problem initially and I am still facing this problem with the solver. Unfortunately, you or Bruno cannot obtain the same error with my scripts and that's a big mystery for me at the moment... Since another person seems to face the same problem (so this is not just a curse), we need to figure out if this is really a problem with the solver or the scripts... I must say that I am confident in saying that this is related to the solver in association with the computer set up since I don't have any problem when I run my test with useSolver=0.

I put here again the script that gives the following behaviors when useSolver=3:

1) when the packing is not imported from a text file:
- works every time with yade installed with sudo apt-get install yade ( VTK file corrupted)
- works randomly with yade installed from git sources if no other simulation is running at the same time
- never works with yade installed from git sources if another simulation is running at the same time

2) when the packing is imported from a text file:
- works every time with yade installed with sudo apt-get install yade (VTK file corrupted)
- works randomly with yade installed from git sources if no other simulation is running at the same time
- never works with yade installed from git sources if another simulation is running at the same time

The script:

from yade import pack, ymport

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

#### bodies

X=Y=Z=1
sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((0.,0.,0.),(X,Y,Z)),spheresInCell=1000,radius=X/40.,returnSpherePack=True)
sp.toSimulation(material=sphereMat)

## if you uncomment, it imports a pre-existing packing
#PACKING='111_10k'
#O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),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

Rmean=0
nbSph=0
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
Rmean=Rmean/nbSph

print "nb spheres=",nbSph, " | mean Diameter=", 2*Rmean, ",X,Y,Z=", X,Y,Z

mn,mx=Vector3(xinf,yinf,zinf),Vector3(xsup,ysup,zsup)
walls=utils.aabbWalls(extrema=[mn,mx],oversizeFactor=1.5,thickness=min(X,Y,Z)/100.,material=sphereMat)
ids=O.bodies.append(walls)
print 'walls ids=',ids

#### engines
flow=FlowEngine(
        isActivated=1
        ,useSolver=0
        ,wallIds=ids
        ,boundaryUseMaxMin=[0,0,0,0,0,0]
        ,permeabilityFactor=1
        ,viscosity=0.001
        #,debug=1
)

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()])
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
 )
 ,flow
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=0.1*utils.PWaveTimeStep())
 ,NewtonIntegrator(damping=0.)
]

#### simulation
## to block the particles (better for permeability test)
for o in O.bodies:
 o.state.blockedDOFs+='xyzXYZ'

## flow along X direction
print 'Flow along X!'
flow.bndCondIsPressure = [1,1,0,0,0,0]
flow.bndCondValue = [1,0,0,0,0,0]
O.step()
## getBoundaryFlux get the total discharge [m3/s]
Qin = flow.getBoundaryFlux(ids[0])
Qout = flow.getBoundaryFlux(ids[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?"
print "Permeability [m2]=",permeability," || Hydraulic conductivity [m/s]=",conductivity
flow.saveVtk() # if you want to see the result in Paraview

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

I'm confident in the fact that something random occurs. This is consistently reported by you (Luc) and Lingran (thanks for detailed report).
This is most likely unrelated to Yade, which is at the same time good news: someone else will fix it hopefully; and bad news: we can't do much except testing various compile/config by trial and error until openblas+cholmod+eigen will work deterministicaly.
Alternatively, a solution could be to switch to a different solver, i.e. useSolver=1 or reimplement another (or many other) linear solvers (useSolver=2,3,4,...).

If I find time I will prepare "hello world" examples for the cholmod-through-eigen method (the current one) and send them to you Luc/Lingran for testing.

I would also ask you Luc/Lingran to test yade compiled with disabled openmp, and possibly to compile openblas on your side with also disabled openmp.

Bruno

Revision history for this message
Lingran Zhang (lingran) said :
#26

Hi, Bruno,

Just let you know, I've installed openblas from git (compiled with 'make FC=gfortran USE_OPENMP=0') as well as installed yade from git (with cmake -DENABLE_OPENMP=0).

However, the situation has not been improved, something random still occurs.

Thank you.
Lingran

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

@Lingran,
Thank you very much for reporting. Further, you could compile with a replacement of openblas (atlas, or legacy blas) since I suspect openblas could be the problem. This may need to trick cmake a little though, and I can't explain now in a few lines how to do that.
Bruno

Revision history for this message
Luc Scholtès (luc) said :
#28

Hi,

Just an additional comment:

Running my simulations with debug=1, I obtain the exact same output when it works and when it does not work.

Even at the very last line of the output (total force), I obtain the same values for each case (which is always !=0).

I haven't gone through the entire code yet but is this "total force" computed from fluid pressure? If yes, would that mean that the pressure is computed correctly and that the problem could be in the "writing" of the pressure in the cells?

Luc

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

The total force is entirely defined by the fluid boundary conditions (not from an algorithmic point of view, but mathematically it is the case).
For instance the total force on a permeameter column is ΔP*A whatever the details of pore pressure inside the column. That's why it looks correct.

debug=1 is helpless because the problem is an internal bug in a third-party library (we'll have to know which one).

Bruno

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

Hi Luc, Lingran,
Did you check if some older versions (that "were" working) would work for you?
I'd like to be sure that it is not a problem introduced recently in the code.
You can go back to every previous version with git.
Bruno

Revision history for this message
Lingran Zhang (lingran) said :
#31

Hi Bruno,

I've just installed yade from git (Branch:2017.01, uncomment the '#define DFNFLOW') and tested the 5 scripts provided above by Luc.
The random thing in Flow Engine seems disappeared:
For all the 5 scripts:
 useSolver=3 or 0, the simulations always give Qin=Qout !=0
 useSolver=1 or 2, the simulations always give Qin>>Qout =0
Besides, running the simulation in full script or line by line gives the same results.

Cheers!
Lingran

Revision history for this message
Luc Scholtès (luc) said :
#32

Hi,

Just to confirm that I don't have the random behavior with yade 2017.10.a as mentioned by Lingran. Why? I have no idea but I'll stick to this version for now on.

However, this version does not allow to run simulations with DFNFlowEngine as it produces a segfault as soon as a crack or fracture is present in the medium. I managed to avoid the segfault by modifying the trickPermeability() function which was at the origin of the segfault. However, now, when a cell is identified as cracked (and so its permeability is modified/increased with respect to the non cracked cells), I get this warning:

CHOLMOD warning: matrix not positive definite
something went wrong in Cholesky factorization, use LDLt as fallback this time1

and the flow computation is wrong.

Any idea where I should look at for this particular issue?

Luc

Revision history for this message
Robert Caulk (rcaulk) said :
#33

Hey Luc,

>However, this version does not allow to run simulations with DFNFlowEngine as it produces a segfault

Yes, the first DFNFlow code I encountered in Yade was 2017.10a. And it segfaulted by default...this is the bug that I reported and fixed (?) [1]. I would be curious to someday see Papachristos' DFNFlow code used to publish [2]. I digress.

>CHOLMOD warning: matrix not positive definite
>something went wrong in Cholesky factorization, use LDLt as fallback this time1

Have you tried using a smaller timestep?

Best,

Robert

[1]https://bugs.launchpad.net/yade/+bug/1666339
[2]Papachristos, E., Scholtès, L., Donzé, F. V., & Chareyre, B. (2017). Intensity and volumetric characterizations of hydraulically driven fractures by hydro-mechanical simulations. International Journal of Rock Mechanics and Mining Sciences, 93(January), 163–178. http://doi.org/10.1016/j.ijrmms.2017.01.011

Revision history for this message
Luc Scholtès (luc) said :
#34

Robert, you made my day!

There was another problem linked to the declaration of the trickPermeability() function but, your bug report was what I was looking for. It seems that I can reproduce previous simulations without any randomness in the computations. I'll need to confirm that on other configurations (e.g. injection) but, ANW, thank you very much!

Regarding the code that Timos used, I think it should be very similar to the one I have now. But only Timos can tell.

The question now is: what is the difference between this 2017.01 version and the latest one that produces the random behavior?

I don't close the thread as the initial problem (randomness with latest version) is still present.

Luc

Revision history for this message
Robert Caulk (rcaulk) said :
#35

Glad to help!

Correct me if I am wrong, but the randomness is caused by DFNFlow, you observe it with the basic FlowEngine as well in your original post. So I do not think it is DFNFlow in particular that is causing the randomness. Instead, it seems to be caused by third party libraries employed by the shared solver.

Either way, we can compare the differences between current and 2017.01.a FlowEngines using Github compare [1]. Good luck!

[1]https://github.com/yade/trunk/compare/2017.01...master

Revision history for this message
Robert Caulk (rcaulk) said :
#36

*randomness is [b]not[/b] caused

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

Hi,
As far as I understand DFN is not the question, I agree with Robert.
Third party libraries may also not be to blame in the end, if a previous revision is ok (I don't expect changes in terms of which solver is used, maybe in some solver settings?).
B

Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#38
Revision history for this message
Luc Scholtès (luc) said :
#39

Yes, DFNFlow is not the problem. The randomness appears with the "whole FlowEngine". Sorry if my message made you think so.

Since it seems that my version built from 2017.01 does not show this randomness, there must be something wrong in the source code rather than in the third party libraries as pointed out by Bruno.

I'll try to identify the origin of the problem but that's not going to be an easy task.

Shall we open a bug report instead of going on with this thread?

Luc

Revision history for this message
Robert Caulk (rcaulk) said :
#40

Hey Bruno,

I am responsible for the second commit you posted [1]. I don't think it is to blame, since Luc's original script here does not use multithread=True, and therefore does not factorize matrices in the background.

It is good that we know the 2017.10.a version doesn't have randomness on anyone's computers. But I'm still very confused by why I (and Bruno?) have trunk installed on my machine and I cannot reproduce this randomness...what does this mean??

Robert

[1]https://github.com/yade/trunk/commit/f6970362d9e6e866e8adbc8cbea18e54f677f785

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

@Robert
I tried to identify the commits which touched relevant code, but I could not spot anything wrong precisely.
Bruno

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

@Luc
Yes, let's have a bug.
B

Revision history for this message
Janek Kozicki (cosurgi) said :
#43

Bruno Chareyre said: (by the date of Wed, 08 Nov 2017 17:57:59 -0000)

> Question #659557 on Yade changed:
> https://answers.launchpad.net/yade/+question/659557
>
> Bruno Chareyre posted a new comment:
> @Robert
> I tried to identify the commits which touched relevant code, but I could not spot anything wrong precisely.

you can always try `git blame ` :)

--
Janek Kozicki http://janek.kozicki.pl/ |

Revision history for this message
Robert Caulk (rcaulk) said :
#44

Looks like Bruno pointed us in the right direction at the end but I disagreed :-/.

Luc, I just reread this thread and noticed you were installing i386 libraries? Is your computer 32 bit? I think Bruno and I could not reproduce the bug because we are on AMD64 AFAIK. Just a guess.

Either way, this was clearly a bug of my doing and I sincerely apologize for all the time you wasted trying to chase it down.

Revision history for this message
Luc Scholtès (luc) said :
#45

No worries. I cannot really blame you for trying to improve the code. That's part of the game I would say.

Regarding my computer, it is a 64 bit... Proc is intel core i7. Maybe that's the origin of the i386 lib...

ANW, I think we can close this thread.

Thank you to everyone.

Luc