Adding some integration tests (caught a few bugs in the process)

openpnm/geometry/GenericGeometry.py

@@ -106,7 +106,10 @@ def _set_locations(self, element, indices=[], mode='add'):
         r"""
         """
         element = self._parse_element(element=element, single=True)
-        indices = self._parse_indices(indices)
+        # Use the network's _parse_indices, since indicies could be 'network'
+        # length boolean masks
+        network = self.simulation.network
+        indices = network._parse_indices(indices)
 
         net = self.simulation.network
         sim = self.simulation

openpnm/physics/GenericPhysics.py

@@ -33,6 +33,7 @@ def __init__(self, network=None, phase=None, geometry=None, name=None):
         super().__init__(name=name, simulation=network.simulation)
         # Create a settings attribute
         self.settings['local_data'] = network.simulation.settings['local_data']
+        self.settings['boss'] = phase.name
         # Initialize a label dictionary in the associated phase and network
         phase['pore.'+self.name] = False
         phase['pore.'+self.name][network.pores(geometry.name)] = True
@@ -45,10 +46,10 @@ def __init__(self, network=None, phase=None, geometry=None, name=None):
 
     def __getitem__(self, key):
         element = key.split('.')[0]
+        boss = self.simulation[self.settings['boss']]
         if key.split('.')[-1] == '_id':
-            phase = self.simulation.find_phase(self)
-            inds = phase._get_indices(element=element, labels=self.name)
-            vals = phase[element+'._id'][inds]
+            inds = boss._get_indices(element=element, labels=self.name)
+            vals = boss[element+'._id'][inds]
         # Convert self.name into 'all'
         elif key.split('.')[-1] in [self.name]:
             vals = self[element+'.all']
@@ -57,23 +58,33 @@ def __getitem__(self, key):
             vals = super(Base, self).__getitem__(key)
         # Otherwise retrieve from network
         else:
-            phase = self.simulation.find_phase(self)
-            inds = phase._get_indices(element=element, labels=self.name)
-            vals = phase[key][inds]
+            inds = boss._get_indices(element=element, labels=self.name)
+            vals = boss[key][inds]
         return vals
 
     def __setitem__(self, key, value):
         if self.settings['local_data']:
             super().__setitem__(key, value)
         else:
-            phase = self.simulation.find_phase(self)
+            boss = self.simulation[self.settings['boss']]
             element = self._parse_element(key.split('.')[0], single=True)
-            inds = self._map(ids=self[element+'._id'], element=element,
+            inds = boss._map(ids=self[element+'._id'], element=element,
                              filtered=True)
-            # If array not in network, create it first
-            if key not in phase.keys():
+            # If array not in phase, create it first
+            if key not in boss.keys():
                 if value.dtype == bool:
-                    phase[key] = False
+                    boss[key] = False
                 else:
-                    phase[key] = sp.zeros_like(value)*sp.nan
-            phase[key][inds] = value
+                    dtype = value.dtype
+                    if dtype.name == 'object':
+                        boss[key] = sp.zeros(1, dtype=object)
+                    else:
+                        Nt = len(boss[element+'.all'])
+                        dim = sp.size(value[0])
+                        if dim > 1:
+                            arr = sp.zeros(dim, dtype=dtype)
+                            temp = sp.tile(arr, reps=(Nt, 1))*sp.nan
+                        else:
+                            temp = sp.zeros(Nt)*sp.nan
+                        boss[key] = temp
+            boss[key][inds] = value

tests/integration/test_drainage_algorithm.py

@@ -0,0 +1,176 @@
+import scipy as sp
+import OpenPNM
+import matplotlib
+
+pn = OpenPNM.Network.Cubic(shape=[25, 25, 25], spacing=0.00005)
+pn['pore.internal'] = True
+pn['throat.internal'] = True
+pn.add_boundary_pores(pores=pn.pores('top'),
+                      offset=[0, 0, 0.0001],
+                      apply_label='boundary_top')
+pn.add_boundary_pores(pores=pn.pores('bottom'),
+                      offset=[0, 0, -0.0001],
+                      apply_label='boundary_bottom')
+geom = OpenPNM.Geometry.Toray090(network=pn,
+                                 pores=pn.pores('internal'),
+                                 throats=pn.throats('internal'))
+boun = OpenPNM.Geometry.Boundary(network=pn,
+                                 pores=pn.pores('internal', mode='not'),
+                                 throats=pn.throats('internal', mode='not'))
+water = OpenPNM.Phases.Water(network=pn)
+air = OpenPNM.Phases.Air(network=pn)
+phys = OpenPNM.Physics.Standard(network=pn,
+                                phase=water,
+                                pores=pn.Ps,
+                                throats=pn.Ts)
+drainage = OpenPNM.Algorithms.Drainage(network=pn)
+
+
+# def test_basic():
+drainage.setup(invading_phase=water, defending_phase=air)
+drainage.set_inlets(pores=pn.pores('boundary_top'))
+drainage.run()
+data = drainage.get_drainage_data()
+
+assert sp.amin(data['invading_phase_saturation']) == 0.0
+assert sp.amax(data['invading_phase_saturation']) == 1.0
+
+
+# def test_residual():
+drainage.setup(invading_phase=water, defending_phase=air)
+drainage.set_inlets(pores=pn.pores('boundary_top'))
+Ps = sp.random.randint(0, pn.Np, 1000)
+Ts = sp.random.randint(0, pn.Nt, 1000)
+drainage.set_residual(pores=Ps, throats=Ts)
+drainage.run()
+data = drainage.get_drainage_data()
+
+assert sp.amin(data['invading_phase_saturation']) > 0
+assert sp.amax(data['invading_phase_saturation']) == 1.0
+
+
+# def test_trapping():
+drainage.setup(invading_phase=water, defending_phase=air, trapping=True)
+drainage.set_inlets(pores=pn.pores('boundary_top'))
+drainage.set_outlets(pores=pn.pores('boundary_bottom')[0:300])
+drainage.run()
+data = drainage.get_drainage_data()
+
+assert sp.amin(data['invading_phase_saturation']) == 0.0
+assert sp.amax(data['invading_phase_saturation']) < 1.0
+
+
+# def test_late_pore_filling():
+phys.models.add(propname='pore.fractional_filling',
+                model=OpenPNM.Physics.models.multiphase.late_pore_filling,
+                Pc=0,
+                Swp_star=0.2,
+                eta=1)
+phys.regenerate()
+drainage.setup(invading_phase=water, defending_phase=air,
+               pore_filling='pore.fractional_filling')
+drainage.set_inlets(pores=pn.pores('boundary_top'))
+drainage.run()
+data = drainage.get_drainage_data()
+assert sp.amin(data['invading_phase_saturation']) == 0.0
+assert sp.amax(data['invading_phase_saturation']) < 1.0
+
+drainage.return_results(Pc=5000)
+assert 'pore.occupancy' in water.keys()
+assert 'pore.partial_occupancy' in water.keys()
+
+
+# def test_late_throat_filling():
+mod = OpenPNM.Physics.models.multiphase.late_throat_filling
+phys.models.add(propname='throat.fractional_filling',
+                model=mod,
+                Pc=0,
+                Swp_star=0.2,
+                eta=1)
+phys.regenerate()
+drainage.setup(invading_phase=water, defending_phase=air,
+               throat_filling='throat.fractional_filling')
+drainage.set_inlets(pores=pn.pores('boundary_top'))
+drainage.run()
+data = drainage.get_drainage_data()
+
+assert sp.amin(data['invading_phase_saturation']) == 0.0
+assert sp.amax(data['invading_phase_saturation']) < 1.0
+
+drainage.return_results(Pc=5000)
+assert 'throat.occupancy' in water.keys()
+assert 'throat.partial_occupancy' in water.keys()
+
+
+# def test_late_pore_and_throat_filling():
+phys.models.add(propname='pore.fractional_filling',
+                model=OpenPNM.Physics.models.multiphase.late_pore_filling,
+                Pc=0,
+                Swp_star=0.2,
+                eta=1)
+mod = OpenPNM.Physics.models.multiphase.late_throat_filling
+phys.models.add(propname='throat.fractional_filling',
+                model=mod,
+                Pc=0,
+                Swp_star=0.2,
+                eta=1)
+phys.regenerate()
+drainage.setup(invading_phase=water, defending_phase=air,
+               pore_filling='pore.fractional_filling',
+               throat_filling='throat.fractional_filling')
+drainage.set_inlets(pores=pn.pores('boundary_top'))
+drainage.run()
+data = drainage.get_drainage_data()
+assert sp.amin(data['invading_phase_saturation']) == 0.0
+assert sp.amax(data['invading_phase_saturation']) < 1.0
+
+drainage.return_results(Pc=5000)
+assert 'pore.occupancy' in water.keys()
+assert 'throat.occupancy' in water.keys()
+assert 'pore.partial_occupancy' in water.keys()
+assert 'throat.partial_occupancy' in water.keys()
+
+
+# def test_ploting():
+drainage.setup(invading_phase=water, defending_phase=air)
+drainage.set_inlets(pores=pn.pores('boundary_top'))
+drainage.run()
+data = drainage.get_drainage_data()
+a = drainage.plot_drainage_curve(data)
+assert isinstance(a, matplotlib.figure.Figure)
+
+
+# def test_residual_and_lpf():
+phys.models.add(propname='pore.fractional_filling',
+                model=OpenPNM.Physics.models.multiphase.late_pore_filling,
+                Pc=0,
+                Swp_star=0.2,
+                eta=1)
+phys.models.add(propname='throat.fractional_filling',
+                model=OpenPNM.Physics.models.multiphase.late_throat_filling,
+                Pc=0,
+                Swp_star=0.2,
+                eta=1)
+phys.regenerate()
+drainage.setup(invading_phase=water, defending_phase=air,
+               pore_filling='pore.fractional_filling',
+               throat_filling='throat.fractional_filling')
+drainage.set_inlets(pores=pn.pores('boundary_top'))
+mask = sp.random.random(len(pn.pores('internal'))) < 0.1
+resPs = pn.pores('internal')[mask]
+mask = sp.random.random(len(pn.throats('internal'))) < 0.1
+resTs = pn.throats('internal')[mask]
+drainage.set_residual(pores=resPs, throats=resTs)
+drainage.run()
+drainage.return_results(Pc=5000)
+data = drainage.get_drainage_data()
+assert sp.all(water["pore.partial_occupancy"][resPs] == 1.0)
+assert sp.all(water["throat.partial_occupancy"][resTs] == 1.0)
+assert sp.amin(data['invading_phase_saturation']) > 0.0
+assert sp.amax(data['invading_phase_saturation']) < 1.0
+assert sp.all(water["pore.occupancy"]+air["pore.occupancy"] == 1.0)
+total_pp = water["pore.partial_occupancy"]+air["pore.partial_occupancy"]
+assert sp.all(total_pp == 1.0)
+assert sp.all(water["throat.occupancy"]+air["throat.occupancy"] == 1.0)
+total_pt = water["throat.partial_occupancy"]+air["throat.partial_occupancy"]
+assert sp.all(total_pt == 1.0)

tests/integration/test_multiple_geoms.py

@@ -0,0 +1,29 @@
+import openpnm as op
+
+pn = op.network.Cubic(shape=[10, 10, 40], spacing=1)
+
+# Generate Ps as boolean mask
+Ps = pn['pore.coords'][:, 2] <= 20
+Ts = pn.find_neighbor_throats(pores=Ps, mode='intersection', flatten=True)
+geom1 = op.geometry.GenericGeometry(network=pn, pores=Ps, throats=Ts)
+
+# Convert Ps to indices
+Ps = pn.toindices(pn['pore.coords'][:, 2] > 20)
+Ts = pn.find_neighbor_throats(pores=Ps, mode='union')
+geom2 = op.geometry.GenericGeometry(network=pn, pores=Ps, throats=Ts)
+
+water = op.phases.Water(network=pn)
+
+phys1 = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom1)
+phys1.add_model(propname='pore.test',
+                model=op.geometry.models.pore_misc.constant,
+                value=1.0)
+phys1.regenerate_models()
+
+phys2 = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom2)
+# Copy models from phys1 to phys2 directly
+phys2.models = phys1.models
+# Ensure the models dicts are equal
+assert phys1.models == phys2.models
+# But they are actually distinct objects
+assert phys1.models is not phys2.models

tests/integration/test_network_methods.py

@@ -0,0 +1,93 @@
+import scipy as sp
+import OpenPNM
+import pytest
+
+
+def test_find_connected_pores():
+    pn = OpenPNM.Network.Cubic(shape=(10, 10, 10))
+    a = pn.find_connected_pores(throats=[0, 1])
+    assert sp.all(a.flatten() == [0, 1, 1, 2])
+    a = pn.find_connected_pores(throats=[0, 1], flatten=True)
+    assert sp.all(a == [0, 1, 2])
+    Tind = sp.zeros((pn.Nt,), dtype=bool)
+    Tind[[0, 1]] = True
+    a = pn.find_connected_pores(throats=Tind, flatten=True)
+    assert sp.all(a == [0, 1, 2])
+    a = pn.find_connected_pores(throats=[], flatten=True)
+    assert sp.shape(a) == (0, )
+
+
+def test_find_neighbor_pores():
+    pn = OpenPNM.Network.Cubic(shape=(10, 10, 10))
+    a = pn.find_neighbor_pores(pores=[])
+    assert sp.size(a) == 0
+    Pind = sp.zeros((pn.Np,), dtype=bool)
+    Pind[[0, 1]] = True
+    a = pn.find_neighbor_pores(pores=Pind)
+    assert sp.all(a == [2, 10, 11, 100, 101])
+    a = pn.find_neighbor_pores(pores=[0, 2], mode='union')
+    assert sp.all(a == [1, 3, 10,  12, 100, 102])
+    a = pn.find_neighbor_pores(pores=[0, 2], mode='intersection')
+    assert sp.all(a == [1])
+    a = pn.find_neighbor_pores(pores=[0, 2], mode='not_intersection')
+    assert sp.all(a == [3, 10, 12, 100, 102])
+    a = pn.find_neighbor_pores(pores=[0, 2],
+                               mode='union',
+                               excl_self=False)
+    assert sp.all(a == [0, 1, 2, 3, 10, 12, 100, 102])
+    a = pn.find_neighbor_pores(pores=[0, 2],
+                               mode='intersection',
+                               excl_self=False)
+    assert sp.all(a == [1])
+    a = pn.find_neighbor_pores(pores=[0, 2],
+                               mode='not_intersection',
+                               excl_self=False)
+    assert sp.all(a == [0, 2, 3, 10, 12, 100, 102])
+
+
+def test_find_neighbor_throats():
+    pn = OpenPNM.Network.Cubic(shape=(10, 10, 10))
+    a = pn.find_neighbor_throats(pores=[])
+    assert sp.size(a) == 0
+    Pind = sp.zeros((pn.Np,), dtype=bool)
+    Pind[[0, 1]] = True
+    a = pn.find_neighbor_throats(pores=Pind)
+    assert sp.all(a == [0, 1, 900, 901, 1800, 1801])
+    a = pn.find_neighbor_throats(pores=[0, 2], mode='union')
+    assert sp.all(a == [0, 1, 2, 900, 902, 1800, 1802])
+    a = pn.find_neighbor_throats(pores=[0, 2], mode='intersection')
+    assert sp.size(a) == 0
+    a = pn.find_neighbor_throats(pores=[0, 2], mode='not_intersection')
+    assert sp.all(a == [0, 1, 2, 900, 902, 1800, 1802])
+
+
+def test_num_neighbors():
+    pn = OpenPNM.Network.Cubic(shape=(10, 10, 10))
+    a = pn.num_neighbors(pores=[])
+    assert sp.size(a) == 0
+    Pind = sp.zeros((pn.Np,), dtype=bool)
+    Pind[0] = True
+    a = pn.num_neighbors(pores=Pind)
+    assert a == 3
+    a = pn.num_neighbors(pores=[0, 2], flatten=True)
+    assert a == 6
+    assert isinstance(a, int)
+    a = pn.num_neighbors(pores=[0, 2], flatten=False)
+    assert sp.all(a == [3, 4])
+    a = pn.num_neighbors(pores=0, flatten=False)
+    assert sp.all(a == [3])
+    assert isinstance(a, sp.ndarray)
+
+
+def test_find_interface_throats():
+    pn = OpenPNM.Network.Cubic(shape=(3, 3, 3))
+    pn['pore.domain1'] = False
+    pn['pore.domain2'] = False
+    pn['pore.domain3'] = False
+    pn['pore.domain1'][[0, 1, 2]] = True
+    pn['pore.domain2'][[5, 6, 7]] = True
+    pn['pore.domain3'][18:26] = True
+    a = pn.find_interface_throats(labels=['domain1', 'domain2'])
+    assert a == [20]
+    a = pn.find_interface_throats(labels=['domain1', 'domain3'])
+    assert sp.size(a) == 0