Ultra-Fast 3D Porous Media Generation: a GPU- Accelerated List-Indexed Explicit Time-Stepping QSGS Algorithm

Efficient generation of high-resolution synthetic microstructures is essential in digital rock physics, yet classical Quartet Structure Generation Set (QSGS) algorithms become prohibitively expensive on large three-dimensional grids. We develop a list-indexed explicit time-stepping (LIETS) formulation of QSGS that restricts stochastic growth operations to an explicit active front instead of the entire voxel grid. The method is implemented in Python using NumPy on CPUs and CuPy on GPUs, and incorporates seed-spacing control via diamond dilation together with a volume-fraction-dependent directional growth probability. For a 400^3 domain, LIETS reduces generation time from tens of minutes for a serial CPU implementation and several minutes for vectorized CPU and GPU QSGS to about 24 s on a consumer-grade RTX 4060, achieving peak throughputs up to 2.7x10^7 nodes/s. A Fontainebleau sandstone benchmark at 500^3 resolution shows that LIETS reproduces the dependence of pore and grain size distributions on seed spacing (optimal s=30 voxels) and yields permeability-porosity trends within the experimental envelope and consistent with previously published Fast-QSGS results.

💡 Research Summary

The paper addresses a critical bottleneck in digital rock physics (DRP): the generation of high‑resolution synthetic microstructures using the Quartet Structure Generation Set (QSGS). Classical QSGS implementations update the entire voxel grid at every growth iteration, drawing random numbers, applying masks, and shifting arrays for each direction regardless of whether a voxel lies near the active growth front. This full‑field approach leads to computational costs that scale with the total number of cells (N_cells) and quickly becomes prohibitive for three‑dimensional domains of several hundred voxels per side.

To overcome this limitation, the authors propose a List‑Indexed Explicit Time‑Stepping (LIETS) formulation. The key idea is to maintain an explicit list A of active sites—solid voxels that have at least one neighboring void voxel. During each global iteration, for every growth direction i, the algorithm gathers the neighboring void candidates of the sites in A, draws a single random number per candidate, and converts the voxel to solid if the number is below the directional growth probability G_i. Newly solidified voxels are added to a temporary set A_new, which after pruning (removing duplicates and sites that have become fully surrounded) replaces A for the next iteration. By restricting all stochastic operations to the active front, the per‑step work scales with the size of the front (N_active) rather than the entire domain, yielding substantial speedups especially when the solid fraction is low.

Implementation is carried out in Python using NumPy for CPU execution and CuPy for GPU acceleration. CuPy enables the same high‑level code to run on consumer‑grade GPUs while handling list management and neighbor indexing directly on the device, thereby minimizing global memory traffic. The authors also incorporate seed‑spacing control via diamond dilation and make the directional growth probabilities volume‑fraction dependent, preserving the original QSGS statistical behavior.

Performance benchmarks demonstrate dramatic improvements. For a 400³ voxel domain, a serial CPU version requires tens of minutes, a vectorized CPU/GPU version needs several minutes, whereas the LIETS GPU implementation on an RTX 4060 completes in approximately 24 seconds, achieving peak throughputs of up to 2.7 × 10⁷ nodes per second. Accuracy is validated using a Fontainebleau sandstone benchmark at 500³ resolution. By varying the seed‑spacing parameter s, the authors show that the pore‑size and grain‑size distributions match those reported for Fast‑QSGS, with an optimal spacing of s = 30 voxels. Subsequent pore‑scale flow simulations (pnflow) produce permeability‑porosity curves that lie within experimental envelopes and agree closely with previously published Fast‑QSGS results, confirming that the algorithmic acceleration does not compromise physical realism.

The paper’s contributions are threefold: (1) a novel algorithmic redesign that reduces QSGS’s computational complexity from O(N_cells · N_dir) to O(N_active · N_dir); (2) a portable, high‑level Python/CuPy code base that can be readily integrated into existing DRP pipelines; and (3) a thorough validation showing that statistical and petrophysical properties are preserved. The authors discuss extensions such as multi‑phase simultaneous growth, anisotropic directional probabilities, and scaling to multi‑GPU clusters. They also suggest coupling LIETS directly with lattice‑Boltzmann solvers for on‑the‑fly flow simulations, opening the door to real‑time uncertainty quantification in large‑scale digital rock studies. In summary, LIETS delivers order‑of‑magnitude speedups for 3D porous media generation while maintaining the stochastic fidelity of QSGS, making high‑resolution synthetic rock modeling feasible on mainstream desktop hardware.