Yade

Pore volume - FlowEngine

Asked by Othman Sh on 2020-10-06

Hi All,

I am trying to obtain the volume of individual pores in my sphere packing. I found this function in Yade documentation:

volume( [ (int)id=0 ] ) → float
Returns the volume of Voronoi’s cell of a sphere

I have few questions here:

1- Is this the correct approach to get the pore volume for each individual cell ID?
2- Does this function gives the volume of the whole tetrahedron element (which includes part of the spheres surrounding the pore as in Fig.3a in [1]) or just the pore (i.e. fluid phase)?
3- In my code below, this function is not working and I get this error: Segmentation fault (core dumped)
I have tried using this function in the terminal and I found that it sometimes give a value for some cell IDs and most of the time it will crush and this error occurs. My code is copied below. The packing txt file link is here [2]

I will appreciate your help in answering my questions.

Thank you,
Othman

[1] https://onlinelibrary.wiley.com/doi/full/10.1002/nag.2198
[2] https://drive.google.com/file/d/1JqFe3W1cVFgKEa9lqU-MsCBxY7RDgY89/view?usp=sharing
---------------------------------------------

from yade import plot,pack, export, ymport
import numpy as np

P=4347 #Pa
visc=1e-3 #Pa.sec
density=1000 #kg/m3
g=9.81 #m/s2
radiuscyl=.05
########## create walls ##########
allx,ally,allz,r=np.genfromtxt('S30.txt', unpack=True) #access the packed cyl file (txt)
mnx=min(allx)*0.999
mny=min(ally)*0.999
mnz=min(allz)*0.999
mxx=max(allx)*1.001
mxy=max(ally)*1.001
mxz=max(allz)*1.001

mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz)
walls=aabbWalls([mn,mx],thickness=0)
wallIds=O.bodies.append(walls)

########## spheres ##########
spheres=O.bodies.append(ymport.text('/home/auser/DEM_codes/Clogging/S30.txt')) #please change this when you download the packing file
yade.qt.View()
Height=max(utils.aabbDim())
dP=P/Height #Pa/m

for i in spheres:
 body= O.bodies[i]
 body.state.blockedDOFs='xyzXYZ'
print ('porosity = ', utils.porosity())

Height=max(utils.aabbDim())

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 FlowEngine(dead=1,label="flow"),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 NewtonIntegrator(damping=0.2)
]

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.dead=0
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=visc
flow.bndCondIsPressure=[1,1,1,1,1,1]
flow.bndCondValue=[0,0,0,0,0,P]
flow.boundaryUseMaxMin=[1,1,1,1,1,1]
O.dt=1e-6
O.dynDt=False
flow.emulateAction()

#########################################################################################################################
################################################# Get pores volume#####################################
##########################################################################################################################
cellsHit = []
numPoints = 20
xs = np.linspace(0.95*mnx,1.05*mxx,numPoints)
ys = np.linspace(0.95*mny,1.05*mxy,numPoints)
zs = np.linspace(0.95*mnz,1.05*mxz,numPoints)

v = np.array([0,0,0])

for x,y,z in itertools.product(xs, ys, zs):
 cellId = flow.getCell(x,y,z) #getCell return cell id for a point position xyz
 if cellId in cellsHit: continue #this prevents counting a cell multiple times
 cellsHit.append(cellId)
 VV=flow.volume(cellId)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Othman Sh
Solved:
2020-10-22
Last query:
2020-10-22
Last reply:
2020-10-22

This question was reopened

Robert Caulk (rcaulk) said : #1

Hello,

> I found this function in Yade documentation:

volume( [ (int)id=0 ] ) → float
Returns the volume of Voronoi’s cell of a sphere

Please link to the documentation!

>My code is copied below. The packing txt file link is here [2]

Please don't link to external sources! Please review our posting guidelines if you wish to receive efficient assitance with your Yade related problems -> -> -> [1] <- <- <-

>1- Is this the correct approach to get the pore volume for each individual cell ID?

No, use getCellVolume [2]

>2- Does this function gives the volume of the whole tetrahedron element (which includes part of the spheres surrounding the pore as in Fig.3a in [1]) or just the pore (i.e. fluid phase)?

getCellVolume [2] is the full volume of the tetrahedra. getInvVoidVolume [3] is the inverse of the void volume.

Thanks again for following our forum guidelines ([1]) for your follow up posts/questions :-)

Cheers,

Robert

[1]https://www.yade-dem.org/wiki/Howtoask
[2]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.getCellVolume
[3]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.getCellInvVoidVolume

Othman Sh (othman-sh) said : #2

Hi Robert,

I rearranged my code to have the packing made during the simulation without the need for external link. Sorry about not following the forum guideline in my previous post.

I was able to get the cell volumes as in the code below, but I actually want to get the volume of individual pores/voids. So I believe the voids volume = 1/getInvVoidVolume, right?

Sorry if this looks like a dumb question, but I couldn't figure out how to use the getInvVoidVolume tool. In [1], the documentation reads getCellInvVoidVolume((FlowEngineT)arg1, (int)arg2)

What is (int) arg2?

Also in the code blow, I tried using flow.getInvVoidVolume((x,y,z)) but I got this error (I am using Yadedaily):

'FlowEngine' object has no attribute 'getInvVoidVolume'

Thanks for helping with this issue.

Othman

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.getCellInvVoidVolume

-------------------------------------------------------------------------
from yade import plot,pack, export, ymport
import numpy as np

young=1e6
compFricDegree = 3
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,0.05,.2)
sp.toSimulation(material='spheres')

yade.qt.View()

print ('porosity = ', utils.porosity())
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 FlowEngine(dead=1,label="flow"),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 NewtonIntegrator(damping=0.2)
]

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.dead=0
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=1e-3
flow.bndCondIsPressure=[1,1,1,1,1,1]
flow.bndCondValue=[0,0,0,0,0,10]
flow.boundaryUseMaxMin=[1,1,1,1,1,1]
O.dt=1e-6
O.dynDt=False
flow.emulateAction()

################## Get pores volume ##################
box=aabbExtrema()
lower_corner=box[0]
upper_corner=box[1]
mnx=lower_corner[0]
mxx=upper_corner[0]
mny=lower_corner[1]
mxy=upper_corner[1]
mnz=lower_corner[2]
mxz=upper_corner[2]
cellsHit = []
numPoints = 20
xs = np.linspace(0.95*mnx,1.05*mxx,numPoints)
ys = np.linspace(0.95*mny,1.05*mxy,numPoints)
zs = np.linspace(0.95*mnz,1.05*mxz,numPoints)
v = np.array([0,0,0])

for x,y,z in itertools.product(xs, ys, zs):
 cellId = flow.getCell(x,y,z)
 if cellId in cellsHit: continue
 vv=flow.getCellVolume((x,y,z))
 print(vv)

Robert Caulk (rcaulk) said : #3

>I was able to get the cell volumes as in the code below, but I actually want to get the volume of individual pores/voids. So I believe the voids volume = 1/getInvVoidVolume, right?

Yup. Except you should use getCellInvVoidVolume since getInvVoidVolume doesn't exist.

>What is (int) arg2?

The id of the cell of interest.

>'FlowEngine' object has no attribute 'getInvVoidVolume'

FlowEngine does not have that attribute. It does, however, have the attribute the getCellInvVoidvolume that you reference by link...

Robert Caulk (rcaulk) said : #4

This function, getCellInvVoidVolume() has a bug, I do not expect it to function correctly until this merge request [1] is pushed to master (likely 1 or 2 days), then it wil take another 1 or 2 days to make it to yadedaily.

[1]https://gitlab.com/yade-dev/trunk/-/merge_requests/528

Othman Sh (othman-sh) said : #5

My bad, looks like I just copied getInvVoidVolume from post #1 and I should have copied it from the reference link. I tried the correct function getCellInvVoidVolume and it gave me all zeros. I will wait for the bug to be fixed.

Thanks for your help. This solved my problem

Othman Sh (othman-sh) said : #6

Hello,

It has been a couple of weeks since the merge request mentioned in reply #4 but the getCellInvVoidVolume function is not functioning and it always return 0 values. I tried it using the code below. I appreciate any help in fixing this issue.

Thank you,
Othman

------------------------------------------------------
from yade import plot,pack, export, ymport
import numpy as np

young=1e6
compFricDegree = 3
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,0.05,.2)
sp.toSimulation(material='spheres')

yade.qt.View()

print ('porosity = ', utils.porosity())
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 FlowEngine(dead=1,label="flow"),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 NewtonIntegrator(damping=0.2)
]

flow.dead=0
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=1e-3
flow.bndCondIsPressure=[1,1,1,1,1,1]
flow.bndCondValue=[0,0,0,0,0,10]
flow.boundaryUseMaxMin=[1,1,1,1,1,1]
O.dt=1e-6
O.dynDt=False
flow.emulateAction()

################## Get pores volume ##################
box=aabbExtrema()
lower_corner=box[0]
upper_corner=box[1]
mnx=lower_corner[0]
mxx=upper_corner[0]
mny=lower_corner[1]
mxy=upper_corner[1]
mnz=lower_corner[2]
mxz=upper_corner[2]
cellsHit = []
numPoints = 20
xs = np.linspace(0.95*mnx,1.05*mxx,numPoints)
ys = np.linspace(0.95*mny,1.05*mxy,numPoints)
zs = np.linspace(0.95*mnz,1.05*mxz,numPoints)
v = np.array([0,0,0])

for x,y,z in itertools.product(xs, ys, zs):
 cellId = flow.getCell(x,y,z)
 if cellId in cellsHit: continue
 vv=flow.getCellInvVoidVolume(cellId)
 print(vv)

Robert Caulk (rcaulk) said : #7

Hello,

Your Yade version/installation method is a necessary piece of information
if you hope for any kind of assistance regarding Yade versions [1].

Assuming you have the correct Yade version: you do not appear to be
following the instructions dictated by the doc [2].

Cheers,

Robert

[1]https://www.yade-dem.org/wiki/Howtoask
[2]
https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.getCellInvVoidVolume

On Thu, Oct 22, 2020 at 1:36 AM Othman Sh <
<email address hidden>> wrote:

> Question #693297 on Yade changed:
> https://answers.launchpad.net/yade/+question/693297
>
> Status: Solved => Open
>
> Othman Sh is still having a problem:
> Hello,
>
> It has been a couple of weeks since the merge request mentioned in reply
> #4 but the getCellInvVoidVolume function is not functioning and it
> always return 0 values. I tried it using the code below. I appreciate
> any help in fixing this issue.
>
> Thank you,
> Othman
>
> ------------------------------------------------------
> from yade import plot,pack, export, ymport
> import numpy as np
>
>
> young=1e6
> compFricDegree = 3
> mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing
>
>
> O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
>
> O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
> walls=aabbWalls([mn,mx],thickness=0,material='walls')
> wallIds=O.bodies.append(walls)
>
>
> sp=pack.SpherePack()
> sp.makeCloud(mn,mx,0.05,.2)
> sp.toSimulation(material='spheres')
>
> yade.qt.View()
>
> print ('porosity = ', utils.porosity())
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_FrictPhys()],
> [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
> ),
> FlowEngine(dead=1,label="flow"),
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
> NewtonIntegrator(damping=0.2)
> ]
>
>
> flow.dead=0
> flow.useSolver=3
> flow.permeabilityFactor=1
> flow.viscosity=1e-3
> flow.bndCondIsPressure=[1,1,1,1,1,1]
> flow.bndCondValue=[0,0,0,0,0,10]
> flow.boundaryUseMaxMin=[1,1,1,1,1,1]
> O.dt=1e-6
> O.dynDt=False
> flow.emulateAction()
>
>
> ################## Get pores volume ##################
> box=aabbExtrema()
> lower_corner=box[0]
> upper_corner=box[1]
> mnx=lower_corner[0]
> mxx=upper_corner[0]
> mny=lower_corner[1]
> mxy=upper_corner[1]
> mnz=lower_corner[2]
> mxz=upper_corner[2]
> cellsHit = []
> numPoints = 20
> xs = np.linspace(0.95*mnx,1.05*mxx,numPoints)
> ys = np.linspace(0.95*mny,1.05*mxy,numPoints)
> zs = np.linspace(0.95*mnz,1.05*mxz,numPoints)
> v = np.array([0,0,0])
>
> for x,y,z in itertools.product(xs, ys, zs):
> cellId = flow.getCell(x,y,z)
> if cellId in cellsHit: continue
> vv=flow.getCellInvVoidVolume(cellId)
> print(vv)
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Othman Sh (othman-sh) said : #8

Hi Robert,

I am using Yadedaily. I have added this "flow.iniVoidVolumes=True" to my code that I posted in #6 and it is working now. Thanks for your help. This solves my problem.