Adjust sa model

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 @@
  %(network)s
  %(Dp)s
  %(Act)s
+ amin : float
+ The lower limit for surface area, useful for preventing
 
  Returns
  -------
@@ -42,15 +49,19 @@
  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 @@
  %(network)s
  %(Dp)s
  %(Act)s
+ amin : float
+ The lower limit for surface area, useful for preventing
 
  Returns
  -------
@@ -73,14 +86,18 @@
  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 @@
  %(network)s
  %(Dp)s
  %(Act)s
+ amin : float
+ The lower limit for surface area, useful for preventing
 
  Returns
  -------
@@ -101,14 +120,18 @@
  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 @@
  %(network)s
  %(Dp)s
  %(Act)s
+ amin : float
+ The lower limit for surface area, useful for preventing
 
  Returns
  -------
@@ -129,5 +154,8 @@
  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

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__':