Merge pull request #2813 from PMEAL/enhance_create_im

Enhanced `create_incidence_matrix` to allow 2*Nt long weights #enh
Fixed bug in pore surface area models due to incorrect use of `flatten` on conns #bug

openpnm/algorithms/_transport.py

@@ -311,6 +311,7 @@ def rate(self, pores=[], throats=[], mode='group'):
         if g.size == self.Nt:
             g = np.tile(g, (2, 1)).T    # Make conductance an Nt by 2 matrix
         # The next line is critical for rates to be correct
+        # We could also do "g.T.flatten()" or "g.flatten('F')"
         g = np.flip(g, axis=1)
         Qt = np.diff(g*X12, axis=1).ravel()
 

openpnm/models/geometry/pore_surface_area/_funcs.py

@@ -40,7 +40,7 @@ def sphere(
     R = network[pore_diameter] / 2
     value = 4 * _np.pi * R**2
     Tca = network[throat_cross_sectional_area]
-    _np.subtract.at(value, network.conns.flatten(), _np.repeat(Tca, repeats=2))
+    _np.subtract.at(value, network.conns.T.flatten(), _np.repeat(Tca, repeats=2))
     return value
 
 
@@ -69,7 +69,7 @@ def circle(
     """
     value = _np.pi * network[pore_diameter]
     Tca = network[throat_cross_sectional_area]
-    _np.subtract.at(value, network.conns.flatten(), _np.repeat(Tca, repeats=2))
+    _np.subtract.at(value, network.conns.T.flatten(), _np.repeat(Tca, repeats=2))
     return value
 
 
@@ -95,7 +95,7 @@ def cube(
     D = network[pore_diameter]
     value = 6.0 * D**2
     Tca = network[throat_cross_sectional_area]
-    _np.subtract.at(value, network.conns.flatten(), _np.repeat(Tca, repeats=2))
+    _np.subtract.at(value, network.conns.T.flatten(), _np.repeat(Tca, repeats=2))
     return value
 
 
@@ -121,5 +121,5 @@ def square(
     D = network[pore_diameter]
     value = 4.0 * D
     Tca = network[throat_cross_sectional_area]
-    _np.subtract.at(value, network.conns.flatten(), _np.repeat(Tca, repeats=2))
+    _np.subtract.at(value, network.conns.T.flatten(), _np.repeat(Tca, repeats=2))
     return value

openpnm/network/_network.py

@@ -359,15 +359,16 @@ def create_incidence_matrix(self, weights=None, fmt='coo',
         # Check if provided data is valid
         if weights is None:
             weights = np.ones((self.Nt,), dtype=int)
-        elif np.shape(weights)[0] != self.Nt:
+        if np.shape(weights)[0] == self.Nt:
+            weights = np.append(weights, weights)
+        elif np.shape(weights)[0] != 2*self.Nt:
             raise Exception('Received dataset of incorrect length')
 
         conn = self['throat.conns']
-        row = conn[:, 0]
-        row = np.append(row, conn[:, 1])
+        row = conn[:, 1]
+        row = np.append(row, conn[:, 0])
         col = np.arange(self.Nt)
         col = np.append(col, col)
-        weights = np.append(weights, weights)
 
         temp = sprs.coo.coo_matrix((weights, (row, col)), (self.Np, self.Nt))
 

tests/unit/models/geometry/PoreSurfaceAreaTest.py

@@ -6,64 +6,66 @@
 
 class PoreSurfaceAreaTest:
     def setup_class(self):
-        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
-        self.net['pore.diameter'] = 1
-        self.net['throat.cross_sectional_area'] = 0.1
+        self.net3D = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
+        self.net3D['pore.diameter'] = 1
+        self.net3D['throat.cross_sectional_area'] = 0.1
+        self.net2D = op.network.Cubic(shape=[5, 5, 1], spacing=1.0)
+        self.net2D['pore.diameter'] = 1
+        self.net2D['throat.cross_sectional_area'] = 0.1
 
     def test_sphere(self):
-        self.net.add_model(propname='pore.surface_area',
-                           model=mods.sphere,
-                           regen_mode='normal')
+        self.net3D.add_model(propname='pore.surface_area',
+                             model=mods.sphere,
+                             regen_mode='normal')
         a = np.array([2.54159265, 2.64159265, 2.74159265, 2.84159265])
-        b = np.unique(self.net['pore.surface_area'])
+        b = np.unique(self.net3D['pore.surface_area'])
         assert_allclose(a, b)
 
     def test_circle(self):
-        self.net.add_model(propname='pore.surface_area',
+        self.net2D.add_model(propname='pore.surface_area',
                            model=mods.circle,
                            regen_mode='normal')
-        a = np.array([2.54159265, 2.64159265, 2.74159265, 2.84159265])
-        b = np.unique(self.net['pore.surface_area'])
+        a = np.array([2.74159265, 2.84159265, 2.94159265])
+        b = np.unique(self.net2D['pore.surface_area'])
         assert_allclose(a, b)
 
     def test_cube(self):
-        self.net.add_model(propname='pore.surface_area',
-                           model=mods.cube,
-                           regen_mode='normal')
+        self.net3D.add_model(propname='pore.surface_area',
+                             model=mods.cube,
+                             regen_mode='normal')
         a = np.array([5.4, 5.5, 5.6, 5.7])
-        b = np.unique(self.net['pore.surface_area'])
+        b = np.unique(self.net3D['pore.surface_area'])
         assert_allclose(a, b)
 
     def test_square(self):
-        self.net.add_model(propname='pore.surface_area',
+        self.net2D.add_model(propname='pore.surface_area',
                            model=mods.square,
                            regen_mode='normal')
-        a = np.array([3.4, 3.5, 3.6, 3.7])
-        b = np.unique(self.net['pore.surface_area'])
+        a = np.array([3.6, 3.7, 3.8])
+        b = np.unique(self.net2D['pore.surface_area'])
         assert_allclose(a, b)
 
     def test_circle_multi_geom(self):
-        self.net = op.network.Cubic(shape=[10, 1, 1], spacing=1.0)
-        self.net['pore.left'] = False
-        self.net['pore.left'][[0, 1, 2]] = True
-        self.net['pore.diameter'] = 1
-        self.net['throat.cross_sectional_area'] = 0.1
-        self.net['throat.cross_sectional_area'][[0, 1, 2]] = 0.3
-        self.net['pore.right'] = ~self.net['pore.left']
-        self.net.add_model(propname='pore.surface_area',
-                           model=mods.circle,
-                           domain='left',
-                           regen_mode='normal')
-        self.net.add_model(propname='pore.surface_area',
-                           model=mods.circle,
-                           domain='right',
-                           regen_mode='normal')
-        a = np.array([2.84159265, 2.54159265, 2.54159265])
-        b = self.net['pore.surface_area@left']
-        c = np.array([2.74159265, 2.94159265, 2.94159265,
-                      2.94159265, 2.94159265, 2.94159265,
-                      3.04159265])
-        d = self.net['pore.surface_area@right']
+        net = op.network.Cubic(shape=[10, 1, 1], spacing=1.0)
+        net['pore.left'] = False
+        net['pore.left'][[0, 1, 2]] = True
+        net['pore.diameter'] = 1
+        net['throat.cross_sectional_area'] = 0.1
+        net['throat.cross_sectional_area'][[0, 1, 2]] = 0.3
+        net['pore.right'] = ~net['pore.left']
+        net.add_model(propname='pore.surface_area',
+                        model=mods.circle,
+                        domain='left',
+                        regen_mode='normal')
+        net.add_model(propname='pore.surface_area',
+                      model=mods.circle,
+                      domain='right',
+                      regen_mode='normal')
+        a = np.array([2.84159265, 2.74159265, 2.74159265])
+        b = net['pore.surface_area@left']
+        c = np.array([2.74159265, 2.74159265, 2.74159265, 2.94159265,
+                      2.94159265, 2.94159265, 3.04159265])
+        d = net['pore.surface_area@right']
         assert_allclose(a, b)
         assert_allclose(c, d)