Repository: PMEAL/OpenPNM
Pull request #1251 (2019-04-21 23:55:57)

Add bundle of tubes
Adding a class to the ``materials`` module for creating a bundle of tubes.  Returns a project object, including a network and a geometry.


Patch:
diff --git a/example_bot.py b/example_bot.py
new file mode 100644
index 000000000..1130d9aac
--- /dev/null
+++ b/example_bot.py
@@ -0,0 +1,18 @@
+import openpnm as op
+ws = op.Workspace()
+ws.settings['loglevel'] = 30
+proj = ws.new_project()
+
+net, geo, phase = op.materials.BundleOfTubes(shape=50, spacing=0.0001, length=0.001,
+                                             psd_params={'distribution': 'weibull',
+                                                         'loc': 0.000001,
+                                                         'scale': 0.00004,
+                                                         'shape': 2.2},
+                                             settings={'adjust_psd': 'normalize',
+                                                       'seed': 0})
+hg = op.phases.Mercury(network=net)
+phys = op.physics.Classic(network=net, phase=hg, geometry=geo)
+mip = op.algorithms.Porosimetry(network=net, phase=hg)
+mip.set_inlets(pores=net.pores('top'))
+mip.run()
+mip.plot_intrusion_curve()
diff --git a/openpnm/__init__.py b/openpnm/__init__.py
index 298f5728c..0ae07c6f0 100644
--- a/openpnm/__init__.py
+++ b/openpnm/__init__.py
@@ -52,7 +52,7 @@ It consists of the following submodules:
 
 """
 
-__version__ = '2.0.3'
+__version__ = '2.0.3-dev'
 
 from . import utils
 from .utils import Workspace, Project
diff --git a/openpnm/materials/BundleOfTubes.py b/openpnm/materials/BundleOfTubes.py
new file mode 100644
index 000000000..b5f0594b6
--- /dev/null
+++ b/openpnm/materials/BundleOfTubes.py
@@ -0,0 +1,184 @@
+import scipy as sp
+import scipy.stats as spst
+from openpnm.utils import logging, Project
+from openpnm.network import Cubic
+from openpnm.geometry import GenericGeometry
+from openpnm.phases import GenericPhase
+from openpnm.topotools import trim
+import openpnm.models as mods
+logger = logging.getLogger(__name__)
+
+defsets = {'adjust_psd': 'clip'}
+
+
+class BundleOfTubes(Project):
+    r"""
+    The materials class creats a network that matches the bundle-of-tubes model
+
+    An OpenPNM project object is returned that contain a network with a
+    bundle-of-tubes topology, and geometry object with the necessary pore
+    size information, and a phase object with pre-defined pore-scale physics
+    models attached.  Note that this phase object does not have any actual
+    thermophysical properties which must be added by the user.
+
+    Parameters
+    ----------
+    shape : array_like or int
+        The number of pores in the X and Y direction of the domain.  It will
+        be 2 pores thick in the Z direction, as required to be a bundle of
+        tubes model.  If an ``int`` is given it will be applied to both the
+        X and Y directions.
+
+    spacing : array_like or float
+        The spacing between the tubes in the X and Y direction.  If a ``float``
+        is given it will be applied to both the X and Y directions.
+
+    length : float
+        The length of the tubes or thickness of the domain in the z-direction.
+
+    psd_params : dictionary
+        The parameters of the statistical distribution of the pore sizes.
+        The dictionary must contain the type of distribution to use (specified
+        as 'distribution'), selected from the ``scipy.stats`` module, and the
+        parameters corresponding to the chosen distribution.
+        The default is: ``{'distribution': 'norm', 'loc': 0.5, 'scale': 0.11}.
+        Note that ``scipy.stats`` uses *loc* and *scale* to be consistent
+        between different distribution types, instead of things like *mean*
+        and *stddev*.
+
+    name : string, optional
+        The name to give the Project
+
+    """
+
+    def __init__(self,
+                 shape,
+                 spacing=1.0,
+                 length=1.0,
+                 psd_params={'distribution': 'norm',
+                             'loc': None,
+                             'scale': None},
+                 name=None,
+                 settings={},
+                 **kwargs):
+        super().__init__(name=name)
+        self.settings.update(defsets)
+        self.settings.update(settings)
+
+        if isinstance(shape, int):
+            shape = sp.array([shape, shape, 2])
+        elif len(shape) == 2:
+            shape = sp.concatenate((sp.array(shape), [2]))
+        else:
+            raise Exception('shape not understood, must be int ' +
+                            ' or list of 2 ints')
+
+        if isinstance(spacing, float) or isinstance(spacing, int):
+            spacing = float(spacing)
+            self.settings['spacing'] = spacing
+            spacing = sp.array([spacing, spacing, length])
+        else:
+            raise Exception('spacing not understood, must be float')
+
+        net = Cubic(shape=shape, spacing=spacing, project=self, **kwargs)
+        Ps_top = net.pores('top')
+        Ps_bot = net.pores('bottom')
+        Ts = net.find_connecting_throat(P1=Ps_top, P2=Ps_bot)
+        Ts = net.tomask(throats=Ts)
+        trim(network=net, throats=~Ts)
+
+        geom = GenericGeometry(network=net, pores=net.Ps, throats=net.Ts)
+
+        geom.add_model(propname='throat.seed',
+                       model=mods.geometry.throat_seed.random)
+        if psd_params['loc'] is None:
+            psd_params['loc'] = spacing[0]/2
+        if psd_params['scale'] is None:
+            psd_params['scale'] = spacing[0]/10
+        if psd_params['distribution'] in ['norm', 'normal', 'gaussian']:
+            geom.add_model(propname='throat.size_distribution',
+                           seeds='throat.seed',
+                           model=mods.geometry.throat_size.normal,
+                           loc=psd_params['loc'], scale=psd_params['scale'])
+        elif psd_params['distribution'] in ['weibull']:
+            geom.add_model(propname='throat.size_distribution',
+                           seeds='throat.seed',
+                           model=mods.geometry.throat_size.weibull,
+                           loc=psd_params['loc'],
+                           scale=psd_params['scale'],
+                           shape=psd_params['shape'])
+        else:
+            temp = psd_params.copy()
+            func = getattr(spst, temp.pop('distribution'))
+            psd = func.freeze(**temp)
+            geom.add_model(propname='throat.size_distribution',
+                           seeds='throat.seed',
+                           model=mods.geometry.throat_size.generic_distribution,
+                           func=psd)
+
+        if sp.any(geom['throat.size_distribution'] < 0):
+            logger.warning('Given size distribution produced negative ' +
+                           'throat diameters...these will be set to 0')
+        geom.add_model(propname='throat.diameter',
+                       model=mods.misc.clip,
+                       prop='throat.size_distribution',
+                       xmin=1e-12, xmax=sp.inf)
+
+        if self.settings['adjust_psd'] is None:
+            if geom['throat.size_distribution'].max() > spacing[0]:
+                logger.warning('Given size distribution produced throats ' +
+                               'larger than the spacing.')
+
+        elif self.settings['adjust_psd'] == 'clip':
+            geom.add_model(propname='throat.diameter',
+                           model=mods.misc.clip,
+                           prop='throat.size_distribution',
+                           xmin=1e-12, xmax=spacing[0])
+            if geom['throat.size_distribution'].max() > spacing[0]:
+                logger.warning('Given size distribution produced throats ' +
+                               'larger than the spacing...tube diameters ' +
+                               'will be clipped between 0 and given spacing')
+
+        elif self.settings['adjust_psd'] == 'normalize':
+            tmin = max(1e-12, geom['throat.size_distribution'].min())
+            geom.add_model(propname='throat.diameter',
+                           model=mods.misc.normalize,
+                           prop='throat.size_distribution',
+                           xmin=tmin, xmax=spacing[0])
+            if geom['throat.size_distribution'].max() > spacing[0]:
+                logger.warning('Given size distribution produced throats ' +
+                               'larger than the spacing...tube diameters ' +
+                               'will be normalized to fit given spacing')
+        else:
+            logger.warning('Settings not understood, ignoring')
+
+        geom.add_model(propname='pore.diameter',
+                       model=mods.geometry.pore_size.from_neighbor_throats,
+                       throat_prop='throat.diameter', mode='max')
+        geom.add_model(propname='pore.diameter',
+                       model=mods.misc.constant, value=0.0)
+        geom.add_model(propname='throat.length',
+                       model=mods.geometry.throat_length.ctc)
+        geom.add_model(propname='throat.area',
+                       model=mods.geometry.throat_area.cylinder)
+        geom.add_model(propname='pore.area',
+                       model=mods.misc.from_neighbor_throats,
+                       throat_prop='throat.area')
+        geom.add_model(propname='pore.volume',
+                       model=mods.misc.constant, value=0.0)
+        geom.add_model(propname='throat.volume',
+                       model=mods.geometry.throat_volume.cylinder)
+
+        geom.regenerate_models()
+
+        # Now create a generic phase with physics models on it
+        phase = GenericPhase(network=net)
+        m = mods.physics.hydraulic_conductance.classic_hagen_poiseuille
+        phase.add_model(propname='throat.hydraulic_conductance',
+                        model=m, regen_mode='deferred')
+        m = mods.physics.diffusive_conductance.classic_ordinary_diffusion
+        phase.add_model(propname='throat.diffusive_conductance',
+                        model=m, regen_mode='deferred')
+        m = mods.physics.diffusive_conductance.classic_ordinary_diffusion
+        phase.add_model(propname='throat.entry_pressure',
+                        model=m, regen_mode='deferred')
diff --git a/openpnm/materials/__init__.py b/openpnm/materials/__init__.py
index 811a14f7b..2eaacb83a 100644
--- a/openpnm/materials/__init__.py
+++ b/openpnm/materials/__init__.py
@@ -26,3 +26,4 @@ The table below gives a list of available Material generators:
 """
 
 from .VoronoiFibers import VoronoiFibers
+from .BundleOfTubes import BundleOfTubes
diff --git a/openpnm/models/geometry/__init__.py b/openpnm/models/geometry/__init__.py
index 48ce5b4e8..c836e9c6a 100644
--- a/openpnm/models/geometry/__init__.py
+++ b/openpnm/models/geometry/__init__.py
@@ -15,6 +15,7 @@ from . import pore_surface_area
 from . import pore_area
 from . import throat_endpoints
 from . import throat_area
+from . import throat_seed
 from . import throat_size
 from . import throat_length
 from . import throat_perimeter
diff --git a/openpnm/models/geometry/pore_size.py b/openpnm/models/geometry/pore_size.py
index fb79bcd91..a223c938d 100644
--- a/openpnm/models/geometry/pore_size.py
+++ b/openpnm/models/geometry/pore_size.py
@@ -51,6 +51,15 @@ def generic_distribution(target, func, seeds='pore.seed'):
 generic_distribution.__doc__ = _misc.generic_distribution.__doc__
 
 
+def from_neighbor_throats(target, throat_prop, mode='max'):
+    return _misc.from_neighbor_throats(target=target,
+                                       throat_prop=throat_prop,
+                                       mode=mode)
+
+
+from_neighbor_throats.__doc__ = _misc.from_neighbor_throats.__doc__
+
+
 def largest_sphere(target, fixed_diameter='pore.fixed_diameter', iters=5):
     r"""
     Finds the maximum diameter pore that can be placed in each location without
diff --git a/openpnm/models/geometry/pore_volume.py b/openpnm/models/geometry/pore_volume.py
index 6963a2cae..9c6cec555 100644
--- a/openpnm/models/geometry/pore_volume.py
+++ b/openpnm/models/geometry/pore_volume.py
@@ -89,3 +89,29 @@ def square(target, pore_diameter='pore.diameter'):
 
     """
     return target[pore_diameter]**2
+
+
+def cylinder(target, pore_diameter='pore.diameter'):
+    r"""
+    Calculate pore volume from diameter assuming a cylindrical pore body
+    with a height equal to its diameter.
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    pore_diameter : string
+        The dictionary key of the pore diameter values
+
+    Returns
+    -------
+    value : NumPy ndarray
+        Array containing pore volume values.
+
+    """
+    diams = target[pore_diameter]
+    value = _pi/4*diams**2*diams
+    return value
diff --git a/openpnm/models/geometry/throat_endpoints.py b/openpnm/models/geometry/throat_endpoints.py
index ad207b79a..0ad31f744 100644
--- a/openpnm/models/geometry/throat_endpoints.py
+++ b/openpnm/models/geometry/throat_endpoints.py
@@ -34,7 +34,7 @@ def cubic_pores(target, pore_diameter='pore.diameter'):
     throats = network.map_throats(throats=target.Ts, origin=target)
     xyz = network['pore.coords']
     cn = network['throat.conns'][throats]
-    L = _ctc(target=target, pore_diameter=pore_diameter) + 1e-15
+    L = _ctc(target=target) + 1e-15
     D1 = network[pore_diameter][cn[:, 0]]
     D2 = network[pore_diameter][cn[:, 1]]
     unit_vec = (xyz[cn[:, 1]] - xyz[cn[:, 0]]) / L[:, None]
@@ -124,7 +124,7 @@ def spherical_pores(target, pore_diameter='pore.diameter',
     throats = network.map_throats(throats=target.Ts, origin=target)
     xyz = network['pore.coords']
     cn = network['throat.conns'][throats]
-    L = _ctc(target=target, pore_diameter=pore_diameter) + 1e-15
+    L = _ctc(target=target) + 1e-15
     Dt = network[throat_diameter][throats]
     D1 = network[pore_diameter][cn[:, 0]]
     D2 = network[pore_diameter][cn[:, 1]]
diff --git a/openpnm/models/geometry/throat_length.py b/openpnm/models/geometry/throat_length.py
index 242d345dc..0116148c8 100644
--- a/openpnm/models/geometry/throat_length.py
+++ b/openpnm/models/geometry/throat_length.py
@@ -11,10 +11,10 @@ from openpnm.utils import logging as _logging
 _logger = _logging.getLogger(__name__)
 
 
-def ctc(target, pore_diameter='pore.diameter'):
+def ctc(target):
     r"""
     Calculate throat length assuming point-like pores, i.e. center-to-center
-    distance between pores. Also, this models assumes that pores and throat
+    distance between pores. Also, this model assumes that pores and throat
     centroids are colinear.
 
     Parameters
@@ -24,26 +24,18 @@ def ctc(target, pore_diameter='pore.diameter'):
         length of the calculated array, and also provides access to other
         necessary properties.
 
-    pore_diameter : string
-        Dictionary key of the pore diameter values
-
     Returns
     -------
     value : NumPy ndarray
         Array containing throat length values.
 
     """
-    _np.warnings.filterwarnings('ignore', category=RuntimeWarning)
-
     network = target.project.network
     throats = network.map_throats(throats=target.Ts, origin=target)
     cn = network['throat.conns'][throats]
     C1 = network['pore.coords'][cn[:, 0]]
     C2 = network['pore.coords'][cn[:, 1]]
     value = _sqrt(((C1 - C2)**2).sum(axis=1))
-
-    _np.warnings.filterwarnings('default', category=RuntimeWarning)
-
     return value
 
 
@@ -79,7 +71,6 @@ def piecewise(target, throat_endpoints='throat.endpoints',
     length. This could be useful for Voronoi or extracted networks.
 
     """
-    _np.warnings.filterwarnings('ignore', category=RuntimeWarning)
     network = target.project.network
     throats = network.map_throats(throats=target.Ts, origin=target)
     # Get throat endpoints
@@ -94,9 +85,6 @@ def piecewise(target, throat_endpoints='throat.endpoints',
             _sqrt(((Ct - EP2)**2).sum(axis=1))
     except KeyError:
         pass
-
-    _np.warnings.filterwarnings('default', category=RuntimeWarning)
-
     return Lt
 
 
@@ -129,8 +117,6 @@ def conduit_lengths(target, throat_endpoints='throat.endpoints',
     keys 'pore1', 'pore2', and 'throat'.
 
     """
-    _np.warnings.filterwarnings('ignore', category=RuntimeWarning)
-
     network = target.project.network
     throats = network.map_throats(throats=target.Ts, origin=target)
     cn = network['throat.conns'][throats]
@@ -149,7 +135,29 @@ def conduit_lengths(target, throat_endpoints='throat.endpoints',
     # Calculate conduit lengths
     L1 = _sqrt(((C1 - EP1)**2).sum(axis=1))
     L2 = _sqrt(((C2 - EP2)**2).sum(axis=1))
+    return {'pore1': L1, 'throat': Lt, 'pore2': L2}
 
-    _np.warnings.filterwarnings('default', category=RuntimeWarning)
 
-    return {'pore1': L1, 'throat': Lt, 'pore2': L2}
+def classic(target, pore_diameter='pore.diameter'):
+    r"""
+    Find throat length as the pore-to-pore center distance, less the radii of
+    each pore.
+
+    Parameters
+    ----------
+    target : OpenPNM object
+        The object which this model is associated with. This controls the
+        length of the calculated array, and also provides access to other
+        necessary properties.
+
+    pore_diameter : string
+        Dictionary key of the pore diameter values
+    """
+    network = target.project.network
+    throats = network.map_throats(throats=target.Ts, origin=target)
+    cn = network['throat.conns'][throats]
+    C1 = network['pore.coords'][cn[:, 0]]
+    C2 = network['pore.coords'][cn[:, 1]]
+    D = _sqrt(((C1 - C2)**2).sum(axis=1))
+    value = D - _np.sum(network[pore_diameter][cn], axis=1)/2
+    return value
diff --git a/openpnm/models/misc/__init__.py b/openpnm/models/misc/__init__.py
index 1799224b1..2091da78b 100644
--- a/openpnm/models/misc/__init__.py
+++ b/openpnm/models/misc/__init__.py
@@ -12,5 +12,5 @@ generating an array of random numbers.
 
 from .statistical_distributions import generic_distribution, normal, weibull, random
 from .simple_equations import linear, polynomial, generic_function
-from .basic_math import constant, product, scaled
+from .basic_math import constant, product, scaled, clip, normalize
 from .neighbor_lookups import from_neighbor_throats, from_neighbor_pores
diff --git a/openpnm/models/misc/basic_math.py b/openpnm/models/misc/basic_math.py
index a8ef1ae5a..6f0a3df3a 100644
--- a/openpnm/models/misc/basic_math.py
+++ b/openpnm/models/misc/basic_math.py
@@ -3,16 +3,63 @@ r"""
 .. autofunction:: openpnm.models.misc.basic_math.constant
 .. autofunction:: openpnm.models.misc.basic_math.product
 .. autofunction:: openpnm.models.misc.basic_math.scaled
+.. autofunction:: openpnm.models.misc.basic_math.clip
+.. autofunction:: openpnm.models.misc.basic_math.normalize
 
 """
-
+import numpy as np
 from openpnm.utils import logging
 logger = logging.getLogger(__name__)
 
 
+def normalize(target, prop, xmin=0, xmax=1):
+    r"""
+    Normalizes the given array between the supplied limits
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object which this model is associated with. This controls the
+        length of the calculated array, and also provides access to other
+        necessary properties.
+
+    xmin : float
+        Lower limit of the re-scaled data
+
+    xmax : float
+        Upper limit of the re-scaled data
+    """
+    vals = target[prop]
+    # Scale to 0 to 1
+    vals = (vals - vals.min())/(vals.max() - vals.min())
+    vals = vals*(xmax - xmin) + xmin
+    return vals
+
+
+def clip(target, prop, xmax, xmin=0):
+    r"""
+    Clips the given array within the supplied limits
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object which this model is associated with. This controls the
+        length of the calculated array, and also provides access to other
+        necessary properties.
+
+    xmin : float
+        Values below this limit will be replaced with ``xmin``.
+
+    xmax : float
+        Values above this limit will be replaced with ``xmax``.
+    """
+    vals = np.clip(target[prop], xmin, xmax)
+    return vals
+
+
 def constant(target, value):
     r"""
-    Places a constant value into the target object
+    Places the given constant value into the target object
 
     Parameters
     ----------


Test patch:
diff --git a/tests/unit/materials/BundleOfTubesTest.py b/tests/unit/materials/BundleOfTubesTest.py
new file mode 100644
index 000000000..a023b6733
--- /dev/null
+++ b/tests/unit/materials/BundleOfTubesTest.py
@@ -0,0 +1,69 @@
+from openpnm.materials import BundleOfTubes
+import scipy as sp
+import pytest
+import matplotlib.pyplot as plt
+
+
+class BundleOfTubesTest:
+
+    def setup_class(self):
+        pass
+
+    def test_instantiate_with_defaults(self):
+        net, geo, phase = BundleOfTubes(shape=[30, 30])
+        assert net.Np == 1800
+        assert net.Nt == (net.Np/2)
+
+    def test_instantiate_custom_shape_and_size(self):
+        proj = BundleOfTubes(shape=30, spacing=0.001, length=0.01,
+                             psd_params={'distribution': 'weibull',
+                                         'scale': 0.0002,
+                                         'shape': 2,
+                                         'loc': 1e-12})
+
+    def test_instantiate_custom_distribution(self):
+        proj = BundleOfTubes(shape=30, spacing=0.001, length=0.01,
+                             psd_params={'distribution': 'lognorm',
+                                         'scale': 0.0002,
+                                         's': 2,
+                                         'loc': 1e-12})
+
+    def test_instantiate_with_exceptions(self):
+        with pytest.raises(Exception):
+            proj = BundleOfTubes(shape=[30, 30, 30], spacing=0.001, length=0.01)
+        with pytest.raises(Exception):
+            proj = BundleOfTubes(shape=30, spacing=[0.001, 0.001], length=0.01)
+
+    def test_instantiate_with_settings(self):
+        net, geo, phase = BundleOfTubes(shape=30, spacing=0.001, length=0.01,
+                                        psd_params={'distribution': 'normal',
+                                                    'loc': 0.02,
+                                                    'scale': 2},
+                                        settings={'adjust_psd': 'clip'})
+        assert geo['throat.size_distribution'].max() > geo['throat.diameter'].max()
+
+        net, geo, phase = BundleOfTubes(shape=30, spacing=0.001, length=0.01,
+                                        psd_params={'distribution': 'norm',
+                                                    'loc': 0.02,
+                                                    'scale': 0.001},
+                                        settings={'adjust_psd': 'clip'})
+        assert geo['throat.size_distribution'].max() > geo['throat.diameter'].max()
+
+        net, geo, phase = BundleOfTubes(shape=30, spacing=1.0, length=0.01,
+                                        psd_params={'distribution': 'normal',
+                                                    'loc': 0.02,
+                                                    'scale': 0.001},
+                                        settings={'adjust_psd': None})
+        assert sp.all(geo['throat.size_distribution'] == geo['throat.diameter'])
+
+
+if __name__ == '__main__':
+
+    t = BundleOfTubesTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
+    self = t