Determining tortuosity using geometric domain decomposition¶
Import necessary packages and functions
import numpy as np
import porespy as ps
import openpnm as op
import matplotlib.pyplot as plt
ps.visualization.set_mpl_style()
[01:08:12] ERROR PARDISO solver not installed, run `pip install pypardiso`. Otherwise, _workspace.py:56 simulations will be slow. Apple M chips not supported.
Generate test image¶
# im = ps.generators.fractal_noise(shape=[100,100,100], seed=1)<0.8
np.random.seed(1)
im = ps.generators.blobs(shape=[100, 100, 100], porosity=0.7)
plt.imshow(ps.visualization.sem(~im))
plt.axis('off');
Run the function on the image¶
from porespy.beta import tortuosity_bt
out = tortuosity_bt(im, block_size=50)
print(out)
[ ] | 0% Completed | 213.49 us
[ ] | 0% Completed | 101.68 ms
[ ] | 0% Completed | 204.38 ms
[ ] | 0% Completed | 305.29 ms
[ ] | 0% Completed | 408.87 ms
[ ] | 0% Completed | 529.60 ms
[ ] | 0% Completed | 658.59 ms
[ ] | 0% Completed | 925.74 ms
[ ] | 0% Completed | 1.03 s
[ ] | 0% Completed | 1.13 s
[ ] | 0% Completed | 1.24 s
[ ] | 0% Completed | 1.34 s
[ ] | 0% Completed | 1.44 s
[ ] | 0% Completed | 1.54 s
[ ] | 0% Completed | 1.64 s
[ ] | 0% Completed | 1.74 s
[ ] | 0% Completed | 1.85 s
[ ] | 0% Completed | 1.95 s
[ ] | 0% Completed | 2.05 s
[############# ] | 33% Completed | 2.18 s
[############# ] | 33% Completed | 2.29 s
[############# ] | 33% Completed | 2.39 s
[############# ] | 33% Completed | 2.56 s
[############# ] | 33% Completed | 2.71 s
[############# ] | 33% Completed | 2.89 s
[############# ] | 33% Completed | 3.02 s
[############# ] | 33% Completed | 3.12 s
[############# ] | 33% Completed | 3.23 s
[############# ] | 33% Completed | 3.33 s
[############# ] | 33% Completed | 3.43 s
[############# ] | 33% Completed | 3.55 s
[####################### ] | 58% Completed | 3.65 s
[########################## ] | 66% Completed | 3.76 s
[########################## ] | 66% Completed | 3.86 s
[########################## ] | 66% Completed | 3.97 s
[########################## ] | 66% Completed | 4.16 s
[########################## ] | 66% Completed | 4.47 s
[########################## ] | 66% Completed | 4.57 s
[########################## ] | 66% Completed | 4.70 s
[########################## ] | 66% Completed | 4.81 s
[########################## ] | 66% Completed | 4.91 s
[########################## ] | 66% Completed | 5.01 s
[#################################### ] | 91% Completed | 5.12 s
[########################################] | 100% Completed | 5.22 s
[1.392510902972281, 1.3999150512069574, 1.3606684244878242]
The three values in the returned list are the tortuosity values in the x, y, and z-direction respectively. However, there is a more useful form of this function.
from porespy.beta import analyze_blocks
out2 = analyze_blocks(im, block_size=50)
out2.iloc[:10,:]
[ ] | 0% Completed | 281.10 us
[ ] | 0% Completed | 102.32 ms
[ ] | 0% Completed | 211.71 ms
[ ] | 0% Completed | 334.44 ms
[ ] | 0% Completed | 517.29 ms
[ ] | 0% Completed | 791.73 ms
[ ] | 0% Completed | 909.18 ms
[ ] | 0% Completed | 1.05 s
[ ] | 0% Completed | 1.15 s
[ ] | 0% Completed | 1.26 s
[ ] | 0% Completed | 1.36 s
[########## ] | 25% Completed | 1.46 s
[############# ] | 33% Completed | 1.57 s
[############# ] | 33% Completed | 1.67 s
[############# ] | 33% Completed | 1.78 s
[############# ] | 33% Completed | 1.92 s
[############# ] | 33% Completed | 2.17 s
[############# ] | 33% Completed | 2.29 s
[############# ] | 33% Completed | 2.39 s
[############# ] | 33% Completed | 2.52 s
[############# ] | 33% Completed | 2.63 s
[############# ] | 33% Completed | 2.73 s
[############# ] | 33% Completed | 2.83 s
[#################### ] | 50% Completed | 2.93 s
[#################### ] | 50% Completed | 3.03 s
[########################## ] | 66% Completed | 3.14 s
[########################## ] | 66% Completed | 3.24 s
[########################## ] | 66% Completed | 3.36 s
[########################## ] | 66% Completed | 3.50 s
[########################## ] | 66% Completed | 3.60 s
[########################## ] | 66% Completed | 3.70 s
[########################## ] | 66% Completed | 3.88 s
[########################## ] | 66% Completed | 3.98 s
[########################## ] | 66% Completed | 4.10 s
[########################## ] | 66% Completed | 4.21 s
[########################## ] | 66% Completed | 4.31 s
[########################## ] | 66% Completed | 4.42 s
[################################# ] | 83% Completed | 4.52 s
[########################################] | 100% Completed | 4.62 s
| eps_orig | eps_perc | g | tau | volume | length | axis | time | |
|---|---|---|---|---|---|---|---|---|
| 0 | 0.705176 | 0.705176 | 25.439633 | 1.414264 | 125000 | 50 | 0 | 1.457192 |
| 1 | 0.721312 | 0.721312 | 26.571451 | 1.385007 | 125000 | 50 | 0 | 1.467454 |
| 2 | 0.705616 | 0.705592 | 26.630394 | 1.351823 | 125000 | 50 | 0 | 1.486075 |
| 3 | 0.651352 | 0.651344 | 21.896334 | 1.517690 | 125000 | 50 | 0 | 1.451312 |
| 4 | 0.706072 | 0.706024 | 25.394748 | 1.418468 | 125000 | 50 | 1 | 1.648835 |
| 5 | 0.717600 | 0.717600 | 27.020315 | 1.354990 | 125000 | 50 | 1 | 1.443812 |
| 6 | 0.691104 | 0.691016 | 24.780591 | 1.422723 | 125000 | 50 | 1 | 1.430844 |
| 7 | 0.679584 | 0.679256 | 22.810414 | 1.519303 | 125000 | 50 | 1 | 1.587346 |
| 8 | 0.695480 | 0.695480 | 25.926987 | 1.368600 | 125000 | 50 | 2 | 1.478277 |
| 9 | 0.683104 | 0.683104 | 23.386926 | 1.490245 | 125000 | 50 | 2 | 1.501230 |
The analyze_blocks function returns a DataFrame containing the tortuosity, diffusive conductance, and porosity values of each block, which can be used to obtain the tortuosity by using OpenPNM to solve the resistor network as follows. Start by assigning the diffusive conductance values to the network, then run the simulation:
net = op.network.Cubic([2, 2, 2])
air = op.phase.Phase(network=net)
air['throat.diffusive_conductance']=np.array(out2.iloc[:,2]).flatten()
fd = op.algorithms.FickianDiffusion(network=net, phase=air)
fd.set_value_BC(pores=net.pores('left'), values=1)
fd.set_value_BC(pores=net.pores('right'), values=0)
fd.run()
rate_inlet = fd.rate(pores=net.pores('left'))[0]
# the length of one block is removed from the total length since the network edge begins
# in the center of the first slice and ends in the center of the last slice, so the image
# length is decreased
L = im.shape[1] - 50
A = im.shape[0] * im.shape[2]
d_eff = rate_inlet * L /(A * (1-0))
e = im.sum() / im.size
D_AB = 1
tau_gdd = e * D_AB / d_eff
tau_gdd
1.3606684244878242
The direct calculation can also be done on the same image, and the results can be compared.
direct = ps.simulations.tortuosity_fd(im, 0)
tau_direct = direct.tortuosity
tau_direct
With similar results, the main benefit to using the “block and tackle” method is the time save on larger images.