Update surface area models

openpnm/__version__.py

@@ -1 +1 @@
-__version__ = '3.4.1.dev0'
+__version__ = '3.4.1.dev5'

openpnm/_skgraph/generators/tools/_funcs.py

@@ -85,7 +85,7 @@ def get_centroid(pts, mode='rigorous'):
     if mode == 'rigorous':
         tri = sptl.Delaunay(pts)
         CoM = center_of_mass(simplices=tri.simplices.astype(np.int_),
-                             points=tri.points.astype(np.float_))
+                             points=tri.points.astype(float))
     elif mode == 'fast':
         CoM = pts.mean(axis=0)
     return CoM
@@ -99,11 +99,11 @@ def get_centroid(pts, mode='rigorous'):
 # @njit
 def center_of_mass(simplices, points):
     A = []
-    centroids = np.zeros((simplices.shape[0], points.shape[1]), dtype=np.float_)
+    centroids = np.zeros((simplices.shape[0], points.shape[1]), dtype=float)
     for i, s in enumerate(simplices):
         xy = points[s]
         cen = np.sum(points[s], axis=0)/simplices.shape[1]
-        centroids[i, :] = cen.astype(np.float_)
+        centroids[i, :] = cen.astype(float)
         if xy.shape[1] == 2:  # In 2D
             # Use Heron's formula to find area of arbitrary triangle
             a = np.sqrt((xy[0, 0] - xy[1, 0])**2 + (xy[0, 1] - xy[1, 1])**2)
@@ -117,9 +117,9 @@ def center_of_mass(simplices, points):
             temp = np.abs(np.linalg.det(d))/6.0
             # temp = 0
         A.append(temp)
-    A = np.array(A, dtype=np.float_)
+    A = np.array(A, dtype=float)
     CoM = np.array([(centroids[:, i]*A/A.sum()*len(A)).mean()
-                    for i in range(centroids.shape[1])], dtype=np.float_)
+                    for i in range(centroids.shape[1])], dtype=float)
     return CoM
 
 

openpnm/algorithms/_reactive_transport.py

@@ -1,14 +1,16 @@
 import logging
 import sys
-
 import numpy as np
 from numpy.linalg import norm
-from scipy.optimize.nonlin import TerminationCondition
+try:  # For scipy < 1.14
+    from scipy.optimize.nonlin import TerminationCondition
+except ImportError:  # For newer Scipy
+    from scipy.optimize._nonlin import TerminationCondition
 from tqdm.auto import tqdm
-
 from openpnm.algorithms import Transport
 from openpnm.utils import Docorator, TypedList
 
+
 __all__ = ["ReactiveTransport"]
 
 

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

openpnm/solvers/_scipy.py

@@ -23,4 +23,7 @@
         if not isinstance(A, (csr_matrix, csc_matrix)):
             A = A.tocsr()
         atol = self._get_atol(b)
-        return cg(A, b, tol=self.tol, atol=atol, **kwargs)
+        try:
+            return cg(A, b, tol=self.tol, atol=atol, **kwargs)
+        except TypeError:
+            return cg(A, b, rtol=self.tol, atol=atol, **kwargs)

openpnm/utils/__init__.py

@@ -17,38 +17,13 @@
 
 def _get_version():
     from openpnm.__version__ import __version__ as ver
+
     suffix = ".dev0"
     if ver.endswith(suffix):
-        ver = ver[:-len(suffix)]
+        ver = ver[: -len(suffix)]
     return ver
 
 
-def _setup_logger():
-    import os
-    import logging
-    # You can add info to the logger message by inserting the desired %(item)
-    # For a list of available items see:
-    # https://docs.python.org/3/library/logging.html#logrecord-attributes
-    # NOTE: If the calling locations appears as 'root' it's because the logger
-    # was not given a name in a file somewhere.  A good option is __name__.
-
-    try:
-        terminal_width = os.get_terminal_size().columns
-    except OSError:
-        terminal_width = 60
-    column_width = terminal_width if terminal_width < 60 else 60
-
-    log_format = \
-        '-' * column_width + '\n\
-        - %(levelname)-7s: %(message)s \n\
-        - SOURCE : %(name)s.%(funcName)s \n\
-        - TIME   : %(asctime)s\
-        \n' + '-' * column_width
-
-    logging.basicConfig(level=logging.WARNING, format=log_format)
-    del log_format
-
-
 def _setup_logger_rich():
     import logging
     from rich.logging import RichHandler
@@ -57,3 +32,6 @@ def _setup_logger_rich():
     logging.basicConfig(
         format=FORMAT, datefmt="[%X]", handlers=[RichHandler(rich_tracebacks=True)]
     )
+
+
+_setup_logger_rich()

openpnm/visualization/_plottools.py

@@ -602,7 +602,7 @@ def plot_tutorial(network,
     fig = plt.gcf()
     fig.tight_layout()
     dims = op.topotools.dimensionality(network)
-    xy_range = network.coords.ptp(axis=0)[dims]
+    xy_range = np.ptp(network.coords, axis=0)[dims]
     aspect_ratio = xy_range[0] / xy_range[1]
     fig.set_size_inches(5, 5 / aspect_ratio)
 
@@ -781,7 +781,7 @@ def _generate_voxel_image(network, pore_shape, throat_shape, max_dim=200):
     # Transform points to satisfy origin at (0, 0, 0)
     xyz0 = xyz.min(axis=0) - delta
     xyz += -xyz0
-    res = (xyz.ptp(axis=0).max() + 2 * delta) / max_dim
+    res = (np.ptp(xyz, axis=0).max() + 2 * delta) / max_dim
     shape = np.rint((xyz.max(axis=0) + delta) / res).astype(int) + 2 * extra_clearance
 
     # Transforming from real coords to matrix coords
@@ -983,7 +983,7 @@ def set_mpl_style():  # pragma: no cover
     image_props = {'interpolation': 'none',
                    'cmap': 'viridis'}
     line_props = {'linewidth': 2,
-                  'markersize': 8,
+                  'markersize': 7,
                   'markerfacecolor': 'w'}
     font_props = {'size': sfont}
     axes_props = {'titlesize': lfont,

setup.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 3.4.1.dev0
+current_version = 3.4.1.dev5
 parse = (?P<major>\d+)\.(?P<minor>\d+)\.(?P<patch>\d+)\.(?P<release>\D+)(?P<build>\d+)?
 serialize = {major}.{minor}.{patch}.{release}{build}
 

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