Repository: PMEAL/OpenPNM

Issue:
pore_area models are needlessly slow
We're using list comprehensions in these models, but could use unbuffered vector operations:

``` python

import numpy as np
import openpnm as op

net = op.network.Cubic(shape=[3, 3])

net['pore.diameter'] = 1
net['throat.diameter'] = np.random.rand(net.Nt)*0.5
net['throat.length'] = 1

net['throat.cross_sectional_area'] =  4/3*np.pi*net['throat.diameter']**2
net['throat.surface_area'] = np.pi*net['throat.diameter']*net['throat.length']
# Find surface are of sphere
temp = 4*np.pi*(0.5*net['pore.diameter'])**2
# Subtract throat openings
a = net['throat.cross_sectional_area']
np.subtract.at(temp, net.conns.flatten(), np.hstack((a, a)))
# Readd throat surface areeas
a = net['throat.surface_area']/2
np.add.at(temp, net.conns.flatten(), np.hstack((a, a)))

```

Discussion:

Patch:
diff --git a/openpnm/models/geometry/pore_surface_area/_funcs.py b/openpnm/models/geometry/pore_surface_area/_funcs.py
index 877347415..b5cc123a2 100644
--- a/openpnm/models/geometry/pore_surface_area/_funcs.py
+++ b/openpnm/models/geometry/pore_surface_area/_funcs.py
@@ -42,11 +42,11 @@ def sphere(
 
     """
     network = target.project.network
-    R = target[pore_diameter] / 2
-    Asurf = 4 * _np.pi * R**2
-    Tn = network.find_neighbor_throats(pores=target.Ps, flatten=False)
-    Tsurf = _np.array([network[throat_cross_sectional_area][Ts].sum() for Ts in Tn])
-    value = Asurf - Tsurf
+    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))
+    value = value[network.pores(target.name)]
     return value
 
 
@@ -74,11 +74,10 @@ def circle(
 
     """
     network = target.project.network
-    R = target[pore_diameter] / 2
-    Asurf = 2 * _np.pi * R
-    Tn = network.find_neighbor_throats(pores=target.Ps, flatten=False)
-    Tsurf = _np.array([network[throat_cross_sectional_area][Ts].sum() for Ts in Tn])
-    value = Asurf - Tsurf
+    value = _np.pi * network[pore_diameter]
+    Tca = network[throat_cross_sectional_area]
+    _np.subtract.at(value, network.conns.flatten(), _np.repeat(Tca, repeats=2))
+    value = value[network.pores(target.name)]
     return value
 
 
@@ -101,10 +100,11 @@ def cube(
 
     """
     network = target.project.network
-    D = target[pore_diameter]
-    Tn = network.find_neighbor_throats(pores=target.Ps, flatten=False)
-    Tsurf = _np.array([network[throat_cross_sectional_area][Ts].sum() for Ts in Tn])
-    value = 6 * D**2 - Tsurf
+    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))
+    value = value[network.pores(target.name)]
     return value
 
 
@@ -127,8 +127,9 @@ def square(
 
     """
     network = target.project.network
-    D = target[pore_diameter]
-    Tn = network.find_neighbor_throats(pores=target.Ps, flatten=False)
-    Tsurf = _np.array([network[throat_cross_sectional_area][Ts].sum() for Ts in Tn])
-    value = 4 * D - Tsurf
+    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))
+    value = value[network.pores(target.name)]
     return value


Test patch:
diff --git a/tests/unit/models/geometry/PoreSurfaceAreaTest.py b/tests/unit/models/geometry/PoreSurfaceAreaTest.py
index e8e877fcd..fd01e628d 100644
--- a/tests/unit/models/geometry/PoreSurfaceAreaTest.py
+++ b/tests/unit/models/geometry/PoreSurfaceAreaTest.py
@@ -49,6 +49,33 @@ class PoreSurfaceAreaTest:
         b = np.unique(self.geo['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.geo1 = op.geometry.GenericGeometry(network=self.net,
+                                                pores=self.net.Ps[0:3],
+                                                throats=self.net.Ts[0:3])
+        self.geo2 = op.geometry.GenericGeometry(network=self.net,
+                                                pores=self.net.Ps[3:],
+                                                throats=self.net.Ts[3:])
+        self.geo1['pore.diameter'] = 1
+        self.geo1['throat.cross_sectional_area'] = 0.3
+        self.geo2['pore.diameter'] = 1
+        self.geo2['throat.cross_sectional_area'] = 0.1
+        self.geo2.add_model(propname='pore.surface_area',
+                            model=mods.circle,
+                            regen_mode='normal')
+        self.geo1.add_model(propname='pore.surface_area',
+                            model=mods.circle,
+                            regen_mode='normal')
+        a = np.array([2.84159265, 2.54159265, 2.54159265])
+        b = self.geo1['pore.surface_area']
+        c = np.array([2.74159265, 2.94159265, 2.94159265,
+                      2.94159265, 2.94159265, 2.94159265,
+                      3.04159265])
+        d = self.geo2['pore.surface_area']
+        assert_allclose(a, b)
+        assert_allclose(c, d)
+
 
 if __name__ == '__main__':