Machine learning picks the fastest linear solvers for porous media simulations

Every step forward in time forces a computer to solve enormous systems of linear equations, and choosing the wrong numerical recipe for that task can waste hours of machine time. The work, published in Computational Science and Engineering, targets a problem that has long been left to expert intuition and trial and error.

The core difficulty is scale. In industry-relevant simulations of porous media, the linear systems generated at each step range from tens of thousands to billions of rows, which makes direct solution methods infeasible and iterative solvers mandatory. Iterative solvers for these coupled multiphysics problems depend critically on a component called a preconditioner, which transforms the system into a form the iterative method can crack quickly. A well-chosen preconditioner with well-tuned parameters delivers good scalability and efficiency; a poorly chosen one can converge slowly or fail outright. Because multiphysics problems are nonlinear, it is hard to know in advance which parameters will work, and solver behavior can shift dramatically as the simulation state evolves.

Preconditioners for such problems are typically assembled from sub-algorithms, each originally designed for a single physical process such as fluid flow, heat transfer or mechanical deformation. The number of ways to combine these building blocks grows exponentially. For the simulations reported in the study, the number of candidate solver configurations was on the order of ten thousand. Collecting enough performance data to identify the optimal configuration across all of them is justified only if the dataset will be reused extensively, since a single data point requires actually solving a linear system with a specific configuration, which can take seconds to minutes.

The new approach, developed by Yury Zabegaev, Inga Berre and Eirik Keilegavlen, sidesteps the need for a massive pre-collected dataset. Instead, it gathers performance data during the target simulation itself and continuously updates a machine learning model that guides solver selection. The pipeline consists of two stages: a classifier that predicts which configurations are likely to succeed and discards the rest, and a regression model that estimates the expected reward, defined as the negative logarithm of the solve time, for the surviving candidates. The configuration with the highest predicted reward is then applied to the next linear system.

The algorithm begins with a short initial exploration phase, limited to 64 linear systems solved with randomly selected configurations. After that, the machine learning pipeline takes over, choosing a solver for each Newton iteration or time step. Feedback from every solve attempt, successful or failed, is accumulated in batches of 64 and used to retrain the models, which are gradient boosting classifiers and regressors from the SciKit-Learn library. Because a simulation requires solving on the order of ten thousand linear systems, and retraining takes seconds while a single solve can take minutes, the authors judged full retraining on each batch to be affordable and more reliable than purely incremental updates.

The team evaluated the method on two time-dependent nonlinear model problems. The first couples fluid flow and heat transfer in a three-dimensional porous medium based on the SPE10 benchmark, a model with a challenging anisotropic and heterogeneous subsurface reservoir structure, with permeability values spanning six orders of magnitude. The second is a thermo-hydromechanics problem in fractured porous media, where fracture deformation is governed by frictional contact mechanics, allowing fractures to stick, slide or open. For the first problem the researchers varied the permeability field across 15 simulations; for the second they changed the locations of injection and production wells across 25 simulations, mimicking the sensitivity studies and optimization workflows common in practice.

The results were striking. Fewer than 3 percent of all available solver configurations were ever used, yet the algorithm identified well-performing ones and stuck with them, raising the success rate above 96 percent for both problem sequences. The mean successful solve time corresponded to configurations within the top 1.3 percent for the flow-heat problem and the top 2.4 percent for the contact-THM problem. The price was modest: 64 random solves during exploration plus a machine learning overhead of roughly 0.006 seconds per linear system for the first problem and 0.010 seconds for the second, values that did not grow with the amount of collected data.

The researchers also benchmarked their approach against two baselines. The first used a fixed solver configuration with default PETSc parameters, representing what a practitioner would run without solver selection. The second was an oracle, an otherwise identical algorithm granted access to all performance data from prior simulations, a setting that is not practical but shows what near-ideal selection could achieve. Both baselines achieved 99 percent success rates, showing that default options are robust for this class of problem, but the histograms revealed that defaults form a robust yet inefficient middle ground. After its necessary exploration phase, the learning algorithm overtook the PETSc-default baseline, and in the fractured-media sequence it even edged past the oracle, whose heavier retraining costs offset its informational advantage.

The study is careful about its scope. The experiments were run on a single CPU core, each sequence run was repeated five times with shuffled simulation orders and different random seeds, and the authors note that the algorithm is non-deterministic due to random exploration and noise in timing measurements. They also emphasize that the initial exploration carries risk: simulation states that appear later in a run may be poorly represented, and some configurations may never be tested. In the flow-heat problem, where 35 percent of configurations always fail, the algorithm needed roughly 200 linear systems to stabilize and avoid bad decisions, and it overtook the default baseline only after about 1500 linear systems.

Beyond raw speed, the collected performance data proved useful for understanding the solvers themselves. Analysis showed that simulation context, such as the average Fourier flux or the size of the time step, measurably affects solver performance, confirming that context features are necessary for accurate predictions. The data also revealed how individual algorithmic choices, such as the choice between the CPR and System-AMG preconditioning strategies or between V- and W-cycles in the mechanics multigrid subsolver, separate the run-time distributions, offering a way to refine the configuration space for future simulations. The authors suggest the approach holds particular value when expert knowledge on linear solver tuning is not readily available, letting simulation engineers and researchers obtain robust, efficient runs without manual tuning. All datasets and code, implemented in Python with PETSc4Py and PorePy, are publicly available on GitHub and Zenodo.

Subject of Research: Data-driven automated selection and tuning of preconditioned iterative linear solvers for multiphysics simulations in porous media

Article Title: Data-driven linear solver selection and performance tuning for multiphysics simulations in porous media

Article References: Zabegaev, Y., Berre, I., & Keilegavlen, E. (2026). Data-driven linear solver selection and performance tuning for multiphysics simulations in porous media. Computational Science and Engineering, 3(1), Article 4. https://doi.org/10.1007/s44207-026-00013-y

Image Credits: AI Generated

DOI: 10.1007/s44207-026-00013-y

Keywords: linear solvers, preconditioners, machine learning, porous media, multiphysics simulation, GMRES, algebraic multigrid, thermo-poromechanics, fractured reservoirs, computational science, PETSc, incremental learning