# [OpenPNM] Part 2 ——Geometric Properties

Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------

# 指定几何属性

• 手动计算孔和流道的几何属性（不推荐，太麻烦）

• 使用库中的孔隙尺度模型（推荐）

• 依赖项处理程序概述

import numpy as np
import matplotlib.pyplot as plt
import openpnm as op
op.visualization.set_mpl_style()
np.random.seed(0)
pn = op.network.Cubic([20, 20, 20], spacing=5e-5)
print(pn)


# 手动计算几何属性

## 从分布中添加孔流道直径

scipy.stats模块定义了许多统计学分布，首先，让我们使用正态分布来生成1到50 um范围内的孔径

np.random.seed(0)
import scipy.stats as spst
psd = spst.norm.rvs(loc=25, scale=6, size=pn.Np)
plt.hist(psd, edgecolor='k')
plt.xlim([0, 50]);


Definition : rvs(*args, **kwds)
Random variates of given type.

Parameters

arg1, arg2, arg3,...array_like
The shape parameter(s) for the distribution (see docstring of the instance object for more information).

loc:array_like, optional
Location parameter (default=0).

scale:array_like, optional
Scale parameter (default=1).

size:int or tuple of ints, optional
Defining number of random variates (default is 1).

random_state:{None, int, numpy.random.Generator,
numpy.random.RandomState}, optional

If seed is None (or np.random), the numpy.random.RandomState singleton is used. If seed is an int, a new RandomState instance is used, seeded with seed. If seed is already a Generator or RandomState instance then that instance is used.

Returns

rvsndarray or scalar
Random variates of given size.

plt.hist()相关说明可以参考：https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.hist.html

#上面的分布看起来不错，但是要确保我们的分布不是太宽，以至于它的值大于 50 um：

print(psd.min())
print(psd.max())
'''
2.5593961722893255
47.80996128980269
'''
#现在，我们将这些值转换单位并分配给网络：

pn['pore.diameter'] = psd*1e-6  # um to m


np.random.seed(0)
tsd = spst.weibull_min.rvs(c=1.5, loc=.5, scale=7.5, size=pn.Nt)
plt.hist(tsd, edgecolor='k')
plt.xlim([0, 50]);


spst.weibull_min.rvs()函数的相关信息可以参考：https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.weibull_min.html

#同样，让我们检查高值和低值：

print(tsd.min())
print(tsd.max())
'''
0.5130345457395142
36.96861960231873
'''
#所以我们可以看到我们的流道直径小到500纳米，大到37微米。这些也可以分配给网络：

pn['throat.diameter'] = tsd*1e-6  # um to m


hits = np.any(pn['pore.diameter'][pn.conns].T < pn['throat.diameter'], axis=0)
hits.sum()
'''
564
'''


np.any()函数的相关信息可以参考：https://numpy.org/doc/stable/reference/generated/numpy.any.html

tsd = np.amin(pn['pore.diameter'][pn.conns], axis=1)
plt.hist(tsd, edgecolor='k', density=True, alpha=0.5)#alpha设置透明度，0为完全透明
plt.hist(pn['pore.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.xlim([0, 50e-6]);


pn['throat.diameter'] = tsd*0.5
plt.hist(pn['throat.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.hist(pn['pore.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.xlim([0, 50e-6]);


https://en.wikipedia.org/wiki/Random_seed
https://www.geeksforgeeks.org/random-seed-in-python/

lo, hi = 0.01, 0.99
pn['pore.seed'] = np.random.rand(pn.Np)*(hi - lo) + lo
#设置种子的和值可以防止位于分布尾部较远的值出现，从而导致过大的孔径或喉道直径

#使用点概率函数返回给定概率下的累积分布值，因此上面计算的种子值可以使用以下代码与大小值相关联
import scipy.stats as spst
psd = spst.norm.ppf(loc=25, scale=6, q=pn['pore.seed'])
pn['pore.diameter'] = psd*1e-6
plt.hist(pn['pore.diameter'], edgecolor='k')
plt.xlim([0, 50e-6]);


Definition : ppf(q, *args, **kwds)
Percent point function (inverse of cdf) at q of the given RV.

Parameters

q:array_like
lower tail probability

arg1, arg2, arg3,...array_like
The shape parameter(s) for the distribution (see docstring of the instance object for more information)

loc:array_like, optional
location parameter (default=0)

scale:array_like, optional
scale parameter (default=1)

Returns

x:array_like
quantile corresponding to the lower tail probability q.

pn['throat.seed'] = np.amin(pn['pore.seed'][pn.conns], axis=1)


tsd = spst.norm.ppf(loc=25, scale=6, q=pn['throat.seed'])
pn['throat.diameter'] = tsd*1e-6
plt.hist(pn['throat.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.hist(pn['pore.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.xlim([0, 50e-6]);


## 计算流道长度

R1, R2 = (pn['pore.diameter'][pn.conns]/2).T
L_total = np.sqrt(np.sum(np.diff(pn.coords[pn.conns], axis=1).squeeze()**2, axis=1))
'''
np.diff():Calculate the n-th discrete difference along the given axis.
np.sum():Sum of array elements over a given axis.
numpy.squeeze() function is used when we want to remove single-dimensional entries from the shape of an array.
'''
Lt = L_total - R1 - R2
print(L_total)


## 计算孔隙表面积

At = pn['throat.diameter']**2
SAp = (pn['pore.diameter']**2)*6
np.subtract.at(SAp, pn.conns[:, 0], At)
np.subtract.at(SAp, pn.conns[:, 1], At)
print(SAp)
'''
[[   0    1]
[   1    2]
[   2    3]
...
[7597 7997]
[7598 7998]
[7599 7999]]
[7.27174552e-09 3.85799662e-09 4.84581355e-09 ... 2.94168620e-09
2.75016941e-09 3.47031746e-09]
'''


1、从数组的第一列中列出的每个孔中减去喉道横截面积。由于一个给定的孔隙与潜在的许多喉道相连，因此必须使用减法函数的特殊版本来完成减法函数，该函数从所指示的位置中减去每一个值，并且以累积的方式这样做。因此，如果孔隙i出现了不止一次，则删除的值就不止一个。该操作在原地执行，因此不会返回任何数组

2、从数组的第二列中列出的每个孔中减去喉道横截面积

# 使用OpenPNM自带的pore-scale模型

np.random.seed(0)
pn = op.network.Cubic([20, 20, 20], spacing=5e-5)

#为每个孔分配一个随机数，使用孔隙尺度模型：np.random.rand

f = op.models.geometry.pore_seed.random
model=f,
seed=None,  # Seed value provided to numpy's random number generator.
num_range=[0.01, 0.99],)

#让我们让每个流道继承其每个相邻毛孔中的最小种子值，OpenPNM 为此提供了一个模型：
f = op.models.geometry.throat_seed.from_neighbor_pores
model=f,
prop='pore.seed',
mode='min')
"""
add_model(propname, model, domain='all', regen_mode='normal', *, map: Mapping[_KT, _VT], iterable: Iterable[Tuple[_KT, _VT]], **kwargs: _VT)

Add a pore-scale model to the object, along with the desired arguments

Parameters
----------
propname : str
The name of the property being computed. E.g. if
propname='pore.diameter' then the computed results will be stored
in obj['pore.diameter'].
model : function handle
The function that produces the values
domain : str
The label indicating the locations in which results generated by
model should be stored. See Notes for more details.
regen_mode : str
How the model should be regenerated. Options are:

============ =====================================================
regen_mode   description
============ =====================================================
normal       (default) The model is run immediately upon being
added, and is also run each time regenerate_models
is called.
deferred     The model is NOT run when added, but is run each time
regenerate_models is called. This is useful for
models that depend on other data that may not exist
yet.
constant     The model is run immediately upon being added, but is
is not run when regenerate_models is called,
effectively turning the property into a constant.
============ =====================================================

kwargs : keyword arguments
All additional keyword arguments are passed on to the model

Notes
-----
The domain argument dictates where the results of model should
be stored. For instance, given propname='pore.diameter' and
domain='left' then when model is run, the results are stored in
in the pores labelled left. Note that if model returns Np
values, then values not belonging to 'pore.left' are discarded.
The following steps outline the process:

1. Find the pore indices:

.. code-block:: python

Ps = obj.pores('left')

2. Run the model:

.. code-block:: python

vals = model(**kwargs)

3. If the model returns a full Np-length array, then extract the
correct values and apply them to the corresponding locations:

.. code-block:: python

if len(vals) == obj.Np:
obj['pore.diameter'][Ps] = vals[Ps]

4. If the model was designed to return only the subset of values then:

.. code-block:: python

if len(vals) == obj.num_pores('left'):
obj['pore.diameter'][Ps] = vals

"""


f = op.models.geometry.pore_size.normal
model=f,
scale=1e-5,
loc=2.5e-5,
seeds='pore.seed')
plt.hist(pn['pore.diameter'], edgecolor='k');


f = op.models.geometry.throat_size.normal
model=f,
scale=1e-5,
loc=2.5e-5)
plt.hist(pn['pore.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.hist(pn['throat.diameter'], edgecolor='k', density=True, alpha=0.5);


f = op.models.geometry.throat_length.spheres_and_cylinders
model=f)
f1 = op.models.geometry.pore_volume.sphere
model=f1)
f2 = op.models.geometry.throat_volume.cylinder
model=f2)

#球形孔与其连接的流道之间存在重叠区域。需要额外的步骤来消除：
f = op.models.geometry.throat_length.spheres_and_cylinders
model=f)
f1 = op.models.geometry.pore_volume.sphere
model=f1)
f2 = op.models.geometry.throat_volume.cylinder
model=f2)
f3 = op.models.geometry.throat_volume.lens
model=f3)
f4 = op.models.misc.difference
model=f4,
props=['throat.total_volume', 'throat.lens_volume'])
#检查一下
print("Volumes of full throats:", pn['throat.total_volume'][:3])
print("Volumes of lenses:", pn['throat.lens_volume'][:3])
print("Actual throat volumes:", pn['throat.volume'][:3])
'''
Volumes of full throats: [2.27203214e-14 2.58686910e-14 2.44028684e-14]
Volumes of lenses: [6.92245838e-15 8.40091924e-15 7.59556144e-15]
Actual throat volumes: [1.57978631e-14 1.74677718e-14 1.68073070e-14]
'''



# 使用预定义的模型集合

pn = op.network.Cubic(shape=[20, 20, 20], spacing=5e-5)
#The models collection is fetched as follows:

from pprint import pprint
mods = op.models.collections.geometry.spheres_and_cylinders
pprint(mods)
'''
{'pore.diameter': {'model': ,
'props': ['pore.max_size', 'pore.seed']},
'pore.max_size': {'iters': 10,
'model': },
'pore.seed': {'element': 'pore',
'model': ,
'num_range': [0.2, 0.7],
'seed': None},
'pore.volume': {'model': ,
'pore_diameter': 'pore.diameter'},
'throat.cross_sectional_area': {'model': ,
'throat_diameter': 'throat.diameter'},
'throat.diameter': {'factor': 0.5,
'model': ,
'prop': 'throat.max_size'},
'throat.diffusive_size_factors': {'model': ,
'pore_diameter': 'pore.diameter',
'throat_diameter': 'throat.diameter'},
'throat.hydraulic_size_factors': {'model': ,
'pore_diameter': 'pore.diameter',
'throat_diameter': 'throat.diameter'},
'throat.length': {'model': ,
'pore_diameter': 'pore.diameter',
'throat_diameter': 'throat.diameter'},
'throat.lens_volume': {'model': ,
'pore_diameter': 'pore.diameter',
'throat_diameter': 'throat.diameter'},
'throat.max_size': {'mode': 'min',
'model': ,
'prop': 'pore.diameter'},
'throat.total_volume': {'model': ,
'throat_diameter': 'throat.diameter',
'throat_length': 'throat.length'},
'throat.volume': {'model': ,
'props': ['throat.total_volume', 'throat.lens_volume']
'''


pn.add_model_collection(mods)
pn.regenerate_models()
print(pn)


# 通过覆盖模型或其参数来自定义模型

OpenPNM 中包含的预先编写的模型和集合很有用，但通常需要调整某些模型以适应特定的应用程序。这可以通过更改现有模型的参数或将它们全部覆盖来完成。

spheres_and_cylinders

print(pn.models['pore.diameter@all'])


f = op.models.geometry.pore_size.normal
model=f,
scale=1e-5,
loc=2.5e-5,
seeds='pore.seed')
print(pn.models['pore.diameter@all'])
plt.hist(pn['pore.diameter'], edgecolor='k')


print(pn.models['pore.seed@all'])
pn.models['pore.seed@all']['num_range'] = [0.01, 0.99]
pn.regenerate_models()
plt.hist(pn['pore.diameter'], edgecolor='k');



# 依赖项处理程序简介

pn = op.network.Cubic(shape=[20, 20, 20], spacing=5e-5)
model=op.models.geometry.pore_seed.random)
model=f,
scale=1e-5,
loc=2.5e-5,
seeds='pore.seed')
model=op.models.geometry.throat_seed.from_neighbor_pores,
prop='pore.seed')
model=f,
scale=1e-5,
loc=2.5e-5,
seeds='throat.seed')
print(pn['pore.seed'])
print(pn['throat.diameter'])
'''
[0.15319619 0.0328116  0.2287953  ... 0.23887453 0.8411705  0.10988224]
[6.59011212e-06 6.59011212e-06 1.75717988e-05 ... 1.13561876e-05
2.66940036e-05 1.27284536e-05]
'''

pn.regenerate_models()
print(pn['pore.seed'])
print(pn['throat.diameter'])
'''
[0.36618201 0.35063831 0.40236746 ... 0.58822388 0.42461491 0.9016327 ]
[2.11640227e-05 2.11640227e-05 2.20171387e-05 ... 2.57748463e-05
2.30989878e-05 3.79091080e-05]
'''


pn.models.dependency_map(style='planar');


Everything not saved will be lost.