Error "Must specify either pores or throats"

[HFenice]

Hi,
I'm trying to compute absolute permeability but I have this error "Must specify either pores or throats". Here is my code:

# Importing library

import math
import openpnm as op
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as spst
from scipy.optimize import curve_fit
from scipy.spatial import distance_matrix
import csv
import pandas as pd
import random
import openpnm.models as mods
from openpnm.geometry import GenericGeometry
from openpnm.io import PerGeos

# Importing data from Avizo

nodes = pd.read_csv("dati_puliti.csv", sep = ";", encoding = "utf-8")
cx = nodes["x"].values.reshape(nodes.shape[0], 1)
cy = nodes["y"].values.reshape(nodes.shape[0], 1)
cz = nodes["z"].values.reshape(nodes.shape[0], 1)
diametro = nodes["d"].values.reshape(nodes.shape[0], 1)

# Network topology

pts = np.hstack([cx,cy,cz])
pn = op.network.Gabriel(points = pts,shape=[5000,5000,5000])  # [4879,5032,4569])
fig = op.topotools.plot_connections(pn)

# Checking network health

h = pn.check_network_health()
print(h)

# My Geometry object by modifing StickAndBall class

class MyGeo(GenericGeometry):
    def __init__(self, random_state = None, **kwargs):
        super().__init__(**kwargs)
        
        self['pore.diameter'] = diametro.ravel()
        
        self.add_model(propname='pore.area',
                       model=mods.geometry.pore_cross_sectional_area.sphere,
                       pore_diameter='pore.diameter') 

        self.add_model(propname='pore.volume',
                       model=mods.geometry.pore_volume.sphere, 
                       pore_diameter='pore.diameter')
        
        self.add_model(propname='throat.max_size',
                       model=mods.misc.from_neighbor_pores,
                       mode='min',
                       prop='pore.diameter')

        self.add_model(propname='throat.diameter',
                       model=mods.misc.scaled,
                       factor=0.5,
                       prop='throat.max_size')
                
        self.add_model(propname='throat.endpoints',
                       model=mods.geometry.throat_endpoints.spherical_pores,
                       pore_diameter='pore.diameter',
                       throat_diameter='throat.diameter')

        self.add_model(propname='throat.length',
                       model=mods.geometry.throat_length.piecewise,          
                       throat_endpoints='throat.endpoints')

        self.add_model(propname='throat.surface_area',
                       model=mods.geometry.throat_surface_area.cylinder, 
                       throat_diameter='throat.diameter',
                       throat_length='throat.length')

        self.add_model(propname='throat.volume',
                       model=mods.geometry.throat_volume.cylinder, 
                       throat_diameter='throat.diameter',
                       throat_length='throat.length')

        self.add_model(propname='throat.area',
                       model=mods.geometry.throat_cross_sectional_area.cylinder, 
                       throat_diameter='throat.diameter')

        self.add_model(propname='throat.conduit_lengths',
                       model=mods.geometry.throat_length.conduit_lengths,  
                       throat_endpoints='throat.endpoints',
                       throat_length='throat.length')

geo = MyGeo(network=pn, pores=pn.Ps, throats=pn.Ts, random_state = 123) 
geo.show_hist('pore.diameter')
geo.show_hist('throat.diameter')

water = op.phases.Water(network=pn) # Phase Object

phys_water =op.physics.Standard(network=pn, phase=water
                                      ,geometry=geo)

# Importing data in Avizo by using PerGeos Class

PerGeos.export_data(network=pn, phases=[water], filename='prova')

# Perform Stokes flow to find Permeability coefficient

sf_z = op.algorithms.StokesFlow(network=pn, phase=water) 
Pin = 0
Pout = 101325
sf_z.set_value_BC(pores=pn.pores('top'), values=Pin)
sf_z.set_value_BC(pores=pn.pores('bottom'), values=Pout)
sf_z.run()
A = (np.pi*((4*10e3)**2)/4)
L = 4*10e3 
mu = water['pore.viscosity'][0] 
Q = sf_z.rate(pores=pn.pores('top'), mode='group') 
Kz = (Q*L*mu/(A*(Pin - Pout))) 
print('The permeability coefficient is:', Kz)

Any ideas?

[ma-sadeghi]

Could you also paste the entire error block you're getting (so we can pinpoint the cause)?

[HFenice]

Yes, the labels 'top' and 'bottom' are empty. How can I fill them?

[jgostick]

You need to identify the pores first. You can use the simple method: Ps = pn.coords[:, 0] < thresh, where thresh is some value that you consider to be on on the left face. You'd then do pn.set_label(pores=Ps, label='left') or pn['pore.left'] = False followed by pn['pore.left][Ps] = True. Or you could try the "find_surface_pores" method which is described here. You can specify where the 'markers' are located so you can use this to 'only' find the left pores if you 'only' pass in one markers. Hopefully the docs and example are sufficient for you. Let me/us know if not.

[HFenice]

Many thanks! I've tried and now it works. Just a question: in this way I' computing the hydraulic conductivity of my network along z direction, is it right?

[HFenice]

Can I use also this?

perm = op.algorithms.StokesFlow(network=pn)
perm.setup(phase=water)
perm.set_value_BC(pores=pn.pores('top'), values=0)
perm.set_value_BC(pores=pn.pores('bottom'), values=101325)
perm.run()
perm.domain_area = (np.pi*((4*10e3)**2)/4)
perm.domain_length = 4000
perm.run()
K = perm.calc_effective_permeability()

[ma-sadeghi]

@HFenice just a quick note: calc_effective_permeability is going to be deprecated in v3 (to be released before 2022, hopefully). It's going to be re-packaged under openpnm.metrics module.