Merge branch 'master' into label_faces

OpenPNM/Algorithms/__InvasionPercolation__.py

@@ -7,6 +7,7 @@
 """
 import heapq as hq
 import scipy as sp
+import numpy as np
 from OpenPNM.Algorithms import GenericAlgorithm
 from OpenPNM.Base import logging
 logger = logging.getLogger(__name__)
@@ -192,3 +193,138 @@ def apply_flow(self, flowrate):
         t = sp.zeros((self.Nt,))
         t[b] = e  # Convert back to original order
         self._phase['throat.invasion_time'] = t
+
+    def apply_trapping(self, outlets):
+        """
+        Apply trapping based on algorithm described by Y. Masson [1].
+        It is applied as a post-process and runs the percolation algorithm in
+        reverse assessing the occupancy of pore neighbors. Consider the
+        following scenario when running standard IP without trapping,
+        3 situations can happen after each invasion step:
+            The number of defending clusters stays the same and clusters can
+            shrink
+            A cluster of size one is suppressed
+            A cluster is split into multiple clusters
+        In reverse the following opposite situations can happen:
+            The number of defending clusters stays the same and clusters can
+            grow
+            A cluster of size one is created
+            Mutliple clusters merge into one cluster
+        With trapping the reversed rules are adjusted so that:
+            Only clusters that do not connect to a sink can grow and merge.
+            At the point that a neighbor connected to a sink is touched the
+            trapped cluster stops growing as this is the point of trapping in
+            forward invasion time.
+
+        Logger info displays the invasion sequence and pore index and a message
+        with condition number based on the modified trapping rules and the
+        assignment of the pore to a given cluster.
+
+        Initially all invaded pores are given cluster label -1
+        Outlets / Sinks are given -2
+        New clusters that grow into fully trapped clusters are either
+        identified at the point of breakthrough or grow from nothing if the
+        full invasion sequence is run, they are assigned numbers from 0 up.
+
+        Ref:
+        [1] Masson, Y., 2016. A fast two-step algorithm for invasion
+        percolation with trapping. Computers & Geosciences, 90, pp.41-48
+
+        Parameters
+        ----------
+        outlets : list or array of pore indices for defending fluid to escape
+        through
+
+        Returns
+        -------
+        Creates a throat array called 'pore.clusters' in the Algorithm
+        dictionary. Any positive number is a trapped cluster
+
+        Also creates 2 boolean arrays Np and Nt long called '<element>.trapped'
+        """
+        # First see if network is fully invaded
+        invaded_ps = self['pore.invasion_sequence'] > -1
+        if ~np.all(invaded_ps):
+            # Put defending phase into clusters
+            clusters = self._net.find_clusters2(~invaded_ps)
+            # Identify clusters that are connected to an outlet and set to -2
+            # -1 is the invaded fluid
+            # -2 is the defender fluid able to escape
+            # All others now trapped clusters which grow as invasion is reversed
+            out_clusters = sp.unique(clusters[outlets])
+            for c in out_clusters:
+                if c >= 0:
+                    clusters[clusters == c] = -2
+        else:
+            # Go from end
+            clusters = np.ones(self._net.Np, dtype=int)*-1
+            clusters[outlets] = -2
+
+        # Turn into a list for indexing
+        inv_seq = np.vstack((self['pore.invasion_sequence'].astype(int),
+                             np.arange(0, self._net.Np, dtype=int))).T
+        # Reverse sort list
+        inv_seq = inv_seq[inv_seq[:, 0].argsort()][::-1]
+        next_cluster_num = np.max(clusters)+1
+        # For all the steps after the inlets are set up to break-through
+        # Reverse the sequence and assess the neighbors cluster state
+        stopped_clusters = np.zeros(self._net.Np, dtype=bool)
+        all_neighbors = self._net.find_neighbor_pores(self._net.pores(),
+                                                      flatten=False)
+        for un_seq, pore in inv_seq:
+            if pore not in outlets:  # Don't bother with outlets
+                nc = clusters[all_neighbors[pore]]  # Neighboring clusters
+                unique_ns = np.unique(nc[nc != -1])  # Unique Neighbors
+                seq_pore = "S:"+str(un_seq)+" P:"+str(pore)
+                if np.all(nc == -1):
+                    # This is the start of a new trapped cluster
+                    clusters[pore] = next_cluster_num
+                    next_cluster_num += 1
+                    msg = (seq_pore+" C:1 new cluster number: " +
+                           str(clusters[pore]))
+                    logger.info(msg)
+                elif len(unique_ns) == 1:
+                    # Grow the only connected neighboring cluster
+                    if not stopped_clusters[unique_ns[0]]:
+                        clusters[pore] = unique_ns[0]
+                        msg = (seq_pore+" C:2 joins cluster number: " +
+                               str(clusters[pore]))
+                        logger.info(msg)
+                    else:
+                        clusters[pore] = -2
+                elif -2 in unique_ns:
+                    # We have reached a sink neighbor, stop growing cluster
+                    msg = (seq_pore+" C:3 joins sink cluster")
+                    logger.info(msg)
+                    clusters[pore] = -2
+                    # Stop growth and merging
+                    stopped_clusters[unique_ns[unique_ns > -1]] = True
+                else:
+                    # We might be able to do some merging
+                    # Check if any stopped clusters are neighbors
+                    if np.any(stopped_clusters[unique_ns]):
+                        msg = (seq_pore+" C:4 joins sink cluster")
+                        logger.info(msg)
+                        clusters[pore] = -2
+                        # Stop growing all neighboring clusters
+                        stopped_clusters[unique_ns] = True
+                    else:
+                        # Merge multiple un-stopped trapped clusters
+                        new_num = unique_ns[0]
+                        clusters[pore] = new_num
+                        for c in unique_ns:
+                            clusters[clusters == c] = new_num
+                            msg = (seq_pore + " C:5 merge clusters: " +
+                                   str(c) + " into "+str(new_num))
+                            logger.info(msg)
+
+        # And now return clusters
+        self['pore.clusters'] = clusters
+        logger.info("Number of trapped clusters" +
+                    str(np.sum(np.unique(clusters) >= 0)))
+        self['pore.trapped'] = self['pore.clusters'] > -1
+        trapped_ts = self._net.find_neighbor_throats(self['pore.trapped'])
+        self['throat.trapped'] = np.zeros([self._net.Nt], dtype=bool)
+        self['throat.trapped'][trapped_ts] = True
+        self['pore.invasion_sequence'][self['pore.trapped']] = np.inf
+        self['throat.invasion_sequence'][self['throat.trapped']] = np.inf

OpenPNM/Geometry/models/throat_diameter.py

@@ -127,6 +127,9 @@ def minpore(network, geometry, factor=1, **kwargs):
     gTs = geometry.throats()
     nTs = geometry.map_throats(network, gTs)
     pDs = network["pore.diameter"][network["throat.conns"][nTs]]
+    if _sp.any(_sp.isnan(pDs)):
+        pDs[_sp.isnan(pDs)] = 0.0
     value = _sp.amin(pDs, axis=1)*factor
     value[value == 0.0] = _sp.amax(pDs, axis=1)[value == 0.0]
+
     return value

OpenPNM/Geometry/models/throat_length.py

@@ -29,7 +29,9 @@ def straight(network, geometry, pore_diameter='pore.diameter',
     C2 = network['pore.coords'][pore2]
     E = _sp.sqrt(_sp.sum((C1-C2)**2, axis=1))  # Euclidean distance between pores
     D1 = network[pore_diameter][pore1]
+    D1[_sp.isnan(D1)] = 0.0
     D2 = network[pore_diameter][pore2]
+    D2[_sp.isnan(D2)] = 0.0
     value = E-(D1+D2)/2.
     value = value[throats]
     if _sp.any(value < 0) and L_negative is not None:

OpenPNM/Utilities/vertexops.py

@@ -286,36 +286,6 @@ def pore2centroid(network):
                     geometry['pore.centroid'][geom_pores[i]]
 
 
-def tortuosity(network=None):
-    r"""
-    Calculate the tortuosity from the angle between throat vectors and principle axes
-    """
-    conns = network['throat.conns']
-    va = network['throat.centroid'] - network['pore.centroid'][conns[:, 0]]
-    vb = network['throat.centroid'] - network['pore.centroid'][conns[:, 1]]
-    x = [1, 0, 0]
-    y = [0, 1, 0]
-    z = [0, 0, 1]
-    f = 180 / np.pi
-    theta_x_a = tr.angle_between_vectors(va, x, directed=False, axis=1)
-    theta_x_b = tr.angle_between_vectors(vb, x, directed=False, axis=1)
-    theta_x = (np.mean(theta_x_a[~np.isnan(theta_x_a)]) +
-               np.mean(theta_x_b[~np.isnan(theta_x_b)])) / 2
-    theta_y_a = tr.angle_between_vectors(va, y, directed=False, axis=1)
-    theta_y_b = tr.angle_between_vectors(vb, y, directed=False, axis=1)
-    theta_y = (np.mean(theta_y_a[~np.isnan(theta_y_a)]) +
-               np.mean(theta_y_b[~np.isnan(theta_y_b)])) / 2
-    theta_z_a = tr.angle_between_vectors(va, z, directed=False, axis=1)
-    theta_z_b = tr.angle_between_vectors(vb, z, directed=False, axis=1)
-    theta_z = (np.mean(theta_z_a[~np.isnan(theta_z_a)]) +
-               np.mean(theta_z_b[~np.isnan(theta_z_b)])) / 2
-    tot_angle = (theta_x+theta_y+theta_z)*f
-    if 180 < tot_angle:
-        logger.error('Something is wrong: ' + str(tot_angle))
-
-    return 1 / np.cos(np.array([theta_x, theta_y, theta_z]))
-
-
 def plot_throat(geometry, throats, fig=None):
     r"""
     Print a given throat or list of throats accepted as [1, 2, 3, ..., n]

docs/userguide/tutorials/advanced_tutorial.rst

@@ -203,7 +203,7 @@ Pore-scale models are written as basic function definitions:
     ...     Dt = network['throat.diameter']
     ...     theta = phase['throat.contact_angle']
     ...     sigma = phase['throat.surface_tension']
-    ...     Pc = -4*sigma*sp.cos(f*sp.deg2rad(theta))/Dt
+    ...     Pc = 4*sigma*sp.cos(f*sp.deg2rad(theta))/Dt
     ...     return Pc[network.throats(physics.name)]
 
 Let's examine the components of above code:

test/unit/Algorithms/InvasionPercolationTest.py

@@ -1,3 +1,76 @@
+import OpenPNM as op
+import numpy as np
+
+
 class InvasionPercolationTest:
-    def test_return_results(self):
-        pass
+
+    def setup_class(self):
+        # Create Topological Network object
+        self.wrk = op.Base.Workspace()
+        self.wrk.loglevel = 50
+        N = 10
+        self.net = op.Network.Cubic(shape=[N, N, 1], spacing=1, name='net')
+        self.net.add_boundaries()
+        self.net.trim(self.net.pores('bottom_boundary'))
+        self.net.trim(self.net.pores('top_boundary'))
+        # Create Geometry object for internal pores
+        Ps = self.net.pores('boundary', mode='not')
+        Ts = self.net.find_neighbor_throats(pores=Ps,
+                                            mode='intersection',
+                                            flatten=True)
+        self.geom = op.Geometry.Stick_and_Ball(network=self.net,
+                                               pores=Ps,
+                                               throats=Ts)
+        # Create Geometry object for boundary pores
+        Ps = self.net.pores('boundary')
+        Ts = self.net.find_neighbor_throats(pores=Ps, mode='not_intersection')
+        self.boun = op.Geometry.Boundary(network=self.net, pores=Ps, throats=Ts)
+        self.phase = op.Phases.Water(network=self.net, name='water')
+        # Create one Physics object for each phase
+        Ps = self.net.pores()
+        Ts = self.net.throats()
+        self.phys = op.Physics.Standard(network=self.net,
+                                        phase=self.phase,
+                                        pores=Ps,
+                                        throats=Ts)
+        # 2. Perform Invasion Percolation
+        step = 1
+        inlets = self.net.pores('front_boundary')
+        ip_inlets = [inlets[x] for x in range(0, len(inlets), step)]
+        self.alg = op.Algorithms.InvasionPercolation(network=self.net)
+        self.alg.run(phase=self.phase, inlets=ip_inlets)
+        self.alg.return_results()
+
+    def _trapping_slow(self, outlets):
+        r"""
+        Implementation of the standard OP trapping logic for every
+        invasion step to benchmark speed
+        """
+        alg = self.alg
+        alg['pore.trapped_slow'] = np.ones([alg.Np, ], dtype=float)*-1
+        for seq in np.sort(alg['pore.invasion_sequence']):
+            invader = alg['pore.invasion_sequence'] <= seq
+            defender = ~invader.copy()
+            clusters = alg._net.find_clusters2(defender)
+            out_clusters = np.unique(clusters[outlets])
+            trapped_pores = ~np.in1d(clusters, out_clusters)
+            trapped_pores[invader] = False
+            if np.sum(trapped_pores) > 0:
+                inds = (alg['pore.trapped_slow'] == -1) * trapped_pores
+                if np.sum(inds) > 0:
+                    alg['pore.trapped_slow'][inds] = seq
+
+    def test_apply_trapping(self):
+        import time
+        t1 = time.time()
+        outlets = self.net.pores(['back_boundary',
+                                  'left_boundary',
+                                  'right_boundary'])
+        self.alg.apply_trapping(outlets)
+        t2 = time.time()
+        self._trapping_slow(outlets)
+        t3 = time.time()
+        assert (t2-t1) < (t3-t2)
+        bulk = self.net.pores('boundary', mode='not')
+        assert np.allclose(self.alg['pore.trapped'][bulk],
+                           (self.alg['pore.trapped_slow'] != -1)[bulk])

test/unit/Utilities/VertexOpsTest.py

@@ -36,10 +36,6 @@ def test_porosity(self):
         por = vo.porosity(self.net)
         assert por < 1.0
 
-    def test_tortuosity(self):
-        tor = vo.tortuosity(self.net)
-        assert sp.all(tor > 1.0)
-
     def test_pore2centroid(self):
         temp_coords = self.net['pore.coords']
         self.geo['pore.centroid'] = sp.ones([self.geo.num_pores(), 3])
@@ -48,16 +44,14 @@ def test_pore2centroid(self):
                       sp.ones([self.geo.num_pores(), 3])) == 0.0
         self.net['pore.coords'] = temp_coords
 
-    # def test_plot_pore(self):
-    #    vo.plot_pore(self.geo, self.geo.pores())
-
-    # def test_plot_throat(self):
-    #    vo.plot_throat(self.geo, [1, 2, 3, 4])
-
     def test_rotate_and_chop(self):
         throat_verts = self.geo["throat.vertices"][0]
         throat_normal = self.geo["throat.normal"][0]
         test = vo.rotate_and_chop(throat_verts, throat_normal, [0, 1, 0])
         r, c = sp.shape(test)
         assert r == len(throat_verts)
         assert c == 2
+
+if __name__ == '__main__':
+    a = VertexOpsTest()
+    a.setup_class()