Merge pull request #2893 from PMEAL/adjust_sa_model

Adjusted surface area model to accept lower limit (amin), defaults to 0 #enh
PMEAL · Jan 27, 2024 · c24a6ae · c24a6ae
2 parents 0273ab4 + 7c45c22
commit c24a6ae
Show file tree

Hide file tree

Showing 2 changed files with 53 additions and 8 deletions.
diff --git a/openpnm/models/geometry/pore_surface_area/_funcs.py b/openpnm/models/geometry/pore_surface_area/_funcs.py
@@ -1,5 +1,9 @@
 import numpy as _np
 from openpnm.models.geometry import _geodocs
+import logging
+
+
+logger = logging.getLogger(__name__)
 
 
 __all__ = ["sphere",
@@ -13,7 +17,8 @@
 def sphere(
     network,
     pore_diameter='pore.diameter',
-    throat_cross_sectional_area='throat.cross_sectional_area'
+    throat_cross_sectional_area='throat.cross_sectional_area',
+    amin=0,
 ):
     r"""
     Calculates internal surface area of pore bodies assuming they are spheres,
@@ -24,6 +29,8 @@ def sphere(
     %(network)s
     %(Dp)s
     %(Act)s
+    amin : float
+        The lower limit for surface area, useful for preventing
 
     Returns
     -------
@@ -42,15 +49,19 @@ def sphere(
     Tca = network[throat_cross_sectional_area]
     _np.subtract.at(value, network.conns[:, 0], Tca)
     _np.subtract.at(value, network.conns[:, 1], Tca)
-    value = _np.clip(value, 1e-12, _np.inf)
+    if amin is not None:
+        value = _np.clip(value, amin, _np.inf)
+    elif _np.any(value < 0):
+        logger.warn('Negative pore surface areas found')
     return value
 
 
 @_geodocs
 def circle(
     network,
     pore_diameter='pore.diameter',
-    throat_cross_sectional_area='throat.cross_sectional_area'
+    throat_cross_sectional_area='throat.cross_sectional_area',
+    amin=0,
 ):
     r"""
     Calculates internal surface area of pore bodies assuming they are
@@ -61,6 +72,8 @@ def circle(
     %(network)s
     %(Dp)s
     %(Act)s
+    amin : float
+        The lower limit for surface area, useful for preventing
 
     Returns
     -------
@@ -73,14 +86,18 @@ def circle(
     Tca = network[throat_cross_sectional_area]
     _np.subtract.at(value, network.conns[:, 0], Tca)
     _np.subtract.at(value, network.conns[:, 1], Tca)
-    value = _np.clip(value, 1e-12, _np.inf)
+    if amin is not None:
+        value = _np.clip(value, amin, _np.inf)
+    elif _np.any(value < 0):
+        logger.warn('Negative pore surface areas found')
     return value
 
 
 def cube(
     network,
     pore_diameter='pore.diameter',
-    throat_cross_sectional_area='throat.cross_sectional_area'
+    throat_cross_sectional_area='throat.cross_sectional_area',
+    amin=0,
 ):
     r"""
     Calculates internal surface area of pore bodies assuming they are cubes
@@ -91,6 +108,8 @@ def cube(
     %(network)s
     %(Dp)s
     %(Act)s
+    amin : float
+        The lower limit for surface area, useful for preventing
 
     Returns
     -------
@@ -101,14 +120,18 @@ def cube(
     Tca = network[throat_cross_sectional_area]
     _np.subtract.at(value, network.conns[:, 0], Tca)
     _np.subtract.at(value, network.conns[:, 1], Tca)
-    value = _np.clip(value, 1e-12, _np.inf)
+    if amin is not None:
+        value = _np.clip(value, amin, _np.inf)
+    elif _np.any(value < 0):
+        logger.warn('Negative pore surface areas found')
     return value
 
 
 def square(
     network,
     pore_diameter='pore.diameter',
-    throat_cross_sectional_area='throat.cross_sectional_area'
+    throat_cross_sectional_area='throat.cross_sectional_area',
+    amin=0,
 ):
     r"""
     Calculates internal surface area of pore bodies assuming they are
@@ -119,6 +142,8 @@ def square(
     %(network)s
     %(Dp)s
     %(Act)s
+    amin : float
+        The lower limit for surface area, useful for preventing
 
     Returns
     -------
@@ -129,5 +154,8 @@ def square(
     Tca = network[throat_cross_sectional_area]
     _np.subtract.at(value, network.conns[:, 0], Tca)
     _np.subtract.at(value, network.conns[:, 1], Tca)
-    value = _np.clip(value, 1e-12, _np.inf)
+    if amin is not None:
+        value = _np.clip(value, amin, _np.inf)
+    elif _np.any(value < 0):
+        logger.warn('Negative pore surface areas found')
     return value
diff --git a/tests/unit/models/geometry/PoreSurfaceAreaTest.py b/tests/unit/models/geometry/PoreSurfaceAreaTest.py
@@ -2,6 +2,7 @@
 import openpnm.models.geometry.pore_surface_area as mods
 import numpy as np
 from numpy.testing import assert_allclose
+import pytest
 
 
 class PoreSurfaceAreaTest:
@@ -69,6 +70,22 @@ def test_circle_multi_geom(self):
         assert_allclose(a, b)
         assert_allclose(c, d)
 
+    def test_negative_surface_area(self):
+        pn = op.network.Cubic([2, 1, 1])
+        pn.add_model_collection(op.models.collections.geometry.circles_and_rectangles)
+        pn.regenerate_models()
+        pn.add_model(propname='pore.surface_area', model=op.models.geometry.pore_surface_area.circle)
+        pn['throat.cross_sectional_area'] = 100
+        pn.regenerate_models('pore.surface_area@all')
+        assert np.all(pn['pore.surface_area'] == 0)
+        with pytest.warns():
+            pn.add_model(
+                propname='pore.surface_area',
+                model=op.models.geometry.pore_surface_area.circle,
+                amin=None,
+            )
+        assert np.all(pn['pore.surface_area'] < 0)
+
 
 if __name__ == '__main__':