Exception: Your network is clustered, making Ax = b ill-conditioned

[JorgeESantos]

Still trying to measure the permeability of some 3d lattices I've been working with. This time, using OpenPNM, v3 I'm getting the following error, after running the flow flow.run().

from functions import Trigo_scaffolds, Voxel_image, TPMS_engine, scatter_3Ddataset
import numpy as np
import porespy as ps
import openpnm as pnm

#Generating Binary image of a porous scaffold
#scaffold, voxel_lenght = Trigo_scaffolds(Geo='b', shape=[2,2,2], elem=27, meshdata=False, min_dimen=13e-3)
scaffold, voxel_lenght = TPMS_engine([2,2,2], elem=40, min_dimen=6.5e-3, iT=0.35, meshdata=False)
#Voxel_image(scaffold, 'oldlace', 'darkgray')

#region Network creations
from skimage.feature import peak_local_max
from functions import region3D_centroids, peaks_MBA_filter, Watershed_regions
from porespy.networks import snow2, regions_to_network
from edt import edt3d
M0 = np.logical_not(scaffold)
peaks = region3D_centroids(peak_local_max(edt3d(M0, black_border=True), indices=False)*1*edt3d(M0, black_border=True))[0]
peaks = peaks_MBA_filter(M0, peaks=peaks, minimum=2.0)[0]
snow_net = snow2(M0, boundary_width=0, accuracy='high',
                 voxel_size=voxel_lenght)
snow_net = ps.networks.label_boundaries(snow_net['network'], labels=[['left', 'right'], ['front', 'back'], ['top', 'bottom']])
Net = regions_to_network(Watershed_regions(M0, peaks, return_img=False), voxel_size=voxel_lenght, accuracy='high')
Net = ps.networks.label_boundaries(Net, labels=[['left', 'right'], ['front', 'back'], ['top', 'bottom']])

#todo tentar perceber as diferenças entre os diâmetros
#O equivalente refere_se ao volume equivalente de uma esfera

#Net = pnm.io.network_from_porespy(Net)
Net = pnm.io.network_from_porespy(snow_net)
print('Network processed')
#endregion

#region Geometrical modifications
Net.add_model_collection(pnm.models.collections.geometry.spheres_and_cylinders)

#Net.add_model(propname='throat.hydraulic_size_factors',
                #model=pnm.models.geometry.hydraulic_size_factors.spheres_and_cylinders)
#todo if I do this here, later raises: KeyError: 'diameter not found'
#todo so is the previous condition doing anything?

Net['pore.diameter'] = Net['pore.equivalent_diameter']
Net['throat.diameter'] = Net['throat.equivalent_diameter']
Net.add_model(propname='throat.length',
              model=pnm.models.geometry.throat_length.spheres_and_cylinders)
Net.add_model(propname='throat.lens_volume',
             model=pnm.models.geometry.throat_volume.lens)
Net.add_model(propname='throat.total_volume',
              model=pnm.models.geometry.throat_volume.cylinder)
Net.add_model(propname='throat.volume',
             model=pnm.models.misc.difference,
             props=['throat.total_volume', 'throat.lens_volume'])
print('Geometry adjusted')
#endregion

#region Phase and physics definition
water = pnm.phase.Water(network=Net, name='water')
water.add_model_collection(pnm.models.collections.geometry.spheres_and_cylinders)
water.add_model_collection(pnm.models.collections.physics.standard)
water['pore.temperature'] = 297
water['throat.temperature'] = 297 #20ºC = 293K
#density
water.add_model(propname='pore.density',
                model=pnm.models.phase.density.water_correlation)
#viscosity
water.add_model(propname='pore.viscosity',
                model=pnm.models.phase.viscosity.water_correlation)
#hydraulic conductance
Net.add_model(propname='throat.hydraulic_size_factors',
                model=pnm.models.geometry.hydraulic_size_factors.spheres_and_cylinders)
#todo I need to do this here in order to prevent errors

water.add_model(propname='throat.hydraulic_conductance',
                model=pnm.models.physics.hydraulic_conductance.generic_hydraulic)
water.regenerate_models()
#endregion

#region Setting Flow Algorithm
flow = pnm.algorithms.StokesFlow(network=Net, phase=water)
flow.set_rate_BC(pores=Net.pores('left'), rates=1/60e-10)
flow.set_value_BC(pores=Net.pores('right'), values=0)

#flow.set_rate_BC(pores=Net.pores('right'), values=1/60e-10)
#todo this raised: "solver did not converge"

flow.run()
print('Flow simulation performed')

This raises: "Exception: Your network is clustered, making Ax = b ill-conditioned"
The error occurs either if I use the SNOW algorithm or not. Is it possible to bypass?

Regards.

[jgostick]

The solution is in the error message. This means you have pore(s) which are not connected to the rest. This makes it impossible to solve the system of equations since they are 'ill conditioned'. If you do Net.check_network_health() what do you get? It should tell you a list of pores which are not connected.

[jgostick]

I'm sorry, I forgot that we removed check_network_health from the network class. You can now get network health informatoin by adding a pore scale model, like this:

pn.add_model(propname='pore.cluster_number',
             model=op.models.network.cluster_number)
pn.add_model(propname='pore.cluster_size',
             model=op.models.network.cluster_size)

These two models will add cluster number and size to the network, so if you have any cluster numbers other than 0, that means your network is flawed. You can decide which pores to trim by check the cluster size...you'll probably have one very large cluster with 99% of the pores, so you'll want to keep that cluster, but trim all pores with cluster numbers other than that. I am realizing that this issue common enough that it is worth a notebook example. I'll add it to the todo list.

[jgostick]

I made an example for this problem because it comes up often enough: https://openpnm.org/examples/reference/networks/managing_clustered_networks.html

[JorgeESantos]

Hello again Professor. By adding both cluster models you suggested and adding a variable health as follows:

`from functions import Trigo_scaffolds, Voxel_image, TPMS_engine, scatter_3Ddataset
import numpy as np
import porespy as ps
import openpnm as pnm

#Generating Binary image of a porous scaffold
#scaffold, voxel_lenght = Trigo_scaffolds(Geo='b', shape=[2,2,2], elem=27, meshdata=False, min_dimen=13e-3)
scaffold, voxel_lenght = TPMS_engine([2,2,2], elem=40, min_dimen=6.5e-3, iT=0.366, meshdata=False)
#Voxel_image(scaffold, 'oldlace', 'darkgray')

#region Network creations
from skimage.feature import peak_local_max
from functions import region3D_centroids, peaks_MBA_filter, Watershed_regions
from porespy.networks import snow2, regions_to_network
from edt import edt3d
M0 = np.logical_not(scaffold)
#SNOW Net
#snow_net = snow2(M0, boundary_width=0, accuracy='high',
#voxel_size=voxel_lenght)
#snow_net = ps.networks.label_boundaries(snow_net['network'], labels=[['left', 'right'], ['front', 'back'], ['top', 'bottom']])
#My NET

peaks = region3D_centroids(peak_local_max(edt3d(M0, black_border=True), indices=False)1edt3d(M0, black_border=True))[0]
peaks = peaks_MBA_filter(M0, peaks=peaks, minimum=2.0)[0]
Net = regions_to_network(Watershed_regions(M0, peaks, return_img=False), voxel_size=voxel_lenght, accuracy='high')
#Net = ps.networks.label_boundaries(Net, labels=[['left', 'right'], ['front', 'back'], ['top', 'bottom']])

Net = pnm.io.network_from_porespy(Net)
#Net = pnm.io.network_from_porespy(snow_net)
Net.add_model_collection(pnm.models.collections.geometry.spheres_and_cylinders)

#Network manipulation
pores_left = Net['pore.coords'][:, 0] <= np.amin([i for i in Net['pore.coords'][:, 0]]) + voxel_lenght
pores_right = Net['pore.coords'][:, 0] >= np.amax([i for i in Net['pore.coords'][:, 0]]) - voxel_lenght
Net.set_label(pores=pores_left, label='left')
Net.set_label(pores=pores_right, label='right')

Net.add_model(propname='pore.cluster_number',
model=pnm.models.network.cluster_number)
Net.add_model(propname='pore.cluster_size',
model=pnm.models.network.cluster_size)
health = pnm.utils.check_network_health(Net)`

The cluster size array returned an array with only two possible values: 177 or 179.

As for the cluster number it alternates between 0 and 1s. The health variable identified 177 disconnected pores. Does this means that my network is divided into two sub-networks? If so, it seems inappropriate to trim one of them due to the fact they have a similar cluster size. Looking forward to ear your advice.

Regards