<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Tutorial 2 of 3
Step1:  Let's generate a cubic network again, but with a different connectivity
Step2:  This Network has pores distributed in a cubic lattice, but connected to diagonal neighbors due to the connectivity being set to 8 (the default is 6 which is orthogonal neighbors).  The various options are outlined in the Cubic class's documentation which can be viewed with the Object Inspector in Spyder.
OpenPNM includes several other classes for generating networks including random topology based on Delaunay tessellations (Delaunay).
It is also possible to import networks &lt;data_io&gt;_ from external sources, such as networks extracted from tomographic images, or that networks generated by external code.
Initialize and Build Multiple Geometry Objects
One of the main functionalities of OpenPNM is the ability to assign drastically different geometrical properties to different regions of the domain to create heterogeneous materials, such as layered structures.  To demonstrate the motivation behind this feature, this tutorial will make a material that has different geometrical properties on the top and bottom surfaces compared to the internal pores.  We need to create one Geometry object to manage the top and bottom pores, and a second to manage the remaining internal pores
Step3:  The above statements result in two distinct Geometry objects, each applying to different regions of the domain.  geom1 applies to only the pores on the top and bottom surfaces (automatically labeled 'top' and 'bottom' during the network generation step), while geom2 applies to the pores 'not' on the top and bottom surfaces.
The assignment of throats is more complicated and illustrates the find_neighbor_throats method, which is one of the more useful topological query methods &lt;topology&gt;_ on the Network class.  In both of these calls, all throats connected to the given set of pores (Ps1 or Ps2) are found; however, the mode argument alters which throats are returned.  The terms 'union' and 'intersection' are used in the "set theory" sense, such that 'union' returns all throats connected to the pores in the supplied list, while 'intersection' returns the throats that are only connected to the supplied pores.  More specifically, if pores 1 and 2 have throats [1, 2] and [2, 3] as neighbors, respectively, then the 'union' mode returns [1, 2, 3] and the 'intersection' mode returns [2].  A detailed description of this behavior is given in 
Step4:  Each of the above lines produced an array of different length, corresponding to the number of pores assigned to each Geometry object.  This is accomplished by the calls to geom1.Np and geom2.Np, which return the number of pores on each object.
Every Core object in OpenPNM possesses the same set of methods for managing their data, such as counting the number of pore and throat values they represent; thus, pn.Np returns 1000 while geom1.Np and geom2.Np return 200 and 800 respectively.
Accessing Data Distributed Between Geometries
The segmentation of the data between separate Geometry objects is essential to the management of pore-scale models, although it does create a complication
Step5:  The following code illustrates the shortcut approach, which accomplishes the same result as above in a single line
Step6:  This shortcut works because the pn dictionary does not contain an array called 'pore.seed', so all associated Geometry objects are then checked for the requested array(s).  If it is found, then OpenPNM essentially performs the interleaving of the data as demonstrated by the manual approach and returns all the values together in a single full-size array.  If it is not found, then a standard KeyError message is received.
This exchange of data between Network and Geometry makes sense if you consider that Network objects act as a sort of master object relative Geometry objects.  Networks apply to all pores and throats in the domain, while Geometries  apply to subsets of the domain, so if the Network needs some values from all pores it has direct access.
Add Pore Size Distribution Models to Each Geometry
Pore-scale models are mathematical functions that are applied to each pore (or throat) in the network to produce some local property value.  Each of the modules in OpenPNM (Network, Geometry, Phase and Physics) have a "library" of pre-written models located under "models" (i.e. Geometry.models).  Below this level, the models are further categorized according to what property they calculate, and there are typical 2-3 models for each.  For instance, under Geometry.models.pore_diameter you will see random, normal and weibull among others.
Pore size distribution models are assigned to each Geometry object as follows
Step7:  Pore-scale models tend to be the most complex (i.e. confusing) aspects of OpenPNM, so it's worth dwelling on the important points of the above two commands
Step8:  Instead of using statistical distribution functions, the above lines use the neighbor model which determines each throat value based on the values found 'pore_prop' from it's neighboring pores.  In this case, each throat is assigned the minimum pore diameter of it's two neighboring pores.  Other options for mode include 'max' and 'mean'.
We'll also need throat length as well as the cross-sectional area of pores and throats, for calculating the hydraulic conductance model later.
Step9:  Create a Phase Object and Assign Thermophysical Property Models
For this tutorial, we will create a generic Phase object for water, then assign some pore-scale models for calculating their properties.  Alternatively, we could use the prewritten Water class included in OpenPNM, which comes complete with the necessary pore-scale models, but this would defeat the purpose of the tutorial.
Step10:  Note that all Phase objects are automatically assigned standard temperature and pressure conditions when created.  This can be adjusted
Step11:  A variety of pore-scale models are available for calculating Phase properties, generally taken from correlations in the literature.  An empirical correlation specifically for the viscosity of water is available
Step12:  Create Physics Objects for Each Geometry
Physics objects are where geometric information and thermophysical properties are combined to produce the pore and throat scale transport parameters.  Thus we need to create one Physics object for EACH Phase and EACH Geometry
Step13:  Next add the Hagan-Poiseuille model to both
Step14:  The same function (mod) was passed as the model argument to both Physics objects.  This means that both objects will calculate the hydraulic conductance using the same function.  A model must be assigned to both objects in order for the 'throat.hydraulic_conductance' property be defined everywhere in the domain since each Physics applies to a unique selection of pores and throats.
The "pore-scale model" mechanism was specifically designed to allow for users to easily create their own custom models.  Creating custom models is outlined in 
Step15:  Each Physics applies to the same subset for pores and throats as the Geometries so its values are distributed spatially, but each Physics is also associated with a single Phase object.  Consequently, only a Phase object can to request all of the values within the domain pertaining to itself.
In other words, a Network object cannot aggregate the Physics data because it doesn't know which Phase is referred to.  For instance, when asking for 'throat.hydraulic_conductance' it could refer to water or air conductivity, so it can only be requested by water or air.
Pore-Scale Models
Step16:  Now, let's alter the Geometry objects by assigning new random seeds, and adjust the temperature of water.
Step17:  So far we have not run the regenerate command on any of these objects, which means that the above changes have not yet been applied to all the dependent properties.  Let's do this and examine what occurs at each step
Step18:  These two lines trigger the re-calculation of all the size related models on each Geometry object.
Step19:  This line causes the viscosity to be recalculated at the new temperature. Let's confirm that the hydraulic conductance has NOT yet changed since we have not yet regenerated the Physics objects' models
Step20:  Finally, if we regenerate phys1 and phys2 we can see that the hydraulic conductance will be updated to reflect the new sizes on the Geometries and the new temperature on the Phase
Step21:  Determine Permeability Tensor by Changing Inlet and Outlet Boundary Conditions
The 
Step22:  Set boundary conditions for flow in the X-direction
Step23:  The resulting pressure field can be seen using Paraview
Step24:  To find K, we need to solve Darcy's law
Step25:  The dimensions of the network can be determined manually from the shape and spacing specified during its generation
Step26:  The pressure drop was specified as 1 atm when setting boundary conditions, so Kxx can be found as
Step27:  We can either create 2 new Algorithm objects to perform the simulations in the other two directions, or reuse alg by adjusting the boundary conditions and re-running it.
Step28:  The first call to set_boundary_conditions used the overwrite mode, which replaces all existing boundary conditions on the alg object with the specified values.  The second call uses the merge mode which adds new boundary conditions to any already present, which is the default behavior.
A new value for the flow rate must be recalculated, but all other parameters are equal to the X-direction
Step29:  The values of Kxx and Kyy should be nearly identical since both these two directions are parallel to the small surface pores.  For the Z-direction
Step30:  The permeability in the Z-direction is about half that in the other two directions due to the constrictions caused by the small surface pores.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import scipy as sp
import openpnm as op
np.random.seed(10)
ws = op.Workspace()
ws.settings["loglevel"] = 40
Explanation: Tutorial 2 of 3: Digging Deeper into OpenPNM
This tutorial will follow the same outline as Getting Started, but will dig a little bit deeper at each step to reveal the important features of OpenPNM that were glossed over previously.
Learning Objectives
Explore different network topologies, and learn some handy topological query methods
Create a heterogeneous domain with different geometrical properties in different regions
Learn about data exchange between objects
Utilize pore-scale models for calculating properties of all types
Propagate changing geometrical and thermo-physical properties to all dependent properties
Calculate the permeability tensor for the stratified media
Use the Workspace Manager to save and load a simulation
Building a Cubic Network
As usual, start by importing the OpenPNM and Scipy packages:
End of explanation
pn = op.network.Cubic(shape=[20, 20, 10], spacing=0.0001, connectivity=8)
Explanation: Let's generate a cubic network again, but with a different connectivity:
End of explanation
Ps1 = pn.pores(['top', 'bottom'])
Ts1 = pn.find_neighbor_throats(pores=Ps1, mode='union')
geom1 = op.geometry.GenericGeometry(network=pn, pores=Ps1, throats=Ts1, name='boundaries')
Ps2 = pn.pores(['top', 'bottom'], mode='not')
Ts2 = pn.find_neighbor_throats(pores=Ps2, mode='xnor')
geom2 = op.geometry.GenericGeometry(network=pn, pores=Ps2, throats=Ts2, name='core')
Explanation: This Network has pores distributed in a cubic lattice, but connected to diagonal neighbors due to the connectivity being set to 8 (the default is 6 which is orthogonal neighbors).  The various options are outlined in the Cubic class's documentation which can be viewed with the Object Inspector in Spyder.
OpenPNM includes several other classes for generating networks including random topology based on Delaunay tessellations (Delaunay).
It is also possible to import networks &lt;data_io&gt;_ from external sources, such as networks extracted from tomographic images, or that networks generated by external code.
Initialize and Build Multiple Geometry Objects
One of the main functionalities of OpenPNM is the ability to assign drastically different geometrical properties to different regions of the domain to create heterogeneous materials, such as layered structures.  To demonstrate the motivation behind this feature, this tutorial will make a material that has different geometrical properties on the top and bottom surfaces compared to the internal pores.  We need to create one Geometry object to manage the top and bottom pores, and a second to manage the remaining internal pores:
End of explanation
geom1['pore.seed'] = np.random.rand(geom1.Np)*0.5 + 0.2
geom2['pore.seed'] = np.random.rand(geom2.Np)*0.5 + 0.2
Explanation: The above statements result in two distinct Geometry objects, each applying to different regions of the domain.  geom1 applies to only the pores on the top and bottom surfaces (automatically labeled 'top' and 'bottom' during the network generation step), while geom2 applies to the pores 'not' on the top and bottom surfaces.
The assignment of throats is more complicated and illustrates the find_neighbor_throats method, which is one of the more useful topological query methods &lt;topology&gt;_ on the Network class.  In both of these calls, all throats connected to the given set of pores (Ps1 or Ps2) are found; however, the mode argument alters which throats are returned.  The terms 'union' and 'intersection' are used in the "set theory" sense, such that 'union' returns all throats connected to the pores in the supplied list, while 'intersection' returns the throats that are only connected to the supplied pores.  More specifically, if pores 1 and 2 have throats [1, 2] and [2, 3] as neighbors, respectively, then the 'union' mode returns [1, 2, 3] and the 'intersection' mode returns [2].  A detailed description of this behavior is given in :ref:topology.
Assign Static Seed values to Each Geometry
In :ref:getting_started we only assigned 'static' values to the Geometry object, which we calculated explicitly.  In this tutorial we will use the pore-scale models that are provided with OpenPNM.  To get started, however, we'll assign static random seed values between 0 and 1 to each pore on both Geometry objects, by assigning random numbers to each Geometry's 'pore.seed' property:
End of explanation
seeds = np.zeros_like(pn.Ps, dtype=float)
seeds[pn.pores(geom1.name)] = geom1['pore.seed']
seeds[pn.pores(geom2.name)] = geom2['pore.seed']
print(np.all(seeds > 0))  # Ensure all zeros are overwritten
Explanation: Each of the above lines produced an array of different length, corresponding to the number of pores assigned to each Geometry object.  This is accomplished by the calls to geom1.Np and geom2.Np, which return the number of pores on each object.
Every Core object in OpenPNM possesses the same set of methods for managing their data, such as counting the number of pore and throat values they represent; thus, pn.Np returns 1000 while geom1.Np and geom2.Np return 200 and 800 respectively.
Accessing Data Distributed Between Geometries
The segmentation of the data between separate Geometry objects is essential to the management of pore-scale models, although it does create a complication: it's not easy to obtain a single array containing all the values of a given property for the whole network.  It is technically possible to piece this data together manually since we know the locations where each Geometry object applies, but this is tedious so OpenPNM provides a shortcut.  First, let's illustrate the manual approach using the 'pore.seed' values we have defined:
End of explanation
seeds = pn['pore.seed']
Explanation: The following code illustrates the shortcut approach, which accomplishes the same result as above in a single line:
End of explanation
geom1.add_model(propname='pore.diameter',
                model=op.models.geometry.pore_size.normal,
                scale=0.00001, loc=0.00005,
                seeds='pore.seed')
geom2.add_model(propname='pore.diameter',
                model=op.models.geometry.pore_size.weibull,
                shape=1.2, scale=0.00001, loc=0.00005,
                seeds='pore.seed')
Explanation: This shortcut works because the pn dictionary does not contain an array called 'pore.seed', so all associated Geometry objects are then checked for the requested array(s).  If it is found, then OpenPNM essentially performs the interleaving of the data as demonstrated by the manual approach and returns all the values together in a single full-size array.  If it is not found, then a standard KeyError message is received.
This exchange of data between Network and Geometry makes sense if you consider that Network objects act as a sort of master object relative Geometry objects.  Networks apply to all pores and throats in the domain, while Geometries  apply to subsets of the domain, so if the Network needs some values from all pores it has direct access.
Add Pore Size Distribution Models to Each Geometry
Pore-scale models are mathematical functions that are applied to each pore (or throat) in the network to produce some local property value.  Each of the modules in OpenPNM (Network, Geometry, Phase and Physics) have a "library" of pre-written models located under "models" (i.e. Geometry.models).  Below this level, the models are further categorized according to what property they calculate, and there are typical 2-3 models for each.  For instance, under Geometry.models.pore_diameter you will see random, normal and weibull among others.
Pore size distribution models are assigned to each Geometry object as follows:
End of explanation
geom1.add_model(propname='throat.diameter',
                model=op.models.misc.from_neighbor_pores,
                pore_prop='pore.diameter',
                mode='min')
geom2.add_model(propname='throat.diameter',
                model=op.models.misc.from_neighbor_pores,
                mode='min')
pn['pore.diameter'][pn['throat.conns']]
Explanation: Pore-scale models tend to be the most complex (i.e. confusing) aspects of OpenPNM, so it's worth dwelling on the important points of the above two commands:
Both geom1 and geom2 have a models attribute where the parameters specified in the add command are stored for future use if/when needed.  The models attribute actually contains a ModelsDict object which is a customized dictionary for storing and managing this type of information.
The propname argument specifies which property the model calculates.  This means that the numerical results of the model calculation will be saved in their respective Geometry objects as geom1['pore.diameter'] and geom2['pore.diameter'].
Each model stores it's result under the same propname but these values do not conflict since each Geometry object presides over a unique subset of pores and throats.
The model argument contains a handle to the desired function, which is extracted from the models library of the relevant Module (Geometry in this case).  Each Geometry object has been assigned a different statistical model, normal and weibull.  This ability to apply different models to different regions of the domain is reason multiple Geometry objects are permitted.  The added complexity is well worth the added flexibility.
The remaining arguments are those required by the chosen model.  In the above cases, these are the parameters that define the statistical distribution.  Note that the mean pore size for geom1 will be 20 um (set by scale) while for geom2 it will be 50 um, thus creating the smaller surface pores as intended.  The pore-scale models are well documented regarding what arguments are required and their meaning; as usual these can be viewed with Object Inspector in Spyder.
Now that we've added pore diameter models the each Geometry we can visualize the network in Paraview to confirm that distinctly different pore sizes on the surface regions:
<img src="http://i.imgur.com/5F70ens.png" style="width: 60%" align="left"/>
Add Additional Pore-Scale Models to Each Geometry
In addition to pore diameter, there are several other geometrical properties needed to perform a permeability simulation.  Let's start with throat diameter:
End of explanation
geom1.add_model(propname='throat.endpoints',
                model=op.models.geometry.throat_endpoints.spherical_pores)
geom2.add_model(propname='throat.endpoints',
                model=op.models.geometry.throat_endpoints.spherical_pores)
geom1.add_model(propname='throat.area',
                model=op.models.geometry.throat_area.cylinder)
geom2.add_model(propname='throat.area',
                model=op.models.geometry.throat_area.cylinder)
geom1.add_model(propname='pore.area',
                model=op.models.geometry.pore_area.sphere)
geom2.add_model(propname='pore.area',
                model=op.models.geometry.pore_area.sphere)
geom1.add_model(propname='throat.conduit_lengths',
                model=op.models.geometry.throat_length.conduit_lengths)
geom2.add_model(propname='throat.conduit_lengths',
                model=op.models.geometry.throat_length.conduit_lengths)
Explanation: Instead of using statistical distribution functions, the above lines use the neighbor model which determines each throat value based on the values found 'pore_prop' from it's neighboring pores.  In this case, each throat is assigned the minimum pore diameter of it's two neighboring pores.  Other options for mode include 'max' and 'mean'.
We'll also need throat length as well as the cross-sectional area of pores and throats, for calculating the hydraulic conductance model later.
End of explanation
water = op.phases.GenericPhase(network=pn)
air = op.phases.GenericPhase(network=pn)
Explanation: Create a Phase Object and Assign Thermophysical Property Models
For this tutorial, we will create a generic Phase object for water, then assign some pore-scale models for calculating their properties.  Alternatively, we could use the prewritten Water class included in OpenPNM, which comes complete with the necessary pore-scale models, but this would defeat the purpose of the tutorial.
End of explanation
water['pore.temperature'] = 353  # K
Explanation: Note that all Phase objects are automatically assigned standard temperature and pressure conditions when created.  This can be adjusted:
End of explanation
water.add_model(propname='pore.viscosity',
                model=op.models.phases.viscosity.water)
Explanation: A variety of pore-scale models are available for calculating Phase properties, generally taken from correlations in the literature.  An empirical correlation specifically for the viscosity of water is available:
End of explanation
phys1 = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom1)
phys2 = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom2)
Explanation: Create Physics Objects for Each Geometry
Physics objects are where geometric information and thermophysical properties are combined to produce the pore and throat scale transport parameters.  Thus we need to create one Physics object for EACH Phase and EACH Geometry:
End of explanation
mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
phys1.add_model(propname='throat.hydraulic_conductance', model=mod)
phys2.add_model(propname='throat.hydraulic_conductance', model=mod)
Explanation: Next add the Hagan-Poiseuille model to both:
End of explanation
g = water['throat.hydraulic_conductance']
Explanation: The same function (mod) was passed as the model argument to both Physics objects.  This means that both objects will calculate the hydraulic conductance using the same function.  A model must be assigned to both objects in order for the 'throat.hydraulic_conductance' property be defined everywhere in the domain since each Physics applies to a unique selection of pores and throats.
The "pore-scale model" mechanism was specifically designed to allow for users to easily create their own custom models.  Creating custom models is outlined in :ref:advanced_usage.
Accessing Data Distributed Between Multiple Physics Objects
Just as Network objects can retrieve data from separate Geometries as a single array with values in the correct locations, Phase objects can retrieve data from Physics objects as follows:
End of explanation
g1 = phys1['throat.hydraulic_conductance']  # Save this for later
g2 = phys2['throat.hydraulic_conductance']  # Save this for later
Explanation: Each Physics applies to the same subset for pores and throats as the Geometries so its values are distributed spatially, but each Physics is also associated with a single Phase object.  Consequently, only a Phase object can to request all of the values within the domain pertaining to itself.
In other words, a Network object cannot aggregate the Physics data because it doesn't know which Phase is referred to.  For instance, when asking for 'throat.hydraulic_conductance' it could refer to water or air conductivity, so it can only be requested by water or air.
Pore-Scale Models: The Big Picture
Having created all the necessary objects with pore-scale models, it is now time to demonstrate why the OpenPNM pore-scale model approach is so powerful.  First, let's inspect the current value of hydraulic conductance in throat 1 on phys1 and phys2:
End of explanation
geom1['pore.seed'] = np.random.rand(geom1.Np)
geom2['pore.seed'] = np.random.rand(geom2.Np)
water['pore.temperature'] = 370  # K
Explanation: Now, let's alter the Geometry objects by assigning new random seeds, and adjust the temperature of water.
End of explanation
geom1.regenerate_models()
geom2.regenerate_models()
Explanation: So far we have not run the regenerate command on any of these objects, which means that the above changes have not yet been applied to all the dependent properties.  Let's do this and examine what occurs at each step:
End of explanation
water.regenerate_models()
Explanation: These two lines trigger the re-calculation of all the size related models on each Geometry object.
End of explanation
print(np.all(phys1['throat.hydraulic_conductance'] == g1))  # g1 was saved above
print(np.all(phys2['throat.hydraulic_conductance'] == g2) ) # g2 was saved above
Explanation: This line causes the viscosity to be recalculated at the new temperature. Let's confirm that the hydraulic conductance has NOT yet changed since we have not yet regenerated the Physics objects' models:
End of explanation
phys1.regenerate_models()
phys2.regenerate_models()
print(np.all(phys1['throat.hydraulic_conductance'] != g1))
print(np.all(phys2['throat.hydraulic_conductance'] != g2))
Explanation: Finally, if we regenerate phys1 and phys2 we can see that the hydraulic conductance will be updated to reflect the new sizes on the Geometries and the new temperature on the Phase:
End of explanation
alg = op.algorithms.StokesFlow(network=pn, phase=water)
Explanation: Determine Permeability Tensor by Changing Inlet and Outlet Boundary Conditions
The :ref:getting started tutorial &lt;getting_started&gt; already demonstrated the process of performing a basic permeability simulation.  In this tutorial, we'll perform the simulation in all three perpendicular dimensions to obtain the permeability tensor of our heterogeneous anisotropic material.
End of explanation
alg.set_value_BC(values=202650, pores=pn.pores('right'))
alg.set_value_BC(values=101325, pores=pn.pores('left'))
alg.run()
Explanation: Set boundary conditions for flow in the X-direction:
End of explanation
Q = alg.rate(pores=pn.pores('right'))
Explanation: The resulting pressure field can be seen using Paraview:
<img src="http://i.imgur.com/ugX0LFG.png" style="width: 60%" align="left"/>
To determine the permeability coefficient we must find the flow rate through the network to use in Darcy's law.  The StokesFlow class (and all analogous transport algorithms) possess a rate method that calculates the net transport through a given set of pores:
End of explanation
mu = np.mean(water['pore.viscosity'])
Explanation: To find K, we need to solve Darcy's law: Q = KA/(mu*L)(P_in - P_out).  This requires knowing the viscosity and macroscopic network dimensions:
End of explanation
L = 20 * 0.0001
A = 20 * 10 * (0.0001**2)
Explanation: The dimensions of the network can be determined manually from the shape and spacing specified during its generation:
End of explanation
Kxx = Q * mu * L / (A * 101325)
Explanation: The pressure drop was specified as 1 atm when setting boundary conditions, so Kxx can be found as:
End of explanation
alg.set_value_BC(values=202650, pores=pn.pores('front'))
alg.set_value_BC(values=101325, pores=pn.pores('back'))
alg.run()
Explanation: We can either create 2 new Algorithm objects to perform the simulations in the other two directions, or reuse alg by adjusting the boundary conditions and re-running it.
End of explanation
Q = alg.rate(pores=pn.pores('front'))
Kyy = Q * mu * L / (A * 101325)
Explanation: The first call to set_boundary_conditions used the overwrite mode, which replaces all existing boundary conditions on the alg object with the specified values.  The second call uses the merge mode which adds new boundary conditions to any already present, which is the default behavior.
A new value for the flow rate must be recalculated, but all other parameters are equal to the X-direction:
End of explanation
alg.set_value_BC(values=202650, pores=pn.pores('top'))
alg.set_value_BC(values=101325, pores=pn.pores('bottom'))
alg.run()
Q = alg.rate(pores=pn.pores('top'))
L = 10 * 0.0001
A = 20 * 20 * (0.0001**2)
Kzz = Q * mu * L / (A * 101325)
Explanation: The values of Kxx and Kyy should be nearly identical since both these two directions are parallel to the small surface pores.  For the Z-direction:
End of explanation
print(Kxx, Kyy, Kzz)
Explanation: The permeability in the Z-direction is about half that in the other two directions due to the constrictions caused by the small surface pores.
End of explanation
<END_TASK>