[Hotfix] Fixing handling of throat.centroid for overlapping but non-a…

…ligned pores
PMEAL · Apr 15, 2020 · 1546fa1 · 1546fa1
2 parents a1ea84e + 35f1f7b
commit 1546fa1
Show file tree

Hide file tree

Showing 9 changed files with 75 additions and 62 deletions.
diff --git a/VERSIONING.md b/VERSIONING.md
@@ -1,7 +1,7 @@
 ## Release Management and Versioning
 
-OpenPNM uses [Semantic Versioning](http://semver.org) (i.e. X.Y.Z) to label releases.  All major and minor versions (X.Y.z) are available on [PyPI](https://pypi.python.org/pypi), but bugfix releases (x.y.Z) are not generally pushed unless the bug is important.  Bubfix releases are available via download of the source code from Github.
+`OpenPNM` uses [Semantic Versioning](http://semver.org) (i.e. X.Y.Z) to label releases.  All major and minor versions (X.Y.z) are available on [PyPI](https://pypi.python.org/pypi), but bugfix releases (x.y.Z) are not generally pushed unless the bug is important.  Bugfix releases are available via download of the source code from Github.
 
-OpenPNM uses the [Github Flow](https://guides.github.com/introduction/flow/) system of Git branching, except instead of merging PRs into *master*, they are merged into a branch called *dev*. Any code added to *dev* is done via Pull Requests (PRs).  When new PRs are merged into the *dev* branch, they are *not* given a new version number. Once enough new features have been added, the *dev* branch is merged into the *master* branch, and the minor release number (x.Y.z) will be incremented. An exception to this rule are bugfixes which may be found on *master*.  In these cases a PR can be merged into *master* and the version number will be incremented (x.y.Z) to indicate the fix.
+`OpenPNM` uses the [Github Flow](https://guides.github.com/introduction/flow/) system of Git branching, except instead of merging PRs into *master*, they are merged into a branch called *dev*. Any code added to *dev* is done via Pull Requests (PRs).  When new PRs are merged into the *dev* branch, they are *not* given a new version number. Once enough new features have been added, the *dev* branch is merged into the *master* branch, and the minor release number (x.Y.z) will be incremented. An exception to this rule are bugfixes which may be found on *master*.  In these cases a PR can be merged into *master* and the version number will be incremented (x.y.Z) to indicate the fix.
 
-OpenPNM depends on several other packages widely known as the [Scipy Stack](https://www.scipy.org/stackspec.html).  It is our policy to always support the latest version of all these packages and their dependencies.
+`OpenPNM` depends on several other packages widely known as the [Scipy Stack](https://www.scipy.org/stackspec.html).  It is our policy to always support the latest version of all these packages and their dependencies.
diff --git a/openpnm/__init__.py b/openpnm/__init__.py
@@ -52,7 +52,7 @@
 
 """
 
-__version__ = '2.3.0'
+__version__ = '2.3.1'
 
 import numpy as np
 np.seterr(divide='ignore', invalid='ignore')

diff --git a/openpnm/algorithms/GenericTransport.py b/openpnm/algorithms/GenericTransport.py
@@ -706,7 +706,7 @@ def _check_for_nans(self):
             if root_props:
                 msg = (
                     f"Found NaNs in A matrix, possibly caused by NaNs in "
-                    f"{', '.join(root_props)} \n{'-' * 80}\nThe issue might get "
+                    f"{', '.join(root_props)}. The issue might get "
                     f"resolved if you call regenerate_models on the following "
                     f"object(s): {', '.join(root_objs)}"
                 )

diff --git a/openpnm/models/geometry/throat_endpoints.py b/openpnm/models/geometry/throat_endpoints.py
@@ -113,11 +113,11 @@ def spherical_pores(target, pore_diameter='pore.diameter',
     Notes
     -----
     (1) This model should not be applied to true 2D networks. Use
-    `circular_pores` model instead.
+    ``circular_pores`` model instead.
 
     (2) By default, this model assumes that throat centroid and pore
     coordinates are colinear. If that's not the case, such as in extracted
-    networks, `throat_centroid` could be passed as an optional argument, and
+    networks, ``throat_centroid`` could be passed as an optional argument, and
     the model takes care of the rest.
 
     """
@@ -153,11 +153,12 @@ def spherical_pores(target, pore_diameter='pore.diameter',
     EP2 = xyz[cn[:, 1]] + L2[:, None] * unit_vec_P2T
     # Handle throats w/ overlapping pores
     L1 = (4*L**2 + D1**2 - D2**2) / (8*L)
-    # L2 = (4*L**2 + D2**2 - D1**2) / (8*L)
+    L2 = (4*L**2 + D2**2 - D1**2) / (8*L)
     h = (2*_np.sqrt(D1**2/4 - L1**2)).real
     overlap = L - 0.5 * (D1+D2) < 0
     mask = overlap & (Dt < h)
-    EP1[mask] = EP2[mask] = (xyz[cn[:, 0]] + L1[:, None] * unit_vec_P1T)[mask]
+    EP1[mask] = (xyz[cn[:, 0]] + L1[:, None] * unit_vec_P1T)[mask]
+    EP2[mask] = (xyz[cn[:, 1]] + L2[:, None] * unit_vec_P2T)[mask]
     return {'head': EP1, 'tail': EP2}
 
 

diff --git a/openpnm/models/geometry/throat_length.py b/openpnm/models/geometry/throat_length.py
@@ -6,7 +6,7 @@
 
 """
 import numpy as _np
-from scipy import sqrt as _sqrt
+from numpy.linalg import norm as _norm
 from openpnm.utils import logging as _logging
 _logger = _logging.getLogger(__name__)
 
@@ -35,7 +35,7 @@ def ctc(target):
     cn = network['throat.conns'][throats]
     C1 = network['pore.coords'][cn[:, 0]]
     C2 = network['pore.coords'][cn[:, 1]]
-    value = _sqrt(((C1 - C2)**2).sum(axis=1))
+    value = _norm(C1 - C2, axis=1)
     return value
 
 
@@ -77,19 +77,19 @@ def piecewise(target, throat_endpoints='throat.endpoints',
     EP1 = network[throat_endpoints + '.head'][throats]
     EP2 = network[throat_endpoints + '.tail'][throats]
     # Calculate throat length
-    Lt = _sqrt(((EP1 - EP2)**2).sum(axis=1))
+    Lt = _norm(EP1 - EP2, axis=1)
     # Handle the case where pores & throat centroids are not colinear
     try:
         Ct = network[throat_centroid][throats]
-        Lt = _sqrt(((Ct - EP1)**2).sum(axis=1)) + \
-            _sqrt(((Ct - EP2)**2).sum(axis=1))
+        Lt = _norm(Ct - EP1, axis=1) + _norm(Ct - EP2, axis=1)
     except KeyError:
         pass
     return Lt
 
 
 def conduit_lengths(target, throat_endpoints='throat.endpoints',
-                    throat_length='throat.length'):
+                    throat_length='throat.length',
+                    throat_centroid='throat.centroid'):
     r"""
     Calculate conduit lengths. A conduit is defined as half pore + throat
     + half pore.
@@ -130,11 +130,11 @@ def conduit_lengths(target, throat_endpoints='throat.endpoints',
         # Look up throat length if given
         Lt = network[throat_length][throats]
     except KeyError:
-        # Calculate throat length otherwise
-        Lt = _sqrt(((EP1 - EP2)**2).sum(axis=1))
-    # Calculate conduit lengths
-    L1 = _sqrt(((C1 - EP1)**2).sum(axis=1))
-    L2 = _sqrt(((C2 - EP2)**2).sum(axis=1))
+        # Calculate throat length otherwise based on piecewise model
+        Lt = piecewise(target, throat_endpoints, throat_centroid)
+    # Calculate conduit lengths for pore 1 and pore 2
+    L1 = _norm(C1 - EP1, axis=1)
+    L2 = _norm(C2 - EP2, axis=1)
     return {'pore1': L1, 'throat': Lt, 'pore2': L2}
 
 
@@ -156,8 +156,6 @@ def classic(target, pore_diameter='pore.diameter'):
     network = target.project.network
     throats = network.map_throats(throats=target.Ts, origin=target)
     cn = network['throat.conns'][throats]
-    C1 = network['pore.coords'][cn[:, 0]]
-    C2 = network['pore.coords'][cn[:, 1]]
-    D = _sqrt(((C1 - C2)**2).sum(axis=1))
-    value = D - _np.sum(network[pore_diameter][cn], axis=1)/2
+    ctc_dist = ctc(target)
+    value = ctc_dist - network[pore_diameter][cn].sum(axis=1)/2
     return value
diff --git a/openpnm/models/physics/hydraulic_conductance.py b/openpnm/models/physics/hydraulic_conductance.py
@@ -7,17 +7,17 @@
 """
 import scipy as _sp
 import numpy as _np
-from numpy import pi as _pi
 
 
 def hagen_poiseuille(
-        target,
-        pore_area="pore.area",
-        throat_area="throat.area",
-        pore_viscosity="pore.viscosity",
-        throat_viscosity="throat.viscosity",
-        conduit_lengths="throat.conduit_lengths",
-        conduit_shape_factors="throat.flow_shape_factors"):
+    target,
+    pore_area="pore.area",
+    throat_area="throat.area",
+    pore_viscosity="pore.viscosity",
+    throat_viscosity="throat.viscosity",
+    conduit_lengths="throat.conduit_lengths",
+    conduit_shape_factors="throat.flow_shape_factors"
+):
     r"""
     Calculate the hydraulic conductance of conduits in network, where a
     conduit is ( 1/2 pore - full throat - 1/2 pore ). See the notes section.
@@ -94,24 +94,14 @@ def hagen_poiseuille(
         SF2 = phase[conduit_shape_factors + ".pore2"][throats]
     except KeyError:
         SF1 = SF2 = SFt = 1.0
-    # Interpolate pore phase property values to throats
-    try:
-        Dt = phase[throat_viscosity][throats]
-    except KeyError:
-        Dt = phase.interpolate_data(propname=pore_viscosity)[throats]
-    try:
-        D1 = phase[pore_viscosity][cn[:, 0]]
-        D2 = phase[pore_viscosity][cn[:, 1]]
-    except KeyError:
-        D1 = phase.interpolate_data(propname=throat_viscosity)[cn[:, 0]]
-        D2 = phase.interpolate_data(propname=throat_viscosity)[cn[:, 1]]
+    Dt = phase[throat_viscosity][throats]
+    D1, D2 = phase[pore_viscosity][cn].T
     # Find g for half of pore 1, throat, and half of pore 2
-    pi = _sp.pi
-    g1[m1] = A1[m1] ** 2 / (8 * pi * D1 * L1)[m1]
-    g2[m2] = A2[m2] ** 2 / (8 * pi * D2 * L2)[m2]
-    gt[mt] = At[mt] ** 2 / (8 * pi * Dt * Lt)[mt]
+    g1[m1] = A1[m1] ** 2 / (8 * _sp.pi * D1 * L1)[m1]
+    g2[m2] = A2[m2] ** 2 / (8 * _sp.pi * D2 * L2)[m2]
+    gt[mt] = At[mt] ** 2 / (8 * _sp.pi * Dt * Lt)[mt]
     # Apply shape factors and calculate the final conductance
-    return (1 / gt / SFt + 1 / g1 / SF1 + 1 / g2 / SF2) ** (-1)
+    return (1/gt/SFt + 1/g1/SF1 + 1/g2/SF2) ** (-1)
 
 
 def hagen_poiseuille_2D(
@@ -191,22 +181,14 @@ def hagen_poiseuille_2D(
     except KeyError:
         SF1 = SF2 = SFt = 1.0
     # Getting viscosity values
-    try:
-        mut = phase[throat_viscosity][throats]
-    except KeyError:
-        mut = phase.interpolate_data(propname=pore_viscosity)[throats]
-    try:
-        mu1 = phase[pore_viscosity][cn[:, 0]]
-        mu2 = phase[pore_viscosity][cn[:, 1]]
-    except KeyError:
-        mu1 = phase.interpolate_data(propname=throat_viscosity)[cn[:, 0]]
-        mu2 = phase.interpolate_data(propname=throat_viscosity)[cn[:, 1]]
+    mut = phase[throat_viscosity][throats]
+    mu1, mu2 = phase[pore_viscosity][cn].T
     # Find g for half of pore 1, throat, and half of pore 2
     g1 = D1 ** 3 / (12 * mu1 * L1)
     g2 = D2 ** 3 / (12 * mu2 * L2)
     gt = Dt ** 3 / (12 * mut * Lt)
 
-    return (1 / gt / SFt + 1 / g1 / SF1 + 1 / g2 / SF2) ** (-1)
+    return (1/gt/SFt + 1/g1/SF1 + 1/g2/SF2) ** (-1)
 
 
 def hagen_poiseuille_power_law(
@@ -418,7 +400,7 @@ def hagen_poiseuille_power_law(
     gt[mt] = At[mt] ** 2 / ((8 * pi * Lt) * mut)[mt]
 
     # Apply shape factors and calculate the final conductance
-    return (1 / gt / SFt + 1 / g1 / SF1 + 1 / g2 / SF2) ** (-1)
+    return (1/gt/SFt + 1/g1/SF1 + 1/g2/SF2) ** (-1)
 
 
 def valvatne_blunt(

diff --git a/tests/unit/core/ModelsTest.py b/tests/unit/core/ModelsTest.py
@@ -22,7 +22,7 @@ def test_models_dict_print(self):
         geo = op.geometry.StickAndBall(network=net, pores=net.Ps,
                                        throats=net.Ts)
         s = geo.models.__str__().split('\n')
-        assert len(s) == 69
+        assert len(s) == 70
         assert s.count('―'*85) == 15
 
     def test_regenerate_models(self):

diff --git a/tests/unit/models/geometry/ThroatEndpointsTest.py b/tests/unit/models/geometry/ThroatEndpointsTest.py
@@ -61,6 +61,23 @@ def test_spherical_pores_with_throat_centroid(self):
         assert_allclose(EP2, desired=EP2d)
         del self.geo['throat.centroid']
 
+    def test_spherical_pores_with_throat_centroid_and_overlap(self):
+        self.geo['throat.centroid'] = np.array([[0, 0.5, 1],
+                                                [0, 1.5, 0]]) + self.base
+        self.geo['pore.diameter'] = [1.5, 1.0, 0.5]
+        self.geo['throat.diameter'] = 0.25
+        self.geo.add_model(propname='throat.endpoints',
+                            model=mods.spherical_pores,
+                            regen_mode='normal')
+        # Only check throat 1->2, 2->3 is already tested (no-overlap)
+        EP12_head = self.geo['throat.endpoints.head'][0]
+        EP12_tail = self.geo['throat.endpoints.tail'][0]
+        EP12_head_desired = np.array([0, 0.29348392, 0.58696784]) + self.base
+        EP12_tail_desired = np.array([0, 0.84627033, 0.30745935]) + self.base
+        assert_allclose(EP12_head, desired=EP12_head_desired)
+        assert_allclose(EP12_tail, desired=EP12_tail_desired)
+        del self.geo["throat.centroid"]
+
     def test_cubic_pores(self):
         self.geo['pore.diameter'] = 0.5
         self.geo['throat.diameter'] = 0.25

diff --git a/tests/unit/models/geometry/ThroatLengthTest.py b/tests/unit/models/geometry/ThroatLengthTest.py
@@ -29,7 +29,8 @@ def test_piecewise_with_centroid(self):
         actual = self.geo['throat.length'][0]
         assert_approx_equal(actual, desired=1.1216117)
 
-    def test_conduit_lengths(self):
+    def test_conduit_lengths_wo_centroid(self):
+        del self.geo["throat.centroid"]
         self.geo['throat.length'] = 0.5
         self.geo.add_model(propname='throat.conduit_lengths',
                            model=mods.conduit_lengths,
@@ -41,6 +42,20 @@ def test_conduit_lengths(self):
         assert_approx_equal(actual=L2, desired=0.3)
         assert_approx_equal(actual=Lt, desired=0.5)
 
+    def test_conduit_lengths_w_centroid(self):
+        # Delete "throat.length" key, otherwise, it won't get re-calculated
+        del self.geo['throat.length']
+        self.geo["throat.centroid"] = np.array([[0, 0.5, 0.5]]) + self.base
+        self.geo.add_model(propname='throat.conduit_lengths',
+                           model=mods.conduit_lengths,
+                           regen_mode='normal')
+        L1 = self.geo['throat.conduit_lengths.pore1'][0]
+        L2 = self.geo['throat.conduit_lengths.pore2'][0]
+        Lt = self.geo['throat.conduit_lengths.throat'][0]
+        assert_approx_equal(actual=L1, desired=0.2)
+        assert_approx_equal(actual=L2, desired=0.3)
+        assert_approx_equal(actual=Lt, desired=1.12161167)
+
 
 if __name__ == '__main__':