Problem:
I'm working on a project where I need to calculate the connectivity and tortuosity features for a reservoir. The reservoir is represented by porosity and permeability arrays, and I want to generate these two features for each grid in the reservoir. However, I'm encountering some TypeError issues with my current code. My porosity and permeability data are 3D with (51 layers in the z direction, 125 grids in the y direction, and 125 grids in the x direction) = (51, 125, 125). I envision the following steps:
For each grid in the reservoir, the algorithm checks its permeability and porosity values to determine if it exceeds the threshold. If it does, the grid is considered connected to other neighboring grids.
Once the connected grids are identified, the algorithm uses Dijkstra's algorithm to calculate the shortest path lengths from the current grid to all other connected grids. Dijkstra's algorithm iteratively explores the graph, finding the shortest path from the current grid to each of the connected grids.
The longest shortest path length from the current grid is then determined by finding the maximum value among all the shortest path lengths calculated in step 2.
Finally, the tortuosity value for the current grid is calculated by dividing the longest shortest path length by the shortest path length from the current grid to itself. This ratio represents the tortuosity of the path from the current grid to any other connected grid.
The Issue/ Error: Taking for ever...
Explanation:
In my code, I have implemented two functions: calculate_tortuosity() and calculate_connectivity().
The calculate_tortuosity() function takes the porosity and permeability arrays, as well as the source and destination coordinates and points. It creates a graph based on the porosity and permeability values, finds the shortest path between the source and destination nodes using Dijkstra's algorithm, and calculates the tortuosity as the ratio of the actual path length to the shortest possible length. One way to do this is to use a Dijkstra algorithm. The Dijkstra algorithm is a graph search algorithm that finds the shortest path between two points in a graph. In this case, the graph would be the porosity and permeability arrays. The Dijkstra algorithm would start at point A and then explore all of the neighboring points that have a porosity and permeability value greater than a certain threshold. The algorithm would continue exploring neighboring points until it reached point B. The points that the algorithm explored would be the points that the flow passed from point A to point B. In summary, the assigned value of tortuosity to a single grid represents how tortuous the path is from that grid to any other connected grid. The distance between grids is measured in terms of the shortest path length, which is the minimum number of grid-to-grid steps required to reach the destination. The tortuosity value captures the extent of path deviation or complexity in the flow path from one grid to another.
The calculate_connectivity() function calculates the connectivity of the pore space based on a given threshold. Two points in the porosity array are considered connected if they are neighbors and both have a porosity value greater than the threshold. The function uses a sparse graph representation and the connected components algorithm to determine the connectivity. For example, in a porosity array of dimension (3, 5, 5) with values between 0 and 0.25, the following points would be considered to be connected:
(0, 0, 0) and (0, 0, 1) (0, 1, 0) and (0, 1, 1) (2, 4, 3) and (2, 4, 4) The following points would not be considered to be connected:
(0, 0, 0) and (1, 0, 0) (0, 0, 0) and (0, 0, 2) (2, 4, 3) and (1, 4, 3) The reason why the first two points are not connected is because they are not neighbors. The reason why the third point is not connected is because the porosity value at (1, 4, 3) is less than the threshold value. Below is my attempt code using random values:
import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import connected_components
from scipy.sparse.csgraph import dijkstra
from scipy.spatial.distance import euclidean
import networkx as nx
import heapq
from scipy.sparse.csgraph import connected_components
connectivity_threshold = 0.125
tortuosity_threshold = 1
def calculate_tortuosity(permeability_field, porosity_field,
threshold=tortuosity_threshold):
# Get the dimensions of the fields
n_layers, n_rows, n_columns = permeability_field.shape
# Create an empty directed graph
graph = nx.DiGraph()
# Iterate over each grid in the fields
for k in range(n_layers):
for i in range(n_rows):
for j in range(n_columns):
# Get the permeability and porosity values of the
#current grid
permeability = permeability_field[k, i, j]
porosity = porosity_field[k, i, j]
# Calculate the effective permeability as the product of
#permeability and porosity
effective_permeability = permeability * porosity
# Check if the effective permeability exceeds the
#threshold
if effective_permeability > threshold:
# Add an edge from (k, i, j) to its neighboring
#grids
if i > 0 and permeability_field[k, i - 1, j] *
porosity_field[k, i - 1, j] > threshold:
graph.add_edge((k, i, j), (k, i - 1, j))
if i < n_rows - 1 and permeability_field[k, i + 1, j] * porosity_field[k, i + 1, j] > threshold:
graph.add_edge((k, i, j), (k, i + 1, j))
if j > 0 and permeability_field[k, i, j - 1] * porosity_field[k, i, j - 1] > threshold:
graph.add_edge((k, i, j), (k, i, j - 1))
if j < n_columns - 1 and permeability_field[k, i, j + 1] * porosity_field[k, i, j + 1] > threshold:
graph.add_edge((k, i, j), (k, i, j + 1))
if k > 0 and permeability_field[k - 1, i, j] * porosity_field[k - 1, i, j] > threshold:
graph.add_edge((k, i, j), (k - 1, i, j))
if k < n_layers - 1 and permeability_field[k + 1, i, j] * porosity_field[k + 1, i, j] > threshold:
graph.add_edge((k, i, j), (k + 1, i, j))
# Calculate the shortest path lengths using Dijkstra's algorithm
shortest_path_lengths =
dict(nx.all_pairs_shortest_path_length(graph))
# Create an empty array to store the tortuosity values
tortuosity_values = np.zeros_like(permeability_field, dtype=float)
# Iterate over each grid in the fields
for k in range(n_layers):
for i in range(n_rows):
for j in range(n_columns):
# Check if the grid is connected to other grids
if (k, i, j) in shortest_path_lengths:
# Calculate the shortest path length from (k, i, j)
#to all other grids
shortest_lengths = shortest_path_lengths[(k, i, j)]
# Calculate the tortuosity as the ratio of actual
#length to shortest length
tortuosity_values[k, i, j] =
np.max(list(shortest_lengths.values())) /
shortest_lengths[(k, i, j)]
if shortest_lengths[(k, i, j)] == 0:
tortuosity_values[k, i, j] = 0
return tortuosity_values
def calculate_connectivity(porosity,
threshold=connectivity_threshold):
n_layers, n_rows, n_columns = porosity.shape
connectivity = np.zeros((n_layers, n_rows, n_columns))
# Create a binary mask indicating porosity values above the
#threshold
mask = porosity > threshold
# Reshape the mask to a 2D array
mask_2d = mask.reshape(n_layers, -1)
# Convert the 2D mask to a CSR sparse matrix
csr_mask = csr_matrix(mask_2d)
# Calculate connected components
n_components, labels = connected_components(csr_mask)
# Reshape the labels back to the original shape
labels_3d = labels.reshape(n_layers, n_rows, n_columns)
# Assign connectivity values based on connected components
for k in range(n_layers):
for i in range(n_rows):
for j in range(n_columns):
if mask[k, i, j]:
connectivity[k, i, j] = labels_3d[k, i, j] + 1
return connectivity
def main():
# Create the porosity and permeability arrays.
porosity = np.random.rand(51, 125, 125)
permeability = np.random.rand(51, 125, 125)
tortuosity = calculate_tortuosity(permeability, porosity,
threshold=tortuosity_threshold)
connectivity = calculate_connectivity(porosity,
threshold=connectivity_threshold)
# Save the tortuosity and connectivity arrays.
np.save("tortuosity.npy", tortuosity)
np.save("connectivity.npy", connectivity)
if __name__ == "__main__":
main()