Merge pull request #2889 from PMEAL/fix_pore_surface_area

Updated pore surface area calculation models #bug
PMEAL · Jan 26, 2024 · 47e0ac7 · 47e0ac7
2 parents 0baf355 + 51c3715
commit 47e0ac7
Show file tree

Hide file tree

Showing 2 changed files with 21 additions and 13 deletions.
diff --git a/openpnm/models/geometry/pore_surface_area/_funcs.py b/openpnm/models/geometry/pore_surface_area/_funcs.py
@@ -40,7 +40,9 @@ 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.T.flatten(), _np.repeat(Tca, repeats=2))
+    _np.subtract.at(value, network.conns[:, 0], Tca)
+    _np.subtract.at(value, network.conns[:, 1], Tca)
+    value = _np.clip(value, 1e-12, _np.inf)
     return value
 
 
@@ -69,7 +71,9 @@ def circle(
     """
     value = _np.pi * network[pore_diameter]
     Tca = network[throat_cross_sectional_area]
-    _np.subtract.at(value, network.conns.T.flatten(), _np.repeat(Tca, repeats=2))
+    _np.subtract.at(value, network.conns[:, 0], Tca)
+    _np.subtract.at(value, network.conns[:, 1], Tca)
+    value = _np.clip(value, 1e-12, _np.inf)
     return value
 
 
@@ -95,7 +99,9 @@ def cube(
     D = network[pore_diameter]
     value = 6.0 * D**2
     Tca = network[throat_cross_sectional_area]
-    _np.subtract.at(value, network.conns.T.flatten(), _np.repeat(Tca, repeats=2))
+    _np.subtract.at(value, network.conns[:, 0], Tca)
+    _np.subtract.at(value, network.conns[:, 1], Tca)
+    value = _np.clip(value, 1e-12, _np.inf)
     return value
 
 
@@ -121,5 +127,7 @@ def square(
     D = network[pore_diameter]
     value = 4.0 * D
     Tca = network[throat_cross_sectional_area]
-    _np.subtract.at(value, network.conns.T.flatten(), _np.repeat(Tca, repeats=2))
+    _np.subtract.at(value, network.conns[:, 0], Tca)
+    _np.subtract.at(value, network.conns[:, 1], Tca)
+    value = _np.clip(value, 1e-12, _np.inf)
     return value
diff --git a/tests/unit/models/geometry/PoreSurfaceAreaTest.py b/tests/unit/models/geometry/PoreSurfaceAreaTest.py
@@ -23,8 +23,8 @@ def test_sphere(self):
 
     def test_circle(self):
         self.net2D.add_model(propname='pore.surface_area',
-                           model=mods.circle,
-                           regen_mode='normal')
+                             model=mods.circle,
+                             regen_mode='normal')
         a = np.array([2.74159265, 2.84159265, 2.94159265])
         b = np.unique(self.net2D['pore.surface_area'])
         assert_allclose(a, b)
@@ -39,8 +39,8 @@ def test_cube(self):
 
     def test_square(self):
         self.net2D.add_model(propname='pore.surface_area',
-                           model=mods.square,
-                           regen_mode='normal')
+                             model=mods.square,
+                             regen_mode='normal')
         a = np.array([3.6, 3.7, 3.8])
         b = np.unique(self.net2D['pore.surface_area'])
         assert_allclose(a, b)
@@ -54,16 +54,16 @@ def test_circle_multi_geom(self):
         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')
+                      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])
+        a = np.array([2.84159265, 2.54159265, 2.54159265])
         b = net['pore.surface_area@left']
-        c = np.array([2.74159265, 2.74159265, 2.74159265, 2.94159265,
+        c = np.array([2.74159265, 2.94159265, 2.94159265, 2.94159265,
                       2.94159265, 2.94159265, 3.04159265])
         d = net['pore.surface_area@right']
         assert_allclose(a, b)