Repository: PMEAL/OpenPNM
Pull request #407 (2015-06-25 05:20:49)

Ip op cleanup
This PR tries to cleanup the internals of IP and OP, as well as organizing their api a bit better.  Attempts will be made to accomodate code that uses the old api still.  Also, have added integrations tests for OP and IP using both the old and new way.


Patch:
diff --git a/OpenPNM/Algorithms/__InvasionPercolation__.py b/OpenPNM/Algorithms/__InvasionPercolation__.py
index 1c4604e41..ca9089ace 100644
--- a/OpenPNM/Algorithms/__InvasionPercolation__.py
+++ b/OpenPNM/Algorithms/__InvasionPercolation__.py
@@ -5,10 +5,8 @@ InvasionPercolationBasic: Simple IP
 ===============================================================================
 
 """
-
+import heapq as hq
 import scipy as sp
-import bisect
-from collections import deque
 from OpenPNM.Algorithms import GenericAlgorithm
 from OpenPNM.Base import logging
 logger = logging.getLogger(__name__)
@@ -32,9 +30,9 @@ class InvasionPercolation(GenericAlgorithm):
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
 
-    def run(self, phase, inlets, throat_prop='throat.capillary_pressure'):
+    def setup(self, phase, throat_prop='throat.capillary_pressure', **kwargs):
         r"""
-        Perform the algorithm
+        Set up the required parameters for the algorithm
 
         Parameters
         ----------
@@ -42,52 +40,85 @@ class InvasionPercolation(GenericAlgorithm):
             The phase to be injected into the Network.  The Phase must have the
             capillary entry pressure values for the system.
 
-        inlets : array_like
-            The list of inlet pores from which the Phase can enter the Network
-
         throat_prop : string
             The name of the throat property containing the capillary entry
             pressure.  The default is 'throat.capillary_pressure'.
 
         """
-        import heapq as hq
-        queue = []
-        hq.heapify(queue)
         self._phase = phase
-        net = self._net
         # Setup arrays and info
-        t_entry = phase[throat_prop]
+        self['throat.entry_pressure'] = phase[throat_prop]
         # Indices into t_entry giving a sorted list
-        t_sorted = sp.argsort(t_entry, axis=0)
-        t_order = sp.zeros_like(t_sorted)
-        t_order[t_sorted] = sp.arange(0, net.Nt)
-        t_inv = -sp.ones_like(net.Ts)
-        p_inv = -sp.ones_like(net.Ps)
-        p_inv[inlets] = 0
+        self['throat.sorted'] = sp.argsort(self['throat.entry_pressure'], axis=0)
+        self['throat.order'] = sp.zeros_like(self['throat.sorted'])
+        self['throat.order'][self['throat.sorted']] = sp.arange(0, self._net.Nt)
+        self['throat.invaded'] = -sp.ones((self._net.Nt,))
+        self['pore.invaded'] = -sp.ones((self._net.Np,))
+        self._tcount = 0
+
+    def set_inlets(self, pores=None, **kwargs):
+        r"""
+
+        Parameters
+        ----------
+        pores : array_like
+            The list of inlet pores from which the Phase can enter the Network
+        """
+        if 'inlets' in kwargs.keys():
+            pores = kwargs['inlets']
+        self['pore.invaded'][pores] = 0
+
         # Perform initial analysis on input pores
-        Ts = net.find_neighbor_throats(pores=inlets)
-        [hq.heappush(queue, T) for T in t_order[Ts]]
-        tcount = 1
-        while len(queue) > 0:
+        Ts = self._net.find_neighbor_throats(pores=pores)
+        self.queue = []
+        [hq.heappush(self.queue, T) for T in self['throat.order'][Ts]]
+
+    def run(self, n_steps=None, **kwargs):
+        r"""
+        Perform the algorithm
+
+        Parameters
+        ----------
+        n_steps : int
+            The number of throats to invaded during this step
+
+        """
+        if 'throat.entry_pressure' not in self.keys():
+            self.setup(**kwargs)
+        if sp.all(self['pore.invaded'] == -1):
+            self.set_inlets(**kwargs)
+
+        if n_steps is None:
+            n_steps = sp.inf
+
+        queue = self.queue
+        t_sorted = self['throat.sorted']
+        t_order = self['throat.order']
+        t_inv = self['throat.invaded']
+        p_inv = self['pore.invaded']
+
+        count = 0
+        while (len(queue) > 0) and (count < n_steps):
             # Find throat at the top of the queue
             t = hq.heappop(queue)
             # Extract actual throat number
             t_next = t_sorted[t]
-            t_inv[t_next] = tcount
+            t_inv[t_next] = self._tcount
             # If throat is duplicated
             while len(queue) > 0 and queue[0] == t:
                 # Note: Preventing duplicate entries below might save some time here
                 t = hq.heappop(queue)
             # Find pores connected to newly invaded throat
-            Ps = net['throat.conns'][t_next]
+            Ps = self._net['throat.conns'][t_next]
             # Remove already invaded pores from Ps
             Ps = Ps[p_inv[Ps] < 0]
             if len(Ps) > 0:
-                p_inv[Ps] = tcount
-                Ts = net.find_neighbor_throats(pores=Ps)  # Find connected throats
+                p_inv[Ps] = self._tcount
+                Ts = self._net.find_neighbor_throats(pores=Ps)
                 Ts = Ts[t_inv[Ts] < 0]  # Remove already invaded throats from Ts
                 [hq.heappush(queue, T) for T in t_order[Ts]]
-            tcount += 1
+            count += 1
+            self._tcount += 1
         self['throat.invasion_sequence'] = t_inv
         self['pore.invasion_sequence'] = p_inv
 
diff --git a/OpenPNM/Algorithms/__OrdinaryPercolation__.py b/OpenPNM/Algorithms/__OrdinaryPercolation__.py
index aa42bd682..d34e782eb 100644
--- a/OpenPNM/Algorithms/__OrdinaryPercolation__.py
+++ b/OpenPNM/Algorithms/__OrdinaryPercolation__.py
@@ -16,7 +16,7 @@ logger = logging.getLogger(__name__)
 
 class OrdinaryPercolation(GenericAlgorithm):
     r"""
-    Simulates a capillary drainage experiment by looping through a list of
+    Simulates a capillary drainage experiment by applying a list of increasing
     capillary pressures.
 
     Parameters
@@ -24,227 +24,215 @@ class OrdinaryPercolation(GenericAlgorithm):
     network : OpenPNM Network Object
         The network upon which the simulation will be run
 
-    invading_phase : OpenPNM Phase Object
-        The phase to be forced into the network at increasingly high pressures
-
-    defending_phase : OpenPNM Phase Object, optional
-        The phase originally residing in the network prior to invasion.  This
-        is only necessary so that the pressure at which the phase is drained
-        can be stored on the phase.
-
     name : string, optional
         The name to assign to the Algorithm Object
 
-    Examples
-    --------
-    >>> import OpenPNM
-    >>> pn = OpenPNM.Network.TestNet()
-    >>> geo = OpenPNM.Geometry.TestGeometry(network=pn, pores=pn.pores(),
-    ...                                     throats=pn.throats())
-    >>> phase1 = OpenPNM.Phases.TestPhase(network=pn)
-    >>> phase2 = OpenPNM.Phases.TestPhase(network=pn)
-    >>> phys1 = OpenPNM.Physics.TestPhysics(network=pn,
-    ...                                     phase=phase1,
-    ...                                     pores=pn.pores(),
-    ...                                     throats=pn.throats())
-    >>> phys2 = OpenPNM.Physics.TestPhysics(network=pn, phase=phase2,
-    ...                                     pores=pn.pores(), throats=pn.throats())
-    >>> OP = OpenPNM.Algorithms.OrdinaryPercolation(network=pn,
-    ...                                             invading_phase=phase1,
-    ...                                             defending_phase=phase2)
-    >>> OP.run(inlets=pn.pores('top'))
-    >>> med_Pc = sp.median(OP['pore.inv_Pc'])
-    >>> OP.return_results(med_Pc)
-    >>> print(len(phase1.pores('occupancy')))
-    71
-
-    To run this algorithm, use 'setup()' to provide the necessary simulation
     """
 
-    def __init__(self, invading_phase=None, defending_phase=None,
-                 residual_pores=None, residual_throats=None, **kwargs):
-        super().__init__(**kwargs)
-        self._phase_inv = invading_phase
-        self._phase_def = defending_phase
-        self._residual_pores = residual_pores
-        self._residual_throats = residual_throats
-
-        logger.debug('Create Drainage Percolation Algorithm Object')
+    def __init__(self, network, name=None, **kwargs):
+        super().__init__(network=network, name=name)
+        if len(kwargs.keys()) > 0:
+            self.setup(**kwargs)
 
-    def run(self, inlets, outlets=[], npts=25, inv_points=None,
-            capillary_pressure='capillary_pressure', access_limited=True,
-            trapping=False, **kwargs):
+    def setup(self,
+              invading_phase,
+              defending_phase=None,
+              t_entry='throat.capillary_pressure',
+              **kwargs):
         r"""
-        Parameters
-        ----------
-        inlets : array_like
-            The list of pores which are the injection sources
+        invading_phase : OpenPNM Phase Object
+            The invading phase to be injected into the Network
 
-        npts : int, optional
-            The number of pressure points to apply.  The list of pressures
-            is logarithmically spaced between the lowest and highest throat
-            entry pressures in the network.
+        p_inlets : array_like
+            The injection points from which the invading phase accesses the
+            Network.  If no inlets are specified then the algorithm assumes
+            no access limitations apply to the invading phase, which is
+            equivalent to performaing a standard bond ordinary percolation.
 
-        inv_points : array_like, optional
-            A list of specific pressure points to apply.
 
-        access_limited : boolean
-            Only pores and throats connected to the inlet sites can be invaded
+        Notes
+        -----
+        The 'inlet' pores are initially filled with invading fluid to start the
+        simulation.  To avoid the capillary pressure curve showing a non-zero
+        starting saturation at low pressures, it is necessary to apply boundary
+        pores that have zero-volume, and set these as the inlets.
+        """
+        self['throat.entry_pressure'] = invading_phase[t_entry]
+        self['pore.inv_Pc'] = sp.inf
+        self['throat.inv_Pc'] = sp.inf
+        self['pore.inv_sat'] = sp.inf
+        self['throat.inv_sat'] = sp.inf
+        self._inv_phase = invading_phase
+        self._def_phase = defending_phase
+        self._trapping = False
+
+    def set_inlets(self, pores):
+        r"""
+        Specify inlet locations
+
+        Parameters
+        ----------
+        pores : array_like
+            The injection points from which the invading phase accesses the
+            Network.  If no inlets are specified then the algorithm assumes
+            no access limitations apply to the invading phase, which is
+            equivalent to performaing a standard bond ordinary percolation.
 
-        trapping : boolean
-            Wetting phase that is cut-off from the outlets becomes immobile.
-            If outlet pores have not been provided then this argument is
-            ignored.
 
         Notes
         -----
         The 'inlet' pores are initially filled with invading fluid to start the
         simulation.  To avoid the capillary pressure curve showing a non-zero
         starting saturation at low pressures, it is necessary to apply boundary
-        pores that have 0 volume, and set these as the inlets.
+        pores that have zero-volume, and set these as the inlets.
+        """
+        Ps = sp.array(pores)
+        if sp.size(Ps) > 0:
+            if Ps.dtype == bool:
+                Ps = self._net.Ps[Ps]
+            self['pore.inlets'] = False
+            self['pore.inlets'][Ps] = True
+
+    def set_outlets(self, pores, defending_phase=None):
+        r"""
+        Specify outlet locations
+
+        Parameters
+        ----------
+        pores : array_like
+            The pores through which the defending phase exits the Network.
 
+        defending_phase : OpenPNM Phase Object
+            The Phase object defining the defending phase.  The defending Phase
+            may be specified during the ``setup`` step, or through this method.
         """
-        self._inv_sites = inlets
-        self._out_sites = outlets
-        self._npts = npts
-        self._p_cap = capillary_pressure  # Name of throat entry pressure prop
+        if defending_phase is not None:
+            self._def_phase = defending_phase
+
+        self._trapping = True
+
+        Ps = sp.array(pores)
+        if sp.size(Ps) > 0:
+            if Ps.dtype == bool:
+                Ps = self._net.Ps[Ps]
+            self['pore.outlets'] = False
+            self['pore.outlets'][Ps] = True
+
+    def run(self, npts=25, inv_points=None, access_limited=True, **kwargs):
+        r"""
+        Parameters
+        ----------
+        npts : int (default = 25)
+            The number of pressure points to apply.  The list of pressures
+            is logarithmically spaced between the lowest and highest throat
+            entry pressures in the network.
+
+        inv_points : array_like, optional
+            A list of specific pressure point(s) to apply.
+
+        """
+        if 'inlets' in kwargs.keys():
+            print('Inlets recieved, passing to set_inlets')
+            self.set_inlets(pores=kwargs['inlets'])
+        if 'outlets' in kwargs.keys():
+            print('Outlets recieved, passing to set_outlets')
+            self.set_outlets(pores=kwargs['outlets'])
         self._AL = access_limited
-        self._TR = trapping
-
-        # Create pore and throat conditions lists to store inv_val at
-        # which each is invaded
-        self._p_inv = sp.zeros((self._net.num_pores(),), dtype=float)
-        self._p_inv.fill(sp.inf)
-        self._p_seq = sp.zeros_like(self._p_inv, dtype=int)
-        self._t_inv = sp.zeros((self._net.num_throats(),), dtype=float)
-        self._t_inv.fill(sp.inf)
-        self._t_seq = sp.zeros_like(self._t_inv, dtype=int)
-        # Determine the invasion pressures to apply
-        try:
-            self._t_cap = self._phase_inv['throat.' + self._p_cap]
-        except:
-            logger.error('Capillary pressure not assigned to invading phase ' +
-                         self._phase_inv.name +
-                         ', check for capillary pressure in defending phase ' +
-                         self._phase_def.name + ' instead')
-            try:
-                self._t_cap = self._phase_def['throat.' + self._p_cap]
-            except:
-                pass
-                logger.error('Capillary pressure neither assigned to defending \
-                              phase ' + self._phase_def.name +
-                             ' nor to invading phase ' + self._phase_inv.name)
         if inv_points is None:
-            min_p = sp.amin(self._t_cap) * 0.98  # nudge min_p down slightly
-            max_p = sp.amax(self._t_cap) * 1.02  # bump max_p up slightly
             logger.info('Generating list of invasion pressures')
-            if min_p == 0:
-                min_p = sp.linspace(min_p, max_p, self._npts)[1]
-            self._inv_points = sp.logspace(sp.log10(min_p),
-                                           sp.log10(max_p),
-                                           self._npts)
-        else:
-            self._inv_points = inv_points
-        self._do_outer_iteration_stage()
+            min_p = sp.amin(self['throat.entry_pressure']) * 0.98  # nudge down
+            max_p = sp.amax(self['throat.entry_pressure']) * 1.02  # bump up
+            inv_points = sp.logspace(sp.log10(min_p),
+                                     sp.log10(max_p),
+                                     npts)
+
+        self._npts = sp.size(inv_points)
+        # Execute calculation
+        self._do_outer_iteration_stage(inv_points)
 
-    def _do_outer_iteration_stage(self):
+    def _do_outer_iteration_stage(self, inv_points):
         # Generate curve from points
-        for inv_val in self._inv_points:
+        for inv_val in inv_points:
             # Apply one applied pressure and determine invaded pores
             logger.info('Applying capillary pressure: ' + str(inv_val))
             self._do_one_inner_iteration(inv_val)
-        # Store results using networks' get/set method
-        self['pore.inv_Pc'] = self._p_inv
-        self['throat.inv_Pc'] = self._t_inv
+
         # Find invasion sequence values (to correspond with IP algorithm)
-        self._p_seq = sp.searchsorted(sp.unique(self._p_inv), self._p_inv)
-        self._t_seq = sp.searchsorted(sp.unique(self._t_inv), self._t_inv)
-        self['pore.inv_seq'] = self._p_seq
-        self['throat.inv_seq'] = self._t_seq
-        # Calculate Saturations
-        v_total = sp.sum(self._net['pore.volume']) + \
-            sp.sum(self._net['throat.volume'])
-        sat = 0.
-        self['pore.inv_sat'] = 1.
-        self['throat.inv_sat'] = 1.
-        for i in range(self._npts):
-            inv_pores = sp.where(self._p_seq == i)[0]
-            inv_throats = sp.where(self._t_seq == i)[0]
-            new_sat = (sum(self._net['pore.volume'][inv_pores]) +
-                       sum(self._net['throat.volume'][inv_throats])) / v_total
-            sat += new_sat
-            self['pore.inv_sat'][inv_pores] = sat
-            self['throat.inv_sat'][inv_throats] = sat
-        if self._TR:
-            logger.info('Evaluating trapping')
-            self.evaluate_trapping(outlets=self._out_sites)
+        self['pore.inv_seq'] = sp.searchsorted(sp.unique(self['pore.inv_Pc']),
+                                               self['pore.inv_Pc'])
+        self['throat.inv_seq'] = sp.searchsorted(sp.unique(self['throat.inv_Pc']),
+                                                 self['throat.inv_Pc'])
+
+        if self._trapping:
+            self.evaluate_trapping(self['pore.outlets'])
 
     def _do_one_inner_iteration(self, inv_val):
         r"""
-        Determine which throats are invaded at a given applied capillary pressure.
-
-        This function uses the scipy.csgraph module for the cluster labeling
-        algorithm (connected_components).
+        Determine which throats are invaded at a given applied capillary
+        pressure.
 
         """
         # Generate a tlist containing boolean values for throat state
-        Tinvaded = self._t_cap <= inv_val
-        # Finding all pores that can be invaded at specified pressure
-        clusters = self._net.find_clusters(Tinvaded)
-        # Find all pores with at least 1 invaded throat (invaded)
-        Pinvaded = sp.zeros_like(clusters, dtype=bool)
-        Ts = self._net.throats()
-        P12 = self._net.find_connected_pores(Ts)
-        temp = P12[Tinvaded]
-        temp = sp.hstack((temp[:, 0], temp[:, 1]))
-        Pinvaded[temp] = True
+        Tinvaded = self['throat.entry_pressure'] <= inv_val
+        # Find all pores that can be invaded at specified pressure
+        [pclusters, tclusters] = self._net.find_clusters2(mask=Tinvaded,
+                                                          t_labels=True)
         if self._AL:
-            # Add injection sites to Pinvaded
-            Pinvaded[self._inv_sites] = True
-            # Clean up clusters (not invaded = -1, invaded >=0)
-            clusters = clusters * (Pinvaded) - (~Pinvaded)
             # Identify clusters connected to invasion sites
-            inv_clusters = sp.unique(clusters[self._inv_sites])
+            inv_clusters = sp.unique(pclusters[self['pore.inlets']])
+            inv_clusters = inv_clusters[inv_clusters >= 0]
         else:
-            # Clean up clusters (not invaded = -1, invaded >=0)
-            clusters = clusters * (Pinvaded) - (~Pinvaded)
             # All clusters are invasion sites
-            inv_clusters = sp.r_[0:self._net.num_pores()]
-        # Store invasion pressure in pores and throats
-        pmask = np.in1d(clusters, inv_clusters)
-        # Store result of invasion step
-        self._p_inv[(self._p_inv == sp.inf) * (pmask)] = inv_val
-        # Determine Pc_invaded for throats as well
-        temp = self._net['throat.conns']
-        tmask = (pmask[temp[:, 0]] + pmask[temp[:, 1]]) * Tinvaded
-        self._t_inv[(self._t_inv == sp.inf) * (tmask)] = inv_val
-
-    def evaluate_trapping(self, outlets):
+            inv_clusters = pclusters
+        # Find pores on the invading clusters
+        pmask = np.in1d(pclusters, inv_clusters)
+        # Store current applied pressure in newly invaded pores
+        pinds = (self['pore.inv_Pc'] == sp.inf) * (pmask)
+        self['pore.inv_Pc'][pinds] = inv_val
+        # Find throats on the invading clusters
+        tmask = np.in1d(tclusters, inv_clusters)
+        # Store current applied pressure in newly invaded throats
+        tinds = (self['throat.inv_Pc'] == sp.inf) * (tmask)
+        self['throat.inv_Pc'][tinds] = inv_val
+        # Store total network saturation
+        tsat = sp.sum(self._net['throat.volume'][self['throat.inv_Pc'] <= inv_val])
+        psat = sp.sum(self._net['pore.volume'][self['pore.inv_Pc'] <= inv_val])
+        total = sp.sum(self._net['throat.volume']) + sp.sum(self._net['pore.volume'])
+        self['pore.inv_sat'][pinds] = (tsat + psat)/total
+        self['throat.inv_sat'][tinds] = (tsat + psat)/total
+
+    def evaluate_trapping(self, p_outlets):
         r"""
         Finds trapped pores and throats after a full ordinary
-        percolation drainage has been run
+        percolation simulation has been run.
 
         Parameters
         ----------
-        outlets : array_like
+        p_outlets : array_like
             A list of pores that define the wetting phase outlets.
             Disconnection from these outlets results in trapping.
 
+        Returns
+        -------
+        It creates arrays called ``pore.trapped`` and ``throat.trapped``, but
+        also adjusts the ``pore.inv_Pc`` and ``throat.inv_Pc`` arrays to set
+        trapped locations to have infinite invasion pressure.
+
         """
-        self._p_trap = sp.zeros_like(self._p_inv, dtype=float)
-        self._t_trap = sp.zeros_like(self._t_inv, dtype=float)
+        self['pore.trapped'] = sp.zeros([self.Np, ], dtype=float)
+        self['throat.trapped'] = sp.zeros([self.Nt, ], dtype=float)
         try:
-            inv_points = sp.unique(self._p_inv)  # Get points used in OP
+            # Get points used in OP
+            inv_points = sp.unique(self['pore.inv_Pc'])
         except:
-            logger.error('Orindary percolation has not been run!')
-            raise Exception('Aborting algorithm')
+            raise Exception('Orindary percolation has not been run!')
         tind = self._net.throats()
         conns = self._net.find_connected_pores(tind)
         for inv_val in inv_points[0:-1]:
             # Find clusters of defender pores
-            Pinvaded = self._p_inv <= inv_val
+            Pinvaded = self['pore.inv_Pc'] <= inv_val
             Cstate = sp.sum(Pinvaded[conns], axis=1)
-            Tinvaded = self._t_inv <= inv_val
+            Tinvaded = self['throat.inv_Pc'] <= inv_val
             # 0 = all open, 1=1 pore filled,
             # 2=2 pores filled 3=2 pores + 1 throat filled
             Cstate = Cstate + Tinvaded
@@ -252,20 +240,23 @@ class OrdinaryPercolation(GenericAlgorithm):
             # Clean up clusters (invaded = -1, defended >=0)
             clusters = clusters * (~Pinvaded) - (Pinvaded)
             # Identify clusters connected to outlet sites
-            out_clusters = sp.unique(clusters[outlets])
+            out_clusters = sp.unique(clusters[p_outlets])
             trapped_pores = ~sp.in1d(clusters, out_clusters)
             trapped_pores[Pinvaded] = False
             if sum(trapped_pores) > 0:
-                self._p_trap[(self._p_trap == 0) * trapped_pores] = inv_val
+                inds = (self['pore.trapped'] == 0) * trapped_pores
+                self['pore.trapped'][inds] = inv_val
                 trapped_throats = self._net.find_neighbor_throats(trapped_pores)
                 trapped_throat_array = np.asarray([False] * len(Cstate))
                 trapped_throat_array[trapped_throats] = True
-                self._t_trap[(self._t_trap == 0) * trapped_throat_array] = inv_val
-                self._t_trap[(self._t_trap == 0) * (Cstate == 2)] = inv_val
-        self._p_inv[self._p_trap > 0] = sp.inf
-        self._t_inv[self._t_trap > 0] = sp.inf
-        self['pore.inv_Pc'] = self._p_inv
-        self['throat.inv_Pc'] = self._t_inv
+                inds = (self['throat.trapped'] == 0) * trapped_throat_array
+                self['throat.trapped'][inds] = inv_val
+                inds = (self['throat.trapped'] == 0) * (Cstate == 2)
+                self['throat.trapped'][inds] = inv_val
+        self['pore.trapped'][self['pore.trapped'] > 0] = sp.inf
+        self['throat.trapped'][self['throat.trapped'] > 0] = sp.inf
+        self['pore.inv_Pc'][self['pore.trapped'] > 0] = sp.inf
+        self['throat.inv_Pc'][self['throat.trapped'] > 0] = sp.inf
 
     def return_results(self, Pc=0, seq=None, sat=None, occupancy='occupancy'):
         r"""
@@ -274,60 +265,60 @@ class OrdinaryPercolation(GenericAlgorithm):
 
         """
         p_inv = self['pore.inv_Pc']
-        self._phase_inv['pore.inv_Pc'] = p_inv
+        self._inv_phase['pore.inv_Pc'] = p_inv
         t_inv = self['throat.inv_Pc']
-        self._phase_inv['throat.inv_Pc'] = t_inv
+        self._inv_phase['throat.inv_Pc'] = t_inv
         # Apply invasion sequence values (to correspond with IP algorithm)
         p_seq = self['pore.inv_seq']
-        self._phase_inv['pore.inv_seq'] = p_seq
+        self._inv_phase['pore.inv_seq'] = p_seq
         t_seq = self['throat.inv_seq']
-        self._phase_inv['throat.inv_seq'] = t_seq
+        self._inv_phase['throat.inv_seq'] = t_seq
         # Apply saturation to pores and throats
-        self._phase_inv['pore.inv_sat'] = self['pore.inv_sat']
-        self._phase_inv['throat.inv_sat'] = self['throat.inv_sat']
+        self._inv_phase['pore.inv_sat'] = self['pore.inv_sat']
+        self._inv_phase['throat.inv_sat'] = self['throat.inv_sat']
 
         if sat is not None:
             p_inv = self['pore.inv_sat'] <= sat
             t_inv = self['throat.inv_sat'] <= sat
             # Apply occupancy to invading phase
             temp = sp.array(p_inv, dtype=sp.float_, ndmin=1)
-            self._phase_inv['pore.' + occupancy] = temp
+            self._inv_phase['pore.' + occupancy] = temp
             temp = sp.array(t_inv, dtype=sp.float_, ndmin=1)
-            self._phase_inv['throat.' + occupancy] = temp
+            self._inv_phase['throat.' + occupancy] = temp
             # Apply occupancy to defending phase
-            if self._phase_def is not None:
+            if self._def_phase is not None:
                 temp = sp.array(~p_inv, dtype=sp.float_, ndmin=1)
-                self._phase_def['pore.' + occupancy] = temp
+                self._def_phase['pore.' + occupancy] = temp
                 temp = sp.array(~t_inv, dtype=sp.float_, ndmin=1)
-                self._phase_def['throat.' + occupancy] = temp
+                self._def_phase['throat.' + occupancy] = temp
         elif seq is not None:
             p_seq = self['pore.inv_seq'] <= seq
             t_seq = self['throat.inv_seq'] <= seq
             # Apply occupancy to invading phase
             temp = sp.array(p_seq, dtype=sp.float_, ndmin=1)
-            self._phase_inv['pore.' + occupancy] = temp
+            self._inv_phase['pore.' + occupancy] = temp
             temp = sp.array(t_seq, dtype=sp.float_, ndmin=1)
-            self._phase_inv['throat.' + occupancy] = temp
+            self._inv_phase['throat.' + occupancy] = temp
             # Apply occupancy to defending phase
-            if self._phase_def is not None:
+            if self._def_phase is not None:
                 temp = sp.array(~p_seq, dtype=sp.float_, ndmin=1)
-                self._phase_def['pore.' + occupancy] = temp
+                self._def_phase['pore.' + occupancy] = temp
                 temp = sp.array(~t_seq, dtype=sp.float_, ndmin=1)
-                self._phase_def['throat.' + occupancy] = temp
+                self._def_phase['throat.' + occupancy] = temp
         else:
             p_inv = self['pore.inv_Pc'] <= Pc
             t_inv = self['throat.inv_Pc'] <= Pc
             # Apply occupancy to invading phase
             temp = sp.array(p_inv, dtype=sp.float_, ndmin=1)
-            self._phase_inv['pore.' + occupancy] = temp
+            self._inv_phase['pore.' + occupancy] = temp
             temp = sp.array(t_inv, dtype=sp.float_, ndmin=1)
-            self._phase_inv['throat.' + occupancy] = temp
+            self._inv_phase['throat.' + occupancy] = temp
             # Apply occupancy to defending phase
-            if self._phase_def is not None:
+            if self._def_phase is not None:
                 temp = sp.array(~p_inv, dtype=sp.float_, ndmin=1)
-                self._phase_def['pore.' + occupancy] = temp
+                self._def_phase['pore.' + occupancy] = temp
                 temp = sp.array(~t_inv, dtype=sp.float_, ndmin=1)
-                self._phase_def['throat.' + occupancy] = temp
+                self._def_phase['throat.' + occupancy] = temp
 
     def plot_drainage_curve(self, pore_volume='volume', throat_volume='volume',
                             pore_label='all', throat_label='all'):
@@ -343,24 +334,30 @@ class OrdinaryPercolation(GenericAlgorithm):
         throats = self._net.throats(labels=throat_label)
         Snwp_t = sp.zeros_like(PcPoints)
         Snwp_p = sp.zeros_like(PcPoints)
+        Snwp_all = sp.zeros_like(PcPoints)
         Pvol = self._net['pore.' + pore_volume]
         Tvol = self._net['throat.' + throat_volume]
-        Pvol_tot = sum(Pvol)
-        Tvol_tot = sum(Tvol)
+        Pvol_tot = sp.sum(Pvol)
+        Tvol_tot = sp.sum(Tvol)
+        vol_tot = Pvol_tot + Tvol_tot
         for i in range(0, sp.size(PcPoints)):
             Pc = PcPoints[i]
-            Snwp_p[i] = sum(Pvol[self._p_inv[pores] <= Pc]) / Pvol_tot
-            Snwp_t[i] = sum(Tvol[self._t_inv[throats] <= Pc]) / Tvol_tot
-        if sp.mean(self._phase_inv['pore.contact_angle']) < 90:
+            Snwp_p[i] = sp.sum(Pvol[self['pore.inv_Pc'][pores] <= Pc]) / vol_tot
+            Snwp_t[i] = sp.sum(Tvol[self['throat.inv_Pc'][throats] <= Pc]) / vol_tot
+            Snwp_all[i] = (sp.sum(Tvol[self['throat.inv_Pc'][throats] <= Pc]) +
+                           sp.sum(Pvol[self['pore.inv_Pc'][pores] <= Pc])) / vol_tot
+        if sp.mean(self._inv_phase['pore.contact_angle']) < 90:
             Snwp_p = 1 - Snwp_p
             Snwp_t = 1 - Snwp_t
             PcPoints *= -1
+        fig = plt.figure()
+        plt.plot(PcPoints, Snwp_all, 'g.-')
         plt.plot(PcPoints, Snwp_p, 'r.-')
         plt.plot(PcPoints, Snwp_t, 'b.-')
         r"""
         TODO: Add legend to distinguish the pore and throat curves
         """
-        plt.show()
+        return fig
 
     def plot_primary_drainage_curve(self, pore_volume='volume',
                                     throat_volume='volume', pore_label='all',
@@ -377,26 +374,29 @@ class OrdinaryPercolation(GenericAlgorithm):
                             simulation has not been run')
         pores = self._net.pores(labels=pore_label)
         throats = self._net.throats(labels=throat_label)
+        p_inv = self['pore.inv_Pc']
+        t_inv = self['throat.inv_Pc']
         Snwp_t = sp.zeros_like(PcPoints)
         Snwp_p = sp.zeros_like(PcPoints)
         Snwp_all = sp.zeros_like(PcPoints)
         Swp_all = sp.zeros_like(PcPoints)
         Pvol = self._net['pore.' + pore_volume]
         Tvol = self._net['throat.' + throat_volume]
-        Pvol_tot = sum(Pvol)
-        Tvol_tot = sum(Tvol)
+        Pvol_tot = sp.sum(Pvol)
+        Tvol_tot = sp.sum(Tvol)
         for i in range(0, sp.size(PcPoints)):
             Pc = PcPoints[i]
-            Snwp_p[i] = sum(Pvol[self._p_inv[pores] <= Pc]) / Pvol_tot
-            Snwp_t[i] = sum(Tvol[self._t_inv[throats] <= Pc]) / Tvol_tot
-            Snwp_all[i] = (sum(Tvol[self._t_inv[throats] <= Pc]) +
-                           sum(Pvol[self._p_inv[pores] <= Pc])) / \
+            Snwp_p[i] = sp.sum(Pvol[p_inv[pores] <= Pc]) / Pvol_tot
+            Snwp_t[i] = sp.sum(Tvol[t_inv[throats] <= Pc]) / Tvol_tot
+            Snwp_all[i] = (sp.sum(Tvol[t_inv[throats] <= Pc]) +
+                           sp.sum(Pvol[p_inv[pores] <= Pc])) / \
                           (Tvol_tot + Pvol_tot)
             Swp_all[i] = 1 - Snwp_all[i]
+        fig = plt.figure()
         plt.plot(Swp_all, PcPoints, 'k.-')
         plt.xlim(xmin=0)
         plt.xlabel('Saturation of wetting phase')
         plt.ylabel('Capillary Pressure [Pa]')
         plt.title('Primay Drainage Curve')
         plt.grid(True)
-        plt.show()
+        return fig


Test patch:
diff --git a/test/integration/script_test.py b/test/integration/script_test.py
index 870e3d3ed..3747d326c 100644
--- a/test/integration/script_test.py
+++ b/test/integration/script_test.py
@@ -3,48 +3,6 @@ import scipy as sp
 import OpenPNM
 import pytest
 
-def test_cubic_standard_call():
-    pn = OpenPNM.Network.Cubic(shape=[3, 4, 5])
-    np.testing.assert_almost_equal(pn['pore.coords'][0], [0.5, 0.5, 0.5])
-
-
-def test_trim_extend():
-    pn = OpenPNM.Network.Cubic(shape=[5, 5, 5])
-    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0), [1, 5, 25]))
-    assert [pn.Np, pn.Nt] == [125, 300]
-    pn.trim(pores=[0])
-    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0), [1, 5, 25]))
-    assert [pn.Np, pn.Nt] == [124, 297]
-    pn.extend(pore_coords=[0, 0, 0], throat_conns=[[124, 0]])
-    assert [pn.Np, pn.Nt] == [125, 298]
-    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0), [1, 5, 25, 124]))
-
-
-def test_stitch():
-    ctrl = OpenPNM.Base.Controller()
-    [Nx, Ny, Nz] = [10, 10, 10]
-    pn = OpenPNM.Network.Cubic(shape=[Nx, Ny, Nz])
-    pn2 = OpenPNM.Network.Cubic(shape=[Nx, Ny, Nz])
-    pn2['pore.coords'][:, 2] += Nz
-    pn.stitch(donor=pn2,
-              P_network=pn.pores('top'),
-              P_donor=pn2.pores('bottom'),
-              len_max=1,
-              method='nearest')
-    assert pn.Np == 2*pn2.Np  # Ensure number of pores doubled
-    assert pn.Nt == (2*pn2.Nt + Nx*Ny)  # Ensure correct number of new throats
-    assert pn2 not in ctrl.values()  # Donor Network is removed from Controller
-    # Reuse the donor Network in another stitch
-    pn2['pore.coords'][:, 2] -= 2*Nz
-    pn.stitch(donor=pn2,
-              P_network=pn.pores('bottom'),
-              P_donor=pn2.pores('top'),
-              len_max=1,
-              method='nearest')
-    assert pn.Np == 3*pn2.Np  # Ensure number of pores increased again
-    assert pn.Nt == (3*pn2.Nt + 2*Nx*Ny)  # Ensure correct number of new throats
-    ctrl.clear()
-
 
 def test_linear_solvers():
     pn = OpenPNM.Network.Cubic([1, 40, 30], spacing=0.0001)
@@ -227,42 +185,6 @@ def test_thermal_conduction():
     assert sp.amax(np.absolute(diff)) < 0.015
 
 
-def test_OP_standard():
-    pn = OpenPNM.Network.Cubic(shape=[10, 10, 10], spacing=0.0001)
-    geom = OpenPNM.Geometry.Toray090(network=pn, pores=pn.Ps, throats=pn.Ts)
-    water = OpenPNM.Phases.Water(network=pn)
-    phys = OpenPNM.Physics.GenericPhysics(network=pn,
-                                          geometry=geom,
-                                          phase=water)
-    f = OpenPNM.Physics.models.capillary_pressure.washburn
-    phys.models.add(propname='throat.capillary_pressure',
-                    model=f)
-    OP = OpenPNM.Algorithms.OrdinaryPercolation(network=pn,
-                                                invading_phase=water)
-    OP.run(inlets=pn.pores('top'))
-    V_inv = sp.sum(pn['pore.volume'][OP['pore.inv_Pc'] < sp.inf])
-    V_tot = sp.sum(pn['pore.volume'])
-    assert V_inv/V_tot == 1.0
-
-
-def test_OP_trapping():
-    pn = OpenPNM.Network.Cubic(shape=[10, 10, 10], spacing=0.0001)
-    geom = OpenPNM.Geometry.Toray090(network=pn, pores=pn.Ps, throats=pn.Ts)
-    water = OpenPNM.Phases.Water(network=pn)
-    phys = OpenPNM.Physics.GenericPhysics(network=pn,
-                                          geometry=geom,
-                                          phase=water)
-    f = OpenPNM.Physics.models.capillary_pressure.washburn
-    phys.models.add(propname='throat.capillary_pressure',
-                    model=f)
-    OP = OpenPNM.Algorithms.OrdinaryPercolation(network=pn,
-                                                invading_phase=water)
-    OP.run(inlets=pn.pores('top'), outlets=pn.pores('bottom'), trapping=True)
-    V_inv = sp.sum(pn['pore.volume'][OP['pore.inv_Pc'] < sp.inf])
-    V_tot = sp.sum(pn['pore.volume'])
-    assert V_inv/V_tot < 0.95
-
-
 def test_Darcy_alg():
     # Generate Network and clean up some of boundaries
     divs = [1, 50, 10]
diff --git a/test/integration/test_percolation_algorithms.py b/test/integration/test_percolation_algorithms.py
new file mode 100644
index 000000000..e509392ae
--- /dev/null
+++ b/test/integration/test_percolation_algorithms.py
@@ -0,0 +1,106 @@
+import scipy as sp
+import OpenPNM
+import pytest
+ctrl = OpenPNM.Base.Controller()
+
+
+def test_IP_old_approach():
+    ctrl.clear()
+    pn = OpenPNM.Network.Cubic(shape=[30, 30, 1], spacing=0.01)
+    geom = OpenPNM.Geometry.Toray090(network=pn, pores=pn.Ps, throats=pn.Ts)
+    water = OpenPNM.Phases.Water(network=pn)
+    phys = OpenPNM.Physics.Standard(network=pn, phase=water, geometry=geom)
+    inlets = pn.pores('left')
+    IP_1 = OpenPNM.Algorithms.InvasionPercolation(network=pn)
+    IP_1.run(phase=water, inlets=inlets)
+    a = ['pore.all', 'pore.invaded', 'pore.invasion_sequence', 'throat.all',
+         'throat.entry_pressure', 'throat.invaded', 'throat.invasion_sequence',
+         'throat.order', 'throat.sorted']
+    assert sorted(list(IP_1.keys())) == a
+    IP_1.return_results()
+    assert 'throat.invasion_sequence' in water.keys()
+    assert 'pore.invasion_sequence' in water.keys()
+    ctrl.clear()
+
+
+def test_IP_new_approach():
+    ctrl.clear()
+    pn = OpenPNM.Network.Cubic(shape=[30, 30, 1], spacing=0.01)
+    geom = OpenPNM.Geometry.Toray090(network=pn, pores=pn.Ps, throats=pn.Ts)
+    water = OpenPNM.Phases.Water(network=pn)
+    phys = OpenPNM.Physics.Standard(network=pn, phase=water, geometry=geom)
+    inlets = pn.pores('left')
+    IP_1 = OpenPNM.Algorithms.InvasionPercolation(network=pn)
+    IP_1.setup(phase=water)
+    IP_1.set_inlets(pores=inlets)
+    IP_1.run()
+    a = ['pore.all', 'pore.invaded', 'pore.invasion_sequence', 'throat.all',
+         'throat.entry_pressure', 'throat.invaded', 'throat.invasion_sequence',
+         'throat.order', 'throat.sorted']
+    assert sorted(list(IP_1.keys())) == a
+    IP_1.return_results()
+    assert 'throat.invasion_sequence' in water.keys()
+    assert 'pore.invasion_sequence' in water.keys()
+    ctrl.clear()
+
+
+def test_OP_old_approach():
+    ctrl.clear()
+    pn = OpenPNM.Network.Cubic(shape=[30, 30, 1], spacing=0.01)
+    geom = OpenPNM.Geometry.Toray090(network=pn, pores=pn.Ps, throats=pn.Ts)
+    water = OpenPNM.Phases.Water(network=pn)
+    phys = OpenPNM.Physics.Standard(network=pn, phase=water, geometry=geom)
+    OP_1 = OpenPNM.Algorithms.OrdinaryPercolation(network=pn,
+                                                  invading_phase=water)
+    Ps = pn.pores(labels=['left'])
+    OP_1.run(inlets=Ps)
+    OP_1.return_results(Pc=7000)
+    a = ['pore.all', 'pore.inlets', 'pore.inv_Pc', 'pore.inv_sat',
+         'pore.inv_seq', 'throat.all', 'throat.entry_pressure',
+         'throat.inv_Pc', 'throat.inv_sat', 'throat.inv_seq']
+    assert sorted(list(OP_1.keys())) == a
+    ctrl.clear()
+
+
+def test_OP_new_approach():
+    ctrl.clear()
+    pn = OpenPNM.Network.Cubic(shape=[30, 30, 1], spacing=0.01)
+    geom = OpenPNM.Geometry.Toray090(network=pn, pores=pn.Ps, throats=pn.Ts)
+    water = OpenPNM.Phases.Water(network=pn)
+    phys = OpenPNM.Physics.Standard(network=pn, phase=water, geometry=geom)
+    inlets = pn.pores('left')
+    OP = OpenPNM.Algorithms.OrdinaryPercolation(network=pn)
+    OP.setup(invading_phase=water)
+    OP.set_inlets(pores=inlets)
+    OP.run(npts=25)
+    a = ['pore.all', 'pore.inlets', 'pore.inv_Pc', 'pore.inv_sat',
+         'pore.inv_seq', 'throat.all', 'throat.entry_pressure',
+         'throat.inv_Pc', 'throat.inv_sat', 'throat.inv_seq']
+    assert sorted(list(OP.keys())) == a
+    V_inv = sp.sum(pn['pore.volume'][OP['pore.inv_Pc'] < sp.inf])
+    V_tot = sp.sum(pn['pore.volume'])
+    assert V_inv/V_tot == 1.0
+    ctrl.clear()
+
+
+def test_OP_trapping():
+    ctrl.clear()
+    pn = OpenPNM.Network.Cubic(shape=[30, 30, 1], spacing=0.01)
+    geom = OpenPNM.Geometry.Toray090(network=pn, pores=pn.Ps, throats=pn.Ts)
+    water = OpenPNM.Phases.Water(network=pn)
+    air = OpenPNM.Phases.Air(network=pn)
+    phys = OpenPNM.Physics.GenericPhysics(network=pn,
+                                          geometry=geom,
+                                          phase=water)
+    f = OpenPNM.Physics.models.capillary_pressure.washburn
+    phys.models.add(propname='throat.capillary_pressure',
+                    model=f)
+
+    OP = OpenPNM.Algorithms.OrdinaryPercolation(network=pn,
+                                                invading_phase=water,
+                                                defending_phase=air)
+    OP.run(inlets=pn.pores('left'), outlets=pn.pores('right'), trapping=True)
+    V_inv = sp.sum(pn['pore.volume'][OP['pore.inv_Pc'] < sp.inf])
+    V_tot = sp.sum(pn['pore.volume'])
+    assert V_inv/V_tot < 1.0
+    ctrl.clear()
diff --git a/test/integration/test_topological_manipulations.py b/test/integration/test_topological_manipulations.py
index 75e1f2c36..ee7b7bf38 100644
--- a/test/integration/test_topological_manipulations.py
+++ b/test/integration/test_topological_manipulations.py
@@ -1,5 +1,6 @@
 import OpenPNM
 import pytest
+import scipy as sp
 from OpenPNM.Utilities import topology
 ctrl = OpenPNM.Base.Controller()
 topo = topology()
@@ -27,3 +28,41 @@ def test_clone_and_trim():
     assert sorted(list(ctrl.keys())) == ['geo1', 'geo1_clone', 'net',
                                          'net_clone']
     topo.trim(network=pn2, pores=pn2.pores('top'))
+
+
+def test_trim_extend():
+    pn = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0), [1, 5, 25]))
+    assert [pn.Np, pn.Nt] == [125, 300]
+    pn.trim(pores=[0])
+    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0), [1, 5, 25]))
+    assert [pn.Np, pn.Nt] == [124, 297]
+    pn.extend(pore_coords=[0, 0, 0], throat_conns=[[124, 0]])
+    assert [pn.Np, pn.Nt] == [125, 298]
+    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0), [1, 5, 25, 124]))
+
+
+def test_stitch():
+    ctrl = OpenPNM.Base.Controller()
+    [Nx, Ny, Nz] = [10, 10, 10]
+    pn = OpenPNM.Network.Cubic(shape=[Nx, Ny, Nz])
+    pn2 = OpenPNM.Network.Cubic(shape=[Nx, Ny, Nz])
+    pn2['pore.coords'][:, 2] += Nz
+    pn.stitch(donor=pn2,
+              P_network=pn.pores('top'),
+              P_donor=pn2.pores('bottom'),
+              len_max=1,
+              method='nearest')
+    assert pn.Np == 2*pn2.Np  # Ensure number of pores doubled
+    assert pn.Nt == (2*pn2.Nt + Nx*Ny)  # Ensure correct number of new throats
+    assert pn2 not in ctrl.values()  # Donor Network is removed from Controller
+    # Reuse the donor Network in another stitch
+    pn2['pore.coords'][:, 2] -= 2*Nz
+    pn.stitch(donor=pn2,
+              P_network=pn.pores('bottom'),
+              P_donor=pn2.pores('top'),
+              len_max=1,
+              method='nearest')
+    assert pn.Np == 3*pn2.Np  # Ensure number of pores increased again
+    assert pn.Nt == (3*pn2.Nt + 2*Nx*Ny)  # Ensure correct num of new throats
+    ctrl.clear()