Tortuosity via LBM diffusion solver¶

This example computes the diffusional tortuosity of a synthetic porous medium using the D3Q7 BGK lattice Boltzmann solver.

In [1]:

import porespy as ps
import matplotlib.pyplot as plt
import poromics

import porespy as ps import matplotlib.pyplot as plt import poromics

Generate a porous medium image¶

In [2]:

im = ps.generators.blobs(shape=[100, 100, 1], porosity=0.6, blobiness=0.5, seed=42)

fig, ax = plt.subplots(figsize=(4, 4))
ax.imshow(im[:, :, 0], cmap="viridis", interpolation="nearest")
ax.set_title(f"Porous medium (porosity = {im.mean():.2f})")
plt.tight_layout()

im = ps.generators.blobs(shape=[100, 100, 1], porosity=0.6, blobiness=0.5, seed=42) fig, ax = plt.subplots(figsize=(4, 4)) ax.imshow(im[:, :, 0], cmap="viridis", interpolation="nearest") ax.set_title(f"Porous medium (porosity = {im.mean():.2f})") plt.tight_layout()

Compute tortuosity¶

The LBM solver accepts physical SI units. Here we use a diffusivity of $10^{-9}$ m²/s (typical for ions in water) and a voxel size of 1 µm.

In [3]:

result = poromics.tortuosity_lbm(im, axis=1, D=1e-9, voxel_size=1e-6)
print(result)
print(f"Tortuosity: {result.tau:.4f}")
print(f"Effective diffusivity (D_eff/D_0): {result.D_eff:.6f}")
print(f"Formation factor: {result.formation_factor:.4f}")

result = poromics.tortuosity_lbm(im, axis=1, D=1e-9, voxel_size=1e-6) print(result) print(f"Tortuosity: {result.tau:.4f}") print(f"Effective diffusivity (D_eff/D_0): {result.D_eff:.6f}") print(f"Formation factor: {result.formation_factor:.4f}")

2026-03-09 01:25:55.53 | WARNING  | poromics._metrics:tortuosity_lbm:555 - Trimmed 972 non-percolating pore voxels from the image.
[I 03/09/26 01:25:55.670 22015075] [shell.py:_shell_pop_print@23] Graphical python shell detected, using wrapped sys.stdout

[Taichi] version 1.7.4, llvm 15.0.7, commit b4b956fd, osx, python 3.12.9

[Taichi] Starting on arch=metal

TortuosityResult(tau=5.1417, D_eff=0.096290, axis=1)
Tortuosity: 5.1417
Effective diffusivity (D_eff/D_0): 0.096290
Formation factor: 10.3853

Visualize the concentration field¶

The steady-state concentration field shows how species diffuses through the tortuous pore network.

In [4]:

result.plot_concentration(z=0)

result.plot_concentration(z=0)

Out[4]:

(<Figure size 750x400 with 4 Axes>,
 array([<Axes: title={'center': 'Pore structure'}>,
        <Axes: title={'center': 'Concentration field'}>], dtype=object))

Using the low-level solver¶

For more control, use TransientDiffusion directly. This lets you access intermediate state, step manually, or customize convergence criteria.

In [5]:

from poromics.simulation import TransientDiffusion

solver = TransientDiffusion(im, axis=1, D=1e-9, voxel_size=1e-6)
solver.run(n_steps=50_000, tol=1e-2)

print(f"Converged: {solver.converged}")
print(f"Iterations: {solver.n_iterations}")
print(f"Physical time step: {solver.dt:.2e} s")
print(f"Midplane flux: {solver.flux(axis=1):.6f}")

from poromics.simulation import TransientDiffusion solver = TransientDiffusion(im, axis=1, D=1e-9, voxel_size=1e-6) solver.run(n_steps=50_000, tol=1e-2) print(f"Converged: {solver.converged}") print(f"Iterations: {solver.n_iterations}") print(f"Physical time step: {solver.dt:.2e} s") print(f"Midplane flux: {solver.flux(axis=1):.6f}")

Converged: True
Iterations: 7001
Physical time step: 2.50e-04 s
Midplane flux: 0.000304