Yade

Pore volume - FlowEngine

Asked by Othman Sh

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:
Last query:
Last reply:

This question was reopened

Revision history for this message
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

Revision history for this message
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)

Revision history for this message
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...

Revision history for this message
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

Revision history for this message
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

Revision history for this message
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)

Revision history for this message
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
>

Revision history for this message
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.