From 8d14176b52142dc1a4b1d4568e55dd72a8177200 Mon Sep 17 00:00:00 2001 From: jgostick Date: Tue, 11 Apr 2023 11:18:10 -0400 Subject: [PATCH 1/3] lotsa tests --- .../physics/DiffusiveConductanceTest.py | 136 ++++++++++++++++++ 1 file changed, 136 insertions(+) diff --git a/tests/unit/models/physics/DiffusiveConductanceTest.py b/tests/unit/models/physics/DiffusiveConductanceTest.py index 5a114bc0db..660c857daa 100644 --- a/tests/unit/models/physics/DiffusiveConductanceTest.py +++ b/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__': From ded4a0f2074b536216781cf8def0af982e998552 Mon Sep 17 00:00:00 2001 From: jgostick Date: Fri, 14 Apr 2023 23:40:36 -0400 Subject: [PATCH 2/3] tidying conduit length formatting --- openpnm/core/_base2.py | 14 +++- .../models/geometry/conduit_lengths/_funcs.py | 72 ++++++++++++++----- 2 files changed, 64 insertions(+), 22 deletions(-) diff --git a/openpnm/core/_base2.py b/openpnm/core/_base2.py index d2970e6fbb..ca6f320888 100644 --- a/openpnm/core/_base2.py +++ b/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] diff --git a/openpnm/models/geometry/conduit_lengths/_funcs.py b/openpnm/models/geometry/conduit_lengths/_funcs.py index 3355ba44fd..f77fe41238 100644 --- a/openpnm/models/geometry/conduit_lengths/_funcs.py +++ b/openpnm/models/geometry/conduit_lengths/_funcs.py @@ -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,15 @@ 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. + """ 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 +69,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 +97,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 +123,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 +139,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 +172,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,6 +187,7 @@ def hybrid_cones_and_cylinders( ---------- %(network)s %(Dp)s + %(Dt)s %(Tcoords)s Returns @@ -179,7 +195,8 @@ def hybrid_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 @@ -200,7 +217,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 +249,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 +281,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 +313,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 +340,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 +367,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 +395,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 +414,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 +432,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 From 831da2dd36132e3f80717d45fe2b4697da238a1c Mon Sep 17 00:00:00 2001 From: jgostick Date: Sun, 16 Apr 2023 22:02:36 -0400 Subject: [PATCH 3/3] docs --- openpnm/models/geometry/conduit_lengths/_funcs.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openpnm/models/geometry/conduit_lengths/_funcs.py b/openpnm/models/geometry/conduit_lengths/_funcs.py index f77fe41238..0a26d4878a 100644 --- a/openpnm/models/geometry/conduit_lengths/_funcs.py +++ b/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", @@ -45,7 +45,8 @@ def spheres_and_cylinders( Notes ----- The cylindrical throat is assumed to overlap the spherical pore bodies creating - a hemispherical cap at the area of intersection. + a hemispherical cap at the area of intersection. The length of the cap is + assigned to the throat. """ L_ctc = _get_L_ctc(network)