Revamping conductance stack

openpnm/core/_base2.py

@@ -477,7 +477,10 @@ def get_conduit_data(self, propname):
         Parameters
         ----------
         propname : str
-            The dictionary key of the property to fetch.
+            The dictionary key of the property to fetch. Can be 'pore.diameter'
+            or just 'diameter'.  To fetch different propnames for pores and throats
+            supply them both in a list, like
+            `['pore.diameter', 'throat.inscribed_diameter']`.
 
         Returns
         -------
@@ -486,8 +489,13 @@ def get_conduit_data(self, propname):
             for pore1, throat, and pore2 respectively.
 
         """
-        poreprop = 'pore.' + propname.split('.', 1)[-1]
-        throatprop = 'throat.' + propname.split('.', 1)[-1]
+        if isinstance(propname, str):
+            poreprop = 'pore.' + propname.split('.', 1)[-1]
+            throatprop = 'throat.' + propname.split('.', 1)[-1]
+        elif isinstance(propname, list):
+            poreprop = propname[0] if propname[0].startswith('pore') else propname[1]
+            throatprop = \
+                propname[0] if propname[0].startswith('throat') else propname[1]
         conns = self.network.conns
         try:
             T = self[throatprop]

openpnm/models/geometry/conduit_lengths/_funcs.py

@@ -1,7 +1,7 @@
 import numpy as np
-
 from openpnm.models.geometry import _geodocs
 
+
 __all__ = [
     "spheres_and_cylinders",
     "circles_and_rectangles",
@@ -29,8 +29,6 @@ def spheres_and_cylinders(
     Calculates conduit lengths in the network assuming pores are spheres
     and throats are cylinders.
 
-    A conduit is defined as ( 1/2 pore - full throat - 1/2 pore ).
-
     Parameters
     ----------
     %(network)s
@@ -44,9 +42,16 @@ def spheres_and_cylinders(
         of the pore-throat-pore conduits. The array is formatted as
         ``[pore1, throat, pore2]``.
 
+    Notes
+    -----
+    The cylindrical throat is assumed to overlap the spherical pore bodies creating
+    a hemispherical cap at the area of intersection. The length of the cap is
+    assigned to the throat.
+
     """
     L_ctc = _get_L_ctc(network)
-    D1, Dt, D2 = network.get_conduit_data(pore_diameter.split(".", 1)[-1]).T
+    D1, Dt, D2 = \
+        network.get_conduit_data(propname=[pore_diameter, throat_diameter]).T
 
     # Handle the case where Dt > Dp
     if (Dt > D1).any() or (Dt > D2).any():
@@ -65,7 +70,9 @@ def spheres_and_cylinders(
 
 @_geodocs
 def circles_and_rectangles(
-    network, pore_diameter="pore.diameter", throat_diameter="throat.diameter"
+    network,
+    pore_diameter="pore.diameter",
+    throat_diameter="throat.diameter",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are circles
@@ -91,13 +98,14 @@ def circles_and_rectangles(
     return spheres_and_cylinders(
         network=network,
         pore_diameter=pore_diameter,
-        throat_diameter=throat_diameter,
     )
 
 
 @_geodocs
 def cones_and_cylinders(
-    network, pore_diameter="pore.diameter", throat_diameter="throat.diameter"
+    network,
+    pore_diameter="pore.diameter",
+    throat_diameter="throat.diameter",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are cones
@@ -116,7 +124,8 @@ def cones_and_cylinders(
 
     """
     L_ctc = _get_L_ctc(network)
-    D1, Dt, D2 = network.get_conduit_data(pore_diameter.split(".", 1)[-1]).T
+    D1, Dt, D2 = \
+        network.get_conduit_data(propname=[pore_diameter, throat_diameter]).T
 
     L1 = D1 / 2
     L2 = D2 / 2
@@ -131,7 +140,11 @@ def cones_and_cylinders(
 
 
 @_geodocs
-def intersecting_cones(network, pore_coords="pore.coords", throat_coords="throat.coords"):
+def intersecting_cones(
+    network,
+    pore_coords="pore.coords",
+    throat_coords="throat.coords",
+):
     r"""
     Calculates conduit lengths in the network assuming pores are cones
     that intersect. Therefore, the throat is the cross sectional plane where
@@ -160,7 +173,10 @@ def intersecting_cones(network, pore_coords="pore.coords", throat_coords="throat
 
 @_geodocs
 def hybrid_cones_and_cylinders(
-    network, pore_diameter="pore.diameter", throat_coords="throat.coords"
+    network,
+    pore_diameter="pore.diameter",
+    throat_diameter="throat.diameter",
+    throat_coords="throat.coords",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are cones
@@ -172,14 +188,16 @@ def hybrid_cones_and_cylinders(
     ----------
     %(network)s
     %(Dp)s
+    %(Dt)s
     %(Tcoords)s
 
     Returns
     -------
 
     """
     L_ctc = _get_L_ctc(network)
-    D1, Dt, D2 = network.get_conduit_data(pore_diameter.split(".", 1)[-1]).T
+    D1, Dt, D2 = \
+        network.get_conduit_data(propname=[pore_diameter, throat_diameter]).T
 
     L1 = D1 / 2
     L2 = D2 / 2
@@ -200,7 +218,9 @@ def hybrid_cones_and_cylinders(
 
 @_geodocs
 def trapezoids_and_rectangles(
-    network, pore_diameter="pore.diameter", throat_diameter="throat.diameter"
+    network,
+    pore_diameter="pore.diameter",
+    throat_diameter="throat.diameter",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are
@@ -230,7 +250,9 @@ def trapezoids_and_rectangles(
 
 @_geodocs
 def intersecting_trapezoids(
-    network, pore_coords="pore.coords", throat_coords="throat.coords"
+    network,
+    pore_coords="pore.coords",
+    throat_coords="throat.coords",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are
@@ -260,7 +282,9 @@ def intersecting_trapezoids(
 
 @_geodocs
 def hybrid_trapezoids_and_rectangles(
-    network, pore_diameter="pore.diameter", throat_coords="throat.coords"
+    network,
+    pore_diameter="pore.diameter",
+    throat_coords="throat.coords",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are
@@ -290,7 +314,9 @@ def hybrid_trapezoids_and_rectangles(
 
 @_geodocs
 def pyramids_and_cuboids(
-    network, pore_diameter="pore.diameter", throat_diameter="throat.diameter"
+    network,
+    pore_diameter="pore.diameter",
+    throat_diameter="throat.diameter",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are truncated
@@ -315,7 +341,9 @@ def pyramids_and_cuboids(
 
 @_geodocs
 def intersecting_pyramids(
-    network, pore_coords="pore.coords", throat_coords="throat.coords"
+    network,
+    pore_coords="pore.coords",
+    throat_coords="throat.coords",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are
@@ -340,7 +368,9 @@ def intersecting_pyramids(
 
 @_geodocs
 def hybrid_pyramids_and_cuboids(
-    network, pore_diameter="pore.diameter", throat_coords="throat.coords"
+    network,
+    pore_diameter="pore.diameter",
+    throat_coords="throat.coords",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are truncated
@@ -366,7 +396,9 @@ def hybrid_pyramids_and_cuboids(
 
 @_geodocs
 def cubes_and_cuboids(
-    network, pore_diameter="pore.diameter", throat_diameter="throat.diameter"
+    network,
+    pore_diameter="pore.diameter",
+    throat_diameter="throat.diameter",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are cubes
@@ -383,7 +415,8 @@ def cubes_and_cuboids(
 
     """
     L_ctc = _get_L_ctc(network)
-    D1, Dt, D2 = network.get_conduit_data(pore_diameter.split(".", 1)[-1]).T
+    D1, Dt, D2 = \
+        network.get_conduit_data(propname=[pore_diameter, throat_diameter]).T
 
     L1 = D1 / 2
     L2 = D2 / 2
@@ -400,7 +433,9 @@ def cubes_and_cuboids(
 
 @_geodocs
 def squares_and_rectangles(
-    network, pore_diameter="pore.diameter", throat_diameter="throat.diameter"
+    network,
+    pore_diameter="pore.diameter",
+    throat_diameter="throat.diameter",
 ):
     r"""
     Calculates conduit lengths in the network assuming pores are squares

tests/unit/models/physics/DiffusiveConductanceTest.py

@@ -120,6 +120,142 @@ def test_taylor_aris_diffusion(self):
         desired = np.array([0.121193, 0.126131, 0.118578])
         assert_allclose(actual, desired, rtol=1e-5)
 
+    def test_spheres_and_cylinders(self):
+        pn = op.network.Cubic(shape=[2, 1, 1])
+        Dp1, Dp2 = 0.5, 0.25
+        Dt = 0.1
+        pn['pore.diameter'] = [Dp1, Dp2]
+        pn['throat.diameter'] = Dt
+
+        f = op.models.geometry.diffusive_size_factors.spheres_and_cylinders
+        pn.add_model(propname='throat.diffusive_size_factors', model=f)
+        Scalc = pn['throat.diffusive_size_factors'].flatten()
+        # Pore 1
+        Lp1 = (Dp1/2)*np.cos(np.arcsin(Dt/Dp1))
+        Sp1 = 1/(np.pi*(Dp1/2))*np.arctanh(Lp1/(Dp1/2))
+        assert_allclose(1/Sp1, Scalc[0], rtol=1e-7)
+        # Pore 2
+        Lp2 = (Dp2/2)*np.cos(np.arcsin(Dt/Dp2))
+        Sp2 = 1/(np.pi*(Dp2/2))*np.arctanh(Lp2/(Dp2/2))
+        assert_allclose(1/Sp2, Scalc[2], rtol=1e-7)
+        # Throat
+        Lt = 1 - Lp1 - Lp2
+        St = Lt/(np.pi*(Dt/2)**2)
+        assert_allclose(1/St, Scalc[1], rtol=1e-6)
+        # Finally check conductance
+        f = op.models.physics.diffusive_conductance.generic_diffusive
+        phase = op.phase.Phase(network=pn)
+        phase['pore.diffusivity'] = 1.0
+        phase['throat.diffusivity'] = 1.0
+        phase.add_model(propname='throat.diffusive_conductance', model=f)
+        Gcalc = phase['throat.diffusive_conductance'][0]
+        G = (Sp1 + St + Sp2)**-1
+        assert_allclose(G, Gcalc, rtol=1e-7)
+
+    def test_spheres_and_cylinders_oversized(self):
+        pn = op.network.Cubic(shape=[2, 1, 1])
+        Dp1, Dp2 = 1.75, 1.25
+        Dt = 0.1
+        pn['pore.diameter'] = [Dp1, Dp2]
+        pn['throat.diameter'] = Dt
+
+        f = op.models.geometry.diffusive_size_factors.spheres_and_cylinders
+        pn.add_model(propname='throat.diffusive_size_factors', model=f)
+        Scalc = pn['throat.diffusive_size_factors'].flatten()
+        # Pore 1
+        Lp1 = (Dp1/2)*np.cos(np.arcsin(Dt/Dp1))
+        Sp1 = 1/(np.pi*(Dp1/2))*np.arctanh(Lp1/(Dp1/2))
+        assert_allclose(1/Sp1, Scalc[0], rtol=1e-7)
+        # Pore 2
+        Lp2 = (Dp2/2)*np.cos(np.arcsin(Dt/Dp2))
+        Sp2 = 1/(np.pi*(Dp2/2))*np.arctanh(Lp2/(Dp2/2))
+        assert_allclose(1/Sp2, Scalc[2], rtol=1e-7)
+        # Throat
+        Lt = 1 - Lp1 - Lp2
+        St = Lt/(np.pi*(Dt/2)**2)
+        assert_allclose(1/St, Scalc[1], rtol=1e-6)
+        # Finally check conductance
+        f = op.models.physics.diffusive_conductance.generic_diffusive
+        phase = op.phase.Phase(network=pn)
+        phase['pore.diffusivity'] = 1.0
+        phase['throat.diffusivity'] = 1.0
+        phase.add_model(propname='throat.diffusive_conductance', model=f)
+        Gcalc = phase['throat.diffusive_conductance'][0]
+        G = (Sp1 + St + Sp2)**-1
+        assert_allclose(G, Gcalc, rtol=1e-7)
+
+    def test_cones_and_cylinders(self):
+        pn = op.network.Cubic(shape=[2, 1, 1])
+        Dp1, Dp2 = 0.5, 0.25
+        Dt = 0.1
+        pn['pore.diameter'] = [Dp1, Dp2]
+        pn['throat.diameter'] = Dt
+
+        f = op.models.geometry.diffusive_size_factors.cones_and_cylinders
+        pn.add_model(propname='throat.diffusive_size_factors', model=f)
+        Scalc = pn['throat.diffusive_size_factors'].flatten()
+        # Pore 1
+        Sp1 = (Dp1/2)/(np.pi*(Dp1/2*Dt/2))
+        assert_allclose(1/Sp1, Scalc[0], rtol=1e-7)
+        # Pore 2
+        Sp2 = (Dp2/2)/(np.pi*(Dp2/2*Dt/2))
+        assert_allclose(1/Sp2, Scalc[2], rtol=1e-7)
+        # Throat
+        Lt = 1 - Dp1/2 - Dp2/2
+        St = Lt/(np.pi*(Dt/2)**2)
+        assert_allclose(1/St, Scalc[1], rtol=1e-6)
+        # Finally check conductance
+        f = op.models.physics.diffusive_conductance.generic_diffusive
+        phase = op.phase.Phase(network=pn)
+        phase['pore.diffusivity'] = 1.0
+        phase['throat.diffusivity'] = 1.0
+        phase.add_model(propname='throat.diffusive_conductance', model=f)
+        Gcalc = phase['throat.diffusive_conductance'][0]
+        G = (Sp1 + St + Sp2)**-1
+        assert_allclose(G, Gcalc, rtol=1e-7)
+
+    def test_intersecting_cones(self):
+        pn = op.network.Cubic(shape=[2, 1, 1])
+        Dp1, Dp2 = 0.5, 0.25
+        Dt = 0.1
+        pn['pore.diameter'] = [Dp1, Dp2]
+        pn['throat.diameter'] = Dt
+        pn['throat.coords'] = [[1.0, 0.5, 0.5]]
+
+        f = op.models.geometry.diffusive_size_factors.intersecting_cones
+        pn.add_model(propname='throat.diffusive_size_factors', model=f)
+        Scalc = pn['throat.diffusive_size_factors'].flatten()
+        # Pore 1
+        Lp1 = 0.5
+        Sp1 = (Lp1)/(np.pi*(Dp1/2*Dt/2))
+        assert_allclose(1/Sp1, Scalc[0], rtol=1e-7)
+        # Pore 2
+        Lp2 = 0.5
+        Sp2 = (Lp2)/(np.pi*(Dp2/2*Dt/2))
+        assert_allclose(1/Sp2, Scalc[2], rtol=1e-7)
+        # Throat
+        St = 0
+        # Finally check conductance
+        f = op.models.physics.diffusive_conductance.generic_diffusive
+        phase = op.phase.Phase(network=pn)
+        phase['pore.diffusivity'] = 1.0
+        phase['throat.diffusivity'] = 1.0
+        phase.add_model(propname='throat.diffusive_conductance', model=f)
+        Gcalc = phase['throat.diffusive_conductance'][0]
+        G = (Sp1 + St + Sp2)**-1
+        assert_allclose(G, Gcalc, rtol=1e-7)
+
+
+
+
+
+
+
+
+
+
+
+
 
 if __name__ == '__main__':