Problem with Domain Area and Length#

In order to find network properties such as permeability using Darcy’s law, it is necessary to know the domain length and area. At first glance this might seem as simple as finding the maxima and minima of the pore coordinates in the direction of interest; however it is a bit more subtle than that as will be explained in this notebook. Consider the simple cubic network:

[1]:
import matplotlib.pyplot as plt
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
import numpy as np
np.random.seed(10)
pn = op.network.Cubic(shape=[4, 4, 1])

Now let’s plot the coordinates and connections between them:

[2]:
# NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots()
op.topotools.plot_coordinates(pn, markersize=500, ax=ax)
op.topotools.plot_connections(pn, ax=ax)
[2]:
<matplotlib.collections.LineCollection at 0x7f87745bf2e0>
../../../_images/examples_reference_uncategorized_the_problem_with_domain_length_and_area_3_1.svg

Let’s say we’d like to know the length of the domain from left to right. The pore coordinates for the leftmost pores are:

[3]:
pn['pore.coords'][pn.pores('left'), :]
[3]:
array([[0.5, 0.5, 0.5],
       [0.5, 1.5, 0.5],
       [0.5, 2.5, 0.5],
       [0.5, 3.5, 0.5]])

Note that the ‘minimum’ value is 0.5 rather than 0.0. This is because the pore coordinates point to the center of the region that belongs to the pore. In the case of the simple cubic with a spacing of 1, pore 0 lies at [0.5, 0.5, 0.5], while the bounding box around the pore goes from [0, 0, 0], to [1, 1, 1].

Therefore, if we use pore coordinates to find the domain length we’d do:

[4]:
L = pn['pore.coords'][:, 0].max() - pn['pore.coords'][:, 0].min()
print(L)
3.0

While in reality to domain is 4 pore long and the correct answer should be 4.0. The extra 0.5 of the left and right cells are not included.

It may seem simple enough to automatically add a full lattice cell to the calculation, but this cannot be assumed since many networks are random. Consider the Delaunay network:

[5]:
dn = op.network.Delaunay(points=20, shape=[1, 1, 0])

Let’s first remove the boundary pores to ensure the network is fully random:

[6]:
op.topotools.trim(network=dn, pores=dn.pores(['left', 'right', 'front', 'back']))
[7]:
# NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots()
op.topotools.plot_coordinates(dn, markersize=500, ax=ax)
op.topotools.plot_connections(dn, ax=ax)
[7]:
<matplotlib.collections.LineCollection at 0x7f872fbaac40>
../../../_images/examples_reference_uncategorized_the_problem_with_domain_length_and_area_12_1.svg

Now it’s quite clear that finding the domain size by assuming a bounding box around the pore coordinates does not tell the whole story. How much ‘extra’ length should be added beyond the extreme pores? This is not possible to know based only on the pore coordinates, which is the only information OpenPNM has.

The end result of this scenario is that when computing anything that requires length and area, these values must be specified by the user. Consider the permeability coefficient:

[8]:
pn = op.network.Cubic(shape=[4, 4, 4])
geo = op.geometry.SpheresAndCylinders(network=pn, pores=pn.Ps, throats=pn.Ts)
air = op.phases.Air(network=pn)
phys = op.physics.Standard(network=pn, phase=air, geometry=geo)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [8], in <cell line: 3>()
      1 pn = op.network.Cubic(shape=[4, 4, 4])
      2 geo = op.geometry.SpheresAndCylinders(network=pn, pores=pn.Ps, throats=pn.Ts)
----> 3 air = op.phases.Air(network=pn)
      4 phys = op.physics.Standard(network=pn, phase=air, geometry=geo)

AttributeError: module 'openpnm' has no attribute 'phases'
[9]:
sf = op.algorithms.StokesFlow(network=pn, phase=air)
sf.set_value_BC(pores=pn.pores('left'), values=200000)
sf.set_value_BC(pores=pn.pores('right'), values=100000)
sf.run()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [9], in <cell line: 1>()
----> 1 sf = op.algorithms.StokesFlow(network=pn, phase=air)
      2 sf.set_value_BC(pores=pn.pores('left'), values=200000)
      3 sf.set_value_BC(pores=pn.pores('right'), values=100000)

NameError: name 'air' is not defined

If you do not specify the domain area and length, OpenPNM will attempt to estimate these values using the faulty logic of minima and maxima of boundary pores, but will issue a warning so you know that errors are present.

[10]:
perm=op.metrics.AbsolutePermeability(network=pn)
K=perm.run()
print(K)
------------------------------------------------------------
    WARNING    : Attempting to estimate inlet area...will be low
    SOURCE     : openpnm.topotools._topotools.get_domain_area
    TIME STAMP : 2022-04-20 03:06:54,964
------------------------------------------------------------
------------------------------------------------------------
    WARNING    : Attempting to estimate domain length...could be low if
             boundary pores were not added
    SOURCE     : openpnm.topotools._topotools.get_domain_length
    TIME STAMP : 2022-04-20 03:06:54,966
------------------------------------------------------------
9.53562132281558e-05

To get correct answers, it is necessary to specify the domain sizes as follows:

[11]:
A = 16
L = 4
perm=op.metrics.AbsolutePermeability(network=pn)
perm.settings._update({
    'area': A,
    'length':L})
K=perm.run()
print(K)
7.151715992111685e-05

Of course, the most correct way to calculate K is to do it manually using Darcy’s law:

[12]:
mu = air['pore.viscosity'].mean()
Q = sf.rate(pores=pn.pores('left'))
dP = 200000 - 100000
K = Q*mu*L/(A * dP)
print(K)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [12], in <cell line: 1>()
----> 1 mu = air['pore.viscosity'].mean()
      2 Q = sf.rate(pores=pn.pores('left'))
      3 dP = 200000 - 100000

NameError: name 'air' is not defined