Open PNM not optimized for pore network with trenches?

[babaganousheggplant]

Good day.

First of all the program is really interesting and does nice work in helping visualize pore networks in a different way. I would just like to ask if there is a way to use the program to describe trenched networks? I have been dabbling with help from ChatGPT on how to do this and ChatGPT suggested to use masks. The problem is that after several iterations, I only end up with a network that has no trench or no network at all. I would just like to ask if there is a better way to do this or if OpenPNM is only optimized for regular networks? Also, I am open to discuss with the developers if need be. Anyhow here is the code just in case you are wondering:

import openpnm as op
import numpy as np
import matplotlib.pyplot as plt

Step 1: Create a cubic network (regular lattice)

shape = [10, 10, 10] # Smaller network size for easier visualization
spacing = 1e-6 # Smaller spacing between pores
pn = op.network.Cubic(shape=shape, spacing=spacing)

Step 2: Define the trench region

For simplicity, let's remove pores in a trench along the x-direction

Here, we will remove pores where the x-coordinate is between 9 and 11

pore_coords = pn['pore.coords'] # Get the coordinates of all pores

Step 3: Define the trench range (in the x-direction, for example)

trench_x_range = (9, 11) # Define the x range for the trench

Step 4: Create a mask that is True for pores not in the trench region

mask = (pore_coords[:, 0] < trench_x_range[0]) | (pore_coords[:, 0] > trench_x_range[1])

Step 5: Get the indices of pores to keep (those outside the trench)

pores_to_keep = pn.pores()[mask] # Select pores not in the trench

Step 6: Remove the unwanted pores

Update the pore coordinates

new_pore_coords = pore_coords[mask]
pn['pore.coords'] = new_pore_coords

Now update the throats: remove connections where at least one pore is removed

throats = pn['throat.conns']
mask_throats = np.isin(throats, pores_to_keep) # Keep only throats connecting to valid pores
valid_throats = np.all(mask_throats, axis=1) # Keep throats where both pores are valid

Update the network to remove invalid throats

pn['throat.conns'] = throats[valid_throats]

Step 7: Visualize the new network

Check the number of remaining pores

print(f"Number of remaining pores: {len(pn.pores())}")

Visualize the network in 3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

Plot all remaining pores (no sampling)

ax.scatter(pn['pore.coords'][:, 0], pn['pore.coords'][:, 1], pn['pore.coords'][:, 2], s=50, c='b', alpha=0.7)

Set axis limits for better visualization

ax.set_xlim([0, shape[0]])
ax.set_ylim([0, shape[1]])
ax.set_zlim([0, shape[2]])

Optionally, set a title

ax.set_title("3D Network with Trench Removed")

Show the plot

plt.show()

[hpc-dave]

If I understand correctly, trenching should work at least for visualization. The solver might have more problems with two independent pore clusters. The problem that you are running into is, that you don't update the pore-IDs stored in the throat array. The way it is now, some throats will point at pores that don't exist anymore.
For manipulating the networks, you should use the topotools module
The algorithms there make sure your network stays in a good shape. You could easily use the trim function:

import openpnm as op
# create your network
# and trim those pores that you want to remove
# ....
x_coords = network['pore.coords'][:, 0]
pores_to_remove =  (x_coords >= trench_x_range[0]) & (x_coords <= trench_x_range[1])
pores_to_remove = np.where(pores_to_remove)[0]   # mainly for cross-checking
op.topotools.trim(network, pores=pores_to_remove)

Hope that helps!