Algebraic Wavefront Parallelization for ILU(0) Smoothing in Reservoir Simulation

Incomplete factorization methods are an important part of the linear solver strategy in reservoir simulation. It has been shown earlier that the inherited pressure-decoupling effect of (block)-ILU(0) plays an important role for the convergence of efficient linear solvers like System-AMG or CPR. From Black-Oil to coupled geomechanics. With these specific linear systems, this decoupling is a by-product of the row-wise ILU-elimination. However, this also makes ILU sequential in nature, which is a problem on parallel compute hardware. The parallelization of ILU methods has been a field of active research - and it still is. Various approaches are reported in the literature. All exploit inherited parallelism in the sparse systems to solve. Either by reordering the system accordingly or by setting synchronization points induced by the underlying structure (so-called wavefronts). All of these approaches have certain disadvantages and advantages regarding parallel efficiency and numerical robustness. It depends on the application which approach is best-suited. Re-ordering approaches affect the elimination order. Hence, they can have significant robustness impacts for AMG in reservoir simulation. Wavefront parallelizations guarantee equivalence to the sequential method. However, they either require the parallelization structure to be induced by the geometry. This may be challenging in unstructured cases and voids a main advantage of AMG. Or they perform a row-wise data-dependency scan, with a resulting amount of blocking communication. In this paper, we are going to present a wavefront parallelization for (block-) ILU(0) that does not perform its dependency scan by considering groups of rows. The resulting wavefront setup works analogously to aggregative AMG setups, just with additional constraints. The outcome is a data-dependency graph where one can control the compromise between the frequency of data-exchange and wait-time. The equivalence to the sequential ILU(0) algorithm is still guaranteed. While this approach can't compete with the parallelizability of methods like Jacobi-relaxations, it can exploit inherited parallelism with ILU(0) for both OpenMP and MPI. It maintains the numerical properties of the original algorithm. We will demonstrate both with test problems as well as with ones from industrial reservoir simulations.